# Approach to asymptotically diffusive behavior for Brownian particles in media with periodic diffusivities

David S. Dean, Thomas Guérin

<figure><img src="image_1.png" /></figure>

## To cite this version:

David S. Dean, Thomas Guérin. Approach to asymptotically diffusive behavior for Brownian particles in media with periodic diffusivities. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2014, 90 (6), pp.062114 (1-6). 10.1103/PhysRevE.90.062114. hal-01119884

HAL Id: hal-01119884

https://hal.science/hal-01119884

Submitted on 24 Feb 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Approach to asymptotically diffusive behavior for Brownian particles in media with periodic diffusivities

David S. Dean¹ and Thomas Guérin¹

¹Université de Bordeaux and CNRS, Laboratoire Ondes et Matière d'Aquitaine (LOMA), UMR 5798, F-33400 Talence, France

We analyze the mean squared displacement of a Brownian particle in a medium with a spatially varying local diffusivity which is assumed to be periodic. When the system is asymptotically diffusive the mean squared displacement, characterizing the dispersion in the system, is, at late times, a linear function of time. A Kubo type formula is given for the mean squared displacement which allows the recovery of some known results for the effective diffusion constant $D_e$ in a direct way, but also allows an understanding of the asymptotic approach to the diffusive limit. In particular, as well as as computing the slope of a linear fit to the late time mean squared displacement, we find a formula for the constant where the fit intersects the y axis.

PACS numbers: 05.10.Gg, 02.50.Ey

## I. INTRODUCTION

Many microscopic processes in solids and liquids can be described at a coarse grained level as Langevin or Brownian processes [1]. The transport of Brownian particles subject to an external potential has been widely studied, the effect of the potential on the dispersion of the particle is to slow it down due to the trapping of the particle in local minima of the potential [2, 3]. Many transport processes are asymptotically diffusive, and the late time diffusion constant can be obtained by analyzing the mean squared displacement (MSD). Extracting information from particle trajectories is a delicate business, and even the optimal way to estimate the diffusion constant for pure Brownian motion is far from trivial [4]. However in systems where the diffusion is not a pure Brownian motion, in most cases the diffusive behavior only sets in at late times and short time transients make the statistical analysis of trajectories much more complicated and may even lead to behavior of the MSD similar to that seen for anomalous diffusion [2]. Understanding finite time corrections to dispersive laws is thus potentially very important to correctly interpret experimental data.

In most studies, the local diffusion constant or diffusivity is assumed to be constant meaning that the local mobility of the particle is assumed to be independent of its environment; only the local drift on the particle, generated by a spatially varying potential varies in space. For periodically varying bounded potentials, the mean squared displacement of the particle becomes linear at late times corresponding to normal diffusion. The late time diffusion constant, when it exists, is lower than the local diffusion constant, which turns out to be the short time diffusion constant, as diffusion dominates over drift at very short times due to the fact that it scales as $\sqrt{t}$, while drift scales as $t$. The effective diffusion constant for these models can be computed exactly in one dimension both without [5] and with an external bias [6]. Recently, it was shown that the asymptotic approach to the final diffusive behavior could be studied in these systems via a Kubo-type formula for the MSD [7]. In this paper we show how one can derive a Kubo formula for a Brownian particle with a spatially varying (but constant in time) diffusion coefficient (which is assumed to be periodic), in the case where there is no external potential acting on the particle. The formula we derive is general and valid in any dimension, however the quantities involved can, in general, only be exactly computed in one dimension.

The probability density function (pdf) of a particle diffusing in a medium with a position dependent (but isotropic) diffusion constant (or diffusivity) $\kappa(\mathbf{x})$ obeys the diffusion equation

$$ \frac{\partial p(\mathbf{x}; t)}{\partial t} = \nabla \cdot \kappa(\mathbf{x}) \nabla p(\mathbf{x}; t). \qquad (1) $$

This equation can be used to model a number of different situations where the local particle mobility, or equivalently diffusion constant, depends on the local environment. An example is diffusion in porous media, where the effective diffusion constant is a function of the local pore size. The most widely studied statistical quantity for such systems, which can be extracted from single particle trajectories, or from following particle ensembles, is the MSD. In systems where the asymptotic transport properties are normal (or classically diffusive) the MSD is given by

$$ \langle (\mathbf{X}_t - \mathbf{X}_0)^2 \rangle = 2dD(t)t, \qquad (2) $$

where $D(t)$ is the effective time dependent diffusion constant and $d$ the spatial dimension. When the system is asymptotically diffusive the effective, or late time, diffusion constant is defined via the limit

$$ D_e = \lim_{t \to \infty} D(t). \qquad (3) $$
Clearly, for pure Brownian motion where $\kappa(\mathbf{x}) = \kappa$ is independent of spatial position, the time dependent diffusion constant is constant in time and given simply by $D(t) = \kappa = D_e$.

## II. A KUBO-TYPE FORMULA FOR THE MSD IN A MEDIUM OF VARYING DIFFUSIVITY

In the Ito prescription of the stochastic calculus, the stochastic differential equation obeyed by particles described the the Fokker-Planck equation Eq. (1) is [8]

$$d\mathbf{X}_t = \sqrt{2\kappa(\mathbf{X}_t)}d\mathbf{B}_t + \nabla\kappa(\mathbf{X}_t)dt \quad (4)$$

where $\langle dB_{it}dB_{jt} \rangle = \delta_{ij}dt$ and all increments $d\mathbf{B}_t$ at different times are independent and of mean zero. Integrating Eq. (4) yields

$$\mathbf{X}_t - \mathbf{X}_0 = \int_0^t ds \sqrt{2\kappa(\mathbf{X}_s)}d\mathbf{B}_s + \int_0^t ds \nabla\kappa(\mathbf{X}_s). \quad (5)$$

The derivation of the Kubo formula is based on a simple trick, one subtracts the second term of the right hand side above from both sides and then squares the resulting equation to obtain:

$$\langle (\mathbf{X}_t - \mathbf{X}_0)^2 - 2(\mathbf{X}_t - \mathbf{X}_0) \cdot \int_0^t ds \nabla\kappa(\mathbf{X}_s) + \int_0^t \int_0^t ds ds' \nabla\kappa(\mathbf{X}_s) \cdot \nabla\kappa(\mathbf{X}_{s'}) \rangle = \langle \int_0^t \int_0^t ds ds' \sqrt{2\kappa(\mathbf{X}_s)} \sqrt{2\kappa(\mathbf{X}_{s'})} d\mathbf{B}_s \cdot d\mathbf{B}_{s'} \rangle. \quad (6)$$

We assume that the diffusion constant $\kappa(\mathbf{x})$ is a periodic function in $d$ dimensions and we assume it is of period $l$ in all cartesian coordinates. Clearly, the methods presented here also apply to spaces divided into an infinitely repeating array of elementary unit cells $\Omega$ over which $\kappa$ is periodic, and which are not necessarily rectangular. We will now derive a formula that relates the MSD of a particle over all space to the properties of the diffusion in each individual cell $\Omega$.

In order to proceed we assume that the system is of overall length $L \gg l$, where $L$ is an integer multiple of $l$ and that the system is periodic in $L$. We assume that $\mathbf{X}_0$ is uniformly distributed with pdf $p_0(\mathbf{x}_0) = 1/L^d$ on the large periodic system. This choice of initial conditions will mean that the process $\mathbf{Y}_s = \mathbf{X}_s \mod l$ will be equilibrated and have initial conditions for $\mathbf{Y}_0$ given by the pdf $p_0(\mathbf{y}) = 1/l^d$. The integrated stochastic differential equation Eq. (6) is clearly not valid for a periodic system, however, if the asymptotic diffusion constant is given by $D_e$, so as long as $\langle (\mathbf{X}_t - \mathbf{X}_0)^2 \rangle \sim 2dD_e t \ll L^2$ the use of Eq. (6) is valid. Furthermore the time scale on which the asymptotic diffusive regime sets in will be determined by the length scale $l$ and given by $\tau \sim l^2/D_e$. We now proceed to evaluate the various terms in Eq. (6). By definition we have that

$$\langle (\mathbf{X}_t - \mathbf{X}_0)^2 \rangle = 2dD(t)t. \quad (7)$$

The second term on the left hand side of Eq. (6) can be evaluated using

$$\langle (\mathbf{X}_t - \mathbf{X}_0) \cdot \int_0^t ds \nabla\kappa(\mathbf{X}_s) \rangle = \int_0^t ds \left\{ \int d\mathbf{x} d\mathbf{y} d\mathbf{x}_0 p(\mathbf{x}, \mathbf{y}; t-s) p(\mathbf{y}, \mathbf{x}_0; s) p_0(\mathbf{x}_0) \mathbf{x} \cdot \nabla\kappa(\mathbf{y}) - \int d\mathbf{y} d\mathbf{x}_0 p(\mathbf{y}, \mathbf{x}_0; s) p_0(\mathbf{x}_0) \mathbf{x}_0 \cdot \nabla\kappa(\mathbf{y}) \right\} \quad (8)$$

where $p(\mathbf{x}, \mathbf{x}_0; s)$ is the transition probability density or propagator for a particle started at time 0 at the point $\mathbf{x}_0$ to be at $\mathbf{x}$ at time $s$. Here it is useful to note that the propagator $p(\mathbf{x}, \mathbf{y}; t)$ obeys the diffusion equation

$$\frac{\partial p}{\partial t} = -Hp \quad (9)$$

where $H$ is the self-adjoint operator defined by its action on a function $f$ by

$$H f(\mathbf{x}) = - \int d\mathbf{y} \nabla \cdot \kappa(\mathbf{x}) \nabla \delta(\mathbf{x} - \mathbf{y}) f(\mathbf{y}) \quad (10)$$

The initial condition for Eq. (9) is given by $p(\mathbf{x}, \mathbf{y}; 0) = \delta(\mathbf{x} - \mathbf{y})$ and with this initial condition the propagator can be written in operator notation as

$$p(\mathbf{x}, \mathbf{y}; t) = \exp(-tH)(\mathbf{x}, \mathbf{y}). \quad (11)$$
As $H$ is self-adjoint, we also note that $\exp(-tH)$ is self adjoint and thus $p(\mathbf{x}, \mathbf{y}; t) = p(\mathbf{y}, \mathbf{x}; t)$. This symmetry of the transition density is crucial for the evaluation of the left hand side of Eq. (8). We first note that as $p_0(\mathbf{x})$ is the equilibrium distribution for the process we have that,

$$ \int d\mathbf{x}_0 p(\mathbf{y}, \mathbf{x}_0; s) p_0(\mathbf{x}_0) = p_0(\mathbf{y}). \quad (12) $$

Using this property in Eq. (8), we obtain:

$$ \begin{aligned} & \langle (\mathbf{X}_t - \mathbf{X}_0) \cdot \int_0^t ds \ \nabla \kappa(\mathbf{X}_s) \rangle \\ &= \int_0^t ds \left\{ \int d\mathbf{x}\,d\mathbf{y} \ p(\mathbf{x}, \mathbf{y}; t-s) p_0(\mathbf{y}) \mathbf{x} \cdot \nabla \kappa(\mathbf{y}) - \int d\mathbf{y}\,d\mathbf{x}_0 \ p(\mathbf{y}, \mathbf{x}_0; s) p_0(\mathbf{x}_0) \mathbf{x}_0 \cdot \nabla \kappa(\mathbf{y}) \right\} \end{aligned} \quad (13) $$

Now, changing the time integration variable in the first term from $t-s$ to $s$ and relabeling the spatial integration variables and using the symmetry of the transition density we see that the term in Eq. (13) is identically zero.

As for the third term on the left hand side of Eq. (6), we use the fact that $\kappa(\mathbf{X}_t) = \kappa(\mathbf{Y}_t)$. The process $\mathbf{Y}_t$ has the same Fokker-Planck equation as $\mathbf{X}_t$ but where the operator $H$ acts on functions with $\mathbf{x} \in \Omega \equiv [0, l]^d$ with periodic boundary conditions. We thus find

$$ \begin{aligned} \langle \int_0^t \int_0^t ds ds' \nabla \kappa(\mathbf{X}_s) \cdot \nabla \kappa(\mathbf{X}_{s'}) \rangle &= 2 \int_0^t ds \int_0^s ds' \int_\Omega d\mathbf{x}\,d\mathbf{y} \, p(\mathbf{x}, \mathbf{y}; s-s') p_0(\mathbf{y}) \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}) \\ &= 2 \int_0^t ds \int_0^s ds' \int_\Omega d\mathbf{x}\,d\mathbf{y} \, \exp(-(s-s')H) (\mathbf{x}, \mathbf{y}) p_0(\mathbf{y}) \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}). \end{aligned} \quad (14) $$

The operator $\exp(-tH)$ on $\Omega$ has an eigenfunction decomposition of the form

$$ \exp(-tH)(\mathbf{x}, \mathbf{y}) = \frac{1}{l^d} + \sum_{\lambda > 0} \exp(-\lambda t) \psi_{\lambda}(\mathbf{x}) \psi_{\lambda}(\mathbf{y}). \quad (15) $$

The first term above corresponds to the equilibrium distribution, and the following sum is over the non-zero eigenvalues $\lambda$ of $H$ with eigenfunctions $\psi_\lambda(\mathbf{x})$. Due to the periodicity of $\kappa$, the term coming from the equilibrium distribution does not contribute and it is now straightforward to see that

$$ \langle \int_0^t \int_0^t ds ds' \nabla \kappa(\mathbf{X}_s) \cdot \nabla \kappa(\mathbf{X}_{s'}) \rangle = \frac{2}{l^d} \int_{\Omega} d\mathbf{x}\,d\mathbf{y} \left[ \sum_{\lambda} \left( \frac{t}{\lambda} - \frac{1 - \exp(-\lambda t)}{\lambda^2} \right) \psi_{\lambda}(\mathbf{x}) \psi_{\lambda}(\mathbf{y}) \right] \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}). \quad (16) $$

The equivalent expression in operator notation is:

$$ \langle \int_0^t \int_0^t ds ds' \nabla \kappa(\mathbf{X}_s) \cdot \nabla \kappa(\mathbf{X}_{s'}) \rangle = \frac{2}{l^d} \int_{\Omega} d\mathbf{x}\,d\mathbf{y} \left[ tH'^{-1}(\mathbf{x}, \mathbf{y}) - H'^{-2}(\mathbf{x}, \mathbf{y}) + H'^{-2} \exp(-tH')(\mathbf{x}, \mathbf{y}) \right] \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}), \quad (17) $$

where $H'$ denotes the operator $H$ acting on the subspace of functions orthogonal to the zero eigenvalue eigenfunction $\psi_0(\mathbf{x}) = 1/\sqrt{l^d}$. The operator $H'^{-1}$ is thus the pseudo-Green's function for $H$.

To complete the computation we evaluate the right hand side of Eq. (6). As we are using the Ito prescription of stochastic calculus, we have that for $s > s'$ that the increment $d\mathbf{B}_s$ is independent of $d\mathbf{B}_{s'}$, $\mathbf{X}_{s'}$ and also $\mathbf{X}_s$ – consequently we have that for $s \neq s'$

$$ \langle \sqrt{2\kappa(\mathbf{X}_s)} \sqrt{2\kappa(\mathbf{X}_{s'})} d\mathbf{B}_s \cdot d\mathbf{B}_{s'} \rangle = 0. \quad (18) $$

However for $s = s'$ we have

$$ \langle 2\kappa(\mathbf{X}_s)[d\mathbf{B}_s]^2 \rangle = 2d\langle\kappa(\mathbf{X}_s)\rangle ds \quad (19) $$

where on the right hand side above the remaining average is with respect to the position $\mathbf{X}_s$. However as $\kappa(\mathbf{X}_s) = \kappa(\mathbf{Y}_s)$ and $\mathbf{Y}_s$ is in equilibrium with a uniform distribution we find that

$$ \langle 2\kappa(\mathbf{X}_s)[d\mathbf{B}_s]^2 \rangle = 2d\bar{\kappa} \, ds, \quad (20) $$
where $\bar{\kappa}$ indicates the spatial average of an $l$ periodic function over the region $\Omega$ and so

$$ \bar{\kappa} = \frac{1}{l^d} \int_{\Omega} d\mathbf{x} \, \kappa(\mathbf{x}). \qquad (21) $$

Finally collecting all these results together we obtain

$$ D(t) = \bar{\kappa} - \frac{1}{l^d d} \int_{\Omega} d\mathbf{x} d\mathbf{y} H'^{-1}(\mathbf{x}, \mathbf{y}) \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}) + \frac{1}{dl^d t} \int_{\Omega} d\mathbf{x} d\mathbf{y} [H'^{-2}(\mathbf{x}, \mathbf{y}) - H'^{-2} \exp(-tH')(\mathbf{x}, \mathbf{y})] \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}). \qquad (22) $$

This formula relates explicitly the diffusion coefficient $D(t)$ in infinite space to the properties of the diffusion on the unit domain $\Omega$.

Without going into explicit computations we can make a few general observations about the time dependent diffusion constant $D(t)$ from Eq. (22). First, at short times it is straight forward to see that

$$ D(t) = \bar{\kappa} + \mathcal{O}(t), \qquad (23) $$

thus the short-time diffusion constant is given by the average diffusivity, this is easy to understand as at any time (given our assumption that $\mathbf{Y}$ starts in equilibrium) the distribution of the process $\mathbf{Y}_s$ is uniform. Recall, that in the late time limit, we define the asymptotic diffusion constant as $D_e = \lim_{t \to \infty} D(t)$ and hence in Eq. (22), only the first two terms contribute to this limit yielding

$$ D_e = \bar{\kappa} - \frac{1}{l^d d} \int_{\Omega} d\mathbf{x} d\mathbf{y} H'^{-1}(\mathbf{x}, \mathbf{y}) \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}). \qquad (24) $$

We note that as $H'$ has positive eigenvalues we must have that $D_e < \bar{\kappa}$, and thus the late time diffusion is slower than one would naively expect. The formula of Eq. (24) agrees with two alternative derivations for $D_e$ in the literature. The first such derivation [3, 9] was obtained by generalizing a method proposed by Derrida [10] for discrete random walks to the continuum case. Another method is based on the Kubo formula corresponding to the case for diffusion in a potential and a random time change argument that relates this problem to the varying spatial diffusion coefficient problem under consideration here [3]. However, none of these approaches has to date been extended to examine finite time correction for the problem considered here.

The last term in Eq. (22) represents the finite time corrections to the asymptotic diffusion constant. However, the second term of this correction will decay exponentially quickly as $\exp(-\lambda_0 t)$ where $\lambda_0$ is the lowest non-zero eigenvalue of $H$. We expect that $\lambda_0 \sim 1/l^2$ and thus for $t > l^2$, the time after which a particle has typically diffused through a few periods, we will find

$$ D(t) = D_e + \frac{C}{t}, \qquad (25) $$

with $C$ given by

$$ C = \frac{1}{l^d d} \int_{\Omega} d\mathbf{x} d\mathbf{y} H'^{-2}(\mathbf{x}, \mathbf{y}) \nabla \kappa(\mathbf{x}) \cdot \nabla \kappa(\mathbf{y}). \qquad (26) $$

This means that the MSD behaves as

$$ \langle (\mathbf{X}_t - \mathbf{X}_0)^2 \rangle = 2dD_e t + 2dC. \qquad (27) $$

Thus the late time part the MSD can be fitted by a straight line, the slope of which is $2dD_e$ and the intercept with the $y$ axis is $2dC$.

Interestingly, the dependence of the operator $H$ on $\kappa$, means that, from Eqs. (24,26), multiplying $\kappa(x)$ by a constant leaves the constant $C$ unchanged, while the asymptotic diffusion constant $D_e$ is multiplied by the same constant. Thus, the constant $C$ is insensitive to the magnitude of the function $\kappa(\mathbf{x})$, but depends only on its relative variations.

The explicitly evaluation of $D_e$ and $C$ can only be carried out in one dimension and this will be done in the following section (III). However, in general we can write the effective diffusion constant as

$$ D_e = \bar{\kappa} - \frac{1}{l^d d} \int_{\Omega} d\mathbf{x} \nabla \kappa(\mathbf{x}) \cdot \mathbf{f}(\mathbf{x}), \qquad (28) $$

where $\mathbf{f}(\mathbf{x})$ is given by

$$ \mathbf{f}(\mathbf{x}) = \int_{\Omega} d\mathbf{y} H'^{-1}(\mathbf{x}, \mathbf{y}) \nabla \kappa(\mathbf{y}), \qquad (29) $$
<figure><img src="image_3.png" /><figcaption>FIG. 1: (color online) $D(t)$ estimated from the MSD in a numerical simulation of $5 \times 10^5$ Langevin particles, in one dimension, with local diffusion coefficient $\kappa(x)$ given by Eq. (45) with $\alpha = 0.8$ and $l = 4\pi l_0$ (continuous) black line (shown in units such that $\kappa_0 = 1$ and $l_0 = 1$). Shown as thick dashed (red -color online) late time prediction of $D(t)$ of the form predicted by Eq. (25) $D(t) = D_e + C/t$, where theory predicts that $D_e = 1$ and $C = 1.28$.</figcaption></figure>

and is thus the solution to the equation

$$ \nabla \cdot \kappa(\mathbf{x}) \nabla f(\mathbf{x}) = -\nabla \kappa(\mathbf{x}). \qquad (30) $$

The vector field $f$ is periodic with period $l$ in all coordinates and in addition is in the space orthogonal to $\psi_0(\mathbf{x}) \propto 1$, which means that

$$ \int_{\Omega} d\mathbf{x} \, f(\mathbf{x}) = 0. \qquad (31) $$

The constant $C$ is then simply given by

$$ C = \frac{1}{l^d} \int_{\Omega} d\mathbf{x} \, f^2(\mathbf{x}). \qquad (32) $$

Most results in the literature on effective diffusion constants are obtained via perturbative approaches [11–13]. Interestingly, the functional renormalisation group [14] can also be used to give approximate results which however agree with all known exact results for this problem. Exact results, where they exist, rely on the physical assumption that the self and collective diffusion constants are the same and thus the effective diffusion constant can be derived from the definition of $D_e$ given by

$$ \bar{\mathbf{j}} = -D_e \overline{\nabla \rho} \qquad (33) $$

relating the spatially averaged current to the average particle gradient. The formula given here agrees with this definition [3] and thus is a direct proof of this reasonable physical hypothesis. Using this definition for $D_e$ means that the one dimensional problem can be easily solved, while in certain cases a clever duality argument can give exact results in two dimensional systems [3, 15]. In principle the formulas given here can be analyzed via perturbation expansions to give results for $D_e$ and $C$.

We emphasize that the relations given here are exact, however, in general, the partial differential equation Eq. (30) in two dimensions and higher is not amenable to analytic solution. However, in principle, Eq. (30) should be easily amenable to numerical solution as the domain over which it holds is compact. The resulting integral expressions for $D_e$ (Eq. (28)) and $C$ (Eq. (32)) would then be straightforward to compute numerically.
### III. EXACT RESULTS IN ONE DIMENSION

In one dimension the function $f$ can be easily computed and one finds the general solution of Eq. (30):

$$f(x) = -x + A + B \int_{0}^{x} \frac{dy}{\kappa(y)}, \qquad (34)$$

where $A$ and $B$ are constants. The condition of periodicity in $l$ then gives

$$B = \overline{\kappa^{-1}}^{-1}, \qquad (35)$$

i.e. $B$ is the Harmonic mean of $\kappa$ under spatial averaging. The constant $A$ is then given by the condition $\int_{0}^{l} dx f(x) = 0$ which gives

$$A = \frac{1}{l} \int_{0}^{l} dx \left[ x - \overline{\kappa^{-1}}^{-1} \int_{0}^{x} dy \kappa^{-1}(y) \right]. \qquad (36)$$

If we define the function

$$r(x) = -x + \overline{\kappa^{-1}}^{-1} \int_{0}^{x} dy \kappa^{-1}(y), \qquad (37)$$

then the function $f(x)$ can be written as

$$f(x) = r(x) - \bar{r}, \qquad (38)$$

where

$$\bar{r} = \frac{1}{l} \int_{0}^{l} dx r(x). \qquad (39)$$

Inserting this expression into Eq. (32), we obtain that the constant $C$ is equal to the variance of the function $r(x)$ with respect to spatial averaging:

$$C = \frac{1}{l} \int_{0}^{l} dx (r(x) - \bar{r})^{2}. \qquad (40)$$

Using the solution for $f$ we find that the effective diffusion constant in one dimension is given by

$$D_e = \overline{\kappa^{-1}}^{-1}, \qquad (41)$$

that is to say the Harmonic mean. This result for the one-dimensional diffusion constant is well known. The result is usually established by examining the steady state current in a system subjected to an applied external gradient, mathematically the diffusion constant is the same as the effective conductivity, effective dielectric permittivity or effective porosity of random media when $\kappa(\mathbf{x})$ represents, respectively, a spatially varying conductivity, dielectric constant or porosity (appearing in Darcy's law for incompressible flow in porous media) [3].

If we make the periodicity of the diffusion constant explicit by writing

$$\kappa(x) = K \left( \frac{2\pi x}{l} \right) \qquad (42)$$

we find that the $l$ dependence of $C$ can be made explicit as

$$C = \frac{l^2}{(2\pi)^3} \int_{0}^{2\pi} dz (R(z) - \bar{R})^2, \qquad (43)$$

where

$$R(z) = -y + \frac{1}{\overline{K^{-1}}} \int_{0}^{y} dy' K^{-1}(y'), \qquad (44)$$
<figure><img src="image_4.png" /><figcaption>FIG. 2: (color online) $D(t) - D_e$ (here $D_e = \kappa_0$) estimated from the MSD in the same numerical simulation as Fig. (1) (open circles) compared at all times with a numerical evaluation of the Kubo formula Eq. (22) (solid line- red online)(again shown in units such that $\kappa_0 = 1$ and $l_0 = 1$). Also shown is the asymptotic analytical late time correction $C/t$ to $D(t) - D_e$ (dashed line- blue online) as well as the short time value $D(0) - D_e = \bar{\kappa} - D_e$ (dotted line- green on line).</figcaption></figure>

where $\bar{K}^{-1} = \int_0^{2\pi} dy K^{-1}(y)/2\pi$ is again the inverse of the harmonic mean of the diffusion coefficient and $\bar{R} = \int_0^{2\pi} dy R(y)/2\pi$. We thus see, that in common with the case of diffusion in a periodic potential [7], the constant $C$ is proportional to $l^2$. Indeed it is not hard to see, by dimensional analysis or rescaling of the diffusion equation, that this must be the case.

As an concrete example, consider the system with local diffusion constant given by

$$ \kappa(x) = \frac{\kappa_0}{1 + \alpha \cos(\frac{2\pi x}{l})}, \qquad (45) $$

where we must have $|\alpha| < 1$. Here we find that

$$ D_e = \kappa_0 \text{ and } C = \frac{\alpha^2 l^2}{8\pi^2}. \qquad (46) $$

As a test of our theory we have numerically simulated an ensemble of $5 \times 10^5$ particles obeying a discretized form of Eq. (4) with $dt = 0.0001$. The particles are started in the interval $[0, l]$, with $l = 4\pi$, with a uniform distribution. The values $\kappa_0 = 1$ and $\alpha = 0.8$ are taken in this numerical example. The estimate of $D(t)$ is then made by estimating the ensemble average in Eq. (2) by averaging over the ensemble of particle trajectories. We see in Fig. (1) that at late times the asymptotic theory developed here predicts the correct behavior of $D(t)$, even for times where the value of $D(t)$ is still significantly larger that its asymptotic late time limit $D_e$.

The Kubo formula Eq (22) in principle gives the behavior of the MSD at all times assuming the initial distribution is the equilibrium one. Even in one dimension the full time dependent problem is not analytically tractable. However numerically one can construct a discrete approximation to the operator $H$ on a lattice subspace of $[0, l]$. The resulting operator can then be numerically analyzed to obtain the eigenfunctions $\psi_\lambda$ and corresponding eigenvalues $\lambda$. The last term in Eq. (22) containing the exponentially suppressed terms can then be evaluated numerically. We have carried out this numerical computation for the system we simulated to obtain Fig. (1) and the comparison is shown in Fig (2). We see that the full Kubo formula agrees perfectly with the results of the numerical simulation.
IV. CONCLUSIONS

We have analyzed the temporal behavior of the MSD of Brownian particles diffusing in a medium with a spatially varying local diffusion constant. The assumption behind the calculation is that at the starting time the system is in equilibrium, with a uniform distribution, and we assume that the time at which the diffusion constant attains its final asymptotic limit is always smaller than the time for the particle to feel the finite size of the system. Clearly by taking large enough systems this limit can always be achieved. At short times, as the position distribution is uniform, the diffusion constant is clearly given by $\bar{\kappa}$. At late times the effective diffusion constant is smaller, $D_e < \bar{\kappa}$. This result is not as easy to understand intuitively as the slowing down of diffusion in an external potential, where trapping in the potential's minima clearly provides the mechanism by which dispersion is slowed down. For a periodic local diffusivity, we have shown that the late time correction to the time dependent diffusion constant, as defined by Eq. (2) decays as $C/t$ and we have given a general expression for the coefficient $C$. In one dimension, the effective diffusion constant only depends on the composition of the local diffusivity via its harmonic mean. The constant $C$ depends however on the spatial distribution of $\kappa$, thus two systems with the same value of $D_e = 1/\kappa^{-1}$ will in general have differing values of $C$. The explicit analytical results given here are restricted to the case of one dimension. However, in two dimensional random systems if we write $\kappa(\mathbf{x}) = \kappa_0 \exp(\phi(\mathbf{x}))$ and if the corresponding field $\phi$ is statistically identical to $-\phi$ then the effective diffusion constant is given by $D_e = \kappa_0$ [3, 15]. This result is derived by a duality argument in two dimensions and it would be interesting to see if this argument could be extended to compute $C$.

In principle the method here could also be extended to analyze more general diffusion equations, for example one could extend the study here, and in Ref. [7], by considering diffusion in a periodic potential with a periodic diffusion constant in such a way that the underlying equilibrium distribution is the Gibbs-Boltzmann one (as detailed balance is obeyed). Finally, in this work, and the study of diffusion in a periodic potential [7], the steady state has an equilibrium distribution without a current. It would be very interesting to examine the temporal dependence of the MSD in systems having a steady state with current, for instance Brownian particles in a tilted potential [6].

[1] S. N. Majumdar, Current Science **89**, 2076 (2005)

[2] J.-P. Bouchaud and A. Georges, Phys. Rep. **195**, 127 (1990)

[3] D.S.Dean, I.T. Drummond and R.R. Horgan, *J. Stat. Mech.* **7** P07013 (2007)

[4] D. Boyer and D.S. Dean, J. Phys. A: Math. Gen. **44**, 335003 (2011); D. Boyer, D.S. Dean, C. Mejía-Monasterio and G. Oshanin, Phys. Rev. E **85**, 031136 (2012); *ibid*, Phys. Rev. E **86**, 060101(R) (2012); *ibid*, J. Stat. Mech. P04017 (2013).

[5] S. Lifson and J. L. Jackson, J. Chem. Phys. **36**, 2410 (1962); P.G. de Gennes, J. Stat. Phys. **12**, 463 (1975); R. Zwanzig, Proc. Nat. Acad. Sci. **85**, 2029 (1988).

[6] P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J.M. Rubi and A. Pérez-Madrid, Phys. Rev. Lett. **87**, 010602 (2001); M. Khoury, A. M. Lacasta, J. M. Sancho and K. Lindenberg, Phys. Rev. Lett. **106**, 090602 (2011); M. S. Simon, J. M. Sancho and K. Lindenberg, Phys. Rev. E **88**, 062105 (2013).

[7] D.S. Dean and G. Oshanin, Phys. Rev. E **90**, 022112 (2014).

[8] B. Oksendal, Stochastic Differential Equations, Springer-Verlag, Heidelberg (2013).

[9] D.S. Dean, I.T. Drummond and R.R. Horgan, J. Phys. A **30**, 385 (1997).

[10] B. Derrida, J. Stat. Phys. **31**, 433 (1983).

[11] I.T. Drummond and R R Horgan, J. Phys. A: Math. Gen. **20**, 4661 (1987).

[12] A. De Wit Phys. Fluids **7**, 2553 (1995).

[13] B. Abramovich and P. Indelman J. Phys. A: Math. Gen. **28**, 693 (1995).

[14] D.S. Dean and C. Touya J. Phys. A: Math. Theor. **41**, 335002 (2008).

[15] A.M. Dykhne Sov. Phys.-JETP **32**, 63 (1971); A.M. Dykhne Sov. Phys.-JETP **32**, 348 (1971).