Spreading in narrow channels

0.47 MB·English

by  C. Dotti

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¨at 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 poten-
tials.
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 proper-
ties 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µ˜astheactualchemicalpo-
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 simula-
tion 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 val-
ues 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 assump-
tions 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 de-
fined 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)