Table Of Content

Spreading in narrow channels C. Dotti,1,2 A. Gambassi,1,2 M. N. Popescu,1,2 and S. Dietrich1,2 1Max-Planck-Institut fu¨r Metallforschung, Heisenbergstr. 3, D-70569 Stuttgart, Germany 2Institut fu¨r Theoretische und Angewandte Physik, 7 Universität Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany 0 Westudyalatticemodelforthespreadingoffluidfilms,whichareafewmolecularlayersthick,in 0 narrowchannelswithinertlateralwalls. Wefocusonsystemsconnectedtotwoparticlereservoirsat 2 differentchemical potentials, considering anattractivesubstratepotential at thebottom, confining n side walls, and hard-core repulsive fluid-fluid interactions. Using kinetic Monte Carlo simulations a we find a diffusive behavior. The corresponding diffusion coefficient depends on the density and is J boundedfrombelowbythefreeone-dimensionaldiffusioncoefficient,validforaninertbottomwall. 9 These numerical results are rationalized within thecorresponding continuumlimit. 2 PACSnumbers: 02.50.-r,05.70.Ln,68.15.+e,81.15.Aa ] h c I. INTRODUCTION 17, 18, 19, 20]. The spreading of such monolayers has e m been studied using a two-dimensional lattice gas Ising model [10, 21, 22] in which a half-space is occupied by a - In recent years, substantial progress has been made t particle reservoir. Recent KMC simulations and a con- a in the development of the ”lab on a chip” concept, i.e., tinuum analysis [23] of that model provided results in t s the integration of many physical and chemical processes good qualitative agreement with available experimental . (e.g.,transportthroughmicro-channels,mixing ofdiffer- t data, and a further extension to the case of chemically a ent fluids, chemical reactions) into a single device; en- m patterned substrate has been proposed [24]. tire laboratory setups, like a gas chromatograph, have Fluids in narrow channels have been investigated the- - beenminiaturizedona singlechip(for areview see,e.g., d oreticallyinthe contextofsingle-filediffusion, i.e., when Ref. [1]). In this context, microfluidics is becoming a n fluid particles cannot overtake each other (see, e.g., standardtoolinmanyapplications,rangingfrombiology o Refs.[25,26,27,28,29,30,31,32,33,34]). Suchsystems c (see, e.g., Ref. [2]) to the handling of toxic or rare sub- show the interesting feature of non-diffusive behavior of [ stances. Furtherscalingdowntonanofluidicsisexpected tracerparticles,whichstimulatedexperimental(see,e.g., to take place in the future [3]. Already now it is possi- 1 Refs. [33, 34]) and numerical [27, 29, 30, 31] interest. v ble to sculpture channels with lateral dimensions of few Here we present a lattice model for ultrathin films in 2 tens of nanometers [4] (for a review on such fabrication which multiple occupancy of a site is allowed (generaliz- 2 processes see Ref. [5]) and carbon nanotubes have been ing the single-occupancymodel of Refs. [10, 21, 22]) and 7 proposed as possible pipes in nanofluidics [6, 7]. Chemi- 1 cally patterned substrates have also been suggested as a in which the substrate-particle attractive interaction is 0 decayingasapowerlaw,whereastheparticle-particlein- solution for directed transport, gating, mixing, or sepa- 7 teractionisassumedtobe hard-corerepulsiveonly. This ration of fluids at the micro- and nano-scale [8]. In this 0 mimics the case in which the fluid-substrate interaction / case the channel consists of a strip of wettable material t strongly dominates over the actual attractive long-range a embedded in a non-wettable substrate so that the fluid partofthe fluid-fluidinteraction. Basedonthe phenom- m flows along the wettable region and is laterally confined ena occurring in this minimalistic model, the extension by the chemical contrast. - to the case in which the attractive part of the fluid-fluid d Ifoneofthedimensionsofafluidfilmiscomparableto interaction is relevant will be presented elsewhere. We n o thesizeofthefluidmolecules,ahydrodynamicaldescrip- shall restrict our analysis to a one-dimensional model, c tion of the film is no longer justified [9, 10, 11]. In this which can effectively describe fluids in extremely narrow : case the discrete nature of the fluid becomes important channels with a width which is less than twice the parti- v i andthefluidcannotbetreatedasacontinuuminthecon- cle diameter. The sidewalls act to confine the particles. X fined direction. In order to investigate such systems one The corrugation of the substrate potential both at the r possible approachis to carryout computer simulationof bottom and at the sides is incorporated effectively by a discrete models, e.g., molecular dynamics, kinetic Monte considering a lattice model for the particles. Due to the Carlo (KMC), or lattice Boltzmann simulations; recent small thickness of the channel the transversal variation workinthisdirectionincludesfluidsincarbonnanotubes of the substrate potential can be ignored. This model [6] and on chemically patterned substrates [12]. is supposed to mimic not only molecular fluids but also With the scaling down of microfluidic devices one has colloidal particles in solution, with the colloidal particle to deal with and may exploit the ultrathin precursor setting the length scale. film which spreads ahead of the bulk fluid. Experi- We discuss both the initial dynamics, in which a fluid mental studies have shown that in some cases such pre- filmfedbyareservoirgraduallyfillsthechannel,andthe cursor films have molecular thickness [13, 14, 15, 16, steady state, in which the fluid film in the channel is in 2 contact with two reservoirs at different chemical potentials. The paper is organized such that in Sec. II we define the model whereas in Sec. III the results of our Monte Carlo simulations are presented. The analyses of the y diffusion-like dynamics and of the steady-state properties are presented in Sec. IV. In Sec. V we discuss the mean-fieldcontinuumlimit ofthe model andrationalize, within this approximation, the results for the diffusion a coefficient presented in Sec. IV. Section VI summarizes the main findings and provides our conclusions. a (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) II. THE MODEL 0 (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1)(cid:0)(cid:1) Before specifying the rules defining the model, we fur- a ther describe the general physical picture of the type of systems we have in mind. 0 x (L+1)a As stated above, the fluid is assumed to be confined to a narrow, effectively one-dimensional channel. The FIG. 1: A typical configuration of the model. The possi- sidewalls are very high compared with the fluid particle ble movesin thebulk are indicated bystraight arrows, while diameter, so that the fluid cannot spill out of the chan- the curved arrows denote reservoirs-system exchanges. The nel. Thechannelwallsactontheparticlessuchthatonly substrate, including the exclusion zone of its top layer, cor- the vertical variation of the substrate potential matters. responds to the hatched area. The grey areas at x 0 and The left and the right end of the channel are connected ≤ x (L+1)a indicate reservoirs and the fluid particles are toafeedingandabsorbingparticlereservoir,respectively, sho≥wn as circles. and the channel is initially empty. The fluid film inside the channel is taken to be compact, i.e., molecules are effective diameter of a fluid particle, which is set by the densely packed to form vertical columns without vacan- hard-core repulsion between the particles. In the follow- cies. Thiscorrespondsto the caseinwhichthe substrate ing, the site indices and the distances will be expressed is strongly attractive and vacancies inside columns are inunits ofa. The twosites indexedwith 0 andL+1are eliminated onatime scale muchshorterthanthe typical theboundariesoftheleftandrightparticlereservoirs,re- time for exchanges of particles between columns. We spectively. Theothersites,[1,...,L](called‘bulk’inthe describe these exchanges in terms of rates, which are following),representthechanneloflengthL(assumedto related to the change of the energy of the system due be long, i.e., L 1). to the corresponding move. Particle exchanges between ≫ At every site x [0,L + 1] an occupation number columns and particle insertions and removals near the n N specifies t∈he number of particles piled in the reservoirsaretheonlyprocessesweconsider. Weassume x 0 ∈ column x (see Fig. 1). Within the column, particles cen- thatneitherevaporationnorcondensationtakesplacein- ters are located at integer y positions. Accounting for sidethechannel. Thisminimalisticmodelaimsatidenti- fluid-substrate hard-core repulsion we consider the posi- fyinggeneralaspectsandthemainqualitativefeaturesof tiony =0 aspassingthroughthe centersofthe particles fluids spreading in a strongly confined geometry, rather forming the top layer of the substrate. If the diameter than providing an accurate description of a particular of the substrate particles differs from that of the fluid physical situation. particles,onemayintroduceanextraparametertochar- Inspired by this picture, in Subsec. IIA we specify the acterize the position of the fluid-substrate contact layer; configuration space and the corresponding energy func- forsimplicity,however,weassumethatthefluidparticles tion. In Subsec. IIB the dynamical rules (i.e., the al- in the first layer are located at y =1 (see Fig. 1). lowed changes of the configurations and the associated The substrate is assumed to be uniform, and, consis- rates) governing the time evolution in the bulk are dis- tent with our one-dimensional model, two-dimensional cussed, whereas Subsec. IIC deals with the definition of semi-infinite in the y < 0 - direction. We denote the the dynamics at the feeding and absorbing boundaries. (attractive)pair interactionbetween a substrate particle (p) andafluidparticlebyU ,resemblingdispersionforces: sf A. Configurations and Hamiltonian U(p)(d)= −wds6f, ford≥1, (1) sf , ford<1, The model is defined on a one-dimensional (D = 1) (cid:26)∞ lattice, with sites [0,...,(L+1)a]. The distance a be- where d is the dimensionless distance in units of a be- tweentwoconsecutivesitesisassumedtobeequaltothe tweenthesubstrateparticlelocatedat(x′, y′ 0),and − ≤ 3 the fluid particle located at (1 x L,y 1). In the B. The rates and the dynamics ≤ ≤ ≥ case of pairwise additive interactions, for a semi-infinite substrate (x′ R, y′ R+) and in the continuum limit Inthissubsectionwedefineanddiscusstherateswhich ∈ ∈ (d 1), this leads to a total substrate potential govern the stochastic dynamics in the bulk, i.e., for x ≫ ∈ [1,L]. The dynamics at the boundaries, x = 0 and x = ∞ ∞ 1 L+1, will be discussed in the following subsection. U (y)= w dy′ dx′ , sf − sfZ0 Z−∞ [(y′+y)2+(x−x′)2(]32) jumWpeinatsosuomneeotfhtahteenaecahrepstarnteiicglehbionrt(hNeNc)ocloulmumnnxsxm+ay1 i.e., or x 1. We introduce the rate rCC′(y,y′), which is the − rate for a particle in column x and at given height y, to jump to the next column x+1 and at height y′. Within w′ Usf(y)= − ys4p, fory ≥1, (3) our aforementioned model assumption this process in- ( , fory <1, volves an instantaneous column height reduction by one ∞ incolumnxandacolumnheightincreaseincolumnx+1. This also means that the jumping particle is considered where w′ = 3π w . Within this ansatz, the particle- sf 32 sf to be able to squeeze into column x+1 at positiony′ by substrateinteractioninEq.(3)dependsontheheightyof pushing the particles above this position up by one unit; the particle only. The energy of the fluid configuration on the other hand if y′ is above the top particle of col- n ,...,n exposed to the substrate potential U is { 1 L} sf umnx+1,itfallsdowninordertoformagainacompact thus given by column. Accordingly the configurations C,C′, L nx C = n1,...,nx,nx+1,...,nL Hsf = Usf(yx), (4) C′ ={n ,...,n 1,n +1},...,n , (6) 1 x x+1 L Xx=1yXx=1 { − } represent the initial and the final configurations,for any where the inner sum is defined to be 0 if n = 0. pairy,y′. Analogousconsiderationscanbecarriedoutfor x Note that, following the discussion at the beginning of moves from x+1 to x, where the above configurations the present section, we assume that columns are always areinterchanged. Therefore,within our model, the rates densely packed, so that configurations are as depicted in rCC′(y,y′) depend only on the initial and final configu- Fig. 1: since the configurations are characterized com- rations C and C′, respectively. Accordingly, the depen- pletely by a succession of numbers n ,...,n , the en- denceony,y′isdropped. Weintroducethedimensionless 1 L { } ergy is a function of these numbers only, as in Eq. (4). rateu˜CC′,whichisalsoassumedtodependonlyonC,C′, The same form [Eq. (1)] of the pair potential is as- rCC′ sumed for the fluid particle - fluid particle interaction, u˜CC′ = , (7) ν where the corresponding interaction strength is denoted 0 by wff. Each pair of particles separated by a distance whereν0 fixesthetime-scaleofthemodelandweassume dff 1 contributes to the particle-particle energy, so ittobeindependentofthesourceandtargetcolumnsfill- ≥ thatthe totalenergyduetoparticle-particleinteractions ing, i.e., the same for any particle in the source column. can be written as (ν canbeinterpretedastherateforaparticletojumpto 0 the NN column if the energy happens to be unchanged 1 L L nx nx′ by the move.) In the following, times are measured in Hff = 2 units of ν0, i.e., the dimensionless simulation time t cor- xX=1xX′=1yXx=1yXx′=1 (5) responds to an actual time ta =t/ν0. Uf(pf) dff = (x−x′)2+(yx−yx′)2 , umWnexchtooocsoelutmhenrxat+es1u˜sufocrhatlhlaptodsseitbalielemdobvaelsanfrcoem col- (cid:16) p (cid:17) twoitbheUzef(rpfo)(i0f)n=x =0 a0nodrtnhxe′s=um0.sTovheertyoxtaalnedneyrxg′yafruentcatkioenn uu˜˜CCC′C′ =e−β∆H[C,C′] (8) is H[C]=Hff+Hsf whereC n1,...,nL character- is obeyed, where ∆H[C,C′] = H[C′] H[C] is the en- izes completely each configurat≡io{n. Note th}at this part ergy difference between the final (C′) a−nd the initial (C) of the Hamiltonian is restricted to the bulk; in general configuration. Detailedbalancehasbeenchoseninorder the reservoir-bulk interactions should also be accounted toensurethatthermalequilibriumisreachedinthelong- for separately. In the special case of the absence of long- time limit, if the two reservoirs of particles at the right rangeparticle-particleinteraction,i.e.,wff =0,boththe and the left of the channel are set to the same chemical bulkandthereservoir-bulkcontributionsvanish,andthe potential. A possible choice that satisfies the detailed energyisH[C]=Hsf. AsmentionedintheIntroduction, balance condition is in the following we shall discuss only this situation; the case wff =0 will be presented elsewhere. u˜CC′ =e−β2∆H[C,C′]. (9) 6 4 The chosen form of the rates [Eq. (9)] includes both Beforepassingto the definitionofthe dynamicsatthe “slow” (∆H > 0) and “fast” (∆H < 0) processes, and boundaries,webriefly commentonsimilarmodels which we implicitly assume that it captures essential features havebeenconsideredintheliterature. InRefs.[35,36]a of the real dynamics. classofdynamicalmodels,towhichourmodelbelongs,is Therateu˜CC′ isthesameforanyparticleinthesource introduced and studied. In this class of models the rates column x, so that the total rate uCC′ for a column to depend on both the source and the target column, they decrease its occupation number by one, while a given do not necessarily satisfy detailed balance, and jumps NNcolumnincreasesits ownoccupationnumberby one, occur not only between NN. (These models are known is in the literature as misanthropic processes.) The main uCC′ =nxu˜CC′ =nxe−β2∆H[C,C′]. (10) rtheseuilntfionfiRteesfsq.u[a3r5e,l3a6t]tiicsetihtaitsupnodsseirbcleerttoaionbtcaoinndaitnioenxsaicnt The rates in Eq. (10) are defined on the space of config- expression for the steady-state distribution. A concise urations specified by occupation numbers only. Detailed summary of these results can be found in Refs. [37, 38], balancestillholdsandthecorrespondingBoltzmannsta- where their relevance for the non-equilibrium dynamics tistical weight is ofinteractingparticleshasbeenstressed. Applied toour case, the results in Refs. [35, 36] recover the equilibrium e−βH BoltzmanndistributioninEq.(11),withtheHamiltonian p (n ,...,n ,...,n ) (11) B 1 x L ∝ n1!...nL! defined in Eq. (4), but do not provide information on the dynamics and steady-state distribution if chemical which accounts for ”particle undistinctness” by dividing drive, caused by different chemical potential for the two the Boltzmann factor by n !,...,n ! where n ! is the 1 L x reservoirsat the boundaries, is applied. number of choices to label the n particles in each x x ∈ [1,...,L] column. In the case of the particle-substrate interaction described by the Hamiltonian in Eq. (4), the C. Dynamics at the boundaries rate in Eq. (10) has the following explicit form: β ws′p ws′p We consider now the dynamics at the boundaries and u(n ,n )=n exp . (12) x x+1 x 2 (n +1)4 − n4 discuss two possible implementations. The first choice (cid:26) (cid:20) x+1 x (cid:21)(cid:27) is to fix the occupation number of the columns 0 and This formula emphasizes that u depends only on the oc- L+1atvaluesn andn ,respectively,andtoimpose, cupation numbers of the initial and the target column. 0 L+1 with some additional assumptions, the same dynamics The notation is such that the first argument stands for as in the “bulk”. The boundary dynamics changes the thesourcecolumn(herelocatedatxwithoccupationn ), x occupation number n of the first column of the system, while the second argument represents the target column 1 according to Eq. (12), while the occupation number n (here at x+1 with occupation n ). 0 x+1 is unchangedby the move,andthe same holds atx=L. Assuming that the dynamics leads to a diffusion-like One can physically motivate such a choice by assuming behavior (as will be discussed in Sec. IV), some quali- that the particle exchanges within the reservoir are so tative features of the diffusion coefficient as a function fast, so that a particle extracted from the reservoir is of the local density can be anticipated from the general immediatelyreplaced. Underthisassumptionthedensity properties of the rates in Eq. (12). First, consider the ofparticlesinthereservoirissimplyn . Whilethischoice situation in which both n = n and n = m are large 0 x x±1 seems to be rather natural, as explained in the following comparedto(βw′ )1/4. ThentheexponentinEq.(12)is sp the equilibrium (i.e., for n = n ) properties display very small and u(n,m) n, so that the model reduces 0 L+1 → an unexpected feature, i.e., a jump discontinuity in the to free particles diffusing in D = 1. The same conclu- density between the reservoirsand the system. sion holds for n = m +1, in which case the exponent The total Hamiltonian in Eq. (4) is a sum of single- is zero leading to u(n,m) = n. In general, u(n,m) ≷ n column terms, so that the equilibrium grand canonical if n ≷ m+1; accordingly, jumps from high columns to distribution factorizes: low columns are faster than in the free case, while the oppositeprocessesareslower. Thismeansthatdiffusion, L 1 whichtendstosmoothoutdensitygradients,isenhanced P ( n ,...,n )= p (n )e−µnk, (13) eq 1 L w k { } Z(w,µ) by the exponential factor in Eq. (12). Since at low den- k=1 Y sities most of the configurations are composed either of where empty columns or of columns occupied by one particle, the most probable rate is w(1,0) = 1, which results in ∞ L free diffusion at low densities. Z(w,µ)= p (n)e−µn (14) w These considerations lead to the conclusion that the ! n=0 diffusion coefficient is expected to exhibit a peak at rel- X atively low densities, because the rates exhibit the max- is the total partition function, µ=βµ˜ is the dimension- imal difference with free diffusion rates if the target col- lesschemicalpotential,withµãstheactualchemicalpo- umn is empty. tential, p is a non-normalized single-column statistical w 5 weight, 1.2 1 p (m)= exp[2wh(m)], (15) w m! w =βws′p/2 is a dimensionless quantity, and n0 ρeq = n0 / 0.6 m 1 eq ρ , forn 1, h(m)= k4 ≥ (16) k=1 0X, forn=0. w = 2 w = 4 The chemical potential µ controls the mean density of 0 2 4 particles in the system, n 0 np (n)zn w ρ (z) Neq = z∂ ln[Z(w,z)]= nX≥1 , eq ≡ L L z p(n)zn FIG.2: Theequilibriumdensityρeq [Eqs.(13)-(17)and(19)], dividedbythereservoiroccupation numbern0,asafunction n≥0 X (17) of n0 for w=2 and w=4, respectively. where N is the mean total number of particles in the eq system,andz =e−µ isthefugacity. Detailedbalancefor two reservoirs at the boundaries have different occupa- the rates at the left reservoir reads tion numbers (n = n ); nevertheless we recover the 0 L+1 6 equilibrium density in the first column. This shows how P ( n ,...,n )u(n ,n )= eq 1 L 0 1 { } (18) it is possible to control the densities at the first (x = 1) =P ( n +1,...,n )u(n +1,n ) eq 1 L 1 0 andlast(x=L)site byvaryingthe occupationnumbers { } n and n . In order to obtain arbitrary densities, it where the probability per unit time of inserting a par- 0 L+1 is necessary to take n and n to be continuous, thus ticle into the column at x = 1 is compared with the 0 L loosingthedirectphysicalinterpretationoftheseparam- corresponding probability per unit time of removing the eters. Moreover, as shown in Fig. 2, ρ drops sharply particle in the same column. A similar condition has to eq for n . 0.8, and in the range 0 . ρ . 0.8 a high hold at the right end x = L of the system. Combining 0 eq numerical accuracy would be required to determine the Eq. (18) with Eqs. (12) and (13) leads to corresponding value of n . 0 These two problems can be solved by generalizing the 1 1 z =n0exp −w n4 + (n +1)4 . (19) dynamics at the boundaries as follows. In Eq. (12) the (cid:26) (cid:20) 0 0 (cid:21)(cid:27) terms depending on n and n , i.e., the properties of 0 L+1 the reservoirs, are replaced by constants α, γ, δ, and κ Equations (17) and (19) give the equilibrium density ρ eq in the following way: in the system as a function of n . As expected, if n is 0 0 large ρeq coincides with n0: in Eq. (19) n0 w1/4 im- w w pexlipes(2zw≈ζ(n40)),/ann!d, isnoEtqh.a(t17f)orn0n≫w1/3mimaxpl(≫iwes1/p3w,(wn1)/4≈) uα(n1)=αexp(cid:20)(n1+1)4(cid:21), uγ(n1)=γn1exp(cid:20)−n41(cid:21), 0 ≫ w w Eq. (17) reduces to u (n )=δexp , u (n )=κn exp , δ L (n +1)4 κ L L −n4 (cid:20) L (cid:21) (cid:20) L(cid:21) ∂ (21) ρ e−zz ez =z n , n max(w1/3,w1/4). eq 0 0 ≈ ∂z ≈ ≫ (20) whereuα isthe rateforparticleinsertionintocolumnn1 In the range of substrate potential strength we investi- fromthe leftreservoir,uγ is the rateforparticle removal gated (0.5 < w < 5) the approximation ρeq n0 is from column n1 into the left reservoirand uδ,uκ are the valid if n > 5. As the substrate potential str≈ength is correspondingratesatx=L. Imposingdetailedbalance 0 increased, this threshold increases,while it tends to zero [Eq. (18)] for the rates defined in Eq. (21) gives for w 0. For densities lower than the threshold, the α δ reservo→ir occupation number n does not coincide with e−µ =z = = . (22) 0 γ κ the density ρ of the equilibrium system as obtained eq from Eqs. (13)-(17) and (19) (see Fig. 2). For example, Inthenon-equilibriumcasethefugacityoftherightreser- in the case n = 3 one has ρ 3.19 for w = 2 and voir, denoted as z = δ/κ, and the one of the left 0 eq L+1 ≃ ρ 3.23 for w = 4. These densities are in very good reservoir, i.e., z = α/γ, are different (z = z ). Us- eq 0 0 L+1 ≃ 6 agreementwith the simulationdata for the density ρ in ing these two fugacities in Eq. (17), the densities of the 1 the first column (see, c.f., Fig. 5(b)). In the simulations twocorrespondingequilibriumsystemsarefound;wede- we investigated a non-equilibrium situation in which the fine them to be the reservoir densities. In simulations 6 we proceed backwards: first we choose ρ = ρ (z ) and A. Spreading 0 eq 0 ρ = ρ (z ), and then we find the corresponding L+1 eq L+1 ratios by inverting Eq. (17). Setting γ = κ = 1 implies For the spreading regime the right reservoir has been α=z , δ =z so that the inversion is simpler. 0 L+1 convertedinto a particle sink by setting n =0, while L+1 n = 11 and L = 1000. In the simulations performed 0 with these values of the parameters the leakage of parti- III. MONTE-CARLO SIMULATIONS cles through the sink is negligible for times t . 104. We studied both the initial-time dynamics by setting τ =0 0 The continuous time dynamics defined by the rules andτ 104 andthe long-time behaviorin which both tot described in Subsecs. IIB and IIC is simulated using a reservoir≤s play a role (see, c.f., Fig 4, w = 0.5); in this Kinetic Monte Carlo (KMC) method [39]. At every step latter case we have chosen τ = 5 104 and τ = 104 tot 0 an increment ∆t for the time variable is drawn from the in order to reduce the CPU and me×mory requirements. distribution In the spreading regime we implemented the simulation 1 averageby drawing different sequences of (pseudo-) ran- P(∆t)= exp[S(n0,...,nL+1)∆t], domnumbers while keeping the initial conditionfixed so S(n ,...,n ) 0 L+1 that in this case is the ensemble average; the typical L h·i numberofrunsweaveragedoveris2000. Westudiedthe S(n ,...,n )= [u(n ,n )+u(n ,n )], 0 L+1 x x+1 x+1 x shape of the density profile ρ(x,t), defined in Eq. (24), xX=0 as a function of the interaction strength w in the range (23) 0.5.w.1.5. where S is the total rate to leave the configuration Since forw =0 the dynamicsreducestofree diffusion, n ,...,n . The move to perform is then chosen ac- it is natural to check if in the general case (i.e., w = 0) {cor0dingtoLt+h1e}weightu/S ofitsrate. Weusedaclassical the profiles show a diffusive scaling. Plotting them6 as N-fold way algorithm[40], which has the advantage that a function of λ = x/√t we indeed obtain a collapse of the selected moves are accepted without rejection. The data measured at different times, as shown in Fig. 3. modeldependsonfourparameters: thesubstrateinterac- The rescaledprofiles are similar to the one for free diffu- tionstrengthw =w′ β/2inunitsofthethermalenergy, sion, except for a small bend at densities around 1 (see sp the two boundary densities n and n , and the length Fig 3(b)), which depends on the interaction strength w. 0 L+1 L of the system. The simulations have been performed The diffusive scaling is confirmed by the time evolution up to a maximum time τ and quantities have been of the total mass [Eq. (25)] shown in Fig. 4, which is ex- tot measured after the initial time τ . The simulations cov- pectedtoevolveasM √t. Weobservedeviationsfrom 0 ∝ eredboth the spreadingandthe steady-stateregimeand thisbehavioronlyatshorttimes,whentheboundarydy- in both cases we sampled the same set of quantities. In namics dominates, and at long times, when the leakage order to keep the notation simple we indicate averages of particles through the sink becomes relevant. always with , but, as described in the following, the h·i meaning of the symbol is different in the two situations considered. B. Steady state The mean density of particles, i.e., the density at site x and time t is defined as The steady state is reached after running the simulation for an initial thermalization time τ (τ 105 for 0 0 ρ(x,t)= n (t) , (24) ≈ x L = 1000, chosen by checking that for t > τ the ob- h i 0 servables are time independent) and saving the configu- while the total mean number of particles (or mean total rations generated every sampling time interval τ , with mass) is s τ = 200 for L = 1000. The average for the observ- s h·i L ables defined above is taken over this set of configura- M(t)= n (t) . (25) tions. The choice for τ is a compromise between speed x s h i xX=1 andhavingassmallcorrelationsbetweentheNsmeasure- ments as possible. We assume that the total simulation The transportpropertieshavebeen studiedusingthe in- time τ =N τ +τ is sufficient to explore a significant tegrated particle current at site x and time t, tot s s 0 part of the phase space, so that the performed average coincides with the average over the (unknown) steady J(x,t,∆t)=∆n (t,t+∆t) ∆n (t,t+∆t), x,x+1 x+1,x − state distribution. Note that this assumption is justified (26) because no signs of dynamical phase transitions (which where∆nx,x′(t,t+∆t)isthenumberofparticlesjumping from site x to site x′ within the time interval ∆t. We wouldintroduceextremelylongtime scales)arefoundin the simulations. define a mean instantaneous current as The density ρ(x) [Eq.(24)]exhibits aprofile smoothly J(x,t,∆t) interpolating between the two reservoirs(see Fig. 5) and j(x,t) = lim h i. (27) h i ∆t→0 ∆t slightly deviating from the corresponding free diffusion 7 w=1.5 w = 1.5 t=5.0 ×103 free diff. 8 1.0 × 104 8 t = 5.0 × 103 ρ(x,t) 22..05 ×× 110044 ρ(x,t) 12..00 ×× 110044 4 4 2.5 × 104 (a) (b) 0 0 500 1000 0 1.6 3.2 x λ = x / t1/2 FIG.3: Time-dependentdensityprofilesforspreadinginasystemwithL=1000,forn0 =11,nL+1 =0,andw=1.5attimes t=5 103 (+), 104 (⊡), 2 104 ((cid:4)), and 2.5 104 ( ) as a function of x (a), and of λ=x/√t (b), respectively. The solid × × × ⊙ line indicates thescaling function for free diffusion (i.e., w=0). to local particle conservation. The instantaneous steady state current j(x) (a dependence on x is indicated to h i 12.9 recall the random fluctuations around the mean value j ) can then be obtained from h i 2 J(x,0,τ ) J(x,0,τ ) 1/ j(x) = tot − 0 . (28) / t 12.6 h i Nsτs M Note that in Eq. (28) the current integrated over the w = 0.5 thermalizationtime, J(x,0,τ ), whichdepends onx and 0 w = 1.0 t, has been subtracted. We have calculated j(x) for w = 1.5 x [1,...,L] leading to j =L−1 L j(x) h. i 12.3 ∈ h i x=1h i 0 2 × 104 4 × 104 P t IV. ANALYSIS OF THE SIMULATION DATA FIG.4: TotalmassM/√tforspreadingasafunctionoftime A. Methods to determine the diffusion coefficient t for different values of the interaction w with the substrate: w = 0.5 (total simulation time τtot = 5 104 and initial sampling time τ0 = 104) ( ), w = 1.00×(τtot = 4 104, Guided by the results of the KMC simulations of the τ0 =0) ((cid:4)), and w=1.50 (τ×tot =2.5 104, τ0 =0) (×) For microscopicmodelwe expectthat acontinuum(in space × ⊙ all symbols L=1000, n0 =11, and nL+1 =0. and time) description for the behavior of the model at long times and large spatial scales is possible, i.e., that the hydrodynamic limit exists and is well defined. A profile which is a straight line. This deviation is consid- rigorous proof has been provided for a small number of eredinmoredetailinAppendix B,whereitsdependence models (see, e.g., Ref. [41]). In the present case such an onboththeinteractionandthereservoirsdensitiesisan- explicit derivation appears to be a difficult task. alyzed. Inthesteadystatebothreservoirsplayaroleand Assuming that the particle density ρ and the current finite-sizeeffectshavetobechecked;itturnsoutthatfor j , defined in Eqs. (24) and (27), are smooth functions L > 200 there is no detectable dependence of the data ohfithepositionxandofthetimet,thelocalconservation onthe particularvalue ofLotherthanatrivialrescaling of particle density in the bulk (which is implicit in the of the density profile. dynamicsofthe model) is expectedto takethe formofa We have also determined the current j defined in continuity equation: h i Eq. (27): in the steady state j does not depend on t, h i sothatJ(x,t,∆t)= j ∆tforanysufficientlylargetime ∂ ρ(x,t)= ∂ j(x,t) . (29) t x h i − h i interval ∆t. J can be obtained by measuring the flux of particles between any two sites x and x+1, because in The results of the simulations strongly indicate a dif- thesteadystatethecurrent j doesnotdependonxdue fusive scaling at long times and large spatial scales, i.e., h i 8 3.4 w = 0.4 3 w = 1.5 3.23 free diff. 6 3.19 3 x) x) 5 10 ρ( ρ( 1.5 3 w=2 w=4 free diff. (a) (b) 0 0 0 500 1000 0 500 1000 x x FIG. 5: Steady-state density profiles in a system of length L=1000 for (a) w=0.4 (⊡) and w=1.5 ( ), n0 =8, nL+1 =0; ⊙ (b)w=2(⊡)andw=4( ),n0 =3,nL+1 =0. Theinsetin(b)isaclose-upviewofthevicinityoftheleftreservoirshowing ⊙ ρeq ρ1 >n0 [ρeq 3.19 for w=2 and ρeq 3.23 for w=4; see Eqs. (17) and (19) and the discussion in the main text]. In ≃ ≃ ≃ both (a) and (b) the free diffusion is indicated bya full line. underestimated diffusion coefficient for small values of ρ. The most severe effect is probably due ρ(x,t)=ρ¯(x/√t),suggestingthatthedynamicsamounts to the derivative d ρ¯(λ), while the errors in the dλ to non-linear diffusion: integral can be partially corrected by using larger lattice sizes. j(x,t) = D(ρ)∂ ρ(x,t) x h i − ⇒ (30) ∂ ρ(x,t)=∂ [D(ρ)∂ ρ(x,t)]. t x x Steady state. In the stationary state the current • and the density are constant with respect to time, Basedonthe density profileρ(x,t) fromthe simulations, so that Eq. (30) leads to itispossibletoextractthefunctionD(ρ)fromthedatain the spreading and the steady state regime, respectively, as follows. j D(ρ)= h i, (33) −ρ′ Spreading. Using the scaling behavior ρ(x,t) = • ρ¯(λ = x/√t) (see the results in Sec. III), Eq. (30) where ρ′(x)= d ρ(x). Note that the computation reduces to dx of D in this regime does not require any assumptions on ρ and ρ ′, as in the previous case, and d d d λ ρ¯(λ)= D(ρ¯(λ)) ρ¯(λ) . (31) therefore no systematic deviations from the actual dλ dλ dλ (cid:20) (cid:21) diffusion coefficient are expected. Assumingthat dρ¯D(ρ(¯λ))andρ¯(λ)vanishforλ , which is sudλpported by the simulation da→ta Steady state in quasi-equilibrium. A steady ∞(L 1 is an approximation for L ), inte- • statedeviatingonlyslightlyfromtheequilibriumis grat≫ing Eq. (31) and inverting ρ¯(λ) in→to∞λ(ρ¯), one realizedbyimposingreservoirdensitieswhichdiffer finds only slightly. Equation (33) leads to d ρ¯ 1 ρ1 D(ρ¯)= λ(ρ¯) dρλ(ρ). (32) j = dρD(ρ), (34) dρ¯ Z0 h i LZρL This method might be inaccurate for small densi- where ρ and ρ are the densities at the first and tiesduetoasystematiceffect. Forx Ltheprofile 1 L ≈ thelastsite,respectively,andListhelengthofthe bends in order to fulfill the condition n =0, so L+1 system. For ρ ρ ρ one has thatits derivative d ρ¯is largerthanthe derivative | 1− L|≪ 1 dλ ofa profileρ¯ in the infinite lattice: d ρ¯ < d ρ¯. ∞ dλ ∞ dλ ρ +ρ L j The integral in Eq. (32) is also underestimated, D 1 L h i . (35) since λ(ρ → 0) 9 ∞. These effects lead to an (cid:18) 2 (cid:19)≈ (ρ1−ρL) 9 w=1.5 steady state steady state: w=0.40 1.8 spreading 0.80 Eq. (56) 1.00 Eq. (56): w=0.40 ρ() ρ) 1.4 0.80 D ( 1.00 1 D (a) 1 (b) 0 4 8 0 4 8 ρ ρ FIG. 6: (a) Nonlinear diffusion coefficient D as a function of the density, obtained from the simulation data in the spreading regime((cid:4))andinthesteadystate(⊡),aswellasthecorrespondinganalyticalresult[fullline,Eq.(56)]. Allthedatacorrespond to w=1.5 and L=1000. (b) Nonlinear diffusion coefficient [Eq. (33)] from steady-state simulation data for w =0.40 (+), w =0.80 (⊡), and w =1.00 ((cid:4)) (L=1000). Analytical results (lines) for D(ρ) are calculated from Eq. (56). B. Results from the spreading regime tainedfromthespreadingdataandfromthesteady-state datais goodforρ&0.5andevenbetter for ρ&2,which showsthatthediffusionpicturedescribedinSubsec.IVA In order to extract D(ρ) from the numerical data for leads to consistent results. given w and L, we have considered all the profiles in Simulations under quasi-equilibrium conditions have the scaling regime. We have binned the ρ axis with been restricted to the case w = 2 because the quantita- ∆ρ = 0.05, averaged all the values λ in each bin, and tive agreement between the results obtained from quasi- evaluated the function λ(ρ) by interpolation of the re- sulting data points, while d ρ¯(λ) has been obtained by equilibrium, using Eq.(35) for D(ρ), and those obtained dλ in the steady state is satisfactory (see Fig. 7). Due finite differences. The resultsforD(ρ) obtainedby using to the small difference in density of the two reservoirs Eq. (32) are shown in Fig. 6(a). While it appears that, (δρ = 0.01) the average current is very small and re- for large values of ρ, D(ρ) 1 as expected from the → quires accurate measurements. The data are obtained corresponding discussion in Subsec. IIB, for ρ 0 the → by averaging over 107 configurations or more (approxi- diffusioncoefficientgoestozeroduetothesystematicer- rorinthederivative d ρ¯(λ),asexplainedinSubsec.IVA. mately 100 times more than for the steady state data), dλ leading to a high precision for the density profile, too. The noise at large values of ρ is due to determining the Theautocorrelationtimefortheaveragedensityhasbeen derivatives numerically, because for large ρ the spatial carefullycheckedandthetime intervalsbetweensamples fluctuations of the density are stronger. have been chosen in order to minimize correlations. The The diffusion coefficient is peaked and the substrate procedure allows one to estimate reliably the statistical potential enhances diffusion (see Subsec. IIB). The po- error for the diffusion coefficient, shown by the error- sition of the peak (ρ 1) cannot be predicted by quali- ≃ bars in Fig. 7(a). Note that the results obtained from tative arguments, but is in the range of low densities, as steady-statesimulationswith biggerreservoirdifferences expected from the discussion in Subsec. IIB. lie within the errorbars. C. Results from steady-state and quasi-equilibrium V. CONTINUUM DESCRIPTION regimes In this section we derive the nonlinear diffusion equa- In order to obtain D(ρ) from the steady-state data tion corresponding to the continuum limit of the model. by using Eq. (33) we have measured the average cur- The equation is derived from the microscopic dynamics rent and the average density profile. The latter has by using severalsimplifying assumptions. We start from been appropriately binned (the density is averaged over themasterequation,whichdescribesexactlythedynam- 5sites), inorderto be able toevaluate the derivativevia ics of the model and in its most general form can be finite differences. The corresponding results are shown written as in Figs. 6(a), 6(b), and 7. We note that D(ρ 0) 1. → → ∂ P (C)= (C,C′)P (C′), (36) The correct behavior at low densities is captured, while t t t M for large ρ the data are still rather noisy, because the XC′ method to extract D relies on determining derivatives where C is a generic configuration, while (C,C′) en- M numerically. The overall agreement between results ob- codesthetransitionsfromtheconfigurationC tothecon- 10 (a) w=2 (b) w=4 steady state 2.6 steady state Eq. (56) 2 quasi equil. Eq. (56) ρ) ρ) ( ( D D 1.8 1 1 1 2 0 1 2 3 ρ ρ FIG. 7: Nonlinear diffusion coefficient D(ρ) for large values of w. The open squares (⊡) are obtained from simulation data in the steady-state regime [Eq. (33)] for w = 2 (a) and w = 4 (b) (L = 1000). The lines correspond to analytical calculations [Eq.(56)]. Thecrosses( )witherrorbarsin(a)areobtainedfromsimulationdataunderquasi-equilibriumconditions[Eq.(35)]. × figuration C′. In our case, the operator can be split 1,n ,n + 1,...,n , n ,...,n ,n + 1,n x x+1 L 1 x−1 x x+1 M } { − into bulk ( ) and boundary ( ) terms, so that 1,...,n , n ,...,n ,n 1,n +1,...,n with b s L 1 x−1 x x+1 L M M } { − } x [2,...,L 1], so that using Eqs. (39), (40), (A1), ∈ − ∂ P (C)= [ (C,C′)+ (C,C′)]P (C′). (37) and (A2) one obtains in the bulk t t s b t M M C′ X ∂ ρ(x,t)= [ j(x+1,t) j(x,t) ], (41) t − h − i The operator describes bulk moves, upon which b M particles are exchanged between columns at sites x where ρ(x,t)= n , x [2,...,L 1], and x ∈ h i ∈ − [1,...,L] and which are associated with the rates defined in Eq. (12). The operator Ms, describes boundary hj(x,t)i=hu(nx−1,nx)−u(nx,nx−1)i (42) moves upon which particles are inserted into or removed is the mean local and instantaneous current in Eq. (27). from the system at the sites x = 1 or x = L, and which Assumingthattheprobabilitydistribution[Eq.(37)]fac- are associated with the rates introduced in Subsec. IIC. torizescompletely,i.e.,withinthe meanfieldapproxima- Explicit expressions for the operators and are Mb Ms tion given in Appendix A. The evolution of the ensemble average of a generic L (time-indipendent) operator can be obtained from P(C,t)= p˜ (n ,t), (43) y y O Eq. (36) as y=1 Y ∂ = (C) (C,C′)P (C′), (38) the average rates in Eq. (42) [for the definition of the t t hOi O M rates see Eq. (12)] reduce to C,C′ X ∞ ∞ where (C) is the value ofthe operator for configura- u(n ,n ) tion CO. Recalling that (C,C) = O (C′,C), h x x+1 i≡ ··· it is straightforwardto oMbtain − C′6=CM nX1=1 nXL=1 P L w w (44) p˜ (n ,t) n exp ∂thOi= K(C)P(C), (39) yY=1 y y (cid:26) x (cid:20)(nx+1+1)4 − n4x(cid:21)(cid:27) XC =f˜(w,t,x)f˜(w,t,x+1), 1 2 where (C) is the jump moment of the operator de- K O where fined as ∞ w (C) [ (C′) (C)] (C′,C). (40) f˜(w,t,x)= p˜ (n ,t)n exp (45a) K ≡ O −O M 1 x x x n4 CX′6=C nXx=1 (cid:18) x(cid:19) We consider now the operator , defined as (C) = and x x nx. Accordingly, the configuraNtions C′, for wNhich the ∞ w jump moments defined inEq.(40) arenonzero,areC′ = f˜(w,t,x)= p˜ (n ,t)exp . (45b) {n1,...,nx−1+1,nx−1,nx+1,...,nL},{n1,...,nx−1− 2 nXx=0 x x (cid:20)−(nx+1)4(cid:21)