format: text

text:
# Continuum-particle hybrid coupling for mass, momentum and energy transfers in unsteady fluid flow

## I Introduction

A wide range of systems with important applications are governed by a fine interplay between the atomistic processes occurring within a small region of the system and the slow dynamics occurring within the bulk. A large list of examples arise in complex flows near interfaces (polymers or colloids near surfaces, wetting, drop formation, melting, crystal growth from a fluid phase, moving interfaces of immiscible fluids or membranes, to name only a few). The computational expense of realistic-size simulations of these problems via standard molecular dynamics (MD) is prohibitive, and such kind of studies require new algorithms which can retain the benefit of the atomistic description of matter where it is really needed, while treating the bulk of the system by much less costly continuum fluid mechanics methods.

Several hybrid algorithms of this sort have been proposed in the recent literature. In general, to couple the particle region P and the continuum region C, such hybrid schemes use an overlapping region comprised of two buffers C$``$P and P$``$C, which account for the two-way transfer of information: from C to P and vice versa (see Fig. 1). While the P$``$C transfer essentially consists of a coarse-graining procedure, at C$``$P one needs to reconstruct the dynamics of a large collection of particles with only the limited prescription from the C region as input. Moreover, in performing this reconstruction, the number of unphysical artifacts added (as Maxwell dæmons) should be minimized far as possible. This task is very complicated and represents in fact the heart of any hybrid scheme.

Hybrid algorithms for fluids are relatively recent. The elegant method introduced by Garcia et al. Gar00  for rarefied gases couples fluxes arising from a direct simulation Monte Carlo (DSMC) scheme to another region described by computational fluid dynamics (CFD). The DSMC is set at the finest grid scale of an adaptive mesh refinement hierarchy, while a CFD semi-implicit solver is used at the upper level scales. In passing, we note that the scheme may, in principle, be implemented using an (MD-Continuum) liquid description, although in this case the C-solver must be completely explicit to avoid having to change the particle’s energy in the iterations of the implicit scheme.

In the case of liquids, the state of the art is relatively less developed due to the complications arising from the interparticle forces. A pioneering work by O’Connel and Thompson Tho95  coupled momentum by imposing the local continuum velocity at C$``$P via a crude constraint Lagrangian dynamics. Hadjiconstantinou and Patera Had97  introduced a reservoir region to impose boundary conditions on the P region (this reservoir being the equivalent of the C$``$P domain defined here). While in residence in the reservoir, particles were given at each time step, a velocity drawn from a Maxwellian distribution with mean and variance consistent with the velocity and temperature of the C-flow. To obtain the boundary condition for the C region the authors used a low order polynomial to smooth the field variables derived from P at the P$``$C region. In order to match the boundary conditions for both the P and C regions, Hadjiconstantinou and Patera Had97  implemented an iterative scheme (based on the Schwarz alternating method) that is suitable for steady incompressible flows. Liao and Yip Liao98  proposed a sophisticated method (called the thermodynamic field estimator) to extract continuum fields from the particle data by means of maximum likehood inference. This idea may be used to ameliorate the P$``$C coupling when the flow presents large gradients, albeit at a rather large computational cost. To transfer momentum on the P region Liao and Yip Liao98  proposed a new Maxwell dæmon, called reflecting particle method. A drawback is that the pressure gradient is then an outcome of the simulation, rather than an input. Finally, Flekkoy et al. Flek00  used the idea of coupling-through-fluxes and also implemented mass transfer. However, energy transfer was still not allowed and only steady flows were considered. The main purpose of the present work is to broaden the scope of such hybrid schemes towards a general description allowing mass, momentum and energy coupling in unsteady flows.

A question of central interest is to decide what kind of information needs to be transfered at C$``$P and P$``$C. There are essentially two possibilities, to transfer either generalised forces (fluxes of conserved quantities) or the local values of the averaged-variables. Both kinds of approaches can be found in the published literature. Here, in the context of energy transfer, we show that under unsteady flows it is not sufficient to impose the local C-quantities at the boundary of P; instead, it is necessary to couple through fluxes. Another possible benefit of flux-coupling was pointed out by Flekkoy et al. Flek00  who stated that this procedure transcends the problem of working with fluids whose constitutive relations may be only partially or incompletely known. Although we agree that the flux-based coupling is the correct matching procedure, we show nevertheless that if the transport coefficients at C and P are disparate enough, the hybrid scheme fails to couple the time evolution of both domains. Hence, in such cases, the evaluation of transport coefficients (using standard microscopic techniques, at least for the range of densities and temperatures under study) is an unavoidable requirement for the correct behaviour of the hybrid scheme.

The rest of the paper proceeds as follows. The equations governing the continuum and particle regions and the averaging procedures are presented in Sec. II. The core of the scheme, describing the C$``$P coupling for momentum, energy and mass fluxes, is presented in Sec III. General requirements on the coarse-graining length and time scales are discussed in Section IV. The unsteady flows under which the scheme has been tested (decay of longitudinal and transversal waves) are presented in Sec. V and in Sec. VI we discuss the results of these tests as far as the main hydrodynamic and thermodynamic variables are concerned. Finally Sec. VII is devoted to conclusions.

## II  Overview and geometry of the hybrid coupling scheme

The domain decomposition of the hybrid scheme is depicted in Fig. 1. Two regions need to be distinguished: the particle region P, and the continuum region C. Region P is composed of an ensemble of particles interacting through prescribed interparticle potentials and evolving in time through Newtonian dynamics. In order to illustrate the coupling procedure a Lennard-Jones (LJ) fluid will be considered. Within P a number $`N(t)`$ of particles, located at $`𝐫=\left\{𝐫_𝐢\right\}`$ (the subscript $`i`$ denoting the $`i^{th}`$ particle) interacts through the LJ potential $`\psi (r)=4ϵ^1\left[(\sigma /r)^{12}(\sigma /r)^6\right]`$. Each particle has a mass $`m`$, velocity $`𝐯_𝐢`$ and energy $`ϵ_i=\frac{1}{2}mv_i^2+\mathrm{\Sigma }_j\psi (r_{ij})`$ ($`𝐫_{ij}=𝐫_j𝐫_i`$). Their equations of motion,

$`\dot{𝐫}_i`$ $`=`$ $`𝐯_i,`$ (1)
$`\dot{𝐯}_i`$ $`=`$ $`𝐟_i/m=\mathrm{\Sigma }_{j=1}^N{\displaystyle \frac{d\psi (r_{ij})}{dr_{ij}}}{\displaystyle \frac{𝐫_{ij}}{r_{ij}}},`$ (2)
are solved via standard molecular dynamics (MD) at time steps $`\mathrm{\Delta }t_P10^3\tau `$, where $`\tau =(m\sigma ^2/ϵ)^{1/2}`$ is the characteristic time of the LJ potential. Throughout the rest of the paper, all quantities will be expressed in reduced units of the LJ potential: $`\tau (=0.45\times 10^{13}`$s), $`\sigma (=3.305\times 10^{12}`$cm), $`ϵ`$, $`m(=6.63\times 10^{23}`$g) and $`ϵ/k_B(=119.18`$K) for time, length, energy, mass and temperature , respectively (the numerical values correspond to argon).

On the other hand, within the C region the relevant variables are the macroscopic local densities associated with the conserved quantities, the number density $`\rho (𝐑,𝐭)`$, the energy density $`\rho e(𝐑,𝐭)`$ and the momentum density $`𝐣(𝐑,𝐭)`$ (related to the local mean velocity $`𝐮`$ by $`𝐣=\rho 𝐮`$). In what follows the spatial coordinates of the macroscopic fields are denoted by capital letters $`𝐑`$, while the the microscopic coordinates are designated by lower case letters. The conservation laws for the local densities are

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\rho 𝐮,`$ (3)
$`{\displaystyle \frac{𝐣}{t}}`$ $`=`$ $`\left(\mathrm{𝐣𝐮}+𝚷\right),`$ (4)
$`{\displaystyle \frac{\rho e}{t}}`$ $`=`$ $`\left(\rho e𝐮+𝚷𝐮+𝐪\right),`$ (5)
where the specific energy $`e=u^2/2+3T/2+\varphi `$, includes the translational energy, the thermal kinetic energy and the potential energy $`\varphi `$. The momentum flow contains contributions from convection $`\mathrm{𝐣𝐮}`$ and the pressure tensor $`𝚷=P\mathbf{\hspace{0.17em}1}+𝝉`$, the latter including the local hydrostatic pressure $`P(𝐑,𝐭)`$ and the viscous stress tensor, which satisfy a Newtonian constitutive relation, as shown by previous previous MD descriptions of the LJ-fluid TRANSP ; Ciccotti ,

$$𝝉=\eta \left(𝐮+(𝐮)^T\frac{2}{3}𝐮\right)\xi 𝐮.$$
(6)
The energy current includes convection $`\rho e𝐮`$, dissipation $`𝚷𝐮`$ and conduction $`𝐪`$, which can be expressed in terms of the local temperature gradients and the thermal conductivity $`\kappa _c`$, through Fourier’s law $`𝐪=\kappa _cT(𝐑,t)`$. In order to close the above equations it is necessary to know the caloric $`e(\rho ,T)`$ and thermal equations of state $`P(\rho ,T)`$ and the constitutive relations for the transport coefficients (shear and bulk viscosities and thermal conductivity; $`\eta `$, $`\xi `$ and $`\kappa _c`$ respectively) in terms of a set of independent thermodynamic variables, such as $`\rho `$ and $`T`$. The equations of state for a LJ-fluid were extracted from Johnson et al. EOS  and the transport coefficients $`\eta `$, $`\kappa _c`$ and $`\xi `$ from Heyes TRANSP  and Borgelt et al. BulkVis . The variables relevant to the C region are the slower ones. Using any standard continuum fluid dynamics solver (e.g. based on a finite volume method) the evolution of the C-variables will be traced at time intervals $`\mathrm{\Delta }t_C>>\mathrm{\Delta }t_P`$ and evaluated within cells of volume $`V_l`$ whose size and location is given by the nodes of a certain mesh, $`\left\{𝐑_l\right\}`$, $`l=\{1,\mathrm{},M_c\}`$. It will be assumed that the size of the C$``$P and P$``$C regions are the same size as those of the cells used in the spatial discretisation of the selected continuum solver, say $`V_l=(\mathrm{\Delta }X)^3`$. In general both $`\mathrm{\Delta }X`$ and $`\mathrm{\Delta }t_C`$ may depend on the type of solver used for the C-region, or on the characteristic length of the particular phenomena under study, Nevertheless, various intrinsic constraints on $`\mathrm{\Delta }X`$ and $`\mathrm{\Delta }t_C`$ will be mentioned in Sec.IV.

#### II.0.1 Averages

Averages are needed in order to transfer information from the faster time and shorter length-scale particle dynamics to the slower and longer coarse-grained description. In order to deal with unsteady, non-equilibrium scenarios, averages need to be local on the slower time scale and in the coarse-grained spatial coordinates. For any particle variable, say $`\mathrm{\Phi }_i`$, we define the following averages:

$`\overline{\mathrm{\Phi }}(𝐑_l,t)`$ $``$ $`{\displaystyle \frac{1}{N_l}}{\displaystyle \underset{iV_l}{\overset{N_l}{}}}\mathrm{\Phi }_i,`$ (7)
$`\overline{\mathrm{\Phi }}(𝐑_l,t_C)`$ $``$ $`{\displaystyle \frac{1}{\mathrm{\Delta }t_C}}{\displaystyle _{t_C}^{t_C+\mathrm{\Delta }t_C}}\overline{\mathrm{\Phi }}(𝐑_l,t)𝑑t`$ (8)
where the summation in Eq. (7) is made over the $`N_l`$ particles inside the cell $`l`$.

The averaging procedure is needed to translate the P and C “languages” to and from each domain. This translation is done within the overlapping region, where the two descriptions of matter coexist (see Fig. 1). In particular, within the P$``$C cells the many degrees of freedom arising from the particle dynamics are coarse-grained to provide boundary conditions at the “upper” C-level. As long as the number of degrees of freedom is very much larger at P than at C, this operation is rather straightforward and is based on the microscopic derivation of continuum fluid dynamics Evans . We adopt the approach advocated by Flekkoy et al. Flek00 , in making the information transferred from P to C to be the coarse-grained particle-fluxes of conserved quantities. These are

$`\rho 𝐮𝐧_{PC}`$ $`=`$ $`{\displaystyle \frac{1}{V}}_{PC}\mathrm{\Sigma }_{i=1}^{N_{PC}}m𝐯_i𝐧_{PC}`$ (9)
$`𝚷𝐧_{PC}`$ $`=`$ $`{\displaystyle \frac{1}{V}}_{PC}\left(\mathrm{\Sigma }_{i=1}^{N_{PC}}m𝐯_i𝐯_i{\displaystyle \frac{1}{2}}\mathrm{\Sigma }_{i,j}^{N_{PC}}𝐫_{ij}𝐅_{ij}\right)𝐧_{PC}`$ (10)
$`𝐪𝐧_{PC}`$ $`=`$ $`{\displaystyle \frac{1}{V}}_{PC}\left(\mathrm{\Sigma }_{i=1}^{N_{PC}}mϵ_i𝐯_i{\displaystyle \frac{1}{2}}\mathrm{\Sigma }_{i,j}^{N_{PC}}𝐫_{ij}𝐯_i𝐅_{ij}\right)𝐧_{PC}`$ (11)
where $`N_{PC}`$ is the number of particles inside the P$``$C cell and $`𝐧_{PC}`$ is the surface vector shown in Fig. 1.

By contrast, within the C$``$ P cells, the particle dynamics must be modified to conform to the averaged-dynamics prescribed by the continuum description. In other words, one needs to construct a sort of “generalised boundary condition” for the particle dynamics. As pointed out in all previous papers on the subject Gar00 ; Tho95 ; Had97 ; Flek00 , this represents the most demanding challenge in that one needs to invent a way to reconstruct the microscopic dynamics of a large number of particles, based on only a few properties of the local continuum variables. Moreover, to ensure that the effect on the inner P-region is minimised, it is crucial to reduce as much as possible the unphysical artifacts, such as Maxwell dæmons, that are added to the particle dynamics at C$``$P. The present work is focused on this problem, which lies at the heart of any hybrid scheme.

## III  The C$``$P coupling

This part of the hybrid scheme can be alternatively stated as the imposition of generalised (mass, momentum and energy) boundary conditions on an MD simulation box. To deal with this task we have coupled the particle region to a collection of flows (with explicit analytical solution), which involve the whole set of conserved quantities exchanged (mass, momentum and energy). In this sense, in the present work our C-solver is not numerical but rather analytical. In particular, we use the initial (non-equilibrium) state imposed at P to calculate the time-dependent analytical solution at C. This C-flow is then imposed on the P region during the rest of the simulation, meaning that (appart from the initial state) the hybrid coupling used in the tests presented here works in one direction only (from C to P).

### III.1  Imposing fluxes under unsteady flows

Following Flekkoy et al. Flek00 , at C$``$ P we shall communicate fluxes of conserved quantities. These fluxes correspond to mass, momentum and energy transfers through the outer interface of the C$``$ P cell (the W-surface in Fig. 1). Flekkoy et al. Flek00  obtained these fluxes from the values of the continuum variables at the centre of the control cell, $`x=x_O`$, instead of at the exact position of the C$``$P interface, $`x=x_W`$. We have found that it is essential to take into account this apparently unimportant technicality when dealing with unsteady scenarios. Let us consider a general conservation equation, with

$$\frac{\varphi }{t}+𝐉_\varphi =S_\varphi ,$$
(12)
where $`J_\varphi `$ is the flux of $`\varphi `$ and the source term vanishes, $`S_\varphi =0`$, as in Eqs. (3)-(5). Integrating over the control cell C$``$P,

$$\frac{}{t}_V\varphi 𝑑V+_S𝐉_\varphi 𝐧𝑑s=0.$$
(13)
For illustration we shall restrict analysis to the one-dimensional (1D) situation depicted in Fig. 1. In this case one obtains,

$$\frac{}{t}_V\varphi 𝑑VA𝐉_{\varphi }^{}{}_{E}{}^{}𝐧_W=A𝐉_{\varphi }^{}{}_{W}{}^{}𝐧_W,$$
(14)
where the subscripts $`E`$ (east) and $`W`$ (west) denote that the variables are measured at $`x=x_E`$ and $`x=x_W`$ respectively. The surface vectors, $`𝐧_W`$ and $`𝐧_E`$ are shown in Fig. 1, and in Eq. (14) use has been made of $`𝐧_W=𝐧_E`$. The R.H.S. of Eq. (14) is the flux current of $`\varphi `$ through the interface $`W`$ of the control cell, which is precisely the (generalised) force we want to introduce on the particles at the C$``$P buffer. We note that only under steady flows does $`𝐉_{\varphi }^{}{}_{W}{}^{}=𝐉_{\varphi }^{}{}_{O}{}^{}`$ (to see this, integrate Eq. (12) from $`x=x_O`$ to $`x=x_W`$): hence only in this case does the evaluation of the fluxes at $`x_W`$ using the continuum variables at $`x_O`$ lead to the same converged steady state as if the variables at $`x_W`$ were used (although the transients may of course differ). It is possible to provide an estimate of the global error arising from evaluating the flux at a position $`x_O`$ shifted $`\delta \mathrm{\Delta }X`$ with respect $`x_W`$, over a certain time interval $`\mathrm{\Delta }t`$. In the case of the momentum equation the deviation of the stress contribution to the momentum flux $`J=𝐉𝐧`$ at any instant would be of order $`\mathrm{\Delta }JJ\delta \mathrm{\Delta }X`$, with $`\delta =|x_Ox_W|/\mathrm{\Delta }X`$ is the distance to the C$``$P interface ($`\delta =0.5`$ in Fig. 1). Assuming that the mean velocity field can be expressed as $`u=𝐮𝐧u^{(k)}\mathrm{exp}(ikx)`$ ($`k`$ being the dominant wavenumber) and using the Newtonian constitutive relation for the viscous tensor in Eq. (6), one obtains $`\mathrm{\Delta }J\eta _Lk^2u^{(k)}\delta \mathrm{\Delta }X`$, where $`\eta _L=4\eta /3+\xi `$ is the longitudinal viscosity. As a particular example we consider a longitudinal wave and evaluate the error along a cycle, $`\mathrm{\Delta }t=2\pi /(kc_s)`$ (where $`c_s`$ is the sound velocity). As a crude estimate, the accumulated error of the cell-averaged momentum $`𝐣𝐧=\rho _e\overline{u}`$ is of order $`\rho _e\mathrm{\Delta }\overline{u}\mathrm{\Delta }J\mathrm{\Delta }t/\mathrm{\Delta }X`$; using $`\rho _e=0.5`$, $`\eta _L1`$ and $`c_s5`$, one obtains $`\mathrm{\Delta }\overline{u}/u^{(k)}2\pi \delta \eta _Lk/(c_s\rho _e)0.3`$, for the typical wavenumbers considered here $`k0.2`$. Simulations carried out with the momentum flux evaluated at $`x_O`$ yield relative errors of the averaged velocity at C$``$P of the same order of magnitude as this estimate (see Sec VI.3.2 and Fig. 6b). We observe that most computational fluid dynamics (CFD) codes provide the continuum variables at the centre of the control cells (i.e., at $`x_O`$), so that in order to evaluate the fluxes pertaining to the C$``$P exchange it would first be necessary to make use of simple interpolation techniques.

The fluxes arising from the continuum equations (on the R.H.S. of Eqs. (15-17)) are imposed on the particle ensemble at the C$``$P cells through expressions involving atomistic variables (those on the R.H.S. of (15-17)).

$`m\mathrm{s}`$ $`=`$ $`A\rho 𝐮𝐧,`$ (15)
$`m\mathrm{s}𝐯^{}+{\displaystyle \underset{i}{\overset{N_{CP}}{}}}𝐅_i^{ext}`$ $`=`$ $`A\left(\rho \mathrm{𝐮𝐮}+𝚷\right)𝐧,`$ (16)
$`m\mathrm{s}ϵ^{}+{\displaystyle \underset{i}{\overset{N_{CP}}{}}}𝐅_i^{ext}𝐯_i𝐉_Q^{ext}𝐧`$ $`=`$ $`A\left(\rho 𝐮e+𝚷\mathbf{}𝒖+𝐪\right)𝐧.`$ (17)
where henceforth, $`𝐧`$ indicates the vector of the outermost interface of the C$``$P cell, pointing towards C. The nomenclature used here follows that of Flekkoy et al. Flek00 : $`\mathrm{s}(t)`$ indicates the number of particles inserted ($`\mathrm{s}>0`$) or removed ($`\mathrm{s}<0`$) from C$``$P per unit of time; the velocity of the inserted/removed particles is $`𝐯^{}`$; $`𝐅_i^{ext}`$ is the external force applied to each particle $`i`$ within the C$``$P cell. The total external force is $`\mathrm{\Sigma }^{N_{CP}}𝐅_i^{ext}`$, where the summation is over the $`N_{CP}(t)`$ particles inside C$``$P. Finally, $`ϵ^{}`$ indicates the energy of the inserted/removed particles and $`𝐉_Q^{ext}`$ refers to an externally imposed heat current.

As mentioned by Flekkoy et al. Flek00 , insertion of Eq. (15) into Eq.(16) shows that the rates of change of momentum due to convection and local stresses are correctly introduced if $`𝐯^{}=𝐮`$ and $`_i^{N_{CP}}𝐅_i^{ext}=A𝚷𝐧=A\left(P𝐧+𝝉𝐧\right)`$, respectively. Note that $`P𝐧`$ is the hydrostatic pressure force (pointing inwards the C$``$P cell), while the viscous contribution $`𝝉𝐧`$ depends on the local velocity gradient.

The balance of the energy flux requires some extra conditions. In Eq. (17) the convection, dissipation and conduction of energy are balanced if $`ϵ^{}=e`$, $`_i^N𝐅_i^{ext}𝐯_i=A𝚷𝐮𝐧`$ and $`𝐉_Q^{ext}𝐧=A𝐪𝐧`$ respectively .

Let us now consider in more detail how the scheme deals with momentum, energy and mass transfer from C to P.

### III.2 Momentum exchange

The condition $`𝐯^{}=𝐮`$ ensures the balance of momentum convection. If the mass flux points towards the P region ($`s>0`$), this condition is fulfilled by choosing the velocity of the inserted particles from a Maxwellian distribution according to the local temperature at the C$``$P cell, $`P(𝐯^{})=(1/2\pi mkT)^{3/2}\mathrm{exp}\left(m(𝐯𝐮)^2/2mkT\right)`$. Concerning particle removal $`(s<0)`$, we note that if the average velocity at the C$``$P cell is equal to the continuum velocity $`\overline{𝐯}=𝐮`$, then the average velocity of the subset of extracted particles would be precisely $`𝐯^{}=\overline{𝐯}=𝐮`$. Hence the condition of velocity continuity at C$``$P is needed to ensure the correct balance of momentum convection.

The change of momentum due to local stresses establishes the overall external force exerted on the particle region. Therefore one needs to determine how the overall external force is distributed on each individual particle. Flekkoy et al. Flek00  distributed the force according to a certain function $`g(x)`$ satisfying $`g(x_W)=\mathrm{}`$, $`g(x_O)=g^{}(x_O)=0`$. Normalisation leads to

$$𝐅_i^{ext}=\left(\frac{g(x_i)}{_i^{N_{CP}}g(x_i)}\right)A𝚷𝐧,$$
(18)
where $`x`$ runs perpendicular to the C$``$P interface and the applied force is made constant along each $`\mathrm{\Delta }t_C`$. As $`g(x)`$ tends to infinity as $`xx_W`$, the applied force diverges as one approaches the C$``$P interface; hence the density nearby $`x=x_W`$ is very small or zero. The function $`g(x)`$ is thus endowed with a two-fold purpose: it ensures a limiting extension to P (as the hydrostatic pressure force always points towards the P region, particles will never cross the C$``$P interface outwards) while also guaranteeing the existence of a small region where particles can be inserted with very low risk of overlapping.

Despite the benefits of the $`g(x)`$ function for distributing the externally imposed momentum, we decided to use $`g(x)=1`$ for all $`x`$ inside the C$``$P cell. The reasons for this choice will become clear when explaining the energy exchange, below. The first implication of $`g(x)=1`$ is that the external force is equally distributed among all the particles within the C$``$P cell. In other words, $`𝐅_i^{ext}`$ no longer depends on the particle label

$$𝐅_i^{ext}=𝐅^{ext}=\left(\frac{1}{N_{CP}}\right)A𝚷𝐧.$$
(19)
The second implication of $`g(x)=1`$ is that the particle density profile near the C$``$P interface no longer vanishes, so one needs an efficient way to resolve the problem of overlap on the insertion of new particles. This task is carried out by the usher algorithm, as explained below. Finally in order to ensure a finite extent of the particle region, if a particle ($`i`$) crosses outwards the C$``$P interface (in Fig. 1, $`x_i=x_W\delta `$, with $`\delta >0`$) with velocity $`𝐯_i`$ it is substituted by another one ($`j`$) with $`y_j=y_i;z_j=z_i`$, $`x_j=x_W+\delta `$ and with $`𝐯_j=𝐯_i`$. In this way, the overall momentum is strictly conserved before and after the particle exchange.

### III.3 Energy exchange

#### III.3.1 Advection

The balance of advected energy requires that $`ϵ^{}=e=u^2/2+3T/2+\varphi `$. Decomposing the particle energy into the kinetic and potential parts $`ϵ^{}=[v^{}]^2/2+\psi ^{}`$, one sees that, since the new particles are drawn from a Maxwell distribution, $`[v^{}]^2/2=u^2/2+3T/2`$. By contrast, the balance of potential energy $`\psi ^{}=\varphi `$ is much less straightforward to implement. This condition is fulfilled by the usher algorithm, described below.

#### III.3.2 Dissipation

One needs to satisfy the following balance of heat dissipation:

$$\underset{i}{\overset{N_{CP}}{}}𝐅_i^{ext}𝐯_i=A𝚷𝐮𝐧.$$
(20)
This condition does not generally hold if the external force is distributed according to an arbitrary $`g(x)`$. Indeed it is not even clear that a function $`g(x)`$ exists satisfying $`g(x_W)\mathrm{}`$ and enabling the heat dissipation balance in Eq. (20). In any case, such a function $`g`$ would depend on the particle’s velocity distribution and then the problem of finding $`g`$ would become a formidable task at each time step.

The advantage of using $`g(x)=1`$ now becomes clear. As long as $`𝐅_i^{ext}`$ does not depend on the particle label, one can greatly simplify the L.H.S. of Eq. (20) equation to obtain,

$$𝐅^{ext}\underset{i}{\overset{N_{CP}}{}}𝐯_i=N_{CP}𝐅^{ext}\overline{𝐯}=A𝚷𝐮𝐧.$$
(21)
The last equality follows from construction of the overall force $`N_{CP}𝐅^{ext}=A𝚷𝐧`$ and from the continuity of velocity $`\overline{𝐯}=𝐮`$.

#### III.3.3  Conduction

The condition $`𝐉_Q^{ext}𝐧=A𝐪𝐧`$ requires the establishment of a heat current along the C$``$P cell representing the conduction of energy. This may be implemented by various means; for instance, following the idea of Evans and coworkers Evans  one may include an extra force which pulls the “hotter” particles towards the direction of the heat flux and conserves the overall momentum. Alternatively, one may try to impose a Chapman-Enskog velocity distribution with the desired heat flux, at some region inside the C$``$P buffer. In this work we have made use of the phenomenological Fourier’s law, $`𝐪=\kappa _cT`$. A temperature gradient is imposed along each C$``$P cell by using a set of Nosé-Hoover thermostats (NHT) placed along the direction of the heat flux. The outer and inner thermostats are located a distance $`d`$ apart, and the temperature difference between both set to $`dT𝐧`$. Typically, at each C$``$P cell we have used a set of two or three NHT’s along a distance of $`3\sigma `$ or $`4\sigma `$. The values of the $`Q`$-parameter appearing in the NHT formulation Frenkel.book  were chosen small enough to minimise unphysical dynamics, i.e., we have chosen $`Q5`$. The main benefit of using the Nosé-Hoover formulation is the small distortion these thermostats introduce to the particle dynamics compared with other ways of implementing thermostatting Frenkel.book .

### III.4 Mass exchange: particle insertion

One important condition on the particle insertion, inherited from the balance of potential energy, is $`\psi ^{}=\varphi `$ (see Sec. III.3). To deal with this task, our strategy has been to place the new particles in positions where $`\psi ^{}\varphi `$. As $`\varphi `$ is, roughly speaking, the energy needed to insert a new particle, we shall insert particles with energies close to the local chemical potential. In any case, this implies that the insertions need to be made in very precise positions which depend on the configuration of the rest of the particles. To this end we have developed an algorithm (usher) that guides each new inserted particle to a position where the potential energy is equal to $`\varphi `$ (up to a pre-specified threshold). A brief explanation is given below (see Delgado-Buscalioni and Coveney INSERT  for further details). We first note that while the usher algorithm guides a new particle to a correct location, the rest of the particles remain frozen in position. usher essentially performs the following steps:

1. Place the new particle ($`i=N+1`$) at an initial position inside C$``$P, $`𝐫^{(0)}`$.
2. Evaluate $`𝐟_{N+1}=_{j=1}^N𝐟_{N+1,j}`$ and

$$\delta t_\sigma =\sqrt{\frac{2\delta r}{|𝐟_{N+1}|}}$$

Typically one can use $`\delta r\sigma `$.
3. Move the new particle according to the update rule

$$𝐫^{(n+1)}=𝐫^{(n)}+\frac{1}{2}𝐟_{N+1}^{(n)}\delta t^2$$

where $`\delta t=\mathrm{min}(\mathrm{\Delta }_t,\delta t_\sigma )`$, with $`\mathrm{\Delta }_t0.05`$ in reduced units
4. Evaluate the relative difference between the specific internal energy of the new particle, $`\psi _{N+1}^{}`$, and that prescribed by the continuum, $`\varphi `$: $`Err=|\psi _{N+1}^{}\varphi |/|\varphi |`$
5. The particle is correctly inserted if $`Err`$ is small enough (typically $`0.1`$).

Let us show how the usher algorithm easily overcomes the problem of possible initial overlap with pre-existing particles. An overlap leads to very large values of the interparticle force, $`f_{N+1}>>1`$, so in this case $`\delta t_\sigma <<1`$ and $`\delta t=\delta t_\sigma `$. But by construction, during the interval $`\delta t_\sigma `$ the new particle moves a distance of the order of the particle size $`\sigma `$ in the direction of minimum energy, just enough to avoid any initial overlap. Then, as the particle steadily moves towards a local minimum of energy, $`f_{N+1}`$ decreases, and $`\delta t_\sigma `$ increases until it becomes larger than $`\mathrm{\Delta }_t`$. Then, $`\delta t=\mathrm{\Delta }_t`$ is fixed. For liquid densities varying between $`\rho =0.50.8`$, the usher algorithm typically needs 15-90 iterations (single-force evaluations) to correctly place a new particle note  . By introducing the particles with $`\psi ^{}=\varphi (1\pm Err)`$, we found that, upon averaging over $`\mathrm{\Delta }t_C`$, the condition $`\psi ^{}=\varphi `$ holds within about $`2\%`$ (even using values of $`Err`$ as large as $`0.5`$).

Until now we have not mentioned any limitation on the sizes of the coarse mesh and time step $`\mathrm{\Delta }X`$ and $`\delta t_C`$. Comments on this topic are very scarce in the previous literature on hybrid methods for fluids. Moreover, since the local averages are made using these spatial and temporal windows, it is also appropriate to formulate any condition on $`\mathrm{\Delta }X`$ and $`\mathrm{\Delta }t_C`$ before presenting the results of the tests carried out for different flows.

## IV  Length and time scale prerequisites

#### IV.0.1 Arising from consideration of the microscopic description

Continuum fluid dynamics rest on the local equilibrium assumption. This means that, in order to define a local thermodynamic and hydrodynamic state characterising the coarse-grained variables at the C$``$P and P$``$C cells, the size of these cells needs to be greater than the mean free path $`\lambda `$ extracted from the particle dynamics. Moreover the local equilibrium within each cell should be attained on times scales smaller than $`\mathrm{\Delta }t_C`$. This means that $`\mathrm{\Delta }t_C`$ has to be larger than the collision time $`\tau _{col}`$. In summary, $`\mathrm{\Delta }X>\lambda `$ and $`\mathrm{\Delta }t_C>\tau _{col}`$. In the case of a LJ-fluid it is possible to use the hard-sphere approximation to make an order-of-magnitude estimate, $`\mathrm{\Delta }X>0.2\rho ^1`$ and $`\mathrm{\Delta }t>0.14\rho ^1T^{1/2}`$. These conditions become less restrictive at larger densities; as an example, for $`T=1`$, $`\rho =0.5`$ and a typical MD time-step $`\mathrm{\Delta }t_P10^3`$, local equilibrium would require $`\mathrm{\Delta }t_C/\mathrm{\Delta }t_P100`$ integration steps.

#### IV.0.2 Arising from consideration of the hydrodynamic description

Conditions on $`\mathrm{\Delta }X`$ and $`\mathrm{\Delta }t_C`$ are firstly imposed by the smallest characteristic length and time scales involved in the process under investigation (say $`2\pi /k_{\mathrm{max}}`$ and $`2\pi /\omega _{\mathrm{max}}`$ respectively). Practically, to correctly recover the smallest spatio-temporal flow pattern one needs at least 8 points per period, so $`\mathrm{\Delta }t_C\pi /4\omega _{\mathrm{max}}`$ and $`\mathrm{\Delta }X\pi /4k_{\mathrm{max}}`$. The numerical stability of the C-solver algorithm may also impose limitations. As mentioned in Sec. I, algorithms with explicit time discretisation are better suited for the C-solver of a hybrid scheme. A necessary condition for their numerical stability is $`\mathrm{C}=U\mathrm{\Delta }t_C/\mathrm{\Delta }X1/2`$, where C is the Courant number and $`U`$ the maximum characteristic flow velocity. The value of $`U`$ depends on the physical process one is dealing with but, to provide numbers in the present discussion, let us assume that we are dealing with low or moderate Reynolds numbers. Then, if the process is a diffusive one, $`U=\nu /\mathrm{\Delta }X`$; alternatively, if sound waves are relevant within the flow, $`U=c_s`$. In summary, the computational window for $`\mathrm{\Delta }t_C`$ should be, $`0.14\rho ^1T^{1/2}<\mathrm{\Delta }t_C\mathrm{\Delta }X/(2U)`$. Using the maximum grid-spacing allowed, $`\mathrm{\Delta }X=\pi /4k_{\mathrm{max}}`$, one obtains the computational windows for $`\mathrm{\Delta }t_C`$ shown in Fig. 2 versus $`\rho `$, for $`k_{\mathrm{max}}=\{0.1,0.2\}`$, $`U=\{\nu /\mathrm{\Delta }X,c_s\}`$ and $`T=2.5`$. As expected, a sound wave requires smaller time steps than a diffusive process. For large enough $`k_{\mathrm{max}}`$, the temporal and spatial computational window may be highly localised; e.g. for $`\rho =0.5`$ one should use $`0.2<\mathrm{\Delta }t_C<0.5`$ if waves with wavelength smaller than $`30\sigma `$ need to be captured by the coarse-grained description. As the density decreases these conditions become much more restrictive, until the acoustic time finally becomes smaller than the collision time (see Fig. 2). Also, in rarefied gases, if $`\rho <k_{max}/4`$, the mean free path becomes larger than the wavelength, but here we shall not be concerned with situations where the Navier-Stokes equations are not appropriate (see Garcia et al. for further discussion Gar00 ).

## V  Tests: Hydrodynamic modes

As already mentioned our hybrid scheme has been tested under stationary and unsteady flows. Typical stationary nonequilibrium states were considered, such as heat conduction profiles Ciccotti  and Couette profiles Flek00 ; Ciccotti . The microscopic reconstruction of these flows has been well studied in the literature that nothing new may added here. In passing, we note the transient times to achieve the steady state from the rest solution were found to be in agreement with the diffusive times $`L_x/\kappa `$ and $`L_x/\nu `$. The rest of the discussion will be focused on our choice of unsteady scenarios, which are described by the decay of transversal and longitudinal waves. These flows are now briefly presented, using standard hydrodynamics.

Consider a fluid at equilibrium characterised by homogeneous mass density, $`\rho ^e`$, specific energy $`e^e`$ and a vanishing mean velocity $`𝐮^e=0`$. Our procedure is to perturb this equilibrium state with different hydrodynamic fields $`\{\rho ^p,𝐮,e^p\}`$ (periodic in the $`x`$ direction, i.e. $`𝐤=k𝐢`$) and then make use of the C$``$P coupling scheme described in Sec. III to verify that, within the particle region, the subsequent evolution towards equilibrium is carried out in a hydrodynamically consistent way. If use is to be made of the linear hydrodynamic theory, the externally induced perturbations should be small enough to guarantee that the relaxation process is always governed by the linearised mass, momentum and energy equations (3)-(5) HanMcD . We defer further discussion of this point to Sec. VI below.

As is customary, to solve the linearised set of equations a Laplace-Fourier transform (LFT) is first performed HanMcD . The LFT of any perturbative variable, say $`\mathrm{\Phi }(𝐫,t)`$, will be denoted as

$`\mathrm{\Phi }(𝐤,t)`$ $``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{\Phi }(𝐫,t)\mathrm{exp}(i𝐤𝐫)𝑑𝐫,`$ (22)
$`\widehat{\mathrm{\Phi }}(𝐤,z)`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑z\mathrm{exp}(izt)\mathrm{\Phi }(𝐤,t)𝑑t.`$ (23)
The LFT of the linearised equation (3)-(5) leads to the following algebraic system for $`\widehat{𝚽}=(\widehat{\rho }^p,\widehat{T}^p,\widehat{j_x^p},\widehat{j_y^p},\widehat{j_z^p})`$ HanMcD ,

$$𝐌\widehat{𝚽}^T(𝐤,z)=𝚽^T(𝐤,0),$$
(24)
where the hydrodynamic matrix is,

$`𝐌`$ $`=`$ $`\left(\begin{array}{ccccc}iz& 0& ik& 0& 0\\ 0& iz+\kappa \gamma k^2& ik\frac{\gamma 1}{\rho ^e\alpha }& 0& 0\\ ikc_s^2/\gamma & ikD& iz+bk^2& 0& 0\\ 0& 0& 0& iz+\nu k^2& 0\\ 0& 0& 0& 0& iz+\nu k^2\end{array}\right).`$ (30)
We note that instead of using $`e^p`$, the energy equation is expressed in terms of temperature fluctuations. Also, for clarity it is better to write the solution $`\widehat{𝚽}`$ in terms of the $`t=0`$ Fourier-transformed perturbative heat density $`Q^p`$ and pressure $`P^p`$. These quantities are related to $`\rho ^p`$ and $`T^p`$ through the relations

$`e^p`$ $`=`$ $`c_vT^p+\left({\displaystyle \frac{e}{\rho }}\right)_{T^e}\rho ^p,`$ (31)
$`Q^p`$ $`=`$ $`\rho ^ec_v\left(T^p{\displaystyle \frac{\gamma 1}{\rho ^e\alpha }}\rho ^p\right),`$ (32)
$`P^p`$ $`=`$ $`{\displaystyle \frac{mc_s^2}{\gamma }}\left(\rho ^e\alpha T^p\rho ^p\right).`$ (33)
Here $`c_p`$ and $`c_v`$ are the specific heat at constant pressure and volume, respectively; $`\gamma =c_p/c_v`$, is the adiabatic coefficient; $`\kappa =\kappa _c/c_p\rho ^e`$ is the thermal diffusivity; $`\alpha =(\rho /T)_P/\rho ^e`$ the thermal expansion coefficient; and the kinematic longitudinal viscosity, $`b=(4\eta /3+\xi )/\rho ^e`$ is related to the sound attenuation coefficient, $`\mathrm{\Gamma }`$, through $`2\mathrm{\Gamma }=b+(\gamma 1)\kappa `$. Finally $`D=(P/T)_\rho `$ and the adiabatic speed of sound is $`c_s^2=\gamma \left(P\rho \right)_T`$.

Provided that, in the hydrodynamic limit, the wavelengths of the perturbations are much larger than the mean interparticle distance, it is sufficient to obtain the solution of Eq. (24) up to $`O(k^2)`$ HanMcD . Putting $`z\omega `$ into Eqs. (24) and (30), leads, after some algebra, to the following identities

$`\widehat{T}(𝐤,\omega )`$ $`=`$ $`{\displaystyle \frac{E_\kappa Q(𝐤,0)}{\rho ^ec_p}}+{\displaystyle \frac{\gamma 1}{\rho ^e\alpha }}E_{SR}P(𝐤,0)+{\displaystyle \frac{\gamma 1}{\rho ^e\alpha c_s}}E_{SI}j_x(𝐤,0),`$ (34)
$`\widehat{\rho }(𝐤,\omega )`$ $`=`$ $`{\displaystyle \frac{\alpha }{c_p}}E_\kappa Q(𝐤,0)+{\displaystyle \frac{1}{mc_s^2}}E_{SR}P(𝐤,0)+{\displaystyle \frac{1}{c_s}}E_{SI}j_x(𝐤,0),`$ (35)
$`\widehat{j}_x(𝐤,\omega )`$ $`=`$ $`{\displaystyle \frac{1}{mc_s}}E_{SI}P(𝐤,0)+E_{SR}j_x(𝐤,0),`$ (36)
$`\widehat{j}_y(𝐤,\omega )`$ $`=`$ $`E_\nu j_y(𝐤,0),`$ (37)
$`\widehat{j}_z(𝐤,\omega )`$ $`=`$ $`E_\nu j_z(𝐤,0),`$ (38)
and, using Eqs. (32) and (33),

$`\widehat{Q}(𝐤,\omega )`$ $`=`$ $`E_\kappa Q(𝐤,0),`$ (39)
$`\widehat{P}(𝐤,\omega )`$ $`=`$ $`E_{SR}P(𝐤,0)+mc_sE_{SI}j_x(𝐤,0),`$ (40)
where the following propagators have been introduced,

$`E_\kappa (k,\omega )`$ $`=`$ $`\mathrm{exp}(\kappa k^2t),`$ (41)
$`E_{SR}(k,\omega )`$ $`=`$ $`\mathrm{exp}(\mathrm{\Gamma }k^2t)\mathrm{cos}(c_skt),`$ (42)
$`E_{SI}(k,\omega )`$ $`=`$ $`i\mathrm{exp}(\mathrm{\Gamma }k^2t)\mathrm{sin}(c_skt),`$ (43)
$`E_\nu (k,\omega )`$ $`=`$ $`\mathrm{exp}(\nu k^2t).`$ (44)
The shear or transverse modes correspond to momentum perturbations along either $`y`$ and/or $`z`$ axis (i.e. perpendicular to the wavevector $`k𝐢`$); from Eqs. (37) and (38) it is clear that they are completely decoupled. The remaining hydrodynamic variables $`\{\rho ^p,T^p,j_x^p\}`$ are coupled and conform to three longitudinal modes which can be divided into two subgroups: two sound modes (involving longitudinal momentum and pressure at constant entropy) and one heat mode (involving heat diffusion). The fact that the heat density is an independent mode is evident from Eq. (39), where it is seen that it relaxes diffusively $`\mathrm{exp}(\kappa k^2t)`$. We note that $`Q^p=T^es^p`$ is essentially the fluctuation of the entropy density $`s^p`$ HanMcD . From Eqs. (36) and (40) it is easily shown that the two sound modes ($`\widehat{P}/c_s\pm \widehat{j}_x`$) decay like $`\mathrm{exp}(ikc_st\mathrm{\Gamma }k^2t)`$.

## VI  Results

### VI.1 Set-up and initial states

The coupling scheme was implemented and tested in the set-up shown in Fig. 3. The system is periodic along $`y`$ and $`z`$ directions and the gradients of the continuum variables are set along the $`x`$ direction. The particle region occupies a region of size $`L_x`$ around $`x=0`$ and of size $`L_y=L_z`$ along the periodic directions. The P region is divided into control cells of size $`\mathrm{\Delta }X`$ where in local averages are taken. The centres of the two C$``$P slabs (the outermost cells) are situated at $`x=\pm |L_x\mathrm{\Delta }X/2|`$. The deviation from the local equilibrium assumption was monitored in terms of the relative difference of the cell-averaged pressure and energy with respect the values given by the equation of state of Johnson et al. Ref EOS . Around a distance $`1.5\sigma `$ away from the C$``$P interface the typical maximum deviations were only about $`6\%`$.

The initial perturbative flow was prepared by first letting the P region relax until a vanishing and homogeneous mean flow was obtained. Then, during a small time interval ($`3\tau `$), the particle velocities were periodically changed according to a Maxwellian distribution with the desired velocity profile and local cell temperature. The resulting initial state was then analysed to extract the Fourier components of the whole set of flow variables ($`𝐯`$, $`\rho `$, $`T`$, $`e`$, P). For the sake of consistency these were extracted by Fourier transforms of the cell-averaged variables

$`\overline{\varphi }_{\mathrm{cos}}^{(n)}(t){\displaystyle \frac{c_n}{M_c}}{\displaystyle \underset{l}{\overset{M_c}{}}}\overline{\varphi }(X_l,t)\mathrm{cos}(k_nX_l),`$ (45)
$`\overline{\varphi }_{\mathrm{sin}}^{(n)}(t){\displaystyle \frac{c_n}{M_c}}{\displaystyle \underset{l}{\overset{M_c}{}}}\overline{\varphi }(X_l,t)\mathrm{sin}(k_nX_l),`$ (46)
where $`k_n=nk`$ ($`n𝒩)`$ and; $`c_n=1`$ for $`n=0`$, $`c_n=2`$ otherwise. In any case, it was checked that the Fourier transform of the microscopic variables $`\varphi ^{(n)}=c_n_i^N\varphi (x_i,t)\mathrm{exp}(ik_nx_i)/N`$ yields essentially the same output as Eqs. (45) and (46).

The initial Fourier transforms calculated from Eqs. (45) were injected into Eqs. (34)-(40) to obtain the time-evolution of the continuum variables. These, in turn, were used to calculate the fluxes imposed on the C$``$P cells over time. The transport coefficients used were those reported in the literature TRANSP ; EOS . As an interesting check, it was found that the coupling scheme failed significantly if the transport coefficients used in the C-region differed by more than about $`15\%`$ from those of the LJ fluid. In particular, the oscillation frequency of the averaged particle velocity differed with respect to that of C, while correlations decayed at a faster rate than those of the C-flow. Therefore, in cases where the constitutive relations are not known, this result means that it is first necessary to measure the transport coefficients from the particle dynamics (using any standard molecular technique), before applying the hybrid scheme, particularly if unsteady flows are to be studied.

The wavelengths of the initial perturbations were chosen to be much larger than the mean free path, i.e. $`2\pi \lambda /k<<1`$, in order to work within the hydrodynamic regime. In other words, the dependence of the transport coefficients on the wavenumber was negligible HanMcD . The amplitudes of the initial perturbation were chosen small enough to ensure that the subsequent relaxation process could be described by linear theory. In particular, if $`\overline{𝐯}^{(1)}(t)`$ is the maximum Fourier amplitude of the velocity, the typical values of the Reynolds number at $`t=0`$ \[Re$`=|\overline{𝐯}^{(1)}|\rho _e/(k\eta )`$\] were Re$`(0)3`$. As $`|\overline{𝐯}^{(1)}(t)|`$ decays exponentially, convection was present only in the first stages of the relaxing flow, but it was not strong enough to produce significant deviations from the linear theory (non-linear effects become dominant for Re$`>O(10)`$ TRITTON ). The maximum Mach number was less than $`0.2`$, and density fluctuations were around $`|\overline{\rho }^{(1)}|/\rho _e\mathrm{Ma}^20.05`$.

In another test the hybrid scheme was applied to a fluid in mechanical ($`𝐮=0`$) and thermodynamical equilibrium ($`\rho =0.5`$, $`T=3.5`$, $`e=2.7`$, $`P=3.2`$) during a longer simulation (50$`\tau `$) to check for any possible spurious drift in the overall momentum and energy (note that Eq. (15) ensures the mass conservation by construction). During this calculation, the total momentum inside the P region was conserved up to $`5\times 10^4`$ and the total energy fluctuated $`5\%`$ around its equilibrium value. The size of these fluctuations is consistent with the system size (which contained $`N=1600`$ particles and a specific heat of $`c_v=1.8`$). Note that the total energy of the system cannot be conserved because a part the system is connected to a thermostat and it also receives mechanical energy from C.

### VI.2 Transversal waves

In order to test the transfer of momentum flux along the direction perpendicular to the C$``$P interface, planar shear waves along $`x`$ were excited in a LJ-fluid with $`\rho =0.5`$, $`T=2.5`$ and $`\eta =0.75\pm 0.05`$ (the error bar comes from Ref.TRANSP ). These waves were created by imposing sinusoidal $`y`$-velocity profiles $`v_y(x)=v_{y,\mathrm{sin}}^{(1)}\mathrm{sin}(kx)`$.

Figure 4 shows results for perturbations with $`v_{y,\mathrm{sin}}^{(1)}=1.0`$ in a system with $`L_x=20\sigma `$, $`L_y=L_z=7\sigma `$ and $`M_c=10`$ cells ($`\mathrm{\Delta }X=2`$). The filled circles in Figure 4a correspond to the main Fourier component of the velocity $`\widehat{v}_{y,\mathrm{sin}}^{(1)}`$ for a calculation using wavenumber $`k=0.31`$. The relaxation process is indeed exponential, until the noise-to-signal ratio becomes large enough (around $`t>10`$) and, the observed decay rate $`(0.14\pm 0.1)\tau ^1`$ agrees perfectly with the theoretical value $`k^2\eta /\rho ^e`$. To explicitly appreciate the effect of the hybrid coupling we show (open circles in Fig. 4a) the outcome of a calculation in which the unsteady shear stress contribution to the external force $`F^{ext}`$ was set to zero, leaving just the contribution of the equilibrium hydrostatic pressure. As expected, if less momentum flux is provided towards P, the decay rate is appreciably slower ($`0.07\tau ^1`$) than the hydrodynamic one ($`0.14\tau ^1`$). Figure 4b shows the autocorrelation function (ACF) of $`\widehat{v}_{y,\mathrm{sin}}^{(1)}(t)`$ for another perturbation with an slightly different wavenumber, $`k=0.35`$. The good agreement with the theoretical decay \[$`\mathrm{exp}(0.17t)`$, in dashed line\] shows that the scheme is able to deal with small variations in the perturbation shape.

### VI.3 Longitudinal perturbations

Longitudinal waves transport mass, momentum and both mechanical and thermal energy; they are therefore perfectly suited for an overall test of the hybrid scheme behaviour. In our set-up the wavevector of these waves was perpendicular to the C$``$P interface and they were generated by imposing mean velocities along the $`x`$ direction as either pure cosinusoidal, sinusoidal profiles, or combinations of both profiles. As explained before, the particle velocities were extracted from Maxwell-Boltzmann distributions at the local mean velocities and at a constant temperature. The mean velocity profile induces pressure, density and energy fluctuations with a periodic pattern in the $`x`$ direction. The temporal behaviour of the main flow variables is now described.

#### VI.3.1 Mass

One of the main problems involved in the hybrid mass transfer at C$``$P is that although the continuum flux is a floating point number, one can only possibly exchange an integer number of particles. In order to adhere closely as possible to the prescribed continuum mass flux, the following procedure is followed during each interval $`t_C<tt_C+\mathrm{\Delta }t_C`$, $`t_C=m\mathrm{\Delta }t_C`$, $`m𝒩`$. First one evaluates the quantity

$$\xi (t_C)_{t_C}^{t_C+\mathrm{\Delta }t_C}s(t)𝑑t.$$
(47)
This floating point number, which represents the number of particles that should cross the C$``$P interface along $`t_C<tt_C+\mathrm{\Delta }t_C`$, is converted into an integer $`\delta N(t)`$ by the following construction:

$`\delta N(t)`$ $`=`$ $`\mathrm{NINT}\left[\xi (t_C)\right]+\delta \xi (tt_k)`$ (49)
$`\delta \xi (tt_k)=\mathrm{INT}\left[{\displaystyle _{t_k}^{t\mathrm{\Delta }t_P}}(\xi (t^{})\delta N(t^{}))𝑑t^{}\right]`$
where $`t_k`$ is such that $`|\delta \xi (tt_k)|1`$ and $`0=t_0<t_k<t_{k+1}`$. The deviation $`\delta \xi (tt_k)`$ assimilates the errors made through successive rounding off ($`\xi \mathrm{NINT}[\xi ]`$). When $`|\delta \xi (tt_k)|`$ becomes larger than one (at $`t=\left\{t_k\right\}`$) a particle is added to (or extracted from) $`\delta N`$ and the corrector $`\delta \xi `$ is then reset to zero. To minimise the effect on the remaining particles over each interval $`\mathrm{\Delta }t_C`$, the particle crossings are regularly separated in time, at a rate as close as possible to $`\delta N(t_C)/\mathrm{\Delta }t_C`$. As illustrated in Fig. 5, this kind of procedure enables us to follow rather closely the desired mass flux.

#### VI.3.2  Momentum

Figures 7, 8 and 9 show the time dependence of the Fourier components of the main hydrodynamic and thermodynamic variables. Dashed lines correspond to the theoretical trends obtained via the inverse Fourier transform of Eqs. (34)-(40).

For the reasons explained in Sec III.1, it is particularly important to ascertain that the averaged velocity at C$``$P correctly follows the desired continuum flow velocity. Figure 6 shows the instantaneous $`\overline{v_x}`$ and time averaged velocity $`\overline{v_x}`$ at both C$``$P cells. We obtained very good agreement, the deviation of $`\overline{v}_x`$ with respect the continuum value being less than about $`5\%`$ along most part of the damped oscillation (see Fig. 6). The importance of evaluating the continuum fluxes at precisely the C$``$P interface ($`x_W`$ in Fig. 1), is illustrated in Fig. 6 by comparison with the outcome of a calculation in which the fluxes are evaluated at $`x=x_0`$. For this calculation, the deviation with respect to the continuum prescription is clear and its magnitude agrees with the estimate made in Sec. III.1, $`30\%`$.

Figure 7a shows the time evolution of the Fourier amplitudes of the velocity in one of the calculations (a pure cosinusoidal perturbation with $`\overline{v}_{x,\mathrm{cos}}^{(1)}=0.6`$). The overall velocity $`\overline{v}_{x,\mathrm{cos}}^{(0)}`$ remains close to zero, confirming that the method conserves the initial total momentum. The $`\overline{v}_{x,\mathrm{sin}}^{(1)}`$ component (initially set to zero) has been also included to display the level of noise. The autocorrelation function (ACF) of $`\overline{v}_x^{(1)}=\overline{v}_{(x,\mathrm{cos})}^{(1)}+i\overline{v}_{(x,\mathrm{sin})}^{(1)}`$ is shown in Fig. 7b for two runs with different initial profiles (for more details see the caption of Fig. 7). These data were fitted to the hydrodynamic expression arising from Eq. (36) and in Fig. 7 these fits are shown with dashed lines. The best fit to the velocity was obtained with $`\mathrm{\Gamma }k^2=0.076`$ and $`c_sk=0.867`$, while for the ACF it yielded $`\mathrm{\Gamma }k^2=0.072`$ and $`c_sk=0.867`$. These values coincide, within the error bars, with those imposed by the C-flow (obtained upon insertion of the transport coefficients reported in the literature TRANSP ) $`0.071`$ and $`0.88`$, respectively.

#### VI.3.3 Thermodynamic variables

We start by observing that, unlike longitudinal momentum and pressure, the relaxation of density, temperature and energy perturbations includes not only an acoustic part, but also an entropic contribution proportional to $`\mathrm{exp}(\kappa k^2t)`$ \[see Eqs. (34)-(44)\]. However, as the initial perturbation considered was a mechanical one ($`T^{(n)}0`$), the entropic contribution is rather small. This can be seen in the theoretical expressions written in Fig. 6: the amplitude at $`t=0`$ of the entropic part of any heat-related variable is nearly six times smaller than its mechanical counterpart. This observation led us to a more careful study of the effect of heat conduction. As explained in Sec III.3.3, heat currents through each C$``$P cell were created by imposing several (typically two) Nosé-Hoover thermostats (NHT’s) placed close to the C$``$P interface over a distance of $`4\sigma `$. We found it very informative to study the effect of the number $`p`$ of NHT’s per C$``$P cell (notation: $`p`$-NHT<sub>CP</sub>) on the collective behaviour of the system. In some calculations, this number was reduced to merely 1-NHT<sub>CP</sub>, in such a way that only the local temperature prescribed by the continuum was imposed (and not the heat flux). The results of this comparison may be seen in Figs. 8 and 9. Calculations using 1-NHT<sub>CP</sub> yielded essentially the same evolution of the velocity and pressure as those with a larger number of NHT’s, although the pressure for the 1-NHT<sub>CP</sub> case showed a slight phase-lag, see Fig. 8d. This is not surprising as $`𝐯`$ and P are governed by acoustic terms and are independent of heat transfers. Using 1-NHT<sub>CP</sub>, deviations on $`\rho ^p`$, $`e^p`$ and $`T^p`$ with respect to the hydrodynamic trend are indeed appreciable, while results with 2-NHT<sub>CP</sub> adhere closely to the analytical curves (see Fig. 8). In any case, the calculation of the entropy production is the best way to highlight the completely different qualitative behaviour of 1-NHT<sub>CP</sub> with respect to two (or more) NHT<sub>CP</sub>.

#### VI.3.4 Entropy

The entropy fluctuations, or more precisely the perturbation of heat $`Q^p=T^es^p`$, was calculated using Eq. (32). We note that, in order to reproduce the steady diffusive heat decay in Eq. (35), an exact cancellation of the acoustic oscillations coming from $`\rho ^p`$ and $`T^p`$ is necessary \[see Eqs. (35)and (34)\]. As stated above, the contribution associated with heat diffusion in Eqs. (35) and (34) is rather small compared to the mechanical one. This observation indicates that small mechanically-driven fluctuations around the local thermodynamic equilibrium may become large enough to alter the purely exponential heat decay. In other words, Eqs. (32) and (39) provide a demanding test of the coupling scheme under the present flow. Figure 9 shows the main Fourier component of the heat perturbation $`Q^p`$ obtained for 1- and 2-NHT<sub>CP</sub>. The dashed lines correspond to the theoretical expectation. It is evident that the 1-NHT<sub>CP</sub> case does not obey the second law of thermodynamics at all. On the contrary, a rather good agreement with the theoretical trend is obtained when using at least 2-NHT<sub>CP</sub>. In Fig. 9, it is seen that the behaviour measured in the 2-NHT<sub>CP</sub> case exhibits fluctuations around the theoretical straight line. Typically, the largest excursions last around $`3\tau `$ from pure exponential decay. As previously stated, they may be due to the weakness of the entropy perturbation, but to confirm this statement we plan in the future to study some “heat-driven” flows (such as a heat pulse with initial zero mean velocity).

## VII  Conclusions

We have presented the core of a hybrid continuum-particle method for moderate-to-large fluids which couples mass, momentum and energy transfers between two regions, C and P, described respectively by continuum fluid dynamics and by discrete particle Newtonian dynamics. Both domains overlap within a coupling region divided into two sub-cells which account for the two-way exchange: C$``$P and P$``$C. While the procedure at the P$``$C cell is simply to average the particle-based (mass, momentum and energy) fluxes in order to supply open boundary conditions to the C domain, the operations at the C$``$P cell are much less straightforward as they need to reconstruct a large number of (particles’) degrees of freedom only from the knowledge of the three fluxes of conserved-quantities arising within C. The present work has been concerned with extending the C$``$P coupling to arbitrary rates of mass, momentum and energy transfer. To this end, the proposed method have been tested under unsteady flows which demand conformance to the whole set of conserved variable densities. In particular, we have considered the set of relaxing flows arising from hydrodynamics, namely longitudinal and transversal waves. We have followed the idea proposed by Flekkoy et al. Flek00 , in the sense that the scheme is explicitly based on direct flux exchange between the C and P regions. In order to deal with unsteady scenarios, we have shown that the fluxes injected into the particle region from the continuum region need to be measured exactly at the C$``$P interface and not at the nodes of the continuum lattice.

The implementation of flux exchanges requires the supply of energy currents to the particle system arising from the C domain due to advection, dissipation and conduction. To inject the correct amount of advected energy, the particle-averaged specific energy at the C$``$P cell needs to be equal to the continuum value. This can only be achieved if the new inserted particles are placed at positions where the (inter-particle) potential energy equals the C-specified internal energy per unit mass. This severe condition has been implemented by the usher algorithm, whose purpose is two-fold: to provide the correct mass transfer rate and to ensure the balance of energy advection. In the proposed scheme, the balance of energy dissipation arises naturally, provided that the cell-averaged velocity and the injected momentum flux equal their continuum counterparts. This is made possible by applying the external force according to a flat distribution, instead of a biased one as used by Flekkoy et al. Flek00  but as a result, the new particles have to be inserted within a non-vanishing density environment. This is sorted out by the usher in a very efficient way. Energy conduction has been implemented by using a set of Nosé-Hoover thermostats adjacent to the C$``$P interface, whose temperature and position are determined through the continuum local temperature gradient. Confirmation of the validity of this procedure is obtained from the correct rate of entropy production computed in our simulations of longitudinal waves. We showed that using only one thermostat per C$``$P cell (i.e. providing only the local value of $`T`$ but not the heat flux) leads to negative entropy production. Therefore, in the context of energy transfer, this result reinforces the central importance of coupling-through-fluxes proposed by Flekkoy et al Flek00 .

Enhancements to the present scheme are under investigation and merit some discussion here. The number density at C$``$P may be controlled by a feedback algorithm which preserves the momentum flux balance. Other kinds of implementations for the energy transport by conduction also deserve to be considered. Finally, we plan to implement the P$``$C coupling in conjunction with a finite volume CFD solver in 3D. In order to extend the coupling scheme to higher dimensions some additional complications will need to be adressed. First, a mass flux will be assigned to each cell within an array of neighbourings C$``$P cells. To adhere to mass continuity, particles will then have to be inserted within precisely defined finite regions and the insertion algorithm may have to pay an extra computational cost for this restriction of the search domain. We have checked, however, that the distance travelled by the usher algorithm from the initial trial position to the final insertion site is rather small INSERT  (typically less than 1$`\sigma `$ and less than 0.5$`\sigma `$ on average) so we do not expect any significant extra cost if the search for insertion sites is done within volumes larger than $`(2\sigma )^3`$. In higher dimensions one may also have to smooth to some extent the variations of the mean mechanical quantities imposed along the C$``$P region. To this end, it may be necessary to interpolate the external force along neighbouring C$``$P cells in such a way that the local momentum flux imposed at each C$``$P cell is still preserved. In the same way, although the Maxwell distribution, used here to choose the velocities of the new particles, proved to be sufficient to ensure momentum continuity for 1D coupling (i.e, with no neighbouring C$``$P cells), it may be convenient to use a Chapman-Enskog distribution in higher dimensions. This latter distribution enables the average velocity of the inserted particles to conform to the velocity gradient along neighbouring C$``$P cells. We hope to report our findings in these areas in future publications.

## VIII Acknowledgements

We gratefully acknowledge fruitful discussions with E. Flekkoy, P. Español, G. Ciccotti and R. Winkler and useful comments from B. Boghosian, A. Ladd and I. Paganobarraga. This work is supported by the European Union through a Marie Curie Fellowship (HPMF-CT-2001-01210) to RD-B.


text:
# 1 (a) Sketch of the system with quantum dots embedded between a top gate and back contact. By application of a voltage, one can load single electrons in the QD and change the Fermi energy in the wetting layer; 𝑙≪𝑑. (b) The band diagram of a quantum dot and wetting layer.

Kondo-excitons and Auger processes in self-assembled quantum dots

A. O. Govorov1,2, K. Karrai3, R. J. Warburton4, and A. V. Kalameitsev2
1Department of Physics and Astronomy, Ohio University, Athens, Ohio 45701, USA
2Institute of Semiconductor Physics, 630090 Novosibirsk, Russia
3Center for NanoScience and Sektion Physik, Ludwig-Maximilians-Universität, 80539 München, Germany
4Department of Physics, Heriot-Watt University, Edinburgh EH14 4AS, UK

ABSTRACT

We describe theoretically novel excitons in self-assembled quantum dots interacting with a two-dimensional (2D) electron gas in the wetting layer. In the presence of the Fermi sea, the optical lines become strongly voltage-dependent. If the electron spin is nonzero, the width of optical lines is given by $`k_BT_K`$, where $`T_K`$ is Kondo temperature. If the spin is zero, the exciton couples with the continuum due to Auger-like processes. This leads to anticrossings in a magnetic field. Such states can be called Kondo-Anderson excitons. Some of the described phenomena are observed in recent experiments.

INTRODUCTION

Many-body phenomena and excitons in quantum dots (QDs) attract presently much interest. The number of electrons in QDs is voltage tunable. This makes it possible to study single electron effects. One example of many-body phenomena is the Kondo effect induced by the non-zero spin. So far, the Kondo effect was studied mostly in relation to transport properties Kondo . In optics, Kondo-type effects were discussed with respect to nonlinear shake-up processes in nanostructures Kondo-optics .

Here we study theoretically novel exciton states in self-assembled QDs coupled with delocalized electron states. If the electron spin of the QD is non-zero, the resulting states can be called Kondo-excitons. The Kondo temperature of these excitons is given by the small dimensions of QDs and can be as high as $`10K`$. The Kondo effect manifests itself as peculiar, temperature-dependent optical lines. If the exciton spin is zero, a hybridization of the final state can occur due to Auger-like processes. Another effect is the transition between different exciton ground states on increasing the Fermi energy. In the transition regime, the optical lines become strongly voltage-dependent. Some of above-mentioned effects have been observed in optical spectra of single InAs QDs conference-paper .

Charged excitons $`X^n`$, observed in experiments conference-paper ; Nature ; All , contain $`n+1`$ electrons and one hole. In the exciton, the hole is optically excited; the electrons are supplied by tunneling from the back contact and by optical excitation (fig.1).

The voltage applied between the top and back contact, $`U_{gate}`$, makes it possible to control the number of electrons in a QD. If $`U_{gate}`$ is sufficiently high, electrons fill the wetting layer, a 2D quantum well. The Fermi energy of a 2D gas in the wetting layer is found from the conditions of equilibrium, $`ϵ_F=(a_0^{}/4d)\mathrm{\Delta }U_{gate}`$, where $`a_0^{}`$ is the Bohr radius and $`\mathrm{\Delta }U_{gate}`$ is the voltage measured from the point at which the wetting layer starts to load APL . In our model, an asymmetric QD has two bound, non-degenerate states with energies $`E_s^e`$ and $`E_p^e`$. The indexes $`s`$ and $`p`$ originate from electronic shells in a symmetric QD. The system is described by the Anderson Hamiltonian,

$`\widehat{H}=\widehat{H}_{sp}+\widehat{H}_{Coul}^{dot}+\widehat{H}_{tun},\widehat{H}_{tun}={\displaystyle \underset{\sigma }{}}V_k[c_{k,\sigma }^+a_{p,\sigma }+a_{p,\sigma }^+c_{k,\sigma }],`$ (1)
where $`\widehat{H}_{sp}`$ is the single-particle energy; $`\widehat{H}_{Coul}^{dot}`$ is the operator for intra-dot Coulomb interactions, and $`\widehat{H}_{tun}`$ is the operator for a tunnel hybridization between the QD and Fermi sea; $`c_{k,\sigma }`$ is the operator of a delocalized electron with a momentum $`k`$ and spin $`\sigma `$; $`ϵ_k`$ denotes the kinetic energy of delocalized electrons; $`a_{p,\sigma }^+`$ is related to the $`p`$-state of QD. Here we assume that the tunnel coupling occurs only between the upper $`p`$-state of the QD and the Fermi sea (fig. 1). In our simplified approach, the $`p`$-electron moves in a potential having a barrier at the edges of the QD. This barrier is formed by the short-range QD potential and Coulomb repulsion from the $`s`$-electrons.

We assume that $`\mathrm{\Delta }E>U_{Coul}>E_{hyb}`$, where $`\mathrm{\Delta }E=E_p^eE_s^e`$ and $`U_{Coul}`$ is the intra-dot Coulomb energy, and $`E_{hyb}`$ is the energy due to the hybridization between the QD and the Fermi sea. Thus, the Coulomb interaction in a QD can be included with perturbation theory.

The most intense lines observed in the PL spectra of QDs Nature ; All  are related to the transitions between the $`s`$-states in the conduction and valence bands. The PL spectrum is given by a Green function:

$`I(\omega )=Re{\displaystyle _0^{\mathrm{}}}dte^{i\omega t}<\widehat{V}_{opt}^+(t)\widehat{V}_{opt}(0)>,\widehat{V}_{opt}=V_{opt}(b_{s,\frac{3}{2}}a_{s,}+b_{s,\frac{3}{2}}a_{s,}+c.c.),`$ (2)
where $`\omega `$ is the frequency of the emitted photon and $`\widehat{V}_{opt}`$ is the electron-photon interaction, involving the $`s`$-states of electrons and heavy holes ($`J=3/2`$). An averaging in eq. 2 involves all initial states at a finite temperature $`T`$.

First we consider the limit $`V_k0`$ and $`T0`$. At small Fermi energies $`ϵ_F`$, the ground state of the system is the exciton $`X^2`$ (fig. 2a, upper part). With increasing $`ϵ_F`$, the exciton changes its charge and the ground state becomes $`X^3`$ (fig. 2a, lower part). The transition Fermi energy $`ϵ_{12}`$ is determined by the equation,

$`E_{X2}+ϵ_{12}=E_{X3},`$ (3)
where $`E_{X2}`$ and $`E_{X3}`$ are intra-dot exciton energies. In fig. 2b we show the PL spectrum as a function of $`ϵ_F`$ in the absence of hybridization ($`V_k=0`$). At $`ϵ_F=ϵ_{12}`$, the spectrum changes abruptly because of the ground-state transition. In the regime $`ϵ_F<ϵ_{12}`$, the spectrum has two lines, $`\mathrm{\Delta }E_{X2,s}`$ and $`\mathrm{\Delta }E_{X2,t}`$, related to the singlet and triplet final states Nature . The splitting between these lines is due to the exchange interaction in the exciton. To calculate the optical energies, we use a model of an anisotropic harmonic oscillator. The oscillator energies are taken as follows: $`\mathrm{}\mathrm{\Omega }_x^e=20meV`$, $`\mathrm{}\mathrm{\Omega }_y^e=25meV`$, and $`\mathrm{\Omega }_{x,y}^h=\mathrm{\Omega }_{x,y}^e/2`$, where the upper indexes denote electrons (holes).

KONDO EXCITON $`X^2`$ WITH $`S_{dot}=1/2`$

Now we focus on the interaction between the QD and 2D electrons, assuming $`V_k0`$ and $`T=0`$. In the exciton $`X^2`$, the intra-dot electron spin is $`S_{dot}=1/2`$. So, the exciton wave function is mixed with the Fermi sea, leading to a Kondo state with $`S_e=0`$, where $`S_e`$ is the total electron spin (fig. 2a). A trial function for the Kondo exciton can be written as book :

$`\mathrm{\Psi }_{initial}=[A_0|\varphi _0>+{\displaystyle \underset{k>k_F}{}}A_k|\varphi _k>]|,0>_h,`$ (4)
$`|\varphi _0>=|,;,;\mathrm{\Omega }>_e,|\varphi _k>=(\widehat{c}_{k,}^+|,;,0;\mathrm{\Omega }>_e\widehat{c}_{k,}^+|,;0,;\mathrm{\Omega }>_e)/\sqrt{2}`$
(fig. 2a). The first four indexes in the electron function $`|\gamma _s,\gamma _s;\gamma _p,\gamma _p;\mathrm{\Omega }>_e`$ represent the occupations of the QD; $`\mathrm{\Omega }`$ denotes all states of the Fermi sea with $`k<k_F`$, where $`k_F`$ is the Fermi momentum. The function $`|,0>_h`$ describes the hole with $`J=+3/2`$. Note that the electronic part of the function (4) is a singlet state with $`S_e=0`$. By using the Hamiltonian (1), we find the unknown coefficients in eq. 4 and the ground state energy $`E_i=E_{X2}+ϵ_F\delta `$. The lowering of energy $`\delta `$ plays the role of a Kondo temperature, $`T_K`$. If $`E_p^e+U_pϵ_F>\delta `$, we obtain

$`k_BT_K=\delta =(Dϵ_F)e^{\frac{\pi (E_p^e+U_pϵ_F)}{2\mathrm{\Delta }}},`$ (5)
book . Here $`\mathrm{\Delta }=\pi V_0^2\rho _0`$ is the broadening and $`\rho _0`$ is the 2D density of states (DOS). $`U_p`$ is the Coulomb-blockade energy for the $`p`$-state; $`D`$ is the cut-off parameter for the tunnel matrix element $`V_k`$: $`V_k=V_0`$ for $`ϵ<D`$ and $`0`$ elsewhere D . For example, $`T_K714K`$ for $`ϵ_F=23.6meV`$, $`\mathrm{\Delta }=1meV`$, and $`D30meV`$.

The PL spectrum can be calculated by using eqs. 2 and 4. In the regime $`ϵ_F<ϵ_{12}`$, the PL spectrum contains three lines. The main lines are singlet and triplet states of $`X^2`$. The third structure is a shake-up line due to $`X^3`$. The most intense, triplet line has an asymmetric shape govorov :

$`I_{X2,t}={\displaystyle \frac{3}{4}}V_{opt}^2{\displaystyle \frac{2\mathrm{\Delta }A_0^2}{\pi }}{\displaystyle _0^{Dϵ_F}}𝑑ϵ{\displaystyle \frac{1}{(ϵ+k_BT_K)^2}}Re{\displaystyle \frac{i}{\mathrm{}\omega \mathrm{\Delta }E_{X2,t}^{}+ϵi\gamma _t}},`$ (6)
where $`\mathrm{\Delta }E_{X2,t}^{}=\mathrm{\Delta }E_{X2,t}k_BT_K`$ is the renormalized exciton energy, $`\mathrm{\Delta }E_{X2,t}`$ is the optical energy in the absence of the 2D gas, and $`\gamma _t`$ describes the phonon-induced broadening of the final state. Typically, the broadening of the triplet intra-dot state $`\gamma _t`$ is small, since the energy relaxation in this case involves a spin-flip process.

If $`KT_K\gamma _t`$, the triplet line is asymmetric and its width is equal to $`KT_K`$ (fig. 3). The shape and width of the line originate from the peak of the spectral density of states (SDOS) at the Fermi level in the initial Kondo state given by eq. 4. The line shape and peak position in fig. 4 depend on the Fermi energy and, therefore, are voltage-dependent. Besides, we think that the line width is strongly temperature-dependent, because the Kondo peak in the SDOS vanishes as the temperature increases book .

EXCITON $`X^3`$ WITH $`S_{dot}=0`$: AUGER PROCESSES

In the case $`ϵ_F>ϵ_{12}`$, the main peak in the spectrum comes from the exciton $`X^3`$ where there are two completely filled electron levels. Its intra-dot electron spin is zero. Thus, the Kondo effect does not occur in the initial state. However, after photon emission, the final state of the exciton $`X^3`$ is excited and can interact with the continuum of states due to the Auger-like process induced by the intra-dot Coulomb interaction (fig. 4a). In such a process, the final state is a linear combination of the wave functions $`|f_0>=|,0;,;\mathrm{\Omega }>`$ and $`|f_{1,k}>=\widehat{c}_k^+|,;0,0;\mathrm{\Omega }>`$, where $`k>k_F`$ eh . The energies of the above states can be written as $`ϵ_0=E_0^{dot}`$ and $`ϵ_{1,k}=E_1^{dot}+ϵ_k`$, where $`E_0^{dot}`$ and $`E_1^{dot}`$ are intra-dot energies and $`ϵ_k`$ is the kinetic energy. The parameter $`\delta ϵ_a=E_0^{dot}E_1^{dot}`$ determines the excess kinetic energy, when one of the QD electrons is excited into the continuum (fig. 4a).

By using non-crossing diagrams (fig. 4b), we write the Green function (2) in the form

$`I(\omega )=V_{opt}^2Re[{\displaystyle \frac{i}{\stackrel{~}{\omega }i0\mathrm{\Sigma }_1}}],\mathrm{\Sigma }_1={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{W_ϵ^2\rho (ϵ)[1f_ϵ]dϵ}{\stackrel{~}{\omega }\delta ϵ_a+ϵi0_0^{\mathrm{}}\frac{\rho (ϵ^{})W_ϵ^{}^2f_ϵ^{}dϵ^{}}{\stackrel{~}{\omega }+ϵϵ^{}i0}}},`$ (7)
where $`\mathrm{\Delta }E_{X3}`$ is the optical-line energy, $`\stackrel{~}{\omega }=\omega \mathrm{\Delta }E_{X3}`$, $`\rho `$ is the DOS and $`f_ϵ`$ is the Fermi-distribution function; the Auger matrix element $`W_{ϵ_k}=<f_0|U_{Coul}|f_{1,k}>`$, where $`U_{Coul}`$ is the Coulomb interaction. For the parameter $`W_ϵ`$, we assume: $`W_ϵ=W_0`$ if $`ϵ<D^{}`$ and $`0`$ elsewhere. In our approach, the final states in the PL process are described by the non-interacting Anderson model, in which the Auger matrix element plays the role of tunnel element $`V_k`$.

(i) Interaction with the continuum: empty wetting layer
We now assume that the QD is sufficiently deep and the wetting layer is empty; $`U_{dot}=79meV`$. In this case, we can put $`f(ϵ)=0`$. In the absence of a magnetic field, the PL peak (eq. 7) is close to a Lorentzian. The broadening is given by $`\mathrm{\Delta }_a=\pi W_0^2\rho _0`$. The estimated Auger-broadening $`\mathrm{\Delta }_a`$ can be as high as 2 $`meV`$ govorov ; Bastard . Note that this mechanism of broadening is found in experimental spectra of $`X^3`$ conference-paper . In the magnetic field $`B`$, the DOS becomes strongly modulated, due to formation of Landau levels,

$`\rho (ϵ)=\rho _0{\displaystyle \frac{\mathrm{}\omega _c}{\sqrt{\pi }\gamma _{LL}}}{\displaystyle \underset{n}{}}e^{\frac{(ϵϵ_n)^2}{\gamma _{LL}^2}},`$ (8)
where $`\gamma _{LL}`$ is the Landau level broadening, $`ϵ_n=(n+1/2)\mathrm{}\omega _c`$ is the Landau level energy, and $`\omega _c`$ is the cyclotron frequency. In this case, the optical energies demonstrate anticrossings between the terms: $`\stackrel{~}{\omega }=0`$ and $`\stackrel{~}{\omega }\delta ϵ_a+(n+1/2)\mathrm{}\omega _c=0`$ (fig. 4c). The anticrossing strength increases with the magnetic field as $`\sqrt{\mathrm{\Delta }_a\mathrm{}\omega _c}`$. This is due an increase of the density of states.

(ii) Spectrum of $`X^3`$ in the presence of a 2d gas
If $`T=0`$, $`B=0`$ and $`\delta ϵ_aϵ_F>\mathrm{\Delta }_a`$, the line is again close to a Lorentzian with a high-energy cut-off at the frequency $`\stackrel{~}{\omega }=\delta ϵ_aϵ_F`$. If $`B0`$, the PL spectrum in the presence of 2D gas is more complex due to electron-hole excitations in the Fermi sea.

TRANSITION FROM $`X^2`$ TO $`X^3`$

When $`ϵ_Fϵ_{12}`$, the ground state $`X^2`$ changes to the state $`X^3`$. In the transition regime, the exciton strongly couples with the Fermi sea. Hence, we can expect that the optical energies should strongly depend on $`ϵ_F`$. In fig. 2b, we show this effect by using a simplified, zero bandwidth model govorov ; model . We can see that the energies $`\mathrm{\Delta }E_{X2,t(s)}`$ decrease with $`ϵ_F`$ and, thus, are voltage-dependent. Such a behavior has been observed in recent experiments conference-paper .

CONCLUSIONS

Novel exciton states can appear in QDs interacting with the Fermi gas. The specific feature of these states is a “cloud” of electrons around the charged exciton. The radius of the Kondo-Anderson excitons is given by the size of the electron cloud and can be greater than the size of QD. We describe several striking manifestations of such excitons in PL spectra. Some of described phenomena are seen in experimental spectra of InAs QDs conference-paper .


text:
# Dynamical symmetry breaking as the origin of the zero-𝑑𝑐-resistance state in an 𝑎𝑐-driven system.

## Abstract

Under a strong $`ac`$ drive the zero-frequency linear response dissipative resistivity $`\rho _d(j=0)`$ of a homogeneous state is allowed to become negative. We show that such a state is absolutely unstable. The only time-independent state of a system with a $`\rho _d(j=0)<0`$ is characterized by a current which almost everywhere has a magnitude $`j_0`$ fixed by the condition that the nonlinear dissipative resistivity $`\rho _d(j_0^2)=0`$. As a result, the dissipative component of the $`dc`$ electric field vanishes. The total current may be varied by rearranging the current pattern appropriately with the dissipative component of the $`dc`$-electric field remaining zero. This result, together with the calculation of Durst *et. al.*, indicating the existence of regimes of applied $`ac`$ microwave field and $`dc`$ magnetic field where $`\rho _d(j=0)<0`$, explains the zero-resistance state observed by Mani *et al.* and Zudov *et al.*.

Recently, two experimental groups Mani ; Zudov  reported observations of a novel zero-resistance state in two-dimensional electron systems subjected to a $`dc`$ magnetic field and to strong microwave radiation. When no microwave power is applied, Refs Mani ; Zudov  observe a longitudinal resistivity only weakly dependent on magnetic field, at least for the relatively small fields (filling factor $`\nu >10`$) applied in the experiment. However, when a high level of microwave power was applied the resistance developed a strong oscillatory dependence on the applied magnetic field, with oscillation period controlled only by the ratio of the microwave frequency $`\omega `$ to the cyclotron frequency $`\omega _c`$. At low $`T`$ and in certain field ranges, the dissipative resistance was observed to vanish within the experimental accuracy.

An important step towards the understanding of these observations was taken in Ref. Yale , which presented a calculation of the effect of microwave radiation on the $`dc`$ linear response conductivity of a two dimensional electron gas. A crucial result of Ref. Yale , see also Ref. Ryzhii  for earlier treatment and Ref. Vavilov  for detailed analysis, is the existence of the regimes of magnetic field and applied microwave power for which the longitudinal linear response conductivity is negative,

$$\sigma _{xx}<0.$$
(1)
However, in the literature so far a precise connection between a negative linear response conductivity and the experimental observations has not been presented.

In this Letter we show that Eq. (1) by itself suffices to explain the *zero-$`dc`$-resistance* state observed in Refs. Mani ; Zudov , independent of the details of the microscopic mechanism which gives rise to Eq. (1). The essence of our result is that a negative linear response conductance implies that the zero current state is intrinsically unstable: the system spontaneously develops a non-vanishing local current density, which almost everywhere has a specific magnitude $`j_0`$ determined by the condition that the component of electric field parallel to the local current vanishes. The existence of this instability is shown, under reasonable assumptions, to imply the observed zero resistance state.

We consider $`dc`$ transport in a two-dimensional electron gas exposed to a static magnetic field and to an $`ac`$ electric field. We assume that the local $`dc`$ electric field $`𝑬`$ is related to the local $`dc`$ -current density $`𝒋`$ via

$$𝑬=𝒋\rho _d\left(𝒋^2\right)+\left[𝒋\times 𝒛\right]\rho _H,$$
(2)
where $`𝒛`$ is the unit vector normal to the plane of the system. The crucial quantity in Eq. (2) is the longitudinal (dissipative) resistivity $`\rho _d\left(𝒋^2\right)`$ whose dependence on current determines the physics we consider. The form of $`\rho _d(j^2)`$ is determined by parameters such as the applied magnetic field $`B_{app}`$ and the frequency $`\omega `$ and power $`𝒫_{ac}`$ of the $`ac`$-field, which we do not explicitly write. Also, to simplify the discussion we do not consider nonlinear effects in the Hall resistivity $`\rho _H`$. This is not crucial for the zero resistance state; effects of including it in the theory are discussed briefly at the end of the paper.

We assume that $`\rho _d\left(𝒋^2\right)`$ is a real, continuous function of $`j^2`$ and that (as found, e.g. in the calculations of Ref Yale ) a range of $`B_{app}`$, $`\omega `$ and $`𝒫_{ac}`$ exists for which a spatially homogeneous zero-current state is characterized by the negative dissipative resistivity

$$\rho _d\left(𝒋^2=0\right)<0.$$
(3)
However, at sufficiently large values of the $`dc`$ current $`\rho _d(j^2)`$ must revert to its dark ($`𝒫_{ac}=0`$) value because in this limit the microwave radiation will be a small perturbation on the steady state electron distribution function. Continuity implies that there is a value $`j=j_0`$ at which

$$\rho _d(j_0^2)=0.$$
(4)
We take $`\rho _d(j^2)`$ to have the form shown in the inset of Fig. 1. The main panel of Fig. 1 shows the current-voltage characteristic following from the assumed form of $`\rho _d(j^2).`$ Such a dependence was obtained analytically in Ref. Vavilov .

A negative dissipative resistivity is allowed under non-equilibrium conditions, if the system is continuously supplied with energy. In the situation considered here energy conservation requires only that $`j^2\rho _d\left(𝒋^2\right)+𝒫_{ac}>0`$. However, a negative resistivity may render the system unstable. Specifically, we now show that in a system described by Eq. (3) with resistivity curve as shown in the inset in Fig. 1.

(i) A homogeneous, time independent state characterized by a current $`j`$ of magnitude  less than the critical value $`j_0`$ defined in Eq. (4) is unstable with respect to inhomogeneous current fluctuations.

(ii) The only possible time independent state is one in which the current $`𝒋`$ has magnitude $`j_0`$ everywhere except at isolated singular points (vortex cores) or lines (domain walls), implying vanishing dissipative electric field, $`𝒋𝑬=0`$.

An immediate consequence of (ii) is that by adjusting the details of the current pattern, any net $`dc`$ current less than a threshold value (which we discuss below) can be sustained at vanishing electric field, so that *any microscopic mechanism of non-equilibrium drive resulting in $`\rho _d(𝐣^2=0)<0`$ leads to the observed Mani ; Zudov  zero dissipative differential resistance*:

$$\frac{dV_x}{dI_{dc}}=0.$$
(5)
(Here $`I_{dc}`$ is a sufficiently weak applied current.) If too large a current is imposed, the current structure will collapse and a non-vanishing resistance will be observed. We emphasize, however, that Eq. (5) is obtained on the assumption that the system is in a steady state. Any current pattern consistent with a boundary condition of small net current implies the existence of singularities (domain walls or vortices) in the current distribution; finite density of these objects may lead to a small dissipative resistivity.

Let us pause to discuss the relation of our arguments with phenomena discussed in the literature. The instability of systemsGunn  with absolute negative conductivity is known since the work of Zakharov.Zakharov  The important new feature of the instability and the domain structure of the present paper is that it occurs at large Hall angle; as the result the domains for the current coincide with the domains of the electric field directed perpendicularly to the current. We also would like to notice certain similarity with the model of photoinduced domains proposed by D’yakonovDyakonov  as an explanation of the experiments on Ruby crystals under the intense laser irradiationRuby .

We now present our specific arguments. We begin by considering the fluctuations $`\delta j`$ about a time-independent homogeneous state of current $`j_i`$. Taking the time derivative of Eq. (2) and using the continuity equation,

$$\frac{n}{t}+\mathbf{}𝒋=0,$$
(6)
and the Poisson equation,

$$𝑬=\mathbf{}\widehat{U}n.$$
(7)
we obtain

$$\mathbf{}\left(\widehat{U}\mathbf{}𝒋\right)=\frac{}{t}\left[𝒋\rho _d\left(𝒋^2\right)+\left[𝒋\times 𝒛\right]\rho _H\right].$$
(8)
Here $`n`$ is the electron charge density and $`\widehat{U}`$ is a nonlocal interaction operator which can be expressed in terms of the Green function of the Laplace equation with appropriate boundary conditions. The crucial point for us is that $`\widehat{U}`$ has non-negative eigenvalues. (We assume that the screening radius is equal to zero, and neglect the difference between the electric and electrochemical potentials. This approximation does not alter our conclusions.)

Writing $`𝒋(r,t)=𝒋_i+\delta 𝒋(r,t)`$, linearizing in $`\delta 𝒋`$, and taking the divergence of both sides of Eq. (8), we find

$$\frac{\mathbf{}\delta 𝒋}{t}=(\mathbf{}(\stackrel{~}{\widehat{\rho }}_d+\widehat{\rho }_H)^1\mathbf{}\widehat{U})\mathbf{}\delta 𝒋,$$
(9)
with $`\widehat{\rho }_H=\rho _H\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]`$ the usual Hall resistivity tensor,

$$\stackrel{~}{\widehat{\rho }}_d=\rho _d\mathrm{𝟏}+\alpha _j𝒋_i𝒋_i,$$
(10)
and

$$\alpha _j=2\frac{d\rho _d(𝒋^2)}{d𝒋^2}|_{𝒋^2=𝒋_i^2}.$$
(11)
The Coulomb interaction operator $`\widehat{U}`$ is positive definite, so the stability is determined by the sign of the operator $`\mathbf{}(\stackrel{~}{\widehat{\rho }}_d+\widehat{\rho }_H)^1\mathbf{}`$ in front of it. Doing Fourier transform of this operator we see that in order for any solution of Eq. (9) not to grow with time the conditions

$`\rho _d(𝒋^2)0,`$ (12a)
$`\rho _d(𝒋^2)+\alpha _j𝒋^20,`$ (12b)
must hold.

We therefore conclude that if at least one of $`\rho _d`$ or $`\rho _d+\alpha _j`$, is negative, i.e. if $`j_i<j_0`$ a homogeneous state of uniform current is unstable. However, we may also show that any state with local current density larger than $`j_0`$ but net current density smaller than $`j_0`$ is necessarily time-dependent. In this case the condition $`\mathbf{}𝒋=0`$ requires the presense of circulating currents. The integral $`𝒥=_C𝑑𝐥𝑬`$ along the current flow lines must vanish because $`\mathbf{}\times 𝑬=0`$. On the other hand, from Eq. (2) and $`\mathbf{}𝒋=0`$, we find $`𝒥=_C𝑑𝐥𝒋\rho _d\left(𝒋^2\right)`$. By construction of the contour $`d𝐥𝒋>0`$. Therefore $`𝒥=0`$ can be satisfied together with the stability condition (12a) only for $`\rho _d(𝒋^2)=0`$, i.e. for $`j=j_0`$.

Finally, we examine the stability of general states with $`|\stackrel{}{j}(r)|=j_0`$. In this case $`\rho _d=0`$ but $`\alpha >0`$; substitution into Eq, (9) and use of Eq. (6) leads to

$$\left\{\frac{}{t}+\left(\mathbf{}\left[𝒋_0\times 𝒛\right]\right)\left(\left[𝒛\times 𝒋_0\right]\mathbf{}\right)\frac{\alpha \widehat{U}}{\rho _H^2}\right\}\delta n=0,$$
(13)
$`\mathbf{}\delta 𝒋=_t\delta n`$. Operator $`\widehat{U}`$ is positive definite, while the operator $`\left(\mathbf{}\left[𝒋_0\times 𝒛\right]\right)\left(\left[𝒛\times 𝒋_0\right]\mathbf{}\right)`$ is Hermitian and is non-negative because it can be presented as $`AA^{}`$. Therefore, the state (4) is not unstable, except possibly at singular points. The investigation of the stability of the current pattern in the vicinity of the singular point requires going beyond the local equations (13) and (2) and will not be done in the present paper.

Moreover, one can see from Eq. (13) that all the perturbations decay in time exponentially with the exception of those for which $`\left[𝒋_0\times 𝒛\right]\mathbf{}\widehat{U}\delta n=0`$, i.e. with the electric field directed along $`𝒋_0`$. The physical meaning of these zero-modes is all the perturbations of the current which keep $`𝒋^2=𝒋_0^2`$ and $`\mathbf{}𝒋=0`$ (most trivial example of such perturbation is a homogeneous rotation of vector $`𝒋_0`$). These perturbations have zero eigenvalue and, analogously to Goldstone modes, are a straightforward consequence of the symmetry breaking induced by the applied nonlinear drive.

We now consider the physical consequences of our results. We found that a non-equilibrium system which has a negative linear response resistivity, is unstable to the formation of a state of non-vanishing local current. Almost everywhere in the sample the current has the magnitude $`j_0`$ at which the dissipative resistivity vanishes, but the direction must vary so that the net current is consistent with boundary conditions. The current distribution must contain singular regions, of negligible volume fraction, at which the current takes values different from $`j_0`$. The arguments relating to time-independent states given above may be viewed as showing that it is impossible to construct a time-independent singularity structure for distributions involving currents of magnitude greater than $`j_0`$, whereas it is possible if in almost all of the sample the current has magnitude $`j_0`$. Just as in the theory of superconductivity the detailed nature and structure of the singular regions (domain walls, vortex cores or other structures) presumably depends both on boundary conditions and on short length scale physics. The question cannot be analyzed within the quasicontinuum/local response function approach used in this paper, and is an important topic for future investigation.

For concreteness of further discussion we will consider the obvious choice of singularity shown in Fig. 2, namely a linear domain wall, separating two domains in which current flows parallel and antiparallel to the domain wall. We believe that a structure involving vortices would lead to essentially identical physics. In the presence of a magnetic field, consideration of the Hall component of the current reveals the existence of a discontinuity in the component of the electric field perpendicular to the boundary. If $`\widehat{𝐧}`$ is the vector perpendicular to the wall, and $`\mathrm{\Delta }𝒋=2j_0`$ is the discontinuity in current across the wall (assumed parallel to the wall direction) then the singularity in the electric field is

$$\mathrm{\Delta }E=2𝐧\rho _Hj_0.$$
(14)
This discontinuity requires a charge accumulation which, in a two dimensional situation, is non-local ($`l_0`$ is a cutoff set by the microscopic structure of the domain)

$$n(𝐫)\rho _Hj_0\mathrm{ln}\left(\frac{\left|𝐫\widehat{𝐧}\right|}{l_0}\right).$$
(15)
This charge distribution may be detectable by local probes. The other possible way to detect dynamical symmetry breaking is to measure spontaneous voltages arising inside the domain, voltage $`V_2`$ in Fig. 2.

Figure 2 presents a very natural (albeit probably oversimplified) picture of the experimental situation studied in Refs Mani ; Zudov . In these experiments the current in one direction (say, ’x’) was fixed by current leads to some value $`I`$, and the current in the transverse direction was set to zero. The longitudinal (x) and transverse (y) voltages were measured. Referring to Fig. 2, we see that any value of net current $`I`$ corresponding to a current density much less than $`j_0`$ can be obtained simply by adjusting the height of the domain wall: if $`d`$ is the position of the domain wall relative to the center of the device, then $`I=2dj_0`$ with $`V_x=0`$. Similarly, the total Hall voltage is the sum of a positive voltage in the upper half of the sample and a negative voltage in the lower half, leading to $`V_y=\rho _H\left[j_0\left(\frac{L_y}{2}d\right)j_0\left(\frac{L_y}{2}+d\right)\right]=\rho _HI`$, which will equal the dark (no microwave) result if $`\rho _H`$ is not much affected by the $`ac`$ field. Notice, that for the Corbino disc geometry the applied voltage, (corresponding to $`V_y`$ of Fig. 2) can be also accomodated by the shift of the domain wall without generation of the dissipative current, resulting in the zero-conductance-state Zudov4 .

The equations analyzed in this paper predict threshold behavior in $`I`$ at low temperature $`T:`$ $`V_x`$ is strictly zero for weak applied currents but if the applied current is large enough that the current density becomes of order $`j_0`$ then the domain wall is swept out of the system and a dissipative state corresponding to current densities greater than $`j_0`$ in some parts of the sample will result. Similarly, our equations predict a critical temperature: at very high temperature, the linear response conductivity will be positive even at non-zero (but fixed) microwave power. As $`T`$ is lowered, $`\sigma _{xx}`$ will decrease and at some temperature pass through $`0`$, upon which dissipationless behavior will result. The sharp thresholds in $`I`$ and $`T`$, which are in aparent contradition with Refs. Mani ; Zudov , may be artifacts of the simple treatment given here, which assumed a static singularity structure and zero screening radius.

We did not consider the dependence of the Hall conductivity on the applied current. It is easy to see that this dependence does not change the condition (4) for the circulating currents in the state because the Hall coefficient does not cause dissipation. Singular dependence of $`j_0`$ on the magnetic field will cause the singular features in the magnetic field dependence of the Hall resistivity near the zero-resistance state, The shape and the value of these singularities, however, have nothing to do with the quantized plateaux in the Quantum Hall effect.

Finally, we note that we have assumed an isotropic $`\rho _d`$. In fact, the presence of an $`ac`$ drive will lead to a quadrupolar anisotropy, see Ref. Vavilov  for microscopic derivation, which for the sake of notational clarity we did not write but which can easily be included if desired. This anisotropy will presumably affect the orientations of domain walls, suggesting that it would be interesting to look for differences in threshold behavior for $`dc`$ current parallel or perpendicular to $`ac`$ current.

To summarize, we have shown that the remarkable zero resistance state found by Refs. Mani ; Zudov  may be understood on very general grounds as a consequence of a negative linear response conductivity. Ref. Yale  has presented a calculation, based on a specific microscopic model, showing that this negative linear response conductivity indeed may occur in the regime of magnetic field and microwave frequency in which the zero resistance state occurs. Taken together, we believe the present work and Ref. Yale  capture the essence of the experimental result.

###### Acknowledgements.

We are grateful to H. Stormer for bringing this problem to our attention and and to N. Read, H. Stormer, R. Willett and M. Zudov for valuable discussions. Two of us (A.A. and I.A.) were supported by the Packard foundation. A.A. was supported by the grant NSF-DMR-9984002 and A.J. M acknowledges NSF-DMR-00081075.