Single-file diffusion in spatially inhomogeneous systems

Transport properties in heterogeneous media can often be quantified on large length and timescales in terms of effective transport coefficients. For example, consider a Brownian particle diffusing in one dimension with a local diffusivity $D(x)$ subjected to an external potential $\phi(x)$, both of which vary periodically with the same period $\ell$. The Fokker-Planck equation (FPE) describing the evolution of the particle’s probability density function $p(x,t)$ is

$$\frac{\partial p(x,t)}{\partial t}=\hat{H}p(x,t),\qquad\hat{H}\equiv\frac{\partial}{\partial x}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}+\frac{\partial}{\partial x}\right)\right],$$

(1)

where $\beta=(k_{\mathrm{B}}T)^{-1}$ is the inverse temperature. The mean-squared displacement (MSD) of the particle at late times behaves as

$$\langle[X(t)-X(0)]^{2}\rangle\simeq 2D_{\mathrm{eff}}t,$$

(2)

where $D_{\mathrm{eff}}$ is called the effective (or, late-time) diffusion constant.

The first derivation of $D_{\mathrm{eff}}$ is attributed to Lifson and Jackson [1]. Given that Eq. (1) cannot be solved analytically in general, their argument was based on calculating the mean first-passage time to some large distance from the starting point. They argued that at large lengthscales and late times, the process is still diffusive, and can be described by a simple Brownian motion with an effective diffusion constant, $D_{\mathrm{eff}}$. This means that the particle obeys the effective FPE

$$\frac{\partial p(x,t)}{\partial t}=\hat{H}_{\mathrm{eff}}p(x,t),\qquad\hat{H}_{\mathrm{eff}}\equiv D_{\mathrm{eff}}\frac{\partial^{2}}{\partial x^{2}}.$$

(3)

Upon computing the mean first-passage time for Eq. (1), a comparison with that of Eq. (3) yields the Lifson-Jackson formula,

$$D_{\mathrm{eff}}=\frac{\ell^{2}}{\left[\int_{0}^{\ell}dxe^{\beta\phi(x)}/D(x)\right]\left[\int_{0}^{\ell}dxe^{-\beta\phi(x)}\right]}.$$

(4)

Various direct computations of $D_{\mathrm{eff}}$ [5, 4, 2, 3] from the definition of the MSD in Eq. (2) reassuringly agree with Eq. (4), as one might expect from physical intuition. The Lifson-Jackson formula has been verified in a number of experimental settings, in particular, and of importance for our current study, the dependence on variations in the potential which can be modulated in the context of optical trapping [6, 7, 8]. It can further be shown that the late-time asymptotic correction to Eq. (2) is a constant [9], and so the decay of the time-dependent diffusion coefficient is algebraic ($\sim 1/t$) rather than exponential. One should bear in mind that the mean first-passage time argument occasionally breaks down as in, e.g., the Sinai model [10] for diffusion in a Brownian random potential, where the tracer’s motion is subdiffusive [11].

At the same time, the problem of single-file diffusion (SFD) has received much attention from both a theoretical [15, 16, 12, 14, 17, 18, 19, 20, 21, 29, 22, 23, 24, 26, 27, 28, 13, 25] and, more recently, experimental point of view [30, 31, 32, 34, 33, 35]. SFD occurs in one-dimensional systems where particles cannot cross each other due to a hard-core repulsion term in their interactions. One way of treating this quantitatively is by ignoring the hard-core interaction, and relabelling the particles when they cross each other. The processes so generated are referred to as ‘elastically-colliding stochastic processes’ [29]. Within this picture, a tracer particle which was initially at the origin must have the same number of particles to its left at all times. If the density profile of the particles at time $t$ is given by $\rho(x,t)$, we get the constraint

$$\int_{-\infty}^{Y(t)}dx\rho(x,t)=\int_{-\infty}^{0}dx\rho(x,0),$$

(5)

where $Y(t)$ is the tracer’s (stochastic) position at time $t$. Using this observation, the MSD of the tracer particle can be computed directly at late times by appealing to the central limit theorem (CLT); a method we will use in this paper. Additional methods, useful for the case of pure Brownian motion (i.e., no interactions except for the hard-core repulsion), include the Bethe ansatz for the full joint probability density function [25], and the macroscopic fluctuation theory [27, 36, 37, 38].

There are two types of averages in play when treating SFD. The first average is over the thermal noise acting upon each of the individual particles, which we will denote by $\langle\cdots\rangle$. The second average is over the initial conditions, that is, the position at which each particle starts, which we will denote by $\overline{\cdots}$. Consequently, there are two possible definitions for the MSD [19] — the so-called annealed variance

	$$\langle Y^{2}(t)\rangle_{\mathrm{ac}}=\overline{\langle Y^{2}(t)\rangle}-\overline{\langle Y(t)\rangle}\;\overline{\langle Y(t)\rangle},$$		(6a)
and the quenched variance
	$$\langle Y^{2}(t)\rangle_{\mathrm{qc}}=\overline{\langle Y^{2}(t)\rangle}-\overline{\langle Y(t)\rangle\;\langle Y(t)\rangle}.$$		(6b)

The initial conditions are assumed to be periodic with $\ell_{0}$ and have an average particle density denoted by ${\overline{\rho}}$.

For free Brownian motion, numerous authors (e.g., Refs. [20, 15, 27]) have shown that

	$$\langle Y^{2}(t)\rangle_{\mathrm{ac}}=\frac{2}{{\overline{\rho}}}\sqrt{\frac{Dt}{\pi}}$$		(7a)
and
	$$\langle Y^{2}(t)\rangle_{\mathrm{qc}}=\frac{\sqrt{2}}{{\overline{\rho}}}\sqrt{\frac{Dt}{\pi}},$$		(7b)

where $D$ is the bare diffusion constant of the particles. The subdiffusive behavior is a consequence of the hard-core interaction. We also see that the annealed MSD is larger than the quenched one as the fluctuations in the initial conditions of the former lead to additional fluctuations in the displacement of the tracer.

We now come to the main question of this paper. When considering systems with a periodically varying potential and diffusivity (instead of pure Brownian diffusion), can one simply replace the bare diffusion coefficient $D$ in Eqs. (7a) and (7b) with $D_{\mathrm{eff}}$ at late times? That is, write

	$$\langle Y^{2}(t)\rangle_{\mathrm{ac}}=\frac{2}{{\overline{\rho}}}\sqrt{\frac{D_{\mathrm{eff}}t}{\pi}}$$		(8a)
and
	$$\langle Y^{2}(t)\rangle_{\mathrm{qc}}=\frac{\sqrt{2}}{{\overline{\rho}}}\sqrt{\frac{D_{\mathrm{eff}}t}{\pi}}.$$		(8b)

Simulations of SFD under a periodic potential [39, 40] support Eq. (8a) (and furthermore this idea also works when a constant applied force is added to that produced by the periodic potential). Within liquid theory, a very similar expression to Eq. (8a) for a tracer’s MSD has been derived for ensembles of SFD particles with pairwise interactions [13]. A general inequality relation with the entropy has been derived recently, becoming particularly tight for SFD [41]. As far as we know there have been no numerical studies of the effect of locally varying diffusivity in SFD.

Here we will show analytically that Eqs. (8a) and (8b) do indeed hold at late times. To do so, we will exploit Eq. (5) — the link between SFD and the effusion problem for non-interacting Brownian particles described by the FPE, Eq. (1). The effusion problem, which counts the number of crossings from left (right) to right (left) at the origin, can be solved by a late-time asymptotic analysis of the FPE for a single particle. The method we use is closely related to the mathematical theory of homogenization [42, 44, 43]. The fact that the transport coefficient we compute is that of a tracer makes our prediction particularly relevant for optical trapping experiments [35] where a spatially varying potential can be generated. For example, theoretical results such as Eq. (4) have been confirmed in experiments on single colloids [45, 46].

The paper is organised as follows: In Sec. II, we recall how the MSD of a tracer particle can be described in terms of an effusion problem. There we give an explicit expression for the MSD in terms of the solution to the single-particle FPE. In Sec. III, we carry out an asymptotic analysis of these results to prove Eqs. (8a) and (8b). In Sec. IV, we confirm our analytical predictions using the results of a numerical simulation. Finally, we conclude and discuss possible applications and extensions of our results in Sec. V.

Here we recall the basic arguments given in Refs. [29, 47], which link the displacement of a single-file tracer with the effusion of non-interacting Brownian processes to the left and right of the origin. The key starting point for our analysis is Eq. (5), which can be rewritten as

$$\int_{0}^{Y(t)}dx\rho(x,t)=\int_{-\infty}^{0}dx\rho(x,0)-\int_{-\infty}^{0}dx\rho(x,t).$$

(9)

We first consider the left-hand side. If we assume that the typical value of $Y(t)$ grows as a function of $t$, then at sufficiently late times, we will have $|Y(t)|\gg\ell$. Thus, the integral on the left-hand side can be approximated, using the CLT, by

$$\int_{0}^{Y(t)}dx\rho(x,t)\simeq{\overline{\rho}}Y(t),$$

(10)

with corrections ${\cal O}(\sqrt{{\overline{\rho}}|Y(t)|})$.

The term on the right-hand side is linked to an effusion problem as follows. Suppose the particles initially to the right and left of the origin are two different particle types (say, ‘left type’ and ‘right type’). The two particle types are non-interacting, and both obey the same FPE. Denote their densities as $\rho_{\mathrm{L}}(x,t)$ and $\rho_{\mathrm{R}}(x,t)$. This implies that, in the beginning, $\rho_{\mathrm{L}}(x>0,0)=0$ and $\rho_{\mathrm{R}}(x<0,0)=0$, while the total density is simply $\rho(x,t)=\rho_{\mathrm{L}}(x,t)+\rho_{\mathrm{R}}(x,t)$. We define $N_{+}(t)$ ($N_{-}(t)$) to be the number of particles that were initially to the left (right) at time $0$ and are on the right (left) at time $t$. Clearly,

	$\displaystyle N_{+}(t)$	$\displaystyle=$	$\displaystyle\int_{0}^{\infty}dx\rho_{\mathrm{L}}(x,t),$		(11a)
	$\displaystyle N_{-}(t)$	$\displaystyle=$	$\displaystyle\int_{-\infty}^{0}dx\rho_{\mathrm{R}}(x,t).$		(11b)

We impose conservation of the left-type particle number through

$$\int_{0}^{\infty}dx\rho_{\mathrm{L}}(x,t)+\int_{-\infty}^{0}dx\rho_{\mathrm{L}}(x,t)=\int_{-\infty}^{0}dx\rho_{\mathrm{L}}(x,0).$$

(12)

Combining Eqs. (11) and (12), we relate the change in density to the number of crossings,

	$\displaystyle N_{+}(t)-N_{-}(t)$	$\displaystyle=$	$\displaystyle\int_{-\infty}^{0}dx\rho_{\mathrm{L}}(x,0)-\int_{-\infty}^{0}dx\rho_{\mathrm{L}}(x,t)-\int_{-\infty}^{0}dx\rho_{\mathrm{R}}(x,t)$		(13)
		$\displaystyle=$	$\displaystyle\int_{-\infty}^{0}dx\rho(x,0)-\int_{-\infty}^{0}dx\rho(x,t).$

Using the asymptotic Eq. (10) and the effusion-related Eq. (13), Eq. (9) may be rewritten as

$$Y(t)=\frac{1}{{\overline{\rho}}}\left[N_{+}(t)-N_{-}(t)\right].$$

(14)

Thus we have established the link between SFD and two independent effusion proccesses. The random variables $N_{+}(t)$ and $N_{-}(t)$ are independent and identically distributed. Both the annealed and quenched variances (which are identical for the first moment) are given by $\langle Y(t)\rangle_{\mathrm{a}}=\langle Y(t)\rangle_{\mathrm{q}}=0$ in the absence of bias. Expanding and using the independence and identical distribution of $N_{+}(t)$ and $N_{-}(t)$, we can write the annealed and quenched variances solely in terms of moments of, say, $N_{+}(t)$,

	$\displaystyle\langle Y^{2}(t)\rangle_{\mathrm{ac}}$	$\displaystyle=$	$\displaystyle\frac{2}{{\overline{\rho}}^{2}}\left(\overline{\langle N_{+}^{2}(t)\rangle}-\overline{\langle N_{+}(t)\rangle}\;\overline{\langle N_{+}(t)\rangle}\right)\equiv\frac{2}{{\overline{\rho}}^{2}}\langle N_{+}^{2}(t)\rangle_{\mathrm{ac}},$		(15a)
	$\displaystyle\langle Y^{2}(t)\rangle_{\mathrm{qc}}$	$\displaystyle=$	$\displaystyle\frac{2}{{\overline{\rho}}^{2}}\left(\overline{\langle N_{+}^{2}(t)\rangle}-\overline{\langle N_{+}(t)\rangle\;\langle N_{+}(t)\rangle}\right)\equiv\frac{2}{{\overline{\rho}}^{2}}\langle N_{+}^{2}(t)\rangle_{\mathrm{qc}}.$		(15b)

Therefore, the MSD is given in terms of a single-sided effusion problem as studied in Ref. [47].

Now we examine the effusion problem more closely. The value of $N_{+}(t)$ can be written as

$$N_{+}(t)=\sum_{m=1}^{M}\Theta(X_{m}(t)),$$

(16)

where $\Theta(x)$ is the Heaviside theta function, $X_{m}(t)$ is the diffusion process described by the FPE (Eq. (1)), and $M$ is the number of particles that started to the left of the origin. We proceed by taking the averages over noise and initial conditions explicitly. If $p(x,t|x_{0},0)$ is the solution of Eq. (1) (for the initial condition $p(x,0|x_{0},0)=\delta(x-x_{0})$), the average value of $N_{+}(t)$ over the thermal noise is

$$\langle N_{+}(t)\rangle=\sum_{m=1}^{M}\int_{0}^{\infty}dxp(x,t|x_{0m};0),$$

(17)

where $\{x_{0m}\}$ is a list of the particles’ initial positions that are to the left of the origin. Suppose the particles’ positions are drawn from a periodic initial distribution $P_{0}(x)$, with a wavelength $\ell_{0}$ which does not have to be the same as that of the potential or diffusivity (being $\ell$). We define it to be normalized within one cell, $\int_{0}^{\ell_{0}}dxP_{0}(x)=1$, so the initial density of particles to the left of the origin is then $\overline{\rho_{\mathrm{L}}(x)}=(\ell_{0}M/L)P_{0}(x)$, where $L$ is the (large) size of the left half of the box. Thus, the average value of $\langle N_{+}(t)\rangle$ over the initial positions is

	$\displaystyle\overline{\langle N_{+}(t)\rangle}$	$\displaystyle=$	$\displaystyle\int_{-L}^{0}dx_{0}\overline{\rho_{\mathrm{L}}(x_{0})}\int_{0}^{\infty}dxp(x,t|x_{0m};0)=M\int_{-L}^{0}dx_{0}\frac{\ell_{0}}{L}P_{0}(x_{0})\int_{0}^{\infty}dxp(x,t|x_{0m};0)$		(18)
		$\displaystyle=$	$\displaystyle{\overline{\rho}}\ell_{0}\int_{-\infty}^{0}dx_{0}P_{0}(x_{0})\int_{0}^{\infty}dxp(x,t|x_{0},0),$

where we have taken the limits $M\to\infty$ and $L\to\infty$, while keeping ${\overline{\rho}}=M/L$ fixed.

Finally, we compute the quenched and annealed variances. Since the particles’ Brownian motions (i.e., the thermal noise values) are independent, we have

	$$\langle N^{2}_{+}(t)\rangle=\sum_{m=1}^{M}\langle\Theta(X_{m}(t))\rangle+\sum_{m=1}^{M}\sum_{n\neq m}^{M}\langle\Theta(X_{m}(t))\rangle\langle\Theta(X_{n}(t))\rangle,$$		(19a)
and thus
	$$\langle N^{2}_{+}(t)\rangle-\langle N_{+}(t)\rangle^{2}=\sum_{m=1}^{M}\langle\Theta(X_{m}(t))\rangle-\sum_{m=1}^{M}\langle\Theta(X_{m}(t))\rangle^{2},$$		(19b)

where we have used $\Theta^{2}=\Theta$. Since the particle positions are drawn from the initial distribution independently (unlike the case we consider in Appendix A; see below), averaging Eq. (19a) over the initial positions yields

$$\overline{\langle N^{2}_{+}(t)\rangle}=\int_{-\infty}^{0}dx_{0}\overline{\rho_{\mathrm{L}}(x)}\int_{0}^{\infty}dxp(x,t|x_{0},0)+\overline{\langle N_{+}(t)\rangle}\;\overline{\langle N_{+}(t)\rangle}.$$

(20)

This leads to the annealed variance

	$$\langle N^{2}_{+}(t)\rangle_{\mathrm{ac}}=\int_{-\infty}^{0}dx_{0}\overline{\rho_{\mathrm{L}}(x)}\int_{0}^{\infty}dxp(x,t|x_{0},0).$$		(21a)
On the other hand, we obtain the quenched variance by averaging Eq. (19b) over the initial conditions
	$$\langle N^{2}_{+}(t)\rangle_{\mathrm{qc}}=\int_{-\infty}^{0}dx_{0}\overline{\rho_{\mathrm{L}}(x)}\int_{0}^{\infty}dxp(x,t|x_{0},0)-\int_{-\infty}^{0}dx_{0}\overline{\rho_{\mathrm{L}}(x)}\int_{0}^{\infty}dxp(x,t|x_{0},0)\int_{0}^{\infty}dx^{\prime}p(x^{\prime},t|x_{0},0).$$		(21b)

The two averaging procedures above apply to the same physical system where all initial positions are independent, yielding two different statistics. However, the quenched average has the same mathematical expression as that for the annealed average, carried out on a system where the initial density fluctuations are suppressed [48, 47], in particular if the initial density is periodic or hyperuniform [49]. We discuss this point Appendix A and use it to simplify the simulation procedure in Sec. IV.

Looking at the expressions in Eqs. (15a), (21a), (15b), and (21b), we see that they depend on the solution $p(x,t|x_{0},0)$ of the FPE, Eq. (1), at late times. However, they only depend on terms of the form (which we define as for later convenience)

$$\frac{f(x_{0},t)}{P_{\mathrm{B}}(x_{0})}=\int_{-\infty}^{\infty}dx\ p(x,t|x_{0},0)\Theta(x),$$

(22)

where $P_{\mathrm{B}}(x)$ is the equilibrium Boltzmann factor,

	$$P_{\mathrm{B}}(x)=Z^{-1}e^{-\beta\phi(x)},$$		(23a)
with the partition function
	$$Z=\int_{0}^{\ell}dxe^{-\beta\phi(x)}.$$		(23b)

It is easy to show that $f(x,t)/P_{\mathrm{B}}(x)$ solves the adjoint FPE (or, the so-called Kolmogorov backward equation), $\partial(f/P_{\mathrm{B}})/\partial t=\hat{H}^{\dagger}(f/P_{\mathrm{B}})$, where $\hat{H}^{\dagger}$ is the adjoint operator of $\hat{H}$, and $\hat{H}$ as written in Eq. (1) [50]. We then see that $f(x,t)$ obeys the same FPE (the Kolmogorov forward equation) as Eq. (1),

	$$\frac{\partial f(x,t)}{\partial t}=\hat{H}f(x,t),$$		(24a)
with the initial condition
	$$f(x,0)=\Theta(x)P_{\mathrm{B}}(x).$$		(24b)

We now take the Laplace transform of Eq. (24a), which yields the equation for $\tilde{f}(x,s)=\int_{0}^{\infty}dt\exp(-st)f(x,t)$,

$$\hat{H}\tilde{f}(x,s)=s\tilde{f}(x,s)-f(x,0).$$

(25)

The late time behavior of $f(x,t)$ can be extracted from the behavior of $\tilde{f}(x,s)$ at small $s$. To proceed, we follow the ideas of the mathematical theory of homogenization [42]. The diffusive motion is governed by some ‘macro-variable’, $y=s^{1/2}x$, while the precise location within a cell is described by a ‘micro-variable’, still denoted by $x$, and restricted to a single cell, $0<x<\ell$. Since we take the late-time (small-$s$) limit, we expect there to be a scale separation between $y$ and $x$. Note that the potential $\phi(x)$ and diffusivity $D(x)$ still depend only on the micro-variable $x$, and both are periodic in it with period $\ell$.

To formulate the above mathematically, we write the following expansion

$$\tilde{f}(x,s)=\frac{1}{s}\sum_{n=0}^{\infty}s^{\frac{n}{2}}F_{n}(x,\sqrt{s}x),$$

(26)

where $\{F_{n}(x,y)\}$ are periodic in $x$ with period $\ell$. Upon the introduction of these two spatial variables which are treated independently, the spatial derivatives in Eq. (25) are modified to $(\partial/\partial x)\to(\partial/\partial x)+s^{1/2}\cdot(\partial/\partial y)$. We now insert the expansion of Eq. (26) in Eq. (25), and equate terms of the same order. To obtain an effective late-time diffusion equation of the form of Eq. (3), it turns out that we need the three lowest orders.

Upon equating the leading order $\mathcal{O}(s^{-1})$ terms, we find

$$\frac{\partial}{\partial x}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}F_{0}(x,y)+\frac{\partial F_{0}(x,y)}{\partial x}\right)\right]=0.$$

(27)

We see that we may assume separation of variables, and get a solution of the form

$$F_{0}(x,y)=P_{\mathrm{B}}(x)K_{0}(y),$$

(28)

where $K_{0}(y)$ is a function of $y$ which we will be able to compute later, and $P_{\mathrm{B}}(x)$ is then the Boltzmann factor. Clearly $F_{0}(x,y)$ is periodic in $x$ with period length $\ell$. (Note that the initial distribution in the SFD process, $P_{0}(x)$, need not be the equilibrium Boltzmann factor, $P_{\mathrm{B}}$, or even have the same periodicity $\ell_{0}\neq\ell$. In fact, for the quenched case, we shall use a Dirac comb for $P_{0}(x)$ in Sec. IV, as explained in Appendix A.)

At the next order, $\mathcal{O}(s^{-1/2})$, we find

$$\frac{\partial}{\partial x}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}F_{1}(x,y)+\frac{\partial F_{1}(x,y)}{\partial x}\right)\right]\\ +\frac{\partial}{\partial y}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}F_{0}(x,y)+\frac{\partial F_{0}(x,y)}{\partial x}\right)\right]+\frac{\partial}{\partial x}\left[D(x)\frac{\partial F_{0}(x,y)}{\partial y}\right]=0.$$

(29)

Now, using Eq. (28), and again making the separation of variables ansatz,

$$F_{1}(x,y)=P_{1}(x)\frac{dK_{0}(y)}{dy},$$

(30)

we find the equation for $P_{1}(x)$,

$$\frac{d}{dx}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}P_{1}(x)+\frac{dP_{1}(x)}{dx}+P_{\mathrm{B}}(x)\right)\right]=0.$$

(31)

The last $\mathcal{O}(1)$ terms give

$$\frac{\partial}{\partial x}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}F_{2}(x,y)+\frac{\partial F_{2}(x,y)}{\partial x}\right)\right]\\ +\frac{\partial}{\partial y}\left[D(x)\left(\beta\frac{d\phi(x)}{dx}F_{1}(x,y)+\frac{\partial F_{1}(x,y)}{\partial x}\right)\right]+\frac{\partial}{\partial x}\left[D(x)\frac{\partial F_{1}(x,y)}{\partial y}\right]\\ +D(x)\frac{\partial^{2}F_{0}(x,y)}{\partial y^{2}}=F_{0}(x,y)-\Theta(y)P_{0}(x).$$

(32)

Here, we replace $\Theta(x)\to\Theta(y)$ since $y$ is the macro-variable. Inserting the ansatze of Eqs. (28) and (30), and utilizing the periodicity of the solutions in $x$ to integrate Eq. (32) between $0$ and $\ell$ (thus canceling complete differentials), we find the equation for $K_{0}(y)$,

$$D_{\mathrm{eff}}\frac{\partial^{2}K_{0}(y)}{\partial y^{2}}=K_{0}(y)-\Theta(y).$$

(33)

Equation (33) is of the same form as Eq. (3) (written in Laplace space), where we identified the effective diffusion constant

$$D_{\mathrm{eff}}=\int_{0}^{\ell}dxD(x)\left[\beta\frac{d\phi(x)}{dx}P_{1}(x)+\frac{dP_{1}(x)}{\partial x}+P_{\mathrm{B}}(x)\right].$$

(34)

The solution to Eq. (33) is, using the required smoothness conditions,

$$K_{0}(y)=\begin{cases}1-\frac{1}{2}\exp\left(-y/\sqrt{D_{\mathrm{eff}}}\right),&y>0,\\ \frac{1}{2}\exp\left(y/\sqrt{D_{\mathrm{eff}}}\right).&y<0.\end{cases}$$

(35)

Prior to proceeding, we first find the explicit formula of $D_{\mathrm{eff}}$ (which can only be explicitly computed in one dimension [42].) Integrating Eq. (31), we find that

$$D(x)\left(\beta\frac{d\phi(x)}{dx}P_{1}(x)+\frac{dP_{1}(x)}{dx}+P_{\mathrm{B}}(x)\right)=C,$$

(36)

where $C$ is a constant. Comparing it with Eq. (34), we identify $D_{\mathrm{eff}}=\ell C$. Solving for $P_{1}(x)$ and using its periodicity gives

$$C=\frac{\ell}{Z\left[\int_{0}^{\ell}dxe^{\beta\phi(x)}/D(x)\right]}.$$

(37)

This then yields the Lifson-Jackson formula for the effective diffusion constant $D_{\mathrm{eff}}$ of a Brownian particle, Eq. (4).

Using these late-time solutions, we first compute the annealed variance. Note that we may rewrite Eq. (21a) more compactly in terms of $f$ as

$$\langle N_{+}^{2}(t)\rangle_{\mathrm{ac}}={\overline{\rho}}\ell_{0}\int_{-\infty}^{\infty}dx_{0}\Theta(-x_{0})P_{0}(x_{0})\frac{f(x_{0},t)}{P_{\mathrm{B}}(x_{0})}.$$

(38)

Upon combining the $y<0$ case of Eq. (35), Eq. (28), and Eq. (26), we find, to leading order in small $s$ (corresponding to late times),

$$\tilde{f}(x,s)\simeq\frac{1}{2s}P_{\mathrm{B}}(x)\exp\left(\sqrt{\frac{s}{D_{\mathrm{eff}}}}x\right),$$

(39)

for $y<0$, and see that it satisfies the initial conditions Eq. (24b). Thus, Eq. (38) becomes,

$$\langle\tilde{N}^{2}_{+}(s)\rangle_{\mathrm{ac}}=\frac{{\overline{\rho}}\ell_{0}}{2s}\int_{-\infty}^{0}dx_{0}P_{0}(x_{0})\exp\left(\sqrt{\frac{s}{D_{\mathrm{eff}}}}x_{0}\right).$$

(40)

Here we need to use the fact that $P_{0}(x)$ is periodic, so we can express it as a Fourier series,

$$P_{0}(x)=\frac{1}{\ell_{0}}\left[1+\sum_{n\neq 0}\exp\left(\frac{2\pi inx}{\ell_{0}}\right)\int_{0}^{\ell}dx^{\prime}P_{0}(x^{\prime})\exp\left(-\frac{2\pi inx^{\prime}}{\ell_{0}}\right)\right].$$

(41)

Note that the terms coming from the non-zero Fourier modes $k_{n}=2\pi n/\ell_{0}$ integrate as

	$$\int_{-\infty}^{0}dx_{0}\exp\left(\sqrt{\frac{s}{D_{\mathrm{eff}}}}x_{0}+ik_{n}x_{0}\right)=\frac{1}{ik_{n}+\sqrt{s/D_{\mathrm{eff}}}},$$		(42a)
while for $n=0$ we find
	$$\int_{-\infty}^{0}dx_{0}\exp\left(\sqrt{\frac{s}{D_{\mathrm{eff}}}}x_{0}\right)=\sqrt{\frac{D_{\mathrm{eff}}}{s}}.$$		(42b)

Therefore, at late times, the dominant term in Eq. (40) is the one corresponding to $n=0$. Inverting the Laplace transform, we find

$$\langle N^{2}_{+}(t)\rangle_{\mathrm{ac}}={\overline{\rho}}\sqrt{\frac{D_{\mathrm{eff}}t}{\pi}}.$$

(43)

Using Eq. (15a), we obtain Eq. (8a) for the annealed version of SFD.

To compute the quenched variance, we write Eq. (21b) as

	$$\langle N_{+}^{2}(t)\rangle_{\mathrm{qc}}=\langle N_{+}^{2}(t)\rangle_{\mathrm{ac}}-\lambda(t),$$		(44a)
where, in terms of $f$,
	$$\lambda(t)={\overline{\rho}}\ell_{0}\int_{-\infty}^{\infty}dx_{0}\Theta(-x_{0})P_{0}(x_{0})\frac{f^{2}(x_{0},t)}{P_{\mathrm{B}}^{2}(x_{0})}.$$		(44b)

Since $\lambda(t)$ is not linear in $f$, it is not convenient to work in terms of its Laplace transform. Upon inserting Eq. (35) in Eq. (28), an inverse Laplace transform of Eq. (26) gives for $x<0$, at late times,

$$f(x,t)\simeq\frac{1}{2}P_{\mathrm{B}}(x){\rm erfc}\left(-\frac{x}{2\sqrt{D_{\mathrm{eff}}t}}\right).$$

(45)

We then obtain

$$\lambda(t)=\frac{{\overline{\rho}}\ell_{0}}{4}\int_{-\infty}^{0}dx_{0}P_{0}(x_{0}){\rm erfc}^{2}\left(-\frac{x_{0}}{2\sqrt{D_{\mathrm{eff}}t}}\right),$$

(46)

For large $t$, following the arguments of Eqs. (41) and (42), we may only keep the non-oscillatory component of $P_{0}(x_{0})$, which is $1/\ell_{0}$. Therefore, we obtain

$$\lambda(t)={\overline{\rho}}\left(1-\frac{1}{\sqrt{2}}\right)\sqrt{\frac{D_{\mathrm{eff}}t}{\pi}}.$$

(47)

We thus find that

$$\langle N_{+}^{2}(t)\rangle_{\mathrm{qc}}={\overline{\rho}}\sqrt{\frac{D_{\mathrm{eff}}t}{2\pi}}.$$

(48)

Using Eq. (15b), we obtain Eq. (8b) for the quenched version of SFD, as well.

The above constitutes a proof that both the annealed and quenched MSDs of a tracer particle in SFD can be derived by considering purely Brownian particles but with an effective diffusion constant $D_{\mathrm{eff}}$, given by the Lifson-Jackson formula. It takes fully into account the effects of periodically-varying potential and diffusivity at the level of large-time asymptotics. We find that both MSDs grow asymptotically as $\mathcal{O}[(D_{\mathrm{eff}}t)^{1/2}/{\overline{\rho}}]$. By keeping the next-order terms in the present Sec. III, one can show that the corrections to the MSD, associated with the late-time approximation, are of order $\mathcal{O}[\ell/{\overline{\rho}}]$. A so-called ‘weak convergence theorem’ was proven mathematically for general two-scale homogenization schemes [42, 44, 43], rigorously justifying the replacement of the true $f(x,t)$ with leading order term in its asymptotic expansion. Note, however, that the homogenization-related late-time corrections are smaller than the corrections associated with the application of the CLT (Eq. (10)), which are of order $\mathcal{O}[(D_{\mathrm{eff}}t)^{1/4}/{\overline{\rho}}^{3/2}]$.

In order to verify our analytical findings, we perform Langevin dynamics simulations in the overdamped limit. The Itô stochastic differential equation corresponding to Eq. (1) is [50, 51]

$$\mathrm{d}X(t)=\left(\left.\frac{dD(x)}{dx}\right|_{X(t)}-\beta D(X(t))\left.\frac{d\phi}{dx}\right|_{X(t)}\right)\mathrm{d}t+\sqrt{2D(X(t))}\mathrm{d}B(t).$$

(49)

When this is simulated we take $B(t)$ to be a discretized Brownian motion, with the discrete form of Itô convention where its increments are Gaussian with $\langle\mathrm{d}B(k\mathrm{d}t)\mathrm{d}B(l\mathrm{d}t)\rangle=\delta_{kl}\mathrm{d}t$ and zero mean.

Below, we will consider the case of periodic potential, $\beta\phi(x)=E[1-\cos(2\pi x/\ell)]$, with constant diffusivity $D=D_{0}$. Equation (4) yields $D_{\mathrm{eff}}(E)=D_{0}/I_{0}^{2}(E)$, where $I_{0}(z)$ is the zeroth modified Bessel function of the first kind. In addition to the parameters $E$, recall that the particle density ${\overline{\rho}}$ (number of particles per period length $\ell$) is another control parameter in SFD problems.Thus, using Eqs. (8a) and (8b), we predict

	$\displaystyle\langle Y^{2}(t;E,\rho)\rangle_{\mathrm{ac}}$	$\displaystyle=$	$\displaystyle\frac{2}{{\overline{\rho}}I_{0}(E)}\sqrt{\frac{D_{0}t}{\pi}},$		(50a)
	$\displaystyle\langle Y^{2}(t;E,\rho)\rangle_{\mathrm{qc}}$	$\displaystyle=$	$\displaystyle\frac{\sqrt{2}}{{\overline{\rho}}I_{0}(E)}\sqrt{\frac{D_{0}t}{\pi}}.$		(50b)

This case was studied numerically also in Refs. [39, 40, 41].

The simulations were carried out as follows. We choose the unitless quantities $\ell=1$, $D_{0}=1/2$, and $\beta=1$ and vary the parameters ${\overline{\rho}}$ and $E$. In all simulations, we spread $M=10^{4}\ell{\overline{\rho}}$ particles within the box according to an initial distribution we explain below. In each timestep, we “sweep” the system $M$ times. In each sweep, we pick a random particle, and move it according to Eq. (49) with periodic boundary conditions in a box of size $10^{4}\ell$. We perform a total of $10^{6}$ timesteps of duration $\mathrm{d}t=10^{-3}\ell^{2}/(2D_{0})$. Each time a particle passes one of its neighbors, we relabel the relevant particles such that the ‘leftmost’ particle will be the one with the smallest numeral index (up to periodic boundary conditions). For example, if the $m$th particle skipped its neighbor from the right, then $Y_{m}(t+\mathrm{d}t)=\min(Y_{m}(t)+\mathrm{d}X(t),Y_{m+1}(t))$ and $Y_{m+1}(t+\mathrm{d}t)=\max(Y_{m}(t)+\mathrm{d}X(t),Y_{m+1}(t))$.

We consider two initial conditions. (A) Uniform distribution, where we first pick $M$ positions $X_{m,\mathrm{a}}(0)\sim\mathrm{U}(0,L)$, which are then stored in increasing order within $\{Y_{m,\mathrm{a}}(0)\}$, i.e., $Y_{m,\mathrm{a}}(0)<Y_{m+1,\mathrm{a}}(0)$. The subscript ‘$\mathrm{a}$’ stands for the fact that it would yield the annealed variance. (B) Crystalline configuration, being just a Dirac comb at positions $Y_{m,\mathrm{q}}(0)=L\cdot(m-1)/M$. The subscript ‘$\mathrm{q}$’ reminds that computing the annealed variance with these initial condition coincides with the quenched average for uncorrelated initial conditions [48, 47]. See below and Appendix A.

For each initial condition and set of parameters ${\overline{\rho}}$ and $E$, we repeat the above simulation $S=10^{3}$ times, labelled by $s\in\{1,\cdots,\ S\}$. We estimate the average displacement of the $m$th particle (the tracer) under relabelling at each time-step from the average

$$\langle\Delta Y_{m,\alpha}(t)\rangle=\frac{1}{S}\sum_{s=1}^{S}\left[Y_{m,\alpha}^{(s)}(t)-Y_{m,\alpha}^{(s)}(0)\right].$$

(51)

where $Y_{m,\alpha}^{(s)}(t)$ is its position at time $t$ for the $s$th simulation of either the uniform (annealed, $\alpha=\mathrm{a}$) or crystalline (quenched, $\alpha=\mathrm{q}$) initial distribution. Its second moment is estimated similarly from

$$\langle\Delta Y_{m,\alpha}^{2}(t)\rangle=\frac{1}{S}\sum_{s=1}^{S}\left[Y_{m,\alpha}^{(s)}(t)-Y_{m,\alpha}^{(s)}(0)\right]^{2}.$$

(52)

The estimator for the tracer’s annealed average is then given by

$$\langle\Delta Y^{2}_{m,\alpha}(t)\rangle_{\mathrm{ac}}=\langle\Delta Y_{m,\alpha}^{2}(t)\rangle-\langle\Delta Y_{m,\alpha}(t)\rangle^{2}.$$

(53)

Thus we averaged over realizations for the tracer and, since we chose a different random seed for every run, also over initial conditions for the uniform (annealed) initial condition. Here lies the convenience in estimating the quenched variance from the annealed averaging with crystalline initial conditions; otherwise, we would have needed to run many noise realizations for the same randomly-chosen initial configuration, in addition to picking the initial configurations.

Now, due to the periodic boundary conditions, the system is completely invariant under rotation, and thus any of the $M$ particles could have been chosen as the tracer. Therefore, so to increase data reliability, we take an ensemble average, and find

	$\displaystyle\langle\Delta Y^{2}(t)\rangle_{\mathrm{ac}}$	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{m=1}^{M}\langle\Delta Y^{2}_{m,\mathrm{a}}(t)\rangle_{\mathrm{ac}},$		(54a)
	$\displaystyle\langle\Delta Y^{2}(t)\rangle_{\mathrm{qc}}$	$\displaystyle=$	$\displaystyle\frac{1}{M}\sum_{m=1}^{M}\langle\Delta Y^{2}_{m,\mathrm{q}}(t)\rangle_{\mathrm{ac}}.$		(54b)

Note that performing the sum over particles already with the sum over realizations would have given unnecessary (arguably correlated, as they come from different particles for the same realization) cross terms. We now compare the numerical results with our analytical predictions, Eq. (50).

In Fig. 1 we plot the annealed and quenched variances for density $\rho=2.0$. Figure 1A is a log-log plot, showing in detail the various regimes in typical SFD in a periodic potential, as obtained from the simulations with potential strength $E=2.0$. At first, indeed the particles diffuse normally (with their bare diffusion coefficient, $D_{0}=1/2$), as they have not encountered their neighbors and typically they have not sampled the potential barriers. Afterwards, we have two simultaneous effects. The more apparent one is the (almost) plateau, which has to do with the particles attempting to overcome the potential barrier, which they succeed doing only with a finite Arrhenius rate (not shown; see Ref. [41]). Beyond this time, the particles reach their homogenized regime with $D_{\mathrm{eff}}$. The second effect is the collisions with neighboring particles, which causes a crossover from normal diffusion to anomalous subdiffusion with exponent $1/2$. Since $\rho$ and $E$ are of the same order (in fact, equal) in Fig. 1A, these two effects are not decoupled and occur together in between $t\in[0.1,10]$. Asymptotically, with the above two taking effect, the exponent-$1/2$ subdiffusion occurs with a $D_{\mathrm{eff}}$ equal to the bare diffusion coefficient. Indeed, this is what is seen in Fig. 1B, where all variances grow as $t^{1/2}$. Here, the points shown are the measured variances from the simulations (Eqs. (54)) for different potential strengths $E$ as indicated, while the lines depict the theoretical predictions (Eqs. (50)).

In Fig. 2 we plot the anomalous late-time diffusion coefficient, $F_{\mathrm{eff}}=\langle\Delta Y^{2}(t)\rangle/t^{1/2}$, which we extracted from a least-squares linear fit in the log-log plots of the variances for points at times $t\in[500,1000]$. We see a good agreement of the simulation with theory overall. The numerical values of $F_{\mathrm{eff}}$ exhibit relative deviations $\mathcal{O}(1\%)$ from the analytical ones, commonly being an underestimation of the latter. Indeed the errors are most significant in the low density region (since less collisions occur, and thus the approach to the subdiffusive behaviour is delayed) and high energy barriers (as it takes longer to overcome the potential barriers, and hence the homogenization regime is reached later). This is supported by the fact that the errors get bigger for all points when we include MSD values from earlier times (namely, $t\in[100,1000]$; not shown).