# Single-file diffusion in spatially inhomogeneous systems

Benjamin Sorkin  School of Chemistry and Center for Physics and Chemistry of
Living Systems, Tel Aviv University, 69978 Tel Aviv, Israel  David S. Dean
<EMAIL_ADDRESS>Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33400,
Talence, France.  Team MONC, INRIA Bordeaux Sud Ouest, CNRS UMR 5251, Bordeaux
INP, Univ. Bordeaux, F-33400, Talence, France.

###### Abstract

We study the effect of spatially-varying potential and diffusivity on the
dispersion of a tracer particle in single-file diffusion. Non-interacting
particles in such a system exhibit normal diffusion at late times, which is
characterised by an effective diffusion constant $D_{\mathrm{eff}}$. Here we
demonstrate the physically appealing result that the dispersion of single-file
tracers in this system has the same long-time behavior as that for Brownian
particles in a spatially-homogeneous system with constant diffusivity
$D_{\mathrm{eff}}$. Our results are based on a late-time analysis of the
Fokker-Planck equation, motivated by the mathematical theory of
homogenization. The findings are confirmed by numerical simulations for both
annealed and quenched initial conditions.

##  I Introduction

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.

##  II Link between single-file diffusion and effusion

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.

##  III Late time analysis of the Fokker-Planck equation

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}]$.

##  IV Numerical simulation

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)).

Figure 1: Numerical Simulation, Periodic Potential. Tracer’s annealed and
quenched variances (mean-squared displacement, MSD) versus time, for density
${\overline{\rho}}=2$. Panel A: A log-log plot for $E=2.0$, obtained from the
simulations. It shows a normal-diffusive regime for short times (no collisions
yet and almost-free diffusion around the potential’s minima), a brief plateau
in between (slow-down due to the potential barriers), and exponent-$1/2$
subdiffusion at late times (after homogenization and sufficient interparticle
collisions). Panel B: MSD plots for potential strengths $E$ as indicated.
Points are the variances obtained from simulations, while the immediately-
adjacent lines are the corresponding 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).

Figure 2: Numerical Simulation, Periodic Potential. Tracer’s annealed and
quenched anomalous late-time diffusion constant versus potential strengths
$E$, for different densities ${\overline{\rho}}$, as indicated. Points are the
values obtained from simulations using a least-squares linear fit to the
variances’ log-log plots (as in Fig. 1A), while the immediately-adjacent lines
are the corresponding theoretical predictions.

To test our prediction further, we also consider the case of periodic
diffusivity, given by $D(x)=D_{0}\exp[a\cos(2\pi x/\ell)]$, in the absence of
potential. The effective Lifson-Jackson diffusivity of Eq. (4) is given by
$D_{\mathrm{eff}}(a)=D_{0}/I_{0}(a)$, and so our predictions for the MSDs read

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

The simulation procedure and the calculation of the MSD are identical to
before. (We rely on fewer repetitions of the simulations, $S=200$.) However,
due to the exponentiation of the diffusivity variability $a$, the diffusivity
is very high at $x=k\ell$ ($k\in\mathbb{Z}$) for big $a$s, and so, for $a=3.0$
we had to reduce the timestep duration to $\Delta t=10^{-4}\ell^{2}/(2D_{0})$.
Otherwise, the MSDs are overestimated (as seen, to some extent, for $a=2.0$).
With this correction, Fig. 3 shows that the MSDs extracted from the simulation
agree well with the theory; so to improve them further, one could consider
lower $\Delta t$s, as well as longer-running simulations. A least-squares
linear fit in the log-log plot of the MSDs yields the late-time anomalous
diffusion constants for this case, which are shown along with the theoretical
prediction in Fig. 4.

Figure 3: Numerical Simulation, Periodic Diffusivity. Tracer’s annealed and
quenched variances (mean-squared displacement, MSD) versus time, for density
${\overline{\rho}}=2$, and diffusivity variabilities $a$ as indicated. Points
are the variances obtained from simulations, while the immediately-adjacent
lines are the corresponding theoretical predictions, Eqs. (55). For $a=3.0$, a
shorter timestep duration was used. Figure 4: Numerical Simulation, Periodic
Diffusivity. Tracer’s annealed and quenched anomalous late-time diffusion
constant versus diffusivity variability $a$, for different densities
${\overline{\rho}}$, as indicated. Points are the values obtained from
simulations using a least-squares linear fit to the variances’ log-log plots,
while the immediately-adjacent lines are the corresponding theoretical
predictions. For $a=3.0$, a shorter timestep duration was used.

##  V Conclusions

We have revisited the question of whether one can use the effective diffusion
constant for a single particle in a medium of periodically-varying potential
and diffusivity [39] to describe the late-time dispersion of single-file
diffusing particles in such media. In agreement with previous simulations [39,
40] and the ones of Sec. IV, we prove that one can indeed replace the bare
diffusivity with the effective diffusion constant. This implies that one may
use a coarse-grained description to understand SFD in inhomogeneous systems.
Although it would have been very surprising were it not so (by sheer
intuition), the demonstration is non-trivial.

In our scheme, we have adapted the ideas behind the mathematical theory of
homogenisation [42, 43, 44], couching them in terms of a late-time asymptotic
expansion (more familiar to physicists), which may be of use in other systems.
In the context of SFD, it would be useful to rigorously justify a similar
utilization of the system-appropriate $D_{\mathrm{eff}}$ in cases such as
particles with pairwise additive hydrodynamic correlations and direct
interactions [13] or in the presence of random diffusivity or external
potential [4]. The method we propose rigorously justifies the coarse-graining
procedure evoked to derive the macroscopic fluctuation theory [36, 37, 38],
whereby short-ranged interactions are incorporated into similar effective
diffusion constants and mobilities in a stochastic-hydrodynamic description of
interacting particle systems.

###### Acknowledgements.

We thank Haim Diamant for extensive and enlightening discussions. B.S.
acknowledges support from the Israel Science Foundation (Grant No. 986/18).
D.S.D. would like to thank the Mortimer and Raymond Sackler Institute of
Advance Studies in Tel Aviv for their support of his visit to Tel Aviv
University where this work was initiated and also acknowledges financial
support from the European Union through the European Research Council under
EMet- Brown (ERC-CoG-101039103) grant.

##  Appendix A Mimicking quenched variance from a periodic initial
distribution

In this appendix, we give an argument why choosing periodic crystalline
initial conditions and computing the resulting annealed variance yields the
quenched one. This point is discussed in Refs. [48, 47]. Denote the random
initial particle positions that are to the left of the origin as
$\\{X_{m}(0)\\}$, and define the initial density field as
$\rho_{\mathrm{L}}(x_{0})=\sum_{m=1}^{M}\delta(x_{0}-X_{m}(0))$, which is
random as well. First, we rewrite Eq. (19a) as

$\langle
N^{2}_{+}(t)\rangle=\sum_{m=1}^{M}\langle\Theta(X_{m}(t))\rangle+\sum_{m=1}^{M}\sum_{n=1}^{M}\langle\Theta(X_{m}(t))\rangle\langle\Theta(X_{n}(t))\rangle-\sum_{m=1}^{M}\langle\Theta(X_{m}(t))\rangle\langle\Theta(X_{m}(t))\rangle.$
(56)

Then, we express Eqs. (17) and (56) in terms of $\rho_{\mathrm{L}}(x_{0})$ as

$\displaystyle\langle N_{+}(t)\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dx\int_{-L}^{0}dx_{0}p(x,t|x_{0},0)\rho_{\mathrm{L}}(x_{0}).$
(57a) $\displaystyle\langle N^{2}_{+}(t)\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dx\int_{-L}^{0}dx_{0}p(x,t|x_{0},0)\rho_{\mathrm{L}}(x_{0})$
(57b)
$\displaystyle+\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-L}^{0}dx_{0}\int_{-L}^{0}dx^{\prime}_{0}p(x,t|x_{0},0)p(x^{\prime},t|x^{\prime}_{0},0)\rho_{\mathrm{L}}(x_{0})\rho_{\mathrm{L}}(x^{\prime}_{0})$
$\displaystyle-\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-L}^{0}dx_{0}p(x,t|x_{0},0)p(x^{\prime},t|x_{0},0)\rho_{\mathrm{L}}(x_{0}),$

where we kept the (half-) box size $L$ finite for the moment. Thus, upon
taking the average over initial conditions, the general expression for the
annealed average in terms of $\rho_{\mathrm{L}}(x_{0})$ is

$\langle
N^{2}_{+}(t)\rangle_{\mathrm{ac}}=\int_{0}^{\infty}dx\int_{-L}^{0}dx_{0}\
p(x,t|x_{0},0)\overline{\rho_{\mathrm{L}}(x_{0})}+\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-L}^{0}dx_{0}\int_{-L}^{0}dx^{\prime}_{0}p(x,t|x_{0},0)p(x^{\prime},t|x^{\prime}_{0},0)\times\\\
\times\left[\overline{\rho_{\mathrm{L}}(x_{0})\;\rho_{\mathrm{L}}(x^{\prime}_{0})}-\overline{\rho_{\mathrm{L}}(x_{0})}\delta(x_{0}-x_{0}^{\prime})-\overline{\rho_{\mathrm{L}}(x_{0})}\;\overline{\rho_{\mathrm{L}}(x^{\prime}_{0})}\right].$
(58)

In the case where the particles are initially independent, only the first term
in Eq. (58) contributes, and we recover the annealed result for independent
particles — Eq. (21a).

We have shown in the arguments leading to Eq. (42b) that quantities of the
form

$\int_{-\infty}^{0}dx_{0}p(x,t|x_{0},0)g(x_{0}),$ (59a) can be evaluated at
late times as $\int_{-\infty}^{0}dx_{0}p(x,t|x_{0},0)\langle
g\rangle_{\mathrm{1p}},$ (59b)

if $g$ is a periodic function and $\langle g\rangle_{\mathrm{1p}}$ indicates
the spatial average of $g$ over one period.

In this appendix, we consider a system where the initial positions of the
particles are placed on a periodic lattice with randomness given only by
global small (a multiple of $\ell_{0}$) translations. The density is thus
random but periodic (i.e., hyperuniform). For Eq. (58) it implies (in the
limit $L\to\infty$) that

$\int_{0}^{\infty}dx\int_{-L}^{0}dx_{0}\
p(x,t|x_{0},0)\overline{\rho_{\mathrm{L}}(x_{0})}\to\int_{0}^{\infty}dx\int_{-\infty}^{0}dx_{0}\
p(x,t|x_{0},0)\overline{\langle\rho_{\mathrm{L}}\rangle_{\mathrm{1p}}}$ (60a)
and
$\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-L}^{0}dx_{0}\int_{-L}^{0}dx^{\prime}_{0}p(x,t|x_{0},0)p(x^{\prime},t|x^{\prime}_{0},0)\overline{\rho_{\mathrm{L}}(x_{0})\;\rho_{\mathrm{L}}(x^{\prime}_{0})}\to\\\
\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-\infty}^{0}dx_{0}\int_{-\infty}^{0}dx^{\prime}_{0}p(x,t|x_{0},0)p(x^{\prime},t|x^{\prime}_{0},0)\overline{\langle\rho_{\mathrm{L}}\rangle_{\mathrm{1p}}^{2}}.$
(60b)

This can be seen by carrying out the integrals over $x_{0}$ and
$x_{0}^{\prime}$ prior to the average over initial conditions. The disorder in
the above, which can be due to a random lattice translation, is thus
irrelevant at late times (as should clearly be the case on physical grounds).
This means that the first and third terms in Eq. (58) cancel, yielding

$\langle
N^{2}_{+}(t)\rangle_{\mathrm{ac}}=\int_{0}^{\infty}dx\int_{-\infty}^{0}dx_{0}\
p(x,t|x_{0},0)\overline{\rho_{\mathrm{L}}(x_{0})}-\int_{0}^{\infty}dx\int_{0}^{\infty}dx^{\prime}\int_{-\infty}^{0}dx_{0}p(x,t|x_{0},0)p(x^{\prime},t|x_{0},0)\overline{\rho_{\mathrm{L}}(x_{0})},$
(61)

which corresponds to the quenched variance for uncorrelated particles, as
given by Eq. (21b).

## References

  * [1] S. Lifson and J.L. Jackson, J. Chem. Phys. 36, 2410 (1962).
  * [2] T. Guérin and D.S. Dean, Phys. Rev. Lett. 15, 020601 (2015).
  * [3] T. Guérin and D.S. Dean, Phys. Rev. E 15, 062103 (2015).
  * [4] D.S. Dean, I.T. Drummond, and R.R. Horgan, J. Stat. Mech., P07013 (2007).
  * [5] D.S. Dean, I.T. Drummond, R.R. Horgan, J. Phys. A: Math. Gen. 27 5135 (1994).
  * [6] R.D.L. Hanes, C. Dalle-Ferrier, M. Schmiedeberg, M.C. Jenkins and S. U. Egelhaaf, Soft Matter. 8, 2714 (2012).
  * [7] F. Evers, C. Zunke, R.D.L. Hanes, J. Bewerunge, I. Ladadwa, A. Heuer and S.U. Egelhaaf, Phys. Rev. E 88, 022125 (2013).
  * [8] C. Zunke, J. Bewerunge, F. Olatten, S. U. Egelhaaf and A. Godec, Sci. Adv. 8, eabk0627 (2022).
  * [9] D. S. Dean and G. Oshanin, Phys. Rev. E 90, 022112 (2014).
  * [10] Ya. G. Sinai, Theory of Prob. and Appl. 27, 247 (1982).
  * [11] A. Comtet and D.S. Dean, J. Phys. A: Math. Gen. 31 8595 (1998).
  * [12] T. Bodineau and B. Derrida, Phys. Rev. Lett. 92, 180601 (2004).
  * [13] M. Kollmann, Phys. Rev. Lett. 90, 180602 (2003).
  * [14] S. Prolhac and K. Mallick, J. Phys. A: Math. Theor. 41 175002 (2008).
  * [15] S. Alexander and P. Pincus, Phys. Rev. B 18 2011 (1978).
  * [16] C. Appert-Rolland, B. Derrida, V. Lecomte, F. Van Wijland, Phys. Rev. E. 78, 021122 (2008).
  * [17] M. Prähofer and H. Spohn, Current fluctuations for the totally asymmetric simple exclusion process, in In and out of equilibrium, ed. V. Sidoravicius, Progress in Probability 51, 185 (Birkhauser Boston, 2002).
  * [18] P. L Krapivsky and B. Meerson, Phys. Rev. E 86, 031106 (2012).
  * [19] B. Derrida and A. Gerschenfeld, J. Stat. Phys. 137, 978 (2009).
  * [20] T. Harris, J. Appl. Probab. 2, 323 (1965).
  * [21] F. Spitzer, Adv. Math. 5, 246 (1970)
  * [22] P. M. Richards, Phys. Rev.B 16, 1393 (1977).
  * [23] R. Arratia, Z. Ann. Probab. 11, 362 (1983).
  * [24] T. M. Liggett, Interacting Particle Systems (Springer-Verlag, New York, 1983)
  * [25] C. Rödenbeck, J. Kärger, and K. Hahn, Phys. Rev. E 57, 4382 (1998).
  * [26] L. Lizana and T. Ambjörnsson†, Phys. Rev. Lett. 100, 200601 (2008).
  * [27] P.L. Krapivsky, K. Mallick and T. Sadhu, Phys. Rev. Lett. 113 078101 (2014).
  * [28] N. Leibovich and E. Barkai, Phys. Rev. E, 88, 032107 (2013).
  * [29] D. Dürr, S. Goldstein, and J. L. Lebowitz, Commun. Pure Appl. Math. 38, 573 (1985).
  * [30] J. Kärger and D. Ruthven, Diffusion in zeolites and other microporous solids (Wiley, 1992).
  * [31] Q.-H. Wei, C. Bechinger and P. Leiderer, Science 287, 625 (2000).
  * [32] C. Lutz, M. Kollmann and C. Bechinger, Phys. Rev. Lett. 93, 026001 (2004).
  * [33] S. Y. Yang, J.-A Yang, E.-S. Kim, G. Jeon, E. J. Oh, K. Y. Choi, S. K. Hahn, and J. K. Kim, ACS Nano, 4, 3817 (2010).
  * [34] G.-W. Li, O. G. Berg, and J. Elf, Nat. Phys., 5, 294 (2009).
  * [35] F. Evers, R. D. L. Hanes, C. Zunke, R. F. Capellmann, J. Bewerunge, C. Dalle-Ferrier, M. C. Jenkins, I. Ladadwa, A. Heuer, R. Castañeda-Priego, and S. U. Egelhaaf, Eur. Phys. J. Special Topics 222, 2995–3009 (2013).
  * [36] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim, Phys. Rev. Lett. 87, 040601 (2001).
  * [37] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Phys. Rev. Lett. 94, 030601 (2005).
  * [38] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, J. Stat. Phys. 123237 (2006).
  * [39] A. Taloni and F. Marchesoni, Phys. Rev. Lett. 96, 020601 (2006).
  * [40] D. Lips, A. Ryabov, and P. Maass, Phys. Rev. E 100, 052121 (2019).
  * [41] B. Sorkin, H. Diamant, and G. Ariel, arXiv:2212.12746 (2022).
  * [42] G. A. Pavliotis and A. Stuart, Multiscale Methods: Averaging and Homogenization (Springer New York, 2008).
  * [43] G. Allaire, SIAM J. Math. Anal., 23, 1482 (1992).
  * [44] G. Nguetseng, SIAM J. Math. Anal., 23, 1482 (1989).
  * [45] X.-G. Ma, P.-Y. Lai, and P. Tong, Soft Matter, 9, 8826 (2013).
  * [46] Y. Su, X.-G. Ma, P.-Y. Lai, and P. Tong, Soft Matter, 13, 4773 (2017).
  * [47] D. S. Dean, S. N. Majumdar, and G. Schehr, arXiv:2303.08961 (2023).
  * [48] T. Banerjee, R. L. Jack, and M. E. Cates, Phys. Rev. E 106, L062101 (2022).
  * [49] S. Torquato, Phys. Rep. 745, 1 (2018).
  * [50] B. Øksendal, Stochastic differential equations: An introduction with applications (Springer Berlin Heidelberg, 2003).
  * [51] A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 76, 011123 (2007).