Exact closure and solution for spatial correlations in single-file diffusion.

In single-file transport particles diffuse in narrow channels while not overtaking each other. it is a fundamental model for the tracer subdiffusion observed in confined systems, such as zeolites or carbon nanotubes. This anomalous behavior originates from strong bath-tracer correlations in one dimension. Despite extensive effort, these remained elusive, because they involve an infinite hierarchy of equations. For the symmetric exclusion process, a paradigmatic model of single-file diffusion, we break the hierarchy to unveil and solve a closed exact equation satisfied by these correlations. Beyond quantifying the correlations, the role of this key equation as a tool for interacting particle systems is further demonstrated by its application to out-of-equilibrium situations, other observables, and other representative single-file systems.

INTRODUCTION

Single-file transport, where particles diffuse in narrow channels with the constraint that they cannot bypass each other, is a fundamental model (, ) for tracer subdiffusion in confined systems. The very fact that the initial order is maintained at all times leads to the subdiffusive behavior of the position of a tracer particle (TP) (), in constrast with the regular diffusion scaling . This theoretical prediction has been experimentally observed by microrheology in zeolites, transport of confined colloidal particles, or dipolar spheres in circular channels (–).

The symmetric exclusion process (SEP) is an essential model of single-file diffusion. Particles, present at a density ρ, perform symmetric continuous-time random walks on a one-dimensional infinite lattice with unit jump rate and with the hard-core constraint that there is at most one particle per site (Fig. 1A). The SEP has become a paradigmatic model of statistical physics, and it has generated a huge number of works in the mathematical and physical literature [see, e.g., (, , , )]. A major recent advance has been achieved with the calculation of the large deviation function of the position X of a tracer in the long time limit (, ). It gives access to all the long-time cumulants of X, which are found to behave anomalously as (, ). Similarly, the cumulants of the time integrated current through the origin Q have been shown to also scale as , and the large deviation function has been determined ().

Fig. 1.

Models and bath-tracer correlations.

(A) The SEP. The tracer is located at position X, and Q is the integrated current of particles through the origin up to time t (i.e., the total flux of particles between sites 0 and 1). The occupation numbers of the sites are denoted η. (B) Here, we put forward bath-tracer correlations <ηeλ > / < eλ> (or bath-current correlations in the case of the integrated current Q) as fundamental quantities to analyze single-file diffusion, as they satisfy the simple closed integral equation (Eq. 6). Besides being key technical tools, (α) they characterize the bath-tracer coupling (see Eq. 4), (β) they quantify the response of the bath of particles to the perturbation induced by the tracer [see item (iv) after Eq. 13], and (γ), in turn. they control the large deviations of the subdiffusive motion of the tracer (see Eqs. 1, 2 and 13). Other representative models of single-file systems: (C) the random average process (RAP)—at exponential times, particles on a line can jump in either direction to a random fraction of the distance to the next particle; (D) the Kipnis-Marchioro-Presutti (KMP) model—each site carries a continuous energy variable. At exponential times, the total energy of two neighboring sites is randomly distributed between them.

This collection of anomalous behaviors in the SEP originates from the strong spatial correlations in the single-file geometry, which makes them determining quantities. Even if this has been recognized qualitatively for long and that the case of dense and dilute limits have been recently studied (), up to now, there is no quantitative determination of the bath-tracer correlations at arbitrary density, despite extensive effort. Although the SEP has been studied for more than 40 years, analytic formulas for these functions are still missing. The calculation of these correlations in the SEP actually constitutes an open many-body problem, which we solve here. More generally, we put forward bath-tracer correlations as fundamental quantities to analyze single-file diffusion, because we show that they satisfy a notably simple exact closed equation. The central role of this key equation as a novel tool for interacting particle systems is further demonstrated by showing that it applies to out-of-equilibrium situations, other observables, and other representative single-file systems.

We consider the SEP of average density ρ, with a tracer, of position X at time t, initially at the origin. The bath particles are described by the set of occupation numbers η(t) of each site r ∈ ℤ of the line at time t, with η(t) = 1 if the site is occupied and η(t) = 0 otherwise (see Fig. 1A). The statistics of the position of the tracer is described by the cumulant-generating function, whose expansion defines the cumulants κ of the positionIts evolution equation is given by (see eq. S20)where the generalized density profile (GDP)–generating function is defined byNote that, besides controlling the time evolution of the cumulant-generating function, w (together with ψ) completely characterizes the joint cumulant-generating function of (η, X) and thus the bath-tracer correlations (). The GDP-generating function is, therefore, a key quantity, and the next step consists in writing its evolution equation from the master equation describing the system. However, similarly to Eq. 2, it involves higher-order correlation functions. We are facing an infinite hierarchy of evolution equations, which is the rule for tracer diffusion (and for other observables such as the integrated current through the origin) in interacting particle systems (, –), and whose closure has remained elusive up to now. We provide below a closed equation that allows the determination of the GDP-generating function in the hydrodynamic limit (large time and large distance).

In this limit, the position of the tracer satisfies a large deviation principle (, , ), which implies that the cumulant-generating function scales as . This anomalous behavior originates from the more general scaling formof the GDP-generating function, where the coefficient Φ gives the large-scale limit of the joint cumulant of the tracer’s position X and the occupation number η measured in its frame of reference (Fig. 1B). In the following, we will drop the argument λ of Φ for convenience.

RESULTS

We report here (see Materials and Methods and section SII.A for details) that the two functions (rescaled derivatives of the profiles)are entirely determined by the closed Wiener-Hopf integral equations () with a Gaussian kernelwhere , and we have analytically continued Ω+ to v < 0 and Ω− to v > 0. The parameter ω is determined by the boundary conditions (see Eq. 5)so that the functions Ω±(v) are parametrized by . At this stage, the expression of has not been determined yet, but it can be obtained in the following way. First, Φ is deduced by integration of Ω±, withby definition, and the boundary conditionsThe resulting Φ are at this stage parametrized by and λ. Then, by using the large time limit of Eq. 2 can be written as a function of λ, and we finally obtain the desired GDP-generating function Φ(λ, v).

DISCUSSION

Several comments are in order. (i) We show in the Supplementary Materials (section SI.E) that the boundary condition Equation 9 is exact; furthermore, we argue below that the bulk Equation 6 and, thus, the obtained GDP-generating function Φ are also exact. (ii) The Wiener-Hopf equations (Eq. 6) can be solved explicitly in terms of the one-sided Fourier transforms ()where(iii) As a by-product, our approach yields the cumulant-generating function (or, equivalently, the large deviation function of the tracer’s position)which is shown in the Supplementary Materials (section SII.E) to be identical to the exact expression obtained using the arsenal of integrable probabilities in (, ). (iv) In addition, we obtain a full characterization of the spatial bath-tracer correlations and particularly analytical expressions of the Φ by using the procedure described above [see section SII.C for explicit expressions, which extend, at arbitrary density, the expressions given in () in the dilute (ρ → 0) and dense (ρ → 1) limits and Fig. 2 (A and B) for comparison with numerical simulations].

Fig. 2.

Results for the SEP.

GDPs of order: (A) n = 1 and n = 3 at density ρ = 0.5. Solid lines: numerical simulations (see Materials and Methods for details on the numerical procedure) at t = 3000. Dashed lines: theoretical predictions (eqs. S80 and S81). (B) n = 2 at densities ρ = 0.25, 0.5, and 0.75. Solid lines: numerical simulations at t = 3000. Dashed lines: theoretical predictions (eq. S80). Note that Φ2 exhibits a minimum at a distance v > 0 for ρ < 1/2, which disappears for ρ > 1/2. (C) GDP-generating function (or equivalently, conditional profiles) obtained from the theoretical prediction (Eq. 6) (dashed) compared to the numerical resolution of the MFT equations (solid line) described in section SIV.B. (D) Nonmonotony of the conditional profiles Φ(0+) in front of the tracer, as a function of the rescaled tracer’s position ξ, for different densities ρ.

(v) Our approach also provides the conditional profiles <η∣X = x> defined as the average of the occupation of the site X + r given that the tracer is at position x. In the hydrodynamic limit, where , λ* is defined by and Φ is the GDP-generating function determined above. While the (unconditional) profiles <η> are flat, the conditional profiles allow to probe the response of the bath of particles to the perturbation created by the displacement ξ of the tracer: in particular, for ξ > 0, it leads to an accumulation of bath particles in front of the tracer and a depletion behind (see Fig. 2C), quantified by the simple conservation relationwhich is a consequence of Equations 11 and 12. Another notable feature is the nonmonotony of the conditional profile Φ(λ*,0+) in front of the tracer as a function of the rescaled position of the tracer ξ(λ*) (see Fig. 2D). Unexpectedly, Φ(λ*,0+) does not saturate to 1 as ξ → + ∞ but, instead, returns to its unperturbed value ρ, which results from a global displacement of bath particles induced by the tracer.

Equation 6 describes several other situations of physical relevance. (i) First, it applies to the out-of-equilibrium situation of an initial step of density ρ+ for x > 0 and ρ− for x < 0, with the tracer initially at the origin. This paradigmatic setup has attracted a lot of attention (–, , ) because it remains transient at all times and never reaches a stationary state. The GDP-generating function Φ is obtained from the solution (Eq. 11) by following the procedure described above, upon only changing the boundary condition (Eq. 8) into Φ( ± ∞ ) = ρ±. Again, we recover the results of (, ) on the cumulant-generating function . In addition, we obtain the complete spatial structure of the bath-tracer correlations (see Fig. 3A and section SII.D.1 for explicit expressions). (ii) Second, and notably, it also gives access to the statistics of other observables, as exemplified by the case of the integrated current through the origin Q (see section SII.D.3 for the application to the generalized current, which is extra observable), defined as the total flux of particles between sites 0 and 1 during a time t. This quantity has been the focus of many studies, both in the context of statistical physics (, , , , ) and mesoscopic transport (–), in particular in the nonequilibrium situation ρ− ≠ ρ+ (, ). Note that while the statistics of tracer diffusion and integrated current are easily related in the case of quenched initial conditions (), the relation is more entangled for the annealed case considered here because of the fluctuations of the initial condition. The quantities introduced previously (Eqs. 1 to 3) on the example of tracer diffusion are naturally adapted by substituting Q for X. The corresponding profiles Φ are then obtained as a particular case of Equation 6 by setting ξ = 0, completed by modified boundary conditions (Eqs. 9 and 10) derived from the microscopic model (see section SII.D.2). In particular, the resulting Equation 13 gives back the exact cumulant-generating function of Q obtained in () by Bethe ansatz, because, in this case, we find thatwhich coincides with the single parameter involved in (, ). In addition, the Φ determined here provides the associated spatial structure (see Fig. 3B and section SII.D.2 for explicit expressions). These profiles have been introduced and studied numerically in () for an infinite system [see also () for a finite system between two reservoirs], but no analytical expressions were available until now. (iii) Last, beyond the SEP, it applies to other representative single-file systems of interacting particles with average density ρ [see also (), which is, however, limited to the calculation of the first-order Φ1]. Such systems can be described at large scale by two quantities: the diffusivity D(ρ) and the mobility σ(ρ) (). The case of the SEP considered above corresponds to D(ρ) = 1/2 and σ(ρ) = ρ(1 − ρ). Equation 6 (with adaptations of Eqs. 9 and 10 given in section SIII.C) more generally applies to single-file systems with D(ρ) = 1/2, σ′′(ρ) constant and σ(0) = 0, by replacing the −ω that multiplies the integral in Eq. 6 by ωσ′′(ρ)/2. Important cases covered by our approach include the following (see Fig. 1 for definitions and sections SIII.D and SIII.E for explicit expressions): (a) the model of hard Brownian particles [σ(ρ) = ρ] for which the GDP-generating function of () is recovered; (b) the Kipnis-Marchioro-Presutti (KMP) model () [σ(ρ) = ρ2; see Fig. 3C], which describes situations as varied as force fluctuations in packs of granular beads (), the formation of clouds and gels, self-assembly of molecules in organic and inorganic materials, and distribution of wealth in a society [see () and references therein]; (c) the random average process (RAP), which appears in a variety of problems such as force propagation in granular media, models of mass transport, or models of voting systems (, –). Although σ′′(ρ) is not constant in this case, the GDP-generating function can be deduced from our results thanks to a mapping between tracer diffusion in the RAP and the integrated current in the KMP model (see Fig. 3D) ().

Fig. 3.

Extensions.

(A) The out-of-equilibrium situation of an initial step density. (B) Another observable, the current through the origin. (C and D) Other representative single-file systems: the KMP model (C) and the RAP (D). The dashed lines correspond to the predictions obtained from the central equation (Eq. 6). (A) GDPs at order n = 1 and n = 2 for the SEP with a step initial density ρ− = 0.7 and ρ+ = 0.2 at t = 1500 with 2000 sites. (B) GDPs for the current in the SEP for a density ρ = 0.5 at t = 900, for orders n = 1 and n = 3 (the profile for n = 2 is zero). (C) KMP model [D(ρ) = 1/2 and σ(ρ) = ρ2] for ρ = 1. (C1 to C3) First three orders of the GDPs at t = 900 and 500 sites. (C4) GDP-generating function (D) RAP for ρ = 1 (for a uniform jump distribution, D(ρ) = 1/(4ρ) and σ(ρ) = 1/ρ). (D1 to D3) First three orders of the GDPs at t = 4000 with 5000 particles (solid lines). (D4) GDP-generating function compared to the numerical resolution of the MFT equations (solid line).

We lastly argue that the central equation (Eq. 6) is exact for the following reasons. (i) We show in the Supplementary Materials (section SIV) that the macroscopic fluctuations theory (MFT) () can be used to determine perturbatively the first coefficients Φ analytically. These coefficients computed up to order n = 5 (which is the highest order for which we managed to determine the integrals involved) coincide with those obtained by our approach. Furthermore, the agreement holds also nonperturbatively in λ, as displayed in Figs. 2C and 3 (C4 and D4), where the numerical solution of the MFT equations is compared to the analytical solution (, ). Moreover, and as mentioned above, the exact expression of the cumulant-generating functions of (ii) the tracer position of () and (iii) the integrated current of () are contained in our approach, including the case of an initial step of density.

Together, we have determined analytically the spatial correlations in the SEP, which allowed us to fully quantify the response of the bath to the perturbation induced by a tracer. Besides being paramount physical observables, these correlations have been shown to be fundamental technical quantities, because they satisfy a notably simple closed equation and control large deviations in single-file diffusion. This very same equation applies to a variety of situations involving single-file transport, which makes it a novel and promising tool to tackle interacting particle systems.

MATERIALS AND METHODS

Analytical calculations for the SEP

Details on analytical calculations are provided in the Supplementary Materials. We sketch here the main steps that led to the closed equation (Eq. 6) for the SEP. The starting point is a master equation describing the time evolution of the complete system (bath and tracer in the SEP), from which we obtain the time evolution of the cumulant-generating function ψ and the GDP-generating function (Eq. 3). The main difficulty is that the latter involves higher-order correlation functions.

The next step consists in using the scaling (Eq. 4) of the GDP-generating function and to derive the hydrodynamic limit of the problem (details given in section SI.E). The obtained bulk equation, valid at arbitrary density, is still not closed. We explain in the Supplementary Materials that a closed equation obeyed by Φ has to satisfy the following constraints: (i) It must reduce to the known equations obtained in the limits of high and low density in (); (ii) it should also reproduce, as a by-product, the cumulants of the tracer’s position derived recently in (, ); (iii) additional constraints concern the way the different parameters appear in the equation (see section SII.A for details); and (iv) lastly, the equation that we write should have a “proper scaling” with time determined in section SI.E.

Following these ideas and constraints, we obtain a first closed equation that holds at lowest orders in λ, see eq. S40, which properly reproduces the known cumulant κ for n ≤ 6. This equation is conveniently rewritten by introducing the new functions Ω±(v) defined by Equation 5. Extension of this equation to arbitrary order in λ and then its resummation yields the closed equations (S48 and S49) in the Supplementary Materials. Several technical steps detailed in section SII.A allow us to lastly transform them into the closed Wiener-Hopf integral equation (Eq. 6), which is our central result.

Numerical simulations for the SEP

The numerical simulations of the SEP are performed on a periodic ring of size N, with M = ρN particles at average density ρ. The particles are initially placed uniformly at random. The jumps of the particles are implemented as follows: One picks a particle uniformly at random, along with one direction (left and right with equal probabilities). If the chosen particle has no neighbor in that direction, then the jump is performed; otherwise, it is rejected. In both cases, the time of the simulation is incremented by a random number picked from an exponential distribution of rate N. We keep track of one particle (the tracer) and compute the moments of its displacement and the GDPs. The averaging is performed over 108 simulation.

Extensions

Analytical (section SIII) and numerical (section SV) extensions (other systems than the SEP, other observables, and nonequilibrium situations), following these lines, are described in the Supplementary Materials.