Transition from single-file to two-dimensional diffusion of interacting particles in a quasi-one-dimensional channel
D. Lucena, D. V. Tkachenko, K. Nelissen, V. R. Misko, W. P. Ferreira, G. A. Farias, F. M. Peeters
aa r X i v : . [ c ond - m a t . s o f t ] M a r Transition from single-file to two-dimensional diffusion of interacting particles in aquasi-one-dimensional channel
D. Lucena , , D. V. Tkachenko , K. Nelissen , , V. R. Misko , W. P. Ferreira , G. A. Farias , and F. M. Peeters , Departamento de F´ısica, Universidade Federal do Cear´a,Caixa Postal 6030, Campus do Pici, 60455-760 Fortaleza, Cear´a, Brazil Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerpen, Belgium (Dated: November 8, 2018)Diffusive properties of a monodisperse system of interacting particles confined to a quasi -one-dimensional (Q1D) channel are studied using molecular dynamics (MD) simulations. We calculatenumerically the mean-squared displacement (MSD) and investigate the influence of the width ofthe channel (or the strength of the confinement potential) on diffusion in finite-size channels ofdifferent shapes (i.e., straight and circular). The transition from single-file diffusion (SFD) to thetwo-dimensional diffusion regime is investigated. This transition (regarding the calculation of thescaling exponent ( α ) of the MSD h ∆ x ( t ) i ∝ t α ) as a function of the width of the channel, is shownto change depending on the channel’s confinement profile. In particular the transition can be eithersmooth (i.e., for a parabolic confinement potential) or rather sharp/stepwise (i.e., for a hard-wallpotential), as distinct from infinite channels where this transition is abrupt. This result can beexplained by qualitatively different distributions of the particle density for the different confinementpotentials. PACS numbers: 05.40.-a, 66.10.C-, 82.70.Dd, 83.10.Rs
I. INTRODUCTION
There is a considerable theoretical and practical inter-est in the dynamics of systems of interacting particlesin confined geometries [1]. Single-file diffusion (SFD)refers to a one-dimensional (1D) process where the mo-tion of particles in a narrow channel (e.g., quasi -1D sys-tems) is limited such that particles are not able to crosseach other. As a consequence, the system diffuses asa whole resulting in anomalous diffusion. The mecha-nism of SFD was first proposed by Hodgkin and Keynes[2] in order to study the passage of molecules throughnarrow pores. Since the order of the particles is con-served over time, this results in unusual dynamics ofthe system [3, 4], different from what is predicted fromdiffusion governed by Fick’s law. The main character-istic of the SFD phenomena is that, in the long-timelimit, the MSD (mean-square displacement, defined as h ∆ x ( t ) i = h P Ni =1 (1 /N )[ x i ( t + ∆ t ) − x i ( t )] i ∆ t ) scaleswith time as h ∆ x ( t ) i ∝ t . . (1)This relation was first obtained analytically in the pi-oneering work of Harris [5]. Recent advances in nan-otechnology have stimulated a growing interest in SFD,in particular, in the study of transport in nanopores [6, 7].Ion channels of biological membranes and carbon nan-otubes [8] are examples of such nanopores. The macro-scopic flux of particles through such nanopores is of greatimportance for many practical applications, e.g., parti-cle transport across membranes is a crucial intermediatestep in almost all biological and chemical engineering pro-cesses. SFD was observed in experiments on diffusion ofmolecules in zeolite molecular sieves [9]. Zeolites withunconnected parallel channels may serve as a good real- ization of the theoretically investigated one-dimensionalsystems. SFD is also related to growth phenomena [10].The theoretical background of SFD was developed inearly studies on transport phenomena in 1D channels [11–13]. It is also interesting to learn how the size of the sys-tem will influence the diffusive properties of the system.SFD in finite size systems has been the focus of increas-ing attention since there are few exact theoretical resultsto date [14–16], which showed the existence of differentregimes of diffusion.Colloidal systems, complex plasmas and vortex mat-ter in type-II superconductors are examples of systemswhere SFD may occur. The use of colloidal particles istechnically interesting since it allows real time and spa-tial direct observation of their position, which is a greatadvantage as compared to atoms or molecules, as shownrecently in, e.g., the experimental study of defect inducedmelting [17]. One typically uses micro-meter size col-loidal particles in narrow channels, as shown in [18, 19].The paramagnetic colloidal spheres of 3.6 µm were con-fined in circular trenches fabricated by photolithographyand their trajectories were followed over long periods oftime. Several other studies have focused on the diffusiveproperties of complex plasmas. A complex plasma con-sists of micrometer-sized (“dust”) particles immersed ina gaseous plasma background. Dust particles typicallyacquire a negative charge of several thousand elementarycharges, and thus they interact with each other throughtheir strong electrostatic repulsion [20].Systems of particles moving in space of reduced dimen-sionality or submitted to an external confinement poten-tial exhibit different behavior from their free-of-bordercounterparts [21]. The combined effect of interaction be-tween the particles and the confinement potential plays acrucial role in their physical and chemical properties [22].In Ref. [23], it was found that SFD depends on the inter-particle interaction and can even be suppressed if the in-teraction is sufficiently strong, resulting in a slower sub-diffusive behavior, where h ∆ x ( t ) i ∝ t α , with α < . h ∆ x ( t ) i ∝ t . ( ∝ t . ). Recall thatthe MSD of a tagged hard-sphere particle in a one dimen-sional infinite system is characterized by two limiting dif-fusion behaviors: for time scales shorter than a certaincrossover time τ c = 1 /Dρ , where D is the diffusion coef-ficient and ρ is the particle concentration, h ∆ x ( t ) i ∝ t . which is referred to as the normal diffusion regime [24].For times larger than τ c , the system exhibits a subdiffu-sive behavior, with the MSD h ∆ x ( t ) i ∝ t . , which char-acterizes the single-file diffusion regime. Between thesetwo regimes, there is a transient regime exhibiting a non-trivial functional form.However, in case of a finite system of diffusing particles(e.g., a circular chain or a straight chain in the presence ofperiodic boundary conditions), the SFD regime (i.e., with h ∆ x ( t ) i ∝ t . ) does not hold for t → ∞ , unlike in aninfinite system. Instead, for sufficiently long times, theSFD regime turns to the regime of collective diffusion,i.e., when the whole system diffuses as a single “parti-cle” with a renormalized mass. This diffusive behaviorhas been revealed in experiments [18, 25] and theoreti-cal studies [23, 26–29]. This collective diffusion regime issimilar to the initial short-time diffusion regime and it ischaracterized by either h ∆ x ( t ) i ∝ t . , for overdampedparticles (see, e.g., [18, 23, 26]) or by h ∆ x ( t ) i ∝ t . (followed by the MSD ∝ t . ), for underdamped systems[27, 29]. Correspondingly, the time interval where theSFD regime is observed becomes finite in finite size sys-tems. It depends on the lenght of the chain of diffusingparticles: the longer the chain the longer the SFD timeinterval. Therefore, in order to observe a clear power-law behavior (i.e., h ∆ x ( t ) i ∝ t α ) one should considersufficiently large systems.Here we focus on this intermediate diffusion regimeand we show that it can be characterized by h ∆ x ( t ) i∝ t α , where 0 . < α < .
0, depending on the width (orthe strenght of the confinement potential) of the chan-nel. We analyze the MSD for two different channel ge-ometries: (i) a linear channel, and (ii) a circular channel.These two systems correspond to different experimentalrealizations of diffusion of charged particles in narrowchannels [18, 30]. The latter one (i.e., a circular chan-nel) has obvious advantages: (i) it allows a long-timeobservation of diffusion using a relatively short circuit,and (ii) it provides constant average particle density andabsence of density gradients (which occur in, e.g., a lin-ear channel due to the entry/exit of particles in/fromthe channel). Thus circular narrow channels were used in diffusion experiments with colloids [18] and metalliccharged particles (balls) [25]. Furthermore, using differ-ent systems allows us to demonstrate that the resultsobtained in our study are generic and do not depend onthe specific experimental set-up.This paper is organized as follows. In Sec. II, we in-troduce the model and numerical approach. In Sec. IIIdiffusion in a system of interacting particles, confined toa straight hard-wall or parabolic channel, is studied asa function of the channel width or confinement strength.In Sec. IV, we discuss the possibility of experimental ob-servation of the studied crossover from the SFD to 2Ddiffusive regime. For that purpose, we analyze diffusionin a realistic experimental set-up, i.e., diffusion of mas-sive metallic balls embedded in a circular channel withparabolic confinement whose strength can be controlledby an applied electric field. The long-time limit is ana-lyzed in Sec. V using a discrete site model. Finally, theconclusions are presented in Sec. VI.
II. MODEL SYSTEM AND NUMERICALAPPROACH
Our model system consists of N identical chargedparticles interacting through a repulsive pair potential V int ( ~r ij ). In this study, we use a screened Coulomb po-tential (Yukawa potential), V int ∝ exp( − r/λ D ) /r . Inthe transverse direction, the motion of the particles isrestricted either by a hard-wall or by a parabolic con-finement potential. Thus the total potential energy ofthe system can be written as: H = N X i =1 V c ( ~r i ) + N X i>j =1 V int ( ~r ij ) . (2)The first term in the right-hand side (r.h.s.) of Eq. (2)represents the confinement potential, where V c ( ~r i ) isgiven by: V c ( ~r i ) = (cid:26) | y i | ≤ R w / ∞ for | y i | > R w / , (3)for the hard-wall confinement, V c ( ~r i ) = 12 mω y i , (4)for parabolic one-dimensional potential (in the y -direction), and by V c ( ~r i ) = β ( r − r i ) , (5)for parabolic circular confinement. Here R w is the widthof the channel (for the hard-wall potential), m is the massof the particles, ω is the strenght of the parabolic 1Dconfining potential, r is the coordinate of the minimumof the potential energy and r i is the displacement of the i th particle from r (for the parabolic circular potential).Note that in case of a circular channel, r = r ch , where r ch is the radius of the channel.The second term in the r.h.s. of Eq. (2) representsthe interaction potential between the particles. For thescreened Couloumb potential, V int ( ~r ij ) = q ǫ e −| ~r i − ~r j | /λ D | ~r i − ~r j | , (6)where q is the charge of each particle, ǫ is the dieletricconstant of the medium, r ij = | ~r i − ~r j | is the distance be-tween i th and j th particles, and λ D is the Debye screen-ing length. Substituting (6) into Eq. (2), we obtain thepotential energy of the system H Y : H Y = N X i =1 V c ( ~r i ) + q ǫ N X i>j =1 e −| ~r i − ~r j | /λ D | ~r i − ~r j | . (7)In order to reveal important parameters which character-ize the system, we rewrite the energy H Y in a dimension-less ( H ′ Y ) form by making use of the following variabletransformations: H Y = ( q /ǫa ) H ′ Y , r = r ′ a , where a is the mean inter-particle distance. The energy of thesystem then becomes H ′ Y = N X i =1 V ′ c ( ~r ′ i ) + N X i>j =1 e − κ | ~r ′ i − ~r ′ j | | ~r ′ i − ~r ′ j | , (8)where κ = a /λ D is the screening parameter of the in-teraction potential. In our simulations in Sec. III, weuse a typical value of κ = 1 . λ D = 10 − m.The hard-wall confinement potential is written as V ′ c ( ~r ′ i ) = (cid:26) | y ′ i | ≤ R ′ w / ∞ for | y ′ i | > R ′ w / , (9)where R ′ w is scaled by the inter-particle distance a . Wealso introduce a dimensionless parameter χ = m ( ω a ) k B T , (10)which is a measure of the strenght of the parabolic 1Dconfinement potential.For colloidal particles moving in a nonmagnetic liq-uid, their motion is overdamped and thus the stochasticLangevin equations of motion can be reduced to thosefor Brownian particles [31]: d~r i dt = D i k B T h − X j = i ~ ∇ i V int ( ~r ij ) − ~ ∇ i V c ( ~r i ) + ~F iT ( t ) i . (11)Note, however, that in Sec. IV we will deal with massivemetallic balls and therefore we will keep the inertial termin the Langevin equations of motion. In Eq. (11), ~r i , D i and m i are the position, the self-diffusion coefficient (measured in m /s) and the mass (inkg) of the i th particle, respectively, t is the time (in sec-onds), k B is the Boltzmann constant, and T is the abso-lute temperature of the system. Finally, ~F iT is a randomlyfluctuating force, which obeys the following conditions: h ~F T i = 0 and h F iT ( t ) F i ′ T ( t ′ ) i = 2 ηk B T δ ii ′ δ ( t − t ′ ), where η is the friction coefficient. Eq. (11) can be written indimensionless form as follows: d~r ′ i dt ′ = D ′ i Γ h − X j = i ~ ∇ ′ i V ′ int ( ~r ′ ij ) − ~ ∇ ′ i V ′ conf ( ~r ′ i ) + ~F ′ iT ( t ′ ) i , (12)where we use the following transformation V int =( q /ǫa ) V ′ int , D ′ i = D i /a , and introduced a coupling pa-rameter Γ, which is the ratio of the average potentialenergy to the average kinetic energy, Γ = h V i / h K i , suchthat Γ = q /k B T ǫa . The time t ′ is expressed in secondsand distances are expressed in units of the interparticledistance a . In what follows, we will abandon the prime( ′ ) notation. We have used a first order finite differencemethod (Euler method) to integrate Eq. (12) numerically.In the case of a straight channel, periodic boundary con-ditions (PBC) were applied in the x -direction while inthe y -direction the system is confined either by a hard-wall or by a parabolic potential. Also, we use a timestep∆ t = 0 .
001 and the coupling parameter is set to Γ = 10.For a circular channel, we use polar coordinates ( r, φ ) andmodel a 2D narrow channel of radius r ch with parabolicpotential-energy profile across the channel, i.e., in the r -direction. III. 1D VERSUS 2D DIFFUSION IN ASTRAIGHT CHANNELA. Mean-square displacement (MSD) calculations
In order to characterize the diffusion of the system, wecalculate the MSD as follows: h ∆ x ( t ) i = D N N X i =1 [ x i ( t + ∆ t ) − x i ( t )] E ∆ t , (13)where N is the total number of particles and h ... i ∆ t rep-resents a time average over the time interval ∆ t . Notethat in the general case (e.g., for small circular channelswith the number of particles N = 20 — see Sec. IV)the calculated MSD was averaged over time and over thenumber of ensembles [32]. However, we found that forlarge N (i.e., several hundred) the calculated MSD forvarious ensemble realizations coincide (with a maximumdeviation within the thickness of the line representing theMSD).To keep the inter-particle distance approximately equalto unity, we defined the total number of particles N fora 1D and Q1D system as N = L p − R w ; R w < , (14)where L is the size of the simulation box (in dimensionlessunits) in the x -direction. In our simulations for a straightchannel geometry, we typically used N = 400 −
900 par-ticles. We study the system for two different types ofconfinement potential: (i) a parabolic 1D potential inthe y -direction, which can be tuned by the confinementstrength χ and (ii) a hard-wall potential, where particlesare confined by two parallel walls separated by a distance R w . The results of calculations of the MSD as a functionof time for different values of the confinement strength χ [Eq. (10)] and the width of the channel R w are presentedin Fig.1(a)-(c) and Fig.2(a)-(c), respectively.Initially, in both cases (i.e., a parabolic and a hard-wall confinement potential), the system exhibits a short-time normal diffusion behavior, where h ∆ x ( t ) i ∝ t . .This is the typical initial “free-particle” diffusion regime.After this initial regime, there is an intermediate sub-diffusive regime (ITR). As discussed in Ref. [36], theITR shows an apparent power-law behavior [37], where0 . < α < .
0, and it was also found previously in dif-ferent diffusion models [38, 39]. In the ITR, we founda SFD regime for either a channel with strong parabolicconfinement [ χ = 3 . R w = 0 .
20 (Fig. 2(a))]. This is due to the factthat for large (small) values of χ ( R w ), the confinementprevents particles from passing each other. The resultsfor α in the ITR are shown as a function of χ and R w inFig. 1(d) and Fig. 2(d), respectively. As can be seen inFig. 1(d) [Fig. 2(d)], α increases with decreasing χ [withincreasing R w ] and thus the SFD condition turns out tobe broken. The values of α presented in these figures cor-respond to the minimum of the effective time dependentexponent α ( t ). Following Ref. [40], α ( t ) is calculated us-ing the “double logarithmic time derivative” α ( t ) = d log h ∆ x ( t ) i d log t , (15)and the results are shown in Fig. 3.The different diffusive regimes, i.e. normal diffusionregime ( α = 1 .
0) and SFD ( α = 0 . α -dependence on both the confinement parameters (i.e., α ( χ ) and α ( R w )) presents a different qualitative behav-ior, namely, the SFD regime is reached after a smoothercrossover in the parabolic confinement case as comparedto the hardwall case. A similar smoother crossover isalso found in the case of a circular channel with parabolicconfinement in the radial direction. A more detailed dis-cussion on these two different types of the behavior of α will be provided in Sec. IV. < (cid:1) x ( t ) > < (cid:1) x ( t ) > < (cid:1) x ( t ) > FIG. 1: (Color online) (a)-(c) Log-log plot of the mean-squaredisplacement (MSD) h ∆ x ( t ) i as a function of time for differ-ent values of χ . Different diffusion regimes can be distin-guished: normal diffusion regime ( α = 1 .
0) and intermediatesubdiffusive regime (ITR, α < . χ = 1.5, there is a normal diffusion regime (i.e. α = 1 .
0) afterthe ITR. The dashed and solid lines in (a)-(c) are a guide tothe eye. Panel (d) shows the dependence of the slope ( α ) ofthe MSD curves (in the ITR, characterized by an apparentpower-law; h ∆ x ( t ) i ∝ t α ) on the confinement strength χ . < (cid:1) x ( t ) > < (cid:1) x ( t ) > < (cid:1) x ( t ) > FIG. 2: (Color online) (a)-(c) Log-log plot of the mean-squaredisplacement (MSD) h ∆ x ( t ) i as a function of time for differ-ent values of R w . Different diffusion regimes can be distin-guished: normal diffusion regime ( α = 1 .
0) and intermediatesubdiffusive regime (ITR, α < . R w = 0.60, there is a normal diffusion regime (i.e. α = 1 . α )of the MSD curves (in the ITR, characterized by an apparentpower-law; h ∆ x ( t ) i ∝ t α ) on the confinement parameter R w . -2 -1 α ( t ) t (s) χ = 3.5 χ = 1.5 (a) -2 -1 α ( t ) t (s) R w = 0.20 R w = 0.80 (b) FIG. 3: (Color online) (a)-(b) Exponent α as a function oftime, calculated from Eq. (15) for different values of the con-finement parameters χ and R w , respectively. B. “Long-time” behavior of the MSD curves andcrossing events C ( t ) For small values of the parabolic confinement (e.g., χ = 1 . h ∆ x ( t ) i ∝ t . ; (ii) a subdiffusive regimewith h ∆ x ( t ) i ∝ t α , where 0 . < α < . h ∆ x ( t ) i ∝ t . . Note that the “long-time” termused here is not to be confused with the long-time usedfor infinite systems, as discussed in the Introduction.However, for large values of the parabolic confinement(e.g., χ = 3 . h ∆ x ( t ) i ∝ t . ) and (ii) a SFD regime (i.e., h ∆ x ( t ) i∝ t . ).One question that arises naturally is whether this nor-mal diffusion regime (i.e., h ∆ x ( t ) i ∝ t . for “long-times”) is an effect of the colletive motion of the system(center-of-mass motion) or an effect of the single-particlejumping process, since the confinement potential χ = 1 . C ( t ) as a func-tion of time and results are shown in Fig. 4(a). We foundthat for small values of the confinement potential (e.g., χ = 1 .
5) the number of crossing events grows linearly intime, i.e., C ( t ) ∝ ω c t , where ω c is the rate of crossingevents. On the other hand, a strong confinement poten-tial (e.g., χ = 3 .
5) prevents particles from bypassing, andthus C ( t ) = 0 during the whole simulation time.Therefore, the “long-time” normal diffusive behavior(i.e., h ∆ x ( t ) i ∝ t . for “long-times”) found in our sim-ulations for the case where the SF (single-file) condition isbroken (e.g., χ = 1 .
5) is not due to a collective (center-of-mass) diffusion. Instead, this normal diffusive behaviouris due to a single-particle jumping process, which hap-pens with a constant rate ω c > χ = 1 .
5) and ω c = 0 (for χ = 3 . N = 400 − N = 80 − h ∆ x ( t ) i ∝ t . regime is recovered in the “long-time”limit. In Sec. V, we will further discuss the long-timelimit using a model of discrete sites.As we demonstrated above, the transition from pure1D diffusion (SFD) characterized by α = 0 . α > .
5) could be either more“smooth” (as in Fig. 1(d), for a parabolic confinement)or more “abrupt” (as in Fig. 2(d), for a hard-wall con-finement). One can intuitively expect that this differencein behavior can manifest itself also in the crossing eventsrate ω c , i.e., that ω c as a function of χ (or R w ) shoulddisplay a clear signature of either “smooth” or “abrupt”behavior.However, the link between the two quantities, i.e.,the exponent, α ( χ/R w ), and the crossing events rate, ω c ( χ/R w ) is not that straightforward. To understandthis, let us refer to the long-time limit (which will be ad-dressed in detail within the discrete-site model in Sec. V).As we show, in the long-time limit the exponent α is de-fined by one of the two conditions: ω c = 0 (then α = 0 . ω c > α = 1) and it does not depend on thespecific value of ω c provided it is nonzero. Therefore, inthe long-time limit the transition between 1D to 2D be-havior is not sensitive to the particular behavior of thefunction ω c ( χ/R w ).Although for “intermediate” times (considered in thissection) the condition ω c = 0 or ω c > ω c ( χ/R w ) strongly influences the behavior of the expo-nent α ( χ/R w ). This is illustrated in Figs. 4(b, c). InFig. 4(b), the function ω c ( χ ) gradually decreases from × × × × t (s) C ( t ) χ = 1.5 χ = 3.5 C(t) ~ ω c t x 10 (c) ω c ( s - ) χ (b) R w - d ω c / d R w (c) - d ω c / d χ χ R w (a) FIG. 4: (Color online) (a) Number of crossing events C ( t ) asa function of time for N = 400 particles, for χ = 1 . χ = 3 . C ( t ). Panels (b) and (c) show therate of the crossing events ω c as a function of the confinementpotential parameters ( χ and R w ). The insets in the panels (b)and (c) show the derivatives, dω c ( χ ) /dχ and dω c ( R w ) /dR w ,correspondingly. χ varying in a broad interval from 1.5 to3 (note that the segment of ω c ( χ ) for 2 . < χ < dω c ( χ ) /dχ ). Correspondingly, thetransition from α = 0 . α ≈ . χ is “smooth” (see Fig. 1(d)). On the other hand, thefunction ω c ( R w ) shown in Fig. 4(c) mainly changes (notethe change of the slope dω c ( R w ) /dR w shown in the insetof Fig. 4(c)) in a narrow interval 0 . < R w < .
6. Re-spectively, the transition for the function α ( R w ) occursin the narrow interval 0 . < R w < . C. Distribution of particles along the y -direction For the ideal 1D case, particles are located on a straightline. Increasing the width R w of the confining channelwill lead to a zig-zag transition [20, 41]. This zig-zag FIG. 5: For the hard-wall confinement case, we show typicaltrajectories of particles (i.e. 10 MD simulation steps) con-fined by the channel of width (a) R w = 0 .
20, (b) R w = 0 . R w = 0 . configuration can be seen as a distorted triangular con-figuration in this transition zone. Further increase of R w brings the system into the 2D regime, where the normaldiffusion behavior is recovered (see Fig. 5).For the parabolic 1D confinement, we can see[Fig. 6(a)] that the distribution of particles P ( y ) alongthe channel is symmetric along the axis y = 0. Also,for large values of χ (e.g., χ = 3 .
5) particles are con-fined in the y -direction and thus can move only in the x -direction, forming a single-chain structure. As the con-finement decreases ( χ → P ( y ) broadens resulting in the crossover from theSFD regime ( χ = 3 .
5) to the 2D normal diffusion regime( χ = 0 . χ (e.g., χ = 0 . P ( y ) in Fig. 6(a)), thus allowing par-ticles to pass each other. IV. DIFFUSION IN A CIRCULAR CHANNEL
In the previous section, we analyzed the transition(crossover) from the SFD regime to 2D diffusion in nar-row channels of increasing width. The analysis was per-formed for a straight channel with either hard-wall orparabolic confinement potential. However, in terms ofpossible experimental verification of the studied effect,one faces an obvious limitation of this model: althougheasy in simulation, it is hard to experimentally fulfill theperiodic boundary conditions at the ends of an open chan-nel. Therefore, in order to avoid this difficulty, in SFDexperiments [18, 25] circular channels were used.In this section, we investigate the transition (crossover)from SFD to 2D-diffusion in a system of interacting par- -0.4 -0.2 0 0.2 0.4 y P ( y ) R w = 0.80 -0.3 -0.15 0 0.15 0.3 y P ( y ) R w = 0.60 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 y P ( y ) R w = 0.50 -0.1 -0.05 0 0.05 0.1 y P ( y ) R w = 0.20 -0.5 0 0.5 y P ( y ) χ = 3.5χ = 3.0χ = 2.5χ = 2.0χ = 1.5χ = 1.0χ = 0.5 (a) (b) FIG. 6: (Color online) Probability distribution of the particledensity P ( y ) along the y -direction are shown for (a) differentvalues of χ (parabolic 1D confinement) and (b) four differentvalues of the width R w of the channel (hard-wall confine-ment). ticles diffusing in a channel of circular shape. In par-ticular, we will study the influence of the strength of theconfinement (i.e., the depth of the potential profile acrossthe channel) on the diffusive behavior. Without loss ofgenerality, we will adhere to the specific conditions andparameters of the experimental set-up used in Ref. [25].An additional advantage of this model is that the motionof the system of charged metallic balls [25] is not over-damped, and we will solve the full Langevin equations ofmotion to study the diffusive behavior of the system.We consider N particles, interacting through a Yukawapotential [Eq. (6)], which are embedded in a ring channelof radius r ch . We define a parabolic confinement poten-tial across the channel in the form (5) where parameter β is chosen as follows: β = V γr , V = q ǫ X i = j exp h − κr ch sin (cid:16) φ i − φ j (cid:17)i r ch sin (cid:16) φ i − φ j (cid:17) , (16)when all the particles are equidistantly distributed alongthe bottom of the circular channel. It should be notedthat in this case, V is approximately equal to V gs dueto the weak Yukawa interaction, which slightly shifts theparticles away from the bottom of the channel. Such achoice of V is related to the fact that we study the influ-ence of the confinement on the diffusion and, therefore,the potential energy of the particles must be of the or-der of the inter-particle interaction energy. Parameter r characterizes the distance where the external poten-tial reaches the value V /γ , and V gs is the energy of theground state of the system of N particles as defined byEq. (7). Parameter γ plays the role of a control param-eter. By changing γ we can manipulate the strength ofthe confinement and, therefore, control the fulfillment ofthe single-file condition. Increase in γ corresponds to adecrease in the depth of the confinement (5) which leadsto the expansion of the area of radial localization of par-ticles. Therefore, an increase of γ results in a similareffect (i.e., spatial delocalization of particles) as an in-crease of temperature, i.e., parameter γ can be consideredas an “effective temperature”. Note that such a choiceof the parameter that controls the confinement strengthis rather realistic. In the experiment of Ref. [25] withmetallic balls, the parabolic confinement was created byan external electric field, and the depth of the potentialwas controlled by tuning the strength of the field.To study diffusion of charged metallic balls, we solvethe Langevin equation of motion in the general form (i.e.,with the inertial term ∝ m ), m d ~r i dt = − η d~r i dt − X j,i = j ~ ∇ V int ( ~r ij ) − ~ ∇ V c ( ~r i ) + ~F iT , (17)where m = 2 . × − kg [25] is the mass of a particle, η isthe friction coefficient (inverse to the mobility). Here allthe parameters of the system were chosen following theexperiment [25], and λ D = 4 . × − m , Γ = 1 (which isa typical experimental value, see, e.g., also [18]). Corre-spondingly, mass is measured in kg, length in m, and timein seconds. Also, following Ref. [25], we took a channelof radius r ch = 9 mm (in the experiment [25], the exter-nal radius of the channel was 10 mm, and the channelwidth 2 mm; note that in our model we do not define thechannel width: the motion of a particle in the transversedirection is only restricted by the parabolic confinementpotential). We also took experimentally relevant num-ber of diffusing particles, N , varying from N = 12 to N = 40 (in the experiment [25], the ring channel con-tained N = 12 or N = 16 diffusing balls).Fig. 7 shows the results of calculations of the trajec-tories of N = 20 particles diffusing in a ring of radius (a) (b) (c)(d) (e) (f) FIG. 7: Trajectories of N = 20 particles diffusing in a ring ofradius r ch = 9 mm for 10 consequent time steps for differentvalues of γ . γ =1 (a), 2 (b), 3 (c), 5 (d), 7 (e), 9 (f). r ch = 9 mm for the first 10 MD steps for various valuesof the parameter γ . As can be seen from the presentedsnapshots, the radial localization of particles weakenswith increasing γ . At a certain value of γ this leads tothe breakdown of the single-file behavior (Figs. 7(c)-(f)). A. Breakdown of SFD
It is convenient to introduce the distribution of theprobability density of particles in the channel P rad alongthe radial direction r . In order to calculate the function P rad ( r ) we divided the circular channel in a number ofcoaxial thin rings. The ratio of the number of obser-vations of particles in a sector of radius r i to the totalnumber of observations during the simulation is definedas the probability density P rad ( r i ). In Fig. 8, the proba-bility density P rad ( r ) is presented for different values of γ . With increasing γ , the distribution of the probabilitydensity P rad ( r ) broaden and the maximum of the func-tion P rad ( r ) shifts away from the center of the channel(see Fig. 8). The latter is explained by the softeningof the localization of particles with increasing γ , whichtend to occupy an area with a larger radius due to therepulsive inter-particle interaction. Simultaneously, thedistribution of the probability density P rad ( r ) acquiresan additional bump indicating the nucleation of a two-channel particle distribution [42]. The observed broaden-ing and deformation of the function P rad ( r ) is indicativeof a gradual increase of the probability of mutual bypassof particles (i.e., the violation of the SF (single-file) con-dition, also called the “overtake probability” [43]) withincreasing γ .Let us now discuss a qualitative criterion for the break-down of SFD, i.e., when the majority of particles leavethe SFD mode. For this purpose, let us consider a parti-cle in the potential created by its close neighbor (which FIG. 8: (Color online) The distribution of the probabil-ity density of particles P rad ( r ) in a circular channel of radius r ch = 9 mm along the radial direction r . The different curvescorrespond to various γ . Increasing γ the width of the dis-tribution P rad ( r ) increases due to a weakening of the confine-ment.FIG. 9: (Color online) Spatial distribution of the poten-tial V int ( r, φ ) created by a particle (red (grey) circle) and thequalitative distribution of the probability density of particlesin circular channel P rad ( r ) (green (light grey) line) along theradial direction r . The function ∆ r determines an approx-imate radial distance between particles when the potentialbarrier U bar becomes “permeable” for given temperature T .The function ∆ r sw characterizes a width of the distribution P rad ( r ) at this temperature T . is justified in case of short-range Yukawa interparticle in-teraction and low density of particles in a channel) shownin Fig. 9. Different lines show the interparticle potential V int as a function of angle φ for different radii r . Forsmall values of γ , the center of the distribution P rad ( r )(see Fig. 7) almost coincides with the center of the chan-nel (i.e., with the minimum of the confinement potentialprofile) and the distribution P rad ( r ) is narrow. There- fore, mutual passage of particles is impossible, i.e., theSF condition is fulfilled. The asymmetric broadening ofthe function P rad ( r ) with increasing γ results in an in-creasing probability of mutual bypass of particles whichhave to overcome a barrier U bar (see Fig. 9). This be-comes possible when U bar . k B T . In other words, thethermal energy k B T determines some minimal width ∆ r between adjacent particles when the breakdown of theSF condition becomes possible.It is clear that “massive” violation of the SF condition(i.e., when the majority of particles bypass each other)occurs when the halfwidth ∆ r sw of the distribution ofthe probability density P rad ( r ) obeys the condition:∆ r sw & ∆ r. (18)The function ∆ r sw is defined by the ratio of the thermalenergy k B T to the external potential U conf ( r ) and is ofthe same order as f ∆ r : V γr · ( f ∆ r/ ≈ k B T. (19)Therefore the criterion (18) can be presented in the form:∆ r sw ≈ f ∆ r & ∆ r. (20)This qualitative analysis of the breakdown of the SFDregime clarifies the role of the width and the shape of thedistribution of the probability density influenced by theasymmetry of the circular channel. B. Diffusion regimes
The MSD h ∆ φ ( t ) i is calculated as a function of time t as: (cid:10) ∆ φ ( t ) (cid:11) = * N par N ens X i,j [∆ φ ij ( τ + t ) − ∆ φ ij ( t )] + t , (21)where N par is the total number of particles of an ensem-ble and N ens is the total number of ensembles. In ourcalculations, the number of ensembles was chosen 100 fora system consisting of 20 particles.The time dependence of the MSD for different valuesof γ is shown in Fig. 10(a)–(c). Initially the system ex-hibits normal diffusion, where h ∆ φ i ∝ t . . This regimeis followed by an intermediate subdiffusive regime, wherethe h ∆ φ i ∝ t α (0 . < α < . h ∆ φ i ∝ t . . As in thecase of straight channel geometry, this second crossover(i.e., from intermediate subdiffusion to “long-time” nor-mal diffusion) can also be due to two other reasons: (i)due to a collective (center-of-mass) diffusion or (ii) due toa single-particle jumping process. However, for the sim-ulations in the case of a circular geometry, the numberof particles is relatively small (taking the fact that this0is a finite-size system), and therefore, the crossover fromsublinear to linear regime is due to a collective (center-of-mass) diffusion. We further address this issue in Sec. V,where we consider a discrete site model and we excludethe center-of-mass motion.Fig. 10(d) shows α as a function of γ . The function α ( γ ) experiences a monotonic gradual crossover from the α = 0 . α . γ (Fig. 10) is related to the presence of, though weakbut nonzero, external confinement in the radial direction.This change of the diffusive behavior is explained by aweakening of the average radial localization of particleswith increase of γ (Fig. 8) and, as a consequence, by anincrease of the probability of mutual bypass of particles.The observed crossover between the 1D single-file and2D diffusive regimes, i.e., α ( γ )-dependence, shows a sig-nificant different qualitative behavior as compared to thecase of a hard-wall confinement potential considered inSec. III, where a rather sharp transition between the tworegimes was found [Fig. 2(d)]. The different behavioris due to the different confinement profiles and can beunderstood from the analysis of the distribution of theprobability density of particles for these two cases. In thecase of a hard-wall channel, the uncompensated (i.e., bythe confinement) interparticle repulsion leads to a higherparticle density near the boundaries rather than near thecenter of the channel (see Fig. 5 and Fig. 6(b)). As aconsequence, the breakdown of the SF condition — withincreasing width of the channel — happens simultane-ously for many particles in the vicinity of the boundaryresulting in a sharp transition (see Fig. 2(d)). On thecontrary, in the case of parabolic confinement, the den-sity distribution function has a maximum — sharp orbroad, depending on the confinement strength — nearthe center of the channel (see Figs. 6(a) and 8). Withincreasing the “width” of the channel (i.e., weakeningits strength), only a small fraction of particles undergoesthe breakdown of the SF condition. This fraction grad-ually increases with decreasing strength of the confine-ment, therefore resulting in a smooth crossover betweenthe two diffusion regimes. V. DISCRETE SITE MODEL: THE LONG-TIMELIMIT
The calculated MSD for different geometries and con-finement potentials allowed us to explain the evolution ofthe subdiffusive regime with varying width of the channel(or potential strength in case of a parabolic potential).However, the obtained results are only valid for the in-termediate regime and therefore they only describe the“onset” of the long-time behavior. The problem of ac-cessing the long-time behavior in a finite chain is relatedto the fact that sooner or later (i.e., depending on thechain length) the interacting system will evolve into acollective, or “single-particle”, diffusion mode which is
FIG. 10: (Color online) (a)–(c): Log-log plot of the mean-squared displacement (MSD) h ∆ φ i as a function of time fordifferent values of the “effective” temperature γ = (a) 1, (b)2, and (c) 3. Here ( N ens = 100 , N par = 20). (d) The diffu-sion exponent α as a function of γ . Increase of the “effec-tive” temperature γ leads to the gradual transformation ofthe single-file regime of diffusion into the diffusion regime offree particles. α = 1 .
0. Thus the question is whetherthe observed behavior holds for the long-time limit, i.e.,is the transition from t . to t . behavior smooth?To answer this question, we considered a simple model,i.e., a linear discrete chain of fixed sites filled with eitherparticles or “holes” (i.e., sites not occupied by particles)(for details, see Ref. [38]; this model was also recentlyused in Ref. [44]). The particles can move along the chainonly due to the exchange with adjacent vacancies (i.e.,with holes). Within this model, the long-time diffusionbehavior was described analytically for an infinite linearchain as well as for a finite cyclic chain [38]. In partic-ular, this model predicts that: (i) If the chain is infinitethen the long-time power law of the diffusion curve α is0 . h ∆ x ( t ) i ∝ t . ); (ii) If the chain is finitethen the subdiffusive regime with α = 0 . α = 1 . α = 0 regime, i.e., the regime of satu-ration (if no cyclic boundary condition is imposed [14]).The latter regime is reached for times longer than the“diffusion time” of a “hole” along the whole chain t chain .Let us now apply this model to a finite-size chain ofparticles. For this purpose, we assume that adjacent par-ticles are able to exchange their positions with some prob-ability P at every time step. For example, probability P = 0 . N s = 150 sitesand N h = 1 hole. Averaging was done over 1000 ensem-bles. The calculation was performed for the following val-ues of the probability: P = 0 , − , − , − , . , . h ∆ x ( t ) i ∝ t . and ∝ t . . The characteristic time t chain shifts towardslower values with increasing P . However this analysis(Fig. 11(a)) does not allow to distinguish the contri-butions to the long-time behavior ( ∝ t . ) due to: (i)the breakdown of single-file condition (i.e., diffusion dueto particle exchanges), and (ii) the “collective” diffusion(chain “rotation”).To overcome this difficulty, we exclude the “collective”diffusion of the system and introduce a modified MSD h ∆ x ( t ) i corr (which is so-called “roughness” of the systemof particles, as discussed in Ref. [28]) as follows: h ∆ x i corr = h ( x − ¯ x ) i , where h ... i is the average over time; ¯ x is the average of anensemble of particles at a given time, or “collective” co-ordinate. It should be noted that h x i 6 = ¯ x . If the systemdoes not experience “collective” diffusion then ¯ x ( t ) = 0and the modified MSD coincides with the conventionalone: h ∆ x i corr = h x i . < x > D < x > D < x > D FIG. 11: (Color online) Log-log plot of the MSD h ∆ x ( t ) i (a)and corrected MSD h ∆ x ( t ) i corr (b) as a function of time fordifferent values of the probability P of bypassing. Averagingwas done over N sim = 1000 ensembles. The diffusion curves calculated by using the modifiedMSD are presented in Fig. 11(b). For P = 0, the diffu-sion curve (shown by black open squares) after the sub-diffusive regime reaches saturation (i.e., h ∆ x ( t ) i corr =const). The observed behavior is similar to that of afinite linear chain with fixed ends (see Ref. [14]). For P = 0, all the diffusion curves in the long-time limit arecharacterized by α = 1 . independent of the value of theprobability P , as seen in Fig. 11(b). In other words, thelong-time diffusion does not depends on the probabilityof mutual exchanges of particles and has the same long-time behavior for any probability P = 0. Here we wouldlike to emphasize again that the long-time behavior ofthe diffusion curves is free from the “collective” diffusioneffect and is only determined by particle jump diffusion.Increasing a number of sites in the model corresponds,in fact, approaching to the model of infinite chain. Wehave found that the increasing a number of sites leadsto growth of the h ∆ x ( t ) i corr limit of saturation, on theone hand, and to a shift of t chain to larger t , on the otherhand. Hence, extrapolating our results to the case of infi-2nite chain, we can conclude that in this case as well as inthe case of finite-size chain, the breakdown of single-filecondition leads to an abrupt transition from subdiffusiveto the normal diffusion regime.The difference in the diffusive curves is just the time τ tran from subdiffusive regime to the normal regime: forlow P it ( τ tran ) is long enough while for high P it ( τ tran )is short. It is easy to see that τ tran ∼ /P (%). Thus, wecan conclude that in the long-time limit the transitionfrom t . to t . behavior is abrupt . Note that our cal-culations performed using the modified MSD h ∆ x ( t ) i corr reproduce the results of Ref. [14] for a closed “box”. Thisis explained by the fact that in the closed “box” geom-etry the center of mass (or collective) diffusion is zero,and it is natural that the roughness (see Ref. [28]) andthe particles diffusion coincide. VI. CONCLUSIONS
We have studied a monodisperse system of interact-ing particles subject to three types of confinement po-tentials: (i) a 1D hardwall potential, (ii) a 1D parabolicconfinement potential which both characterize a quasi -1D system, and (iii) a circular confining potential, whichmodels a finite size system. In order to study the diffusiveproperties of the system, we have calculated the mean-squared displacement (MSD) numerically through molec-ular dynamics (MD) simulations. For the case whereparticles diffuse in a straight line in a Q1D channel, dif-ferent diffusion regimes were found for different valuesof the parameters of the confining potential ( χ or R w ).We have found that the normal diffusion is suppressedif the channel width R w is between 0 .
20 and 0 .
50 (or by2 . < χ < .
5, for the case of parabolic 1D confinement),leading the system to a SFD regime for intermediate timescales. For values of R w > .
56, particles will be able tocross each other and the SFD regime will be no longerpresent.The case of a circular channel corresponds to, e.g., the set-up used in experiments with sub-millimetric metallicmassive balls diffusing in a ring with a parabolic potentialprofile created by an external electric field. The strengthof the potential (which determines the effective “width”of the channel) can be tuned by the field strength. Con-trary to the case of hard-wall confinement, where thetransition (regarding the calculation of the scaling expo-nent ( α ) of the MSD h ∆ x ( t ) i∝ t α ) is sharp, a smoothcrossover between the 1D single-file and the 2D diffusiveregimes was observed. This behavior is explained by dif-ferent profiles for the distribution of the particle densityfor the hard-wall and parabolic confinement profiles. Inthe former case, the particle density reaches its maximumnear the boundaries of the channel resulting in a massivebreakdown of the SF condition and thus in a sharp transi-tion between the different diffusive regimes. In the lattercase, on the contrary, the density distribution functionhas a maximum near the center which broadens with de-creasing strength of the confinement. This results in asmooth crossover between the two diffusion regimes, i.e.,SFD and 2D regime. The analysis of the crossing events,i.e., the rate of the crossing events ω c as a function ofthe confinement parameter χ or R w , supports these re-sults: the function ω c ( χ/R w ) displays a clear signatureof either “smooth” or “abrupt” behavior.We also addressed the case of a finite discrete chainof diffusing particles. It was shown that in this casethe breakdown of the single-file condition (i.e., when theprobability P of particles bypassing each other is non-zero) leads to an abrupt transition from a subdiffusiveregime to the normal diffusion regime. Acknowledgments
This work was supported by CNPq, FUNCAP (Pronexgrant), the “Odysseus” program of the Flemish Govern-ment, the Flemish Science Foundation (FWO-Vl), thebilateral program between Flanders and Brazil, and thecollaborative program CNPq - FWO-Vl. [1] C. Cottin-Bizonne, J.-L. Barrat, L. Bocquet, and E.Charlaix, Nature Mater. 2 (2003).[2] L. A. Hodgkin and D. R. Keynes, J. Physiol. , 61(1955).[3] J¨org K¨arger, Phys. Rev. A , 4173 (1992).[4] F. Martin, R. Walczak, A. Boiarski, M. Cohen, T. West,C. Cosentino, J. Shapiro, and M. Ferrari, J. Control. Re-lease , 123 (2005).[5] T. E. Harris, J. Appl. Probab. , 323 (1965).[6] P. Barrozo, A. A. Moreira, J. A. Aguiar, and J. S. An-drade Jr., Phys. Rev. B , 104513 (2009).[7] J. S. Andrade, Jr., G. F. T. da Silva, A. A. Moreira, F.D. Nobre, and E. M. F. Curado, Phys. Rev. Lett. ,260601 (2010).[8] G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature(London) , 188 (2001). [9] M. W. Meier and H. D. Olsen (Editors), Atlas of Zeo-lite Structure Types (Butterworths-Heinemann, London1992).[10] T. Halpin-Healy and Y. C. Zhang, Phys. Rep. , 215(1995).[11] L. J. Lebowitz and K. J. Percus, Phys. Rev. E , 122(1967).[12] D. G. Levitt, Phys. Rev. A , 3050 (1973).[13] P. M. Richards, Phys. Rev. B , 1393 (1977).[14] L. Lizana and T. Ambj¨ornsson, Phys. Rev. E , 051103(2009).[15] V. N. Kharkyanen and Semen O. Yesylevskyy, Phys. Rev.E , 031118 (2009).[16] E. Barkai and R. Silbey, Phys. Rev. Lett. , 050602(2009).[17] A. M. Alsayed, M. F. Islam, J. Zhang, P. J. Collings, and A. D. Yodh, Science , 1207 (2005).[18] H.-Q. Wei, C. Bechinger, and P. Leiderer, Science ,625 (2000).[19] C. Lutz, M. Kollmann, and C. Bechinger, Phys. Rev.Lett. , 026001 (2004).[20] G. Piacente, F. M. Peeters, and J. J. Betouras, Phys.Rev. E , 036406 (2004).[21] W. P. Ferreira, J. C. N. Carvalho, P. W. S. Oliveira, G.A. Farias, and F. M. Peeters, Phys. Rev. B , 014112(2008).[22] W. Yang, K. Nelissen, M. Kong, Z. Zeng, and F. M.Peeters, Phys. Rev. E , 041406 (2009).[23] K. Nelissen, V. R. Misko, and F. M. Peeters, Europhys.Lett. , 56004 (2007).[24] For underdamped systems, the initial fast growth of theMSD (i.e., ballistic regime) is characterized by α = 2.[25] G. Coupier, M. SaintJean, and C. Guthmann, Phys. Rev.E , 031112 (2006).[26] H. L. Tepper, J. P. Hoogenboom, N. F. A. van der Vegt,and W. J. Briels, J. Chem. Phys. , 11511 (1999).[27] J. B. Delfau, C. Coste, and M. Saint-Jean, Phys. Rev. E , 011101 (2011).[28] P. M. Centres and S. Bustingorry, Phys. Rev. E ,061101 (2010).[29] D. V. Tkachenko, V. R. Misko, and F. M. Peeters, Phys.Rev. E , 051102 (2010).[30] B. Lin, M. Meron, B. Cui, S. A. Rice, and H. Diamant,Phys. Rev. Lett. , 216001 (2005).[31] D. L. Ermak and J. A. McCammon, J. Chem. Phys. ,1352 (1978).[32] As recently shown when studying single-trajectory av-erages [33], time average (TA) may differ from ensemble average (EA) MSD not only for nonergodic processes (forexample, for anomalous diffusion described by continuoustime random walks (CTRWs) models [34, 35]), but alsofor some ergodic processes in small complex systems.[33] J-H. Jeon and R. Metzler, Phys. Rev. E , 021147(2012).[34] Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev.Lett. , 058101 (2008).[35] A. Lubelski, I. M. Sokolov, and J. Klafter, Phys. Rev.Lett. , 250602 (2008).[36] R. Kutner, H. van Beijeren, and K. W. Kehr, Phys. Rev.B , 4382 (1984).[37] This apparent power-law behavior is characterized byMSD h ∆ x ( t ) i∝ t α , and it is an intermediate phenom-ena due to the interplay between the crossover from theMSD h ∆ x ( t ) i∝ t . regime to h ∆ x ( t ) i∝ t . regime.For details, see Ref. [36].[38] H. van Beijeren, K. W. Kehr, and R. Kutner, Phys. Rev.B , 5711 (1983).[39] B. J. Alder and W. E. Alley, J. Stat. Phys. , 341 (1978).[40] Y-L. Chou, M. Pleimling, and R. K. P. Zia, Phys. Rev.E , 061602 (2009).[41] G. Piacente, G. Q. Hai, and F. M. Peeters, Phys. Rev. B , 024108 (2010).[42] Note that this is not a zig-zag transition: the function P rad ( r ) still has a single maximum which is shifted fromthe center of the channel.[43] T. Ambj¨ornsson and R. J. Silbey, J. Chem. Phys. ,165103 (2008).[44] S. Savel’ev, F. Marchesoni, A. Taloni, and F. Nori, Phys.Rev. E74