# Kubo formulas for dispersion in heterogeneous periodic nonequilibrium systems

T. Guérin, David S. Dean

► To cite this version:

T. Guérin, David S. Dean. Kubo formulas for dispersion in heterogeneous periodic nonequilibrium systems. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2015, 92 (6), pp.062103 (1-11). 10.1103/PhysRevE.92.062103. hal-01252573

HAL Id: hal-01252573

https://hal.science/hal-01252573

Submitted on 7 Jan 2016

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_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Kubo formulas for dispersion in heterogeneous periodic non-equilibrium systems

T. Guérin¹ and D. S. Dean¹

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

We consider the dispersion properties of tracer particles moving in non-equilibrium heterogeneous periodic media. The tracer motion is described by a Fokker-Planck equation with arbitrary spatially periodic (but constant in time) local diffusion tensors and drift, eventually with the presence of obstacles. We derive a Kubo-like formula for the time dependent effective diffusion tensor valid in any dimension. From this general formula, we derive expressions for the late time effective diffusion tensor and drift in these systems. In addition, we find an explicit formula for the late finite time corrections to these transport coefficients. In one dimension, we give closed analytical formula for the transport coefficients. The formulas derived here are very general and provide a straightforward method to compute the dispersion properties in arbitrary non-equilibrium periodic advection-diffusion systems.

PACS numbers: 05.60.Cd, 02.50.Ey, 05.10.Gg, 05.40.a

# I. INTRODUCTION

Fokker-Planck (FP) equations, and their associated stochastic differential equations (SDE) arise in a huge variety of physical systems [1–3]. In fluid mechanics and hydrology, such equations describe the transport of passive scalars advected by a fluid and locally dispersed by a microscopic molecular diffusivity. They also describe the diffusion of colloids and polymers in soft matter physics, as well as charge and energy transport in solid state physics. In many systems, damping due to the environment is very strong and inertial effects are thus negligible. In such systems, the evolution of the probability density function (pdf) $p(\mathbf{x}, t)$ of a tracer particle at position $\mathbf{x}$ and time $t$ obeys the FP equation

$$ \partial_t p = \partial_{x_i} \{\partial_{x_j} [\kappa_{ij}(\mathbf{x})p] - u_i(\mathbf{x})p\} = -H_{\mathbf{x}}p, \quad (1) $$

where the Einstein convention on indexes summation is used, and the indexes $i, j$ run from 1 to the spatial dimension $d$. The time independent tensor $\kappa_{ij}(\mathbf{x})$ is the local diffusivity tensor and can, in general, vary in space in media of non-uniform composition, since diffusivity is depends on the local material properties of the system. The local diffusivity may also vary in space due to temperature inhomogeneities, or in the case of particle diffusion in viscous fluid flows close to solid boundaries [4]. The term $u_i(\mathbf{x})$ represents a drift term generated by forces such as buoyancy forces or electromagnetic forces acting on colloids, fluid flows, as well as thermodynamic effects, such as thermodiffusion in the presence of temperature gradients [5]. The presence of impenetrable obstacles can also be treated by introducing surfaces at which there is no normal flux,

$$ n_i(\mathbf{x})\{\partial_{x_j} [\kappa_{ij}(\mathbf{x})p] - u_i(\mathbf{x})p\} = 0, \quad \mathbf{x} \in S, \quad (2) $$

where $n_i$ represents the normal vector to the surface $S$ of the obstacles.

Equation (1) describes the motion of tracer particles at small scales in a wide range of heterogeneous media. At a

coarse-grained level, the motion is described by effective transport coefficients characterizing the average flux of particles and the relative spreading with time of initially close particles. From a theoretical point of view, the dispersion of tracer particles is defined in terms of two quantities, the mean displacement $\mu_i(t)$ during a time $t$ in the direction $i$

$$ \mu_i(t) = \langle X_i(t) - X_i(0) \rangle, \quad (3) $$

(where $\mathbf{X}(t)$ denotes the position of a tracer particle and $\langle \cdots \rangle$ denotes ensemble averaging), and the correlation of the dispersion in the directions $i$ and $j$ defined by

$$ \sigma_{ij}(t) = \langle [X_i(t) - X_i(0)][X_j(t) - X_j(0)] \rangle - \mu_i(t)\mu_j(t). \quad (4) $$

The quantity $\sigma_{ij}(t)$ thus characterizes the dispersion of a cloud of particles about its mean position. Of particular interest is the large time behavior of the dispersion properties,

$$ \mu_i(t) \approx V_i t, \quad (t \to \infty) \quad (5) $$

$$ \sigma_{ij}(t) \approx 2 t (D_{ij} + C_{ij}/t), \quad (t \to \infty), \quad (6) $$

where $V_i$ is the average particle drift, $D_{ij}$ is the late time effective diffusion tensor, and $C_{ij}$ gives the leading order corrections and provides information on how fast the system approaches the diffusive limit. We will see that the asymptotic form (6) is valid if the medium is spatially periodic. The quantities $V_i$ and $D_{ij}$ characterize the late time, large scale transport properties of the system and are routinely measured in experimental systems by single particle tracking methods [6] or by ensemble based measurements such as fluorescence recovery after photo-bleaching (FRAP) [7]. Effective transport coefficient are also important for estimating the spread of pollutants, filtration, mixing and chemical reaction times [8].

Transport and dispersion properties in non-equilibrium periodic systems have been investigated at length in various contexts. In a pioneering, classical work, Taylor [9] considered the dispersion of particles moving in viscous
flows in cylindrical containers, and showed that the effective diffusion coefficient could be orders of magnitudes larger than the molecular diffusivity. Effective transport in incompressible fluid flows were also investigated at length, such as in the case of diffusion in Rayleigh-Bénard convection cells [10–12], in frozen turbulent flows [13] or in porous media [14–20] using homogenization theory. However, these approaches have been restricted so far to incompressible flows and it is not clear how to generalize them. In another context, results were also derived, with the methods of statistical physics, for particles diffusing in periodic potentials when the diffusivity is constant [21–25]; in here the potential barrier slows down the diffusion as the particle tends to become trapped in local mining. The non-equilibrium problem of transport in one-dimensional potential in the presence of external constant force (still with uniform microscopic diffusivity) was also considered [26–29]; this system is predicted to exhibit a huge increase of the effective diffusion near a critical tilting force. Results have also been obtained in the more general case where the noise amplitude is a periodic function of position [30–32], notably in the context of the modeling of transport in periodic channels within the Fick-Jacobs approximation (where the full Fokker-Planck equation is approximated by an effective one dimensional one). In many of the above cases, results only exist in one dimension, however their generalization to higher dimensions is implicit in the general results we will present here.

In this paper, we present a general formalism to compute the time-dependent dispersion tensor $\sigma_{ij}(t)$ (section III), the late time diffusion tensor $D_{ij}$ (section IV) and the corrections $C_{ij}$ describing the approach to the diffusive limit (section V). A summary of the results is given in section II for the reader who might not be interested in the technical aspects of the calculations. The results derived are valid for any system described by Eq. (1) for periodic $\kappa_{ij}(\mathbf{x})$ and $u_i(\mathbf{x})$; they encapsulate all known results for this problem and also generalize them to more complex situations. We will explicitly compare them to the expressions of the literature in various cases (section VI), and we will also give an explicit fully analytical solution for the one dimensional problem (section VII). A brief derivation of our main results has already appeared in Ref. [33], where they were used to analyze diffusion in a system with a spatially varying diffusivity in the presence of a uniform applied force. The present paper extends this recently developed approach, new features are the computation of the time-dependent dispersion properties (instead of only the late time diffusive limit) and that the incorporation of the effect of impenetrable obstacles with no-flux boundary conditions. Furthermore, the derivation presented here is different from that appearing in Ref. [33], as it is not based on the SDE formulation used there. Finally, the present paper contains explicit expressions for the first temporal corrections $C_{ij}$ to the late time diffusion tensor. We will show that the $1/t$ decay of the late time correction [Eq. (6)] occurs in all

dimensions. This universality can be understood as, for periodic systems, the Fokker-Planck equation restricted to the basic cell $\Omega$ will in general have a gap and not a continuum (whose density would be dependent on the spatial dimension) of eigenvalues near zero.

## II. NOTATIONS AND SUMMARY OF RESULTS

We assume that the space can be divided into individual cells $\Omega$, and we that the diffusivity tensor $\kappa_{ij}(\mathbf{x})$ and the drift $u_i(\mathbf{x})$ are periodic in space, and thus satisfy

$$ \kappa_{ij}(\mathbf{x}) = \kappa_{ij}(\mathbf{x} + \mathbf{k}), \quad u_i(\mathbf{x}) = u_i(\mathbf{x} + \mathbf{k}), \quad (7) $$

for all vectors $\mathbf{k}$ joining the centers of the cells composing the periodic structure. Denoting by $\mathbf{X}(t)$ the position of an individual particle at time $t$, we define the propagator $p(\mathbf{x}, t|\mathbf{y})$ of the stochastic process $\mathbf{X}(t)$ in infinite space; $p(\mathbf{x}, t|\mathbf{y})$ is thus the probability density to observe a particle at position $\mathbf{x}$ at time $t$ starting from $\mathbf{y}$ at time 0. Mathematically, $p$ is the solution of the equation in infinite space

$$ (\partial_t + H_{\mathbf{x}})p(\mathbf{x}, t|\mathbf{y}) = 0; \quad p(\mathbf{x}, 0|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}). \quad (8) $$

We also introduce the process $\tilde{\mathbf{X}}(t)$ as the position of the particle $\mathbf{X}$ "modulo $\Omega$", $\tilde{\mathbf{X}}$ is defined as the position of the particle in a reference cell $\Omega_0$ obtained after an integer number of translations of $\mathbf{X}$ along the vectors joining the centers of the cells composing the periodic structure. We introduce the propagator $P(\mathbf{x}, t|\mathbf{y})$ for the process modulo $\Omega$, $P$ thus represents the probability density, starting at position $\mathbf{y}$ at time 0, to be at time $t$ at position $\mathbf{x}$ modulo $\Omega$. $P$ is the solution of the FP equation (1) with periodic boundary conditions on the boundaries of $\Omega$ and initial value $P(\mathbf{x}, 0|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y})$. There is an obvious relation between these propagators,

$$ P(\mathbf{x}, t|\mathbf{y}) = \sum_{\mathbf{k}} p(\mathbf{x} + \mathbf{k}, t|\mathbf{y}), \quad (9) $$

where $\mathbf{k}$ represents the vectors joining the center of the cell containing $\mathbf{x}$ to the centers of all the other cells of the periodic structure. The value of $P$ at late times is denoted by $P_0$ and is the steady-state distribution of particles in one period (modulo translations along the lattice vectors). In contrast, the propagator in infinite space $p$ does not tend to a steady state, but becomes Gaussian at long times, with a variance growing linearly as $t$, representing the spreading of the particle density. In what follows, we always assume that at $t = 0$ we are in the steady state, that is to say that the process $\tilde{\mathbf{X}}(t = 0)$ modulo $\Omega$ has the distribution $P_0$.

Before entering the details of the calculations, let us give briefly our main results. First, we will show that the temporal Laplace transform $\hat{\sigma}_{ij}(s) = \int_0^\infty \sigma_{ij}(t)e^{-st}dt$ of
the dispersion tensor is given by

$$ \hat{\sigma}_{ij}(s) = \frac{2}{s^2} \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \\ - \frac{1}{s^2} \iint_{\Omega} d\mathbf{x} d\mathbf{y} [\tilde{u}_i(\mathbf{x}) \tilde{u}_j^*(\mathbf{y}) + \tilde{u}_j(\mathbf{x}) \tilde{u}_i^*(\mathbf{y})] \hat{P}'(\mathbf{x}, s|\mathbf{y}) P_0(\mathbf{y}). \quad (10) $$

In this relatively compact formula, $\hat{P}'(\mathbf{x}, s|\mathbf{y})$ represents the Laplace transform of $P(\mathbf{x}, t|\mathbf{y}) - P_0(\mathbf{x})$. We have also defined an effective velocity field $\tilde{u}$ that takes into account the presence of obstacles (if present)

$$ \tilde{u}_i(\mathbf{x}) = u_i(\mathbf{x}) - n_j(\mathbf{x})\kappa_{ij}(\mathbf{x})\delta_S(\mathbf{x}), \quad (11) $$

where $\delta_S$ is a surface delta function; to be clear with the notations we precise that

$$ \int_{\Omega} d\mathbf{x} n_j(\mathbf{x}) \kappa_{ij}(\mathbf{x}) \delta_S(\mathbf{x}) \psi(\mathbf{x}) = \int_{S_{\Omega}} dS_j \kappa_{ij}(\mathbf{x}) \psi(\mathbf{x}) \quad (12) $$

for any function $\psi$, with the infinitesimal normal surface vector $dS_i$ and the unit normal vector $n_i$ oriented towards the interior of the obstacles. Finally, we have also introduced the drift field $u_i^*$ which can be interpreted as the drift field for the particles upon time reversal, and is given by

$$ u_i^*(\mathbf{x}) = u_i(\mathbf{x}) - 2 \frac{J_i(\mathbf{x})}{P_0(\mathbf{x})}, \quad (13) $$

where we have introduced the current in the steady state $J_i$

$$ J_i(\mathbf{x}) = u_i(\mathbf{x}) P_0(\mathbf{x}) - \partial_{x_k} [\kappa_{ik}(\mathbf{x}) P_0(\mathbf{x})]. \quad (14) $$

The notation $\tilde{u}_i^*(\mathbf{x})$ denotes that we have added surface drift terms which prevent the particles from entering the obstacles,

$$ \tilde{u}_i^*(\mathbf{x}) = u_i^*(\mathbf{x}) - n_j(\mathbf{x})\kappa_{ij}(\mathbf{x})\delta_S(\mathbf{x}). \quad (15) $$

Equation (10) states clearly that the large scale dispersion properties can be obtained from the properties of the propagator (with periodic boundary conditions), which is used to average the drift fields with and without time-reversal; this is the first result of the paper.

The second main result is an expression for the late time diffusion tensor $D_{ij}$. Taking the large time limit of Eq. (10) leads to

$$ D_{ij} = \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \\ + \frac{1}{2} \int_{\Omega} d\mathbf{x} \, [u_i(\mathbf{x}) f_j(\mathbf{x}) + u_j(\mathbf{x}) f_i(\mathbf{x})] \\ - \frac{1}{2} \int_{S_{\Omega}} dS_l(\mathbf{x}) [\kappa_{il}(\mathbf{x}) f_j(\mathbf{x}) + \kappa_{jl}(\mathbf{x}) f_i(\mathbf{x})], \quad (16) $$

where the function $f_j$ is defined by

$$ f_j(\mathbf{x}) = - \int_{\Omega} d\mathbf{y} \, P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \tilde{u}_j^*(\mathbf{y}) \quad (17) $$

and

$$ G(\mathbf{x}|\mathbf{y}) = \int_{0}^{\infty} dt [P(\mathbf{x}, t|\mathbf{y}) - P_0(\mathbf{x})]. \quad (18) $$

$G$ is the pseudo-Green function [34] of the operator $H$. For computational purposes, it is useful to note that $f_i$ can be obtained by solving the non-homogeneous linear partial differential equation

$$ H_{\mathbf{x}} f_i(\mathbf{x}) = - u_i^*(\mathbf{x}) P_0(\mathbf{x}) + P_0(\mathbf{x}) \int_{\Omega} d\mathbf{y} \, u_i^*(\mathbf{y}) P_0(\mathbf{y}) \\ - P_0(\mathbf{x}) \int_{S_{\Omega}} dS_l(\mathbf{y}) P_0(\mathbf{y}) \kappa_{li}(\mathbf{y}) \quad (19) $$

with periodic boundary conditions and also with the condition

$$ n_j\{u_j(\mathbf{x})f_i(\mathbf{x}) - \partial_{x_k}[\kappa_{kj}(\mathbf{x})f_i(\mathbf{x})]\} = -n_j\kappa_{ji}(\mathbf{x})P_0(\mathbf{x}) \quad (20) $$

at the surface of the obstacles, together with the orthogonality condition

$$ \int_{\Omega} d\mathbf{x} \, f_i(\mathbf{x}) = 0. \quad (21) $$

The above expressions are exact, general and can be used to compute the dispersion properties in any advection-diffusion problem with periodic diffusivity tensor and drift fields, in any dimension, in the presence of obstacles with reflecting boundaries. We will show that our expression for $D_{ij}$ contains existing expressions for the effective diffusion tensor in porous media for particles transported by incrompressible fluids [14, 35], diffusion in periodic potentials [36], and diffusion in non-equilibrium one dimensional systems [26, 28, 32].

One may also derive an alternative expression for $D_{ij}$ by defining a function $f_i^*$ as

$$ f_i^*(\mathbf{y}) = - \int_{\Omega} d\mathbf{x} \, P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \tilde{u}_i(\mathbf{x}), \quad (22) $$

in terms of which the effective diffusion tensor reads

$$ D_{ij} = \int_{\Omega} d\mathbf{y} \, \kappa_{ij}(\mathbf{y}) P_0(\mathbf{y}) \\ + \frac{1}{2} \int_{\Omega} d\mathbf{y} \, [u_i^*(\mathbf{y}) f_j^*(\mathbf{y}) + u_j^*(\mathbf{y}) f_i^*(\mathbf{y})] \\ - \frac{1}{2} \int_{S_{\Omega}} dS_l(\mathbf{y}) [\kappa_{il}(\mathbf{y}) f_j^*(\mathbf{x}) + \kappa_{jl}(\mathbf{y}) f_i^*(\mathbf{y})]. \quad (23) $$

This function $f_i^*$ can also be obtained as the solution of the partial differential equations

$$ H_{\mathbf{y}}^* f_i^*(\mathbf{y}) = - u_i(\mathbf{y}) P_0(\mathbf{y}) + P_0(\mathbf{y}) \int_{\Omega} d\mathbf{x} \, [u_i(\mathbf{x}) P_0(\mathbf{x})] \\ - P_0(\mathbf{y}) \int_{S_{\Omega}} dS_l(\mathbf{x}) P_0(\mathbf{x}) \kappa_{li}(\mathbf{x}). \quad (24) $$
<figure><img src="image_2.png" /><figcaption>FIG. 1: (a) Examples of stochastic trajectories of two tracer particles diffusing in a heterogeneous medium (here a 2D array of circular obstacles) under a force F (arrow) that biases the motion to the right. (b) Elements of the dispersion tensor for the system described in (a), calculated from our theoretical formulas Eqs. (16,19) (full and dashed lines) vs results of stochastic simulations (symbols). Parameters: disk radius R = 0.25L with L the length of the side of the periodic cell. Results for $D_{ij}$ are expressed in units of the (constant) diffusion coefficient of the tracers outside the obstacles.</figcaption></figure>

where we have denoted $H^*$ the transport operator for the time-reversed process,

$$H_{\mathbf{y}}^{*}(\cdot) = \partial_{y_i} [u_i^{*}(\mathbf{y})(\cdot)] - \partial_{y_i} \partial_{y_k} [\kappa_{ij}(\mathbf{y})(\cdot)]. \quad (25)$$

The function $f_i^*$ has periodic boundary conditions, and satisfies at the surface of the obstacles

$$n_j(\mathbf{y})\{u_j^*(\mathbf{y})f_i^*(\mathbf{y}) - \partial_{y_k}[\kappa_{kj}(\mathbf{y})f_i^*(\mathbf{y})]\} = -n_j\kappa_{ji}(\mathbf{y})P_0(\mathbf{y}), \quad (26)$$

and also satisfies the orthogonality condition

$$\int_{\Omega} dy f_i^*(y) = 0. \quad (27)$$

Thus, we have two formulas for the effective diffusion tensor, which can be straightforwardly evaluated numerically by solving a set of elliptic partial differential equation on the domain $\Omega$.

Finally, we will show that the tensor $C_{ij}$ describing the approach to the diffusive limit can be obtained from

$$C_{ij} = \frac{1}{2} \int_{\Omega} dx \frac{f_i(\mathbf{x}) f_j^*(\mathbf{x}) + f_j(\mathbf{x}) f_i^*(\mathbf{x})}{P_0(\mathbf{x})}. \quad (28)$$

Although numerous expressions exist in the literature for $D_{ij}$ for particular cases, the correction tensor $C_{ij}$ has to our knowledge only been studied for diffusion in a periodic potential [37] or with periodic diffusivity [38], the above expressions enable its computation in a systematic way and at a computational cost equivalent to that required to calculate $D_{ij}$.

Our results can be applied to study tracer particles diffusing in a 2D array of circular impenetrable obstacles (hard disks) [Fig.1(a)] and subjected to an external force. For this system, the elements of the effective diffusion tensor match perfectly with the results of stochastic simulations (the algorithm is presented in appendix A).

### III. DERIVATION OF THE KUBO FORMULA FOR THE TIME DEPENDENT DISPERSION

Here we describe the derivation of the Kubo formula for the dispersion tensor. We first briefly compute the average displacement in the direction $i$ during a time $t$, defined in Eq. (3), which is given by

$$\mu_i(t) = \int_{\mathbb{R}^d} dx \int_{\Omega} dy P_0(\mathbf{y}) p(\mathbf{x}, t|\mathbf{y})(x_i - y_i). \quad (29)$$

Note that the integral over the whole space $\mathbb{R}^d$ denotes the integral over the volume not occupied by obstacles (where $p$ does not vanish). Taking the time derivative of the above expression, and using Eq. (1), we obtain

$$\partial_t \mu_i(t) = \int_{\mathbb{R}^d} dx \int_{\Omega} dy P_0(\mathbf{y}) \{-\partial_{x_j} [u_j(\mathbf{x}) p(\mathbf{x}, t|\mathbf{y})] + \partial_{x_j} \partial_{x_k} [\kappa_{kj}(\mathbf{x}) p(\mathbf{x}, t|\mathbf{y})]\} (x_i - y_i). \quad (30)$$

Integrating by parts over $\mathbf{x}$ and applying the divergence theorem yields

$$\partial_t \mu_i(t) = \int_{\mathbb{R}^d} dx \int_{\Omega} dy P_0(\mathbf{y}) \{u_i(\mathbf{x}) p(\mathbf{x}, t|\mathbf{y}) - \partial_{x_k} [\kappa_{ik}(\mathbf{x}) p(\mathbf{x}, t|\mathbf{y})]\}, \quad (31)$$

where we have taken into account the reflecting boundary conditions at the obstacles, so that no surface integral has appeared. Decomposing the integral over $\mathbf{x}$ in integrals over each unit cell, using the periodicity property (7) of $u_i$ and $\kappa_{ij}$ and the definition (9), we get

$$\partial_t \mu_i(t) = \int_{\Omega} dx \int_{\Omega} dy P_0(\mathbf{y}) \{u_i(\mathbf{x}) P(\mathbf{x}, t|\mathbf{y}) - \partial_{x_k} [\kappa_{ik}(\mathbf{x}) P(\mathbf{x}, t|\mathbf{y})]\}. \quad (32)$$

We recall that, since $P_0$ is the steady state particle distribution, it satisfies

$$\int_{\Omega} dy P_0(\mathbf{y}) P(\mathbf{x}, t|\mathbf{y}) = P_0(\mathbf{x}). \quad (33)$$
Using this property, and again applying the divergence theorem, (32) becomes

$$ \partial_t \mu_i(t) = \int_{\Omega} d\mathbf{x} \, u_i(\mathbf{x}) P_0(\mathbf{x}) \\ - \int_{S_{\Omega}} dS_j(\mathbf{x}) \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \} = V_i, \quad (34) $$

where $S_Ω$ represents the surface of the obstacles inside the cell $Ω$, and the infinitesimal surface vector $dS_j$ is oriented towards the interior of the obstacles. Eq. (34) shows that, when the initial conditions are those of the steady state, the average drift $V_i$ is constant in time. This is in agreement with the classic result of Stratonovich [39].

Let us now derive a Kubo formula for the time dependent dispersion tensor $\sigma_{ij}$ [defined in Eq. (4)]. Define $\psi_{ij}$ to be

$$ \sigma_{ij}(t) = \psi_{ij}(t) - \mu_i(t)\mu_j(t). \quad (35) $$

By definition, the expression for $\psi_{ij}$ is

$$ \psi_{ij}(t) = \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y) p(\mathbf{x}, t|\mathbf{y})(x_i - y_i)(x_j - y_j). \quad (36) $$

We denote by $\hat{\phi}(s) = \int_0^\infty dt \phi(t) e^{-st}$ the temporal Laplace transform of any function $\phi$. The Laplace transform of the FP equation (1) is given by

$$ s\hat{p}(\mathbf{x}, s|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}) \\ + \partial_{x_i} \{\partial_{x_j}[\kappa_{ij}(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y})] - u_i(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y})\}. \quad (37) $$

Using the above equality, the expression for $\psi_{ij}$ becomes

$$ \hat{\psi}_{ij}(s) = \frac{1}{s} \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y)(x_i - y_i)(x_j - y_j) \\ \partial_{x_k} \{\partial_{x_l}[\kappa_{kl}(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y})] - u_k(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y})\}. \quad (38) $$

We remark that

$$ \partial_{x_k} [(x_i - y_i)(x_j - y_j)] = \delta_{ki}(x_j - y_j) \\ + \delta_{kj}(x_i - y_i), \quad (39) $$

so that the integration by parts over $\mathbf{x}$ in Eq. (38) leads to an expression of the form

$$ \hat{\psi}_{ij}(s) = \hat{E}_{ij}(s) + \hat{E}_{ji}(s) \quad (40) $$

with

$$ \hat{E}_{ij}(s) = \frac{1}{s} \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y)(x_j - y_j) \\ \{u_i(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y}) - \partial_{x_k}[\kappa_{ki}(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y})]\}. \quad (41) $$

Integrating by parts the term containing $\kappa_{ik}$ leads to

$$ \hat{E}_{ij}(s) = \frac{1}{s} \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y)(x_j - y_j)u_i(\mathbf{x})\hat{p}(\mathbf{x}, s|\mathbf{y}) \\ - \frac{1}{s} \int_S dS_l(\mathbf{x}) \int_{\Omega} dy P_0(y)\kappa_{il}(\mathbf{x})(x_j - y_j)\hat{p}(\mathbf{x}, s|\mathbf{y}) \\ + \frac{1}{s} \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y)\hat{p}(\mathbf{x}, s|\mathbf{y})\kappa_{ij}(\mathbf{x}). \quad (42) $$

Using the periodicity property of $\kappa$ and the Laplace transform of the steady state property (33), we obtain

$$ \int_{\mathbb{R}^d} d\mathbf{x} \int_{\Omega} dy P_0(y) \hat{p}(\mathbf{x}, s|\mathbf{y}) \kappa_{il}(\mathbf{x}) = s^{-1} \int_{\Omega} d\mathbf{x} P_0(\mathbf{x}) \kappa_{il}(\mathbf{x}) = s^{-1} \langle \kappa_{ij} \rangle_0, \quad (43) $$

where $\langle \cdots \rangle_0$ represents the average over $P_0$. Decomposing the integral appearing in (42) over $\mathbf{x}$ on all the individual cells of the periodic structure, and changing of variable $\mathbf{x} \to \mathbf{x} + \mathbf{k}$, $y \to y + \mathbf{k}$ in each of them (where $\mathbf{k}$ is the lattice vector such that $\mathbf{x} + \mathbf{k}$ is inside the cell $Ω$), and summing over all lattice vectors $\mathbf{k}$ again, we see that we can exchange the integration domains between $\mathbf{x}$ and $\mathbf{y}$, leading to

$$ \hat{E}_{ij}(s) = \frac{1}{s} \int_{\Omega} d\mathbf{x} u_i(\mathbf{x}) \hat{B}_j(\mathbf{x}, s) \\ - \frac{1}{s} \int_{S_{\Omega}} dS_l(\mathbf{x}) \kappa_{il}(\mathbf{x}) \hat{B}_j(\mathbf{x}, s) + \frac{\langle \kappa_{ij} \rangle_0}{s^2}, \quad (44) $$

where we have defined

$$ \hat{B}_j(\mathbf{x}, s) = \int_{\mathbb{R}^d} dy \hat{p}(\mathbf{x}, s|\mathbf{y}) P_0(y)(x_j - y_j). \quad (45) $$

Now comes the key point of our derivation. In order to calculate $B_j$, we first consider another probability distribution $q$ defined as

$$ q(\mathbf{y}, t|\mathbf{x}) = \frac{p(\mathbf{x}, t|\mathbf{y}) P_0(\mathbf{y})}{P_0(\mathbf{x})}. \quad (46) $$

Using Bayes' theorem, we see that $q(\mathbf{y}, t|\mathbf{x})$ is the probability density at the position $\mathbf{y}$ at time 0 given that $\mathbf{x}$ is the position at the later time $t$; $q$ is thus the propagator of the tracer particles under time reversal. The evolution of $q$ satisfies

$$ \partial_t q(\mathbf{y}, t|\mathbf{x}) P_0(\mathbf{x}) = - \partial_{x_i} [u_i(\mathbf{x}) P_0(\mathbf{x}) q(\mathbf{y}, t|\mathbf{x})] \\ + \partial_{x_i} \partial_{x_j} [\kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) q(\mathbf{y}, t|\mathbf{x})]. \quad (47) $$

Expanding all the derivatives and using the fact that $HP_0 = 0$, we find that $q$ satisfies

$$ \partial_t q(\mathbf{y}, t|\mathbf{x}) = [u_i^*(\mathbf{x})\partial_{x_i} + \kappa_{ij}(\mathbf{x})\partial_{x_i}\partial_{x_j}]q(\mathbf{y}, t|\mathbf{x}), \quad (48) $$

where

$$ \begin{aligned} u_i^*(\mathbf{x}) &= 2 \frac{\partial_{x_j}[\kappa_{ij}(\mathbf{x}) P_0(\mathbf{x})]}{P_0(\mathbf{x})} - u_i(\mathbf{x}), \\ &= u_i(\mathbf{x}) - 2 \frac{J_{0i}(\mathbf{x})}{P_0(\mathbf{x})}, \end{aligned} \quad (49) $$

is the drift of the time reversed process and where we have introduced the current in the steady state $J_i$

$$ J_i(\mathbf{x}) = u_i(\mathbf{x}) P_0(\mathbf{x}) - \partial_{x_k} [\kappa_{ik}(\mathbf{x}) P_0(\mathbf{x})]. \quad (50) $$
If we interpret (48) as a backward Fokker-Planck equation [3], we deduce that $q$ also satisfies the associated forward Fokker-Planck equation

$$ \partial_t q(\mathbf{y}, t|\mathbf{x}) = -\partial_{y_i}[u_i^*(\mathbf{y})q(\mathbf{y}, t|\mathbf{x})] + \partial_{y_i}\partial_{y_j}[\kappa_{ij}(\mathbf{y})q(\mathbf{y}, t|\mathbf{x})]. \quad (51) $$

Therefore $q$ is the propagator of fictive particles, moving in an effective drift field $u_i^*$ instead of $u_i$. Let us say a little bit more about the properties of $q$. From its definition (46), we immediately see that the initial condition for $q$ is

$$ q(\mathbf{y}, 0|\mathbf{x}) = \delta(\mathbf{x} - \mathbf{y}). \quad (52) $$

Given that no particles can flow in or out of the obstacles, the time reversed process must also have no-flux boundary conditions at the surface of obstacles. That is to say that

$$ \begin{aligned} 0 &= n_i\{u_i^*(\mathbf{y})q(\mathbf{y}, t|\mathbf{x}) - \partial_{y_j}[\kappa_{ij}(\mathbf{y})q(\mathbf{y}, t|\mathbf{x})]\} \\ &= \frac{n_i}{P_0(\mathbf{x})}\{\{u_i^*(\mathbf{y})P_0(\mathbf{y}) - \partial_{y_j}[\kappa_{ij}(\mathbf{y})P_0(\mathbf{y})]\}p(\mathbf{x}, t|\mathbf{y}) \\ &\quad - \kappa_{ij}(\mathbf{y})P_0(\mathbf{y})\partial_{y_j}p(\mathbf{x}, t|\mathbf{y})\}, \end{aligned} \quad (53) $$

and as a consequence we see that $n_i\kappa_{ij}(\mathbf{y})\partial_{y_j}p(\mathbf{x}, t|\mathbf{y}) = 0$, which recovers an established boundary condition for the no-flux boundary condition in terms of the starting coordinate $\mathbf{y}$ [3]. It is interesting to note that the steady state current of the time reversed process is given by $-J_i$, i.e. exactly the opposite of the current of the original process. In addition, we see that for currentless steady states, where $J_i = 0$, we have that $u_i^* = u_i$ and thus the original and time reversed processes are statistically identical in that $p(\mathbf{x}, t|\mathbf{y}) = q(\mathbf{x}, t|\mathbf{y})$. Interestingly, for advection of a particle with constant molecular diffusivity by an incompressible flow, one can easily show that $P_0 = 1/|\Omega|$ where $|\Omega|$ is the volume of the unit cell and that $J_i$ is non-zero and given by $J_i = u_i/|\Omega|$, consequently $u_i^* = -u_i$, i.e. the time reversed process has the opposite flow field to the original process.

Eq. (45) may be rewritten in terms of $q$ as

$$ \hat{B}_j(\mathbf{x}, s) = \int_{\mathbb{R}^d} dy \, \hat{q}(\mathbf{y}, s|\mathbf{x}) P_0(\mathbf{x})(x_j - y_j). \quad (54) $$

Laplace transforming the Fokker-Planck equation (51) for $q$ yields

$$ \begin{aligned} s \hat{q}(\mathbf{y}, s|\mathbf{x}) = \delta(\mathbf{x} - \mathbf{y}) \\ - \partial_{y_i}[u_i^*(\mathbf{y})\hat{q}(\mathbf{y}, s|\mathbf{x})] + \partial_{y_i}\partial_{y_j}[\kappa_{ij}(\mathbf{y})\hat{q}(\mathbf{y}, s|\mathbf{x})]. \end{aligned} \quad (55) $$

Inserting the above equality into Eq. (54) and integrating by parts over $\mathbf{y}$, we obtain

$$ \hat{B}_j(\mathbf{x}, s) = \frac{1}{s} \int_{\mathbb{R}^d} dy \, P_0(\mathbf{x}) \{-u_j^*(\mathbf{y})q(\mathbf{y}, s|\mathbf{x}) \\ + \partial_{y_k}[\kappa_{kj}(\mathbf{y})\hat{q}(\mathbf{y}, s|\mathbf{x})]\}. \quad (56) $$

where we have used the no-flux condition (this time of the time reversed process) at the obstacles boundaries. Integrating by parts the term containing $\kappa_{kj}$, and then switching back to the propagator $p$ instead of $q$ and invoking the periodicity property, we obtain

$$ \begin{aligned} \hat{B}_j(\mathbf{x}, s) = & -\frac{1}{s} \int_{\Omega} dy \, P_0(\mathbf{y}) \hat{P}(\mathbf{x}, s|\mathbf{y}) u_j^*(\mathbf{y}) \\ & + \frac{1}{s} \int_{S_{\Omega}} dS_l(\mathbf{y}) P_0(\mathbf{y}) \hat{P}(\mathbf{x}, s|\mathbf{y}) \kappa_{lj}(\mathbf{y}). \end{aligned} \quad (57) $$

Inserting the above expression into Eqs. (44),(36) gives a formal expression for the tensor $\hat{\psi}_{ij}$. Using the fact that $\hat{P}(\mathbf{x}, s|\mathbf{y}) \approx P_0(\mathbf{x})/s$ for $s \to 0$, one can check that, for large times, $\psi_{ij} \approx \mu_i(t)\mu_j(t)$. It is now useful to introduce the difference between the propagator and its stationary value

$$ P'(\mathbf{x}, t|\mathbf{y}) = P(\mathbf{x}, t|\mathbf{y}) - P_0(\mathbf{x}), \quad (58) $$

and thus has Laplace transform

$$ \hat{P}'(\mathbf{x}, s|\mathbf{y}) = \hat{P}(\mathbf{x}, s|\mathbf{y}) - \frac{1}{s} P_0(\mathbf{x}). \quad (59) $$

In terms of $\hat{P}'$, we obtain from Eqs. (35),(36),(44), and (57) the following final expression for the temporal Laplace transform of the dispersion tensor $\sigma_{ij}(t)$

$$ \begin{aligned} \hat{\sigma}_{ij}(s) = & \frac{2}{s^2} \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \\ & + \frac{1}{s^2} \int_{\Omega} d\mathbf{x} [u_i(\mathbf{x}) K_j(\mathbf{x}, s) + u_j(\mathbf{x}) K_i(\mathbf{x}, s)] \\ & - \frac{1}{s^2} \int_{S_{\Omega}} dS_l(\mathbf{x}) [\kappa_{il}(\mathbf{x}) K_j(\mathbf{x}, s) + \kappa_{jl}(\mathbf{x}) K_i(\mathbf{x}, s)], \end{aligned} \quad (60) $$

with

$$ \begin{aligned} K_j(\mathbf{x}, s) = & -\int_{\Omega} dy \, P_0(\mathbf{y}) \hat{P}'(\mathbf{x}, s|\mathbf{y}) u_j^*(\mathbf{y}) \\ & + \int_{S_{\Omega}} dS_l(\mathbf{y}) P_0(\mathbf{y}) \hat{P}'(\mathbf{x}, s|\mathbf{y}) \kappa_{lj}(\mathbf{y}). \end{aligned} \quad (61) $$

and we recall that $u_i^*$ is the effective drift field after time reversal symmetry defined by Eq. (49). Equations (49),(60) and (61) give a closed form expression for the Laplace transform of the dispersion tensor $\sigma_{ij}(t)$ at all times and is the key result of the paper.

## IV. THE EFFECTIVE DIFFUSION TENSOR

Here we extract the late time effective diffusion tensor $D_{ij}$. The Eqs. (60) and (61) are written in Laplace space and the large time limit is thus extracted from the small $s$ behavior, for which

$$ \hat{\sigma}_{ij}(s) \approx \frac{2 D_{ij}}{s^2} \quad (s \to 0). \quad (62) $$
By definition,

$$ \hat{P}'(\mathbf{x}, s|\mathbf{y}) = \int_{0}^{\infty} dt \, e^{-st} [P(\mathbf{x}, t|\mathbf{y}) - P_{0}(\mathbf{x})]. \quad (63) $$

Taking $s = 0$ in the above expression, we find that $\hat{P}'$ is finite at $s = 0$ (as $P(\mathbf{x}, t|\mathbf{y}) \to P_0(\mathbf{x})$ as $t \to \infty$). Therefore, we can define a function $G$ such that

$$ \hat{P}(\mathbf{x}, 0|\mathbf{y}) = \int_{0}^{\infty} dt \, [P(\mathbf{x}, t|\mathbf{y}) - P_{0}(\mathbf{x})] \equiv G(\mathbf{x}|\mathbf{y}). \quad (64) $$

Applying the operator $H$ to the above equality, it is clear that $G$ obeys

$$ H_{\mathbf{x}}G(\mathbf{x}|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}) - P_0(\mathbf{x}), \quad (65) $$

and also from the conservation of probability $G$ satisfies

$$ \int_{\Omega} d\mathbf{x} \, G(\mathbf{x}|\mathbf{y}) = 0. \quad (66) $$

The function $G$ is therefore the pseudo-Green's function [34] of the operator $H$ (that is, the inverse of $H$ in the subspace orthogonal to the uniform function).

Since $\hat{P}'(\mathbf{x}, s = 0|\mathbf{y})$ is finite, we see that the inverse Laplace transform at late times can be extracted from Eqs. (60) and (61) upon setting $s = 0$ in $\hat{P}'$. This immediately yields

$$ \begin{aligned} D_{ij} = & \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \\ & + \frac{1}{2} \int_{\Omega} d\mathbf{x} \, [u_i(\mathbf{x}) f_j(\mathbf{x}) + u_j(\mathbf{x}) f_i(\mathbf{x})] \\ & - \frac{1}{2} \int_{S_{\Omega}} dS_l(\mathbf{y}) [\kappa_{il}(\mathbf{x}) f_j(\mathbf{x}) + \kappa_{jl}(\mathbf{x}) f_i(\mathbf{x})], \end{aligned} \quad (67) $$

with

$$ \begin{aligned} f_j(\mathbf{x}) = & - \int_{\Omega} dy \, P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) u_j^*(\mathbf{y}) \\ & + \int_{S_{\Omega}} dS_l(\mathbf{y}) P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \kappa_{lj}(\mathbf{y}). \end{aligned} \quad (68) $$

For computational purposes, it is useful to find the partial differential equation satisfied by $f_i$. Acting with $H$ on the above expression and using Eq. (65) leads to

$$ \begin{aligned} H_{\mathbf{x}} f_i(\mathbf{x}) = & - u_i^*(\mathbf{x}) P_0(\mathbf{x}) + P_0(\mathbf{x}) \int_{\Omega} dy \, [u_i^*(\mathbf{y}) P_0(\mathbf{y})] \\ & - P_0(\mathbf{x}) \int_{S_{\Omega}} dS_l(\mathbf{y}) P_0(\mathbf{y}) \kappa_{li}(\mathbf{y}), \end{aligned} \quad (69) $$

which is valid for $\mathbf{x} \in \Omega$ (outside the obstacles). An equivalent equation for $f_i$ is

$$ \begin{aligned} H_{\mathbf{x}} f_i(\mathbf{x}) = & - u_i^*(\mathbf{x}) P_0(\mathbf{x}) \\ & + P_0(\mathbf{x}) \int_{\Omega} dy \, \{ u_i^*(\mathbf{y}) P_0(\mathbf{y}) - \partial_{y_k} [P_0(\mathbf{y}) \kappa_{kj}(\mathbf{y})] \}. \end{aligned} \quad (70) $$

This equation for $f_i(\mathbf{x})$ (for $\mathbf{x}$ outside the obstacles) is supplemented by the following conditions. First, because of the conservation of probability, it follows from Eq. (66) that

$$ \int_{\Omega} d\mathbf{x} \, f_i(\mathbf{x}) = 0. \quad (71) $$

Second, $f_i(\mathbf{x})$ is periodic on $\Omega$. Third, at the surface obstacles, $f_i$ satisfies

$$ n_j\{u_j(\mathbf{x})f_i(\mathbf{x}) - \partial_{x_k}[\kappa_{kj}(\mathbf{x})f_i(\mathbf{x})]\} = -n_j\kappa_{ji}(\mathbf{x})P_0(\mathbf{x}) \quad (72) $$

It is actually not obvious to derive the boundary conditions (72) from Eq. (68). The validity of these boundary conditions is checked explicitly in Appendix B, where we demonstrate that Eq. (68) is the actual solution of the Eqs. (70,71,72), proving that our formulation is correct.

One may alternatively express the diffusion tensor by defining a function $f_i^*$ as

$$ \begin{aligned} f_i^*(\mathbf{y}) = & - \int_{\Omega} d\mathbf{x} \, P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) u_i(\mathbf{x}) \\ & + \int_{S_{\Omega}} dS_l(\mathbf{x}) P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \kappa_{li}(\mathbf{x}), \end{aligned} \quad (73) $$

in terms of which the effective diffusion tensor reads

$$ \begin{aligned} D_{ij} = & \int_{\Omega} dy \, \kappa_{ij}(\mathbf{y}) P_0(\mathbf{y}) \\ & + \frac{1}{2} \int_{\Omega} dy \, [u_i^*(\mathbf{y}) f_j^*(\mathbf{y}) + u_j^*(\mathbf{y}) f_i^*(\mathbf{y})] \\ & - \frac{1}{2} \int_{S_{\Omega}} dS_l(\mathbf{y}) [\kappa_{il}(\mathbf{y}) f_j^*(\mathbf{y}) + \kappa_{jl}(\mathbf{y}) f_i^*(\mathbf{y})]. \end{aligned} \quad (74) $$

Obviously,

$$ \begin{aligned} f_i^*(\mathbf{y}) = & - \int_{\Omega} d\mathbf{x} \, P_0(\mathbf{x}) G^*(\mathbf{y}|\mathbf{x}) u_i(\mathbf{x}) \\ & + \int_{S_{\Omega}} dS_l(\mathbf{x}) P_0(\mathbf{x}) G(\mathbf{y}|\mathbf{x}) \kappa_{li}(\mathbf{x}), \end{aligned} \quad (75) $$

where $G^*$ is the pseudo-Green function associated for the operator $H^*$ of the time reversed process, with the drift field $u^*$, that reads

$$ H_{\mathbf{y}}^*(\cdot) = \partial_{y_i} [u_i^*(\mathbf{y})(\cdot)] - \partial_{y_k} [\kappa_{ij}(\mathbf{y})(\cdot)] \quad (76) $$

Hence, $f_i^*$ satisfies the same equations as $f_i$ if one replaces $u_i$ by $u_i^*$ and vice-versa. As a consequence, $f_i^*$ is the solution of

$$ \begin{aligned} H_{\mathbf{y}}^* f_i^*(\mathbf{y}) = & - u_i(\mathbf{y}) P_0(\mathbf{y}) + P_0(\mathbf{y}) \int_{\Omega} d\mathbf{x} \, [u_i(\mathbf{x}) P_0(\mathbf{x})] \\ & - P_0(\mathbf{y}) \int_{S_{\Omega}} dS_l(\mathbf{x}) P_0(\mathbf{x}) \kappa_{li}(\mathbf{x}) \end{aligned} \quad (77) $$

with periodic boundary conditions and also with the condition

$$ n_j\{u_j^*(\mathbf{y})f_i^*(\mathbf{y}) - \partial_{y_k}[\kappa_{kj}(\mathbf{y})f_i^*(\mathbf{y})]\} = -n_j\kappa_{ji}(\mathbf{y})P_0(\mathbf{y}) \quad (78) $$
at the surface of the obstacles, together with the orthog-
onality condition

$$
\int_{\Omega} dy f_{i}^{*}(y) = 0. \qquad (79)
$$

Thus, we have two formulas for the computation of the
effective diffusion tensor.

Having two formulations for the Kubo formula is a
useful check on numerical solutions as the two resulting
results for the diffusion tensor can be compared. Fur-
thermore, we will now see that both the fields fᵢ and
fᵢ* are needed to compute the leading order finite time
correction to the time dependent diffusion tensor.

V. LATE TIME CORRECTIONS TO THE
EFFECTIVE DIFFUSION TENSOR

In section IV we have extracted, from the full time
dependent Kubo formula, the effective late time limit of
the diffusion tensor $D_{ij}$. While many examples of such
formulas exist for special cases, little is known about how
the time dependent diffusion tensor relaxes to its asymp-
totic limit. To our knowledge the only examples known
are for diffusion in a periodic potential [37] and diffusion
in a system with periodic diffusivity [38]. Both of these
examples have steady states with zero current, here we
will extend these results to the most general cases both
with and without current. Here we extract the late time
correction tensor $C_{ij}$ that describes the approach to the
diffusive limit of the system.

Expanding Eq. (63) in powers of s, we obtain

$$
\hat{P}'(\mathbf{x}, s|\mathbf{y}) \approx G(\mathbf{x}|\mathbf{y}) + s G^{(2)}(\mathbf{x}|\mathbf{y}) + \dots \quad (80)
$$

with

$$
G^{(2)}(\mathbf{x}|\mathbf{y}) = - \int_{0}^{\infty} dt \, t [P(\mathbf{x}, t|\mathbf{y}) - P_{0}(\mathbf{x})]. \quad (81)
$$

Since the motion of the tracer particles is Markovian
(memoryless), we can write the equality

$$
P(\mathbf{x}, t|\mathbf{y}) = \int_{\Omega} dz \, P(\mathbf{x}, t - t_1|\mathbf{z}) P(\mathbf{z}, t_1|\mathbf{y}), \quad (82)
$$

which holds for $0 < t_1 < t$. Integrating this relation over $t_1$ leads to

$$
t P(\mathbf{x}, t|\mathbf{y}) = \int_{0}^{t} dt_{1} \int_{\Omega} dz P(\mathbf{x}, t - t_{1}|\mathbf{z}) P(\mathbf{z}, t_{1}|\mathbf{y}). \quad (83)
$$

Substracting $tP_0(\mathbf{x})$ and taking the Laplace transform for small values of $s$, we obtain

$$
G^{(2)}(\mathbf{x}|\mathbf{y}) = - \int_{\Omega} dz \, G(\mathbf{x}|\mathbf{z})G(\mathbf{z}|\mathbf{y}). \quad (84)
$$

The correction tensor $C_{ij}$ is identified from the relation
for small $s$

$$
\hat{\sigma}_{ij}(s) \approx \frac{2}{s^2} (D_{ij} + s C_{ij} + \dots) \quad (s \to 0). \quad (85)
$$

(80) and (84) leads to

$$
C_{ij} = -\frac{1}{2} \int_{\Omega} dz \left\{ \int_{\Omega} dx [u_i(\mathbf{x}) f_j(\mathbf{z}) + u_j(\mathbf{x}) f_i(\mathbf{z})] G(\mathbf{x}|\mathbf{z}) \right. \\
\left. - \int_{S_{\Omega}} dS_l(\mathbf{x}) [\kappa_{il}(\mathbf{x}) f_j(\mathbf{z}) + \kappa_{jl}(\mathbf{x}) f_i(\mathbf{z})] G(\mathbf{x}|\mathbf{z}) \right\}, \quad (86)
$$

where we have used the definition of $f$ in Eq. (68) to per-
form the integration over $\mathbf{y}$. Using the definition of $f_i^*$ in
Eq. (75) and the relation $G^*(\mathbf{y}|\mathbf{x})P_0(\mathbf{x}) = G(\mathbf{x}|\mathbf{y})P_0(\mathbf{y})$,
we finally obtain

$$
C_{ij} = \frac{1}{2} \int_{\Omega} dz \frac{f_i(\mathbf{z}) f_j^*(\mathbf{z}) + f_j(\mathbf{z}) f_i^*(\mathbf{z})}{P_0(\mathbf{z})}. \quad (87)
$$

This relation is a compact and explicit form for the cor-
rection tensor $C_{ij}$ and is the main result of this section.
One can also derive the above result by decomposing $P'$
in terms of left and right eigenvector of the Fokker-Planck
operator $H$. In doing this is straightforward to see
that the temporal corrections at the next order decay as
$\exp(-\lambda_1 t)$, where $\lambda_1$ is the lowest positive eigenvalue of
$H$ which will be strictly positive given that the domain
$\Omega$ is taken to be finite.

VI. SPECIAL CASES AND COMPARISON
WITH EXISTING RESULTS

Here we examine the form the Kubo formula takes for a
number of systems and compare them to existing results
in the literature.

A. Flow in frozen incompressible velocity fields
with isotropic constant diffusion tensor

The Fokker-Planck transport operator for a tracer ad-
vected by an incompressible velocity field **v** with constant
isotropic diffusivity $\kappa_{ij}(\mathbf{x}) = \kappa_0 \delta_{ij}$ is defined via

$$
H_{\mathbf{x}} = -\kappa_0 \partial_{x_i} \partial_{x_i} + v_i(\mathbf{x}) \partial_{x_i}. \quad (88)
$$

Here the drift field is thus $u_i = v_i$. From the incomp-
ressibility condition $\partial_{x_i} v_i = 0$, we see that the steady
state distribution on the unit cell $\Omega$ is uniform and thus
given by $P_0(\mathbf{x}) = 1/|\Omega|$, where $|\Omega|$ denotes the available
volume of the unit cell. The steady state current is thus
given by $J_i(\mathbf{x}) = v_i(\mathbf{x})/|\Omega|$. We also find that $u_i^* = -u_i$.
Applying Eq. (67) leads to

$$
D_{ij} = \delta_{ij}\kappa_0 + \frac{1}{2} \int_{\Omega} dz [f_i(\mathbf{x})v_j(\mathbf{x}) + f_j(\mathbf{x})v_i(\mathbf{x})] \\
- \frac{\kappa_0}{2} \left\{ \int_{S_{\Omega}} dS_i(\mathbf{x})f_j(\mathbf{x}) + \int_{S_{\Omega}} dS_j(\mathbf{x})f_i(\mathbf{x}) \right\}. \quad (89)
$$
The equation satisfied by $f_i$ is a simplified form of Eq. (70):

$$ -\kappa_0 \partial_{x_j}^2 f_i + v_j \partial_{x_j} f_i = \frac{1}{|\Omega|} v_i(\mathbf{x}) - \frac{1}{|\Omega|^2} \int_{\Omega} dy \, v_i(\mathbf{y}). \quad (90) $$

The boundary conditions on the surface of obstacles are [see Eq. (72)]

$$ n_j \partial_{x_j} f_i = n_i / |\Omega|, \quad (91) $$

(where we have assumed that no fluid enters the obstacles, $v_i n_i = 0$). It is then obvious that the last four expressions recover those in the hydrodynamics-homogenization literature (see e.g. Refs. [14, 35]), if we identify the function $f_i$ of this paper with $-f_i/|\Omega|$ to match with the notations of Ref. [35].

## B. Systems with no steady state current

A large class of models studied in statistical mechanics, such as diffusion in a periodic potential or diffusion in a medium having periodic diffusivity, have no current in the steady state. In these cases the original and its time-reversed processes have the same Fokker-Planck evolution operator, and consequently the same pseudo-Green's function. However in general we have that

$$ G(\mathbf{x}|\mathbf{y})P_0(\mathbf{y}) = P_0(\mathbf{x})G^*(\mathbf{y}|\mathbf{x}), \quad (92) $$

where $G^*(\mathbf{y}|\mathbf{x})$ denotes the pseudo-Green's function for $H^*$. In the case where there is no current we know however that $G^*(\mathbf{y}|\mathbf{x}) = G(\mathbf{x}|\mathbf{y})$ and thus we find that the operator

$$ K(\mathbf{x}, \mathbf{y}) = G(\mathbf{x}|\mathbf{y})P_0(\mathbf{y}) = P_0(\mathbf{x})G(\mathbf{y}|\mathbf{x}) \quad (93) $$

is symmetric. In the currentless case we also have

$$ D_{ij} = \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) P_0(\mathbf{x}) \\ - \frac{1}{2} \iint_{\Omega} d\mathbf{x} d\mathbf{y} [\tilde{u}_i(\mathbf{x})\tilde{u}_j(\mathbf{y}) + \tilde{u}_j(\mathbf{x})\tilde{u}_i(\mathbf{y})] K(\mathbf{x}, \mathbf{y}), \quad (94) $$

from Eq. (10). The operator $K$, being symmetric, must have real eigenvalues, however $K$ has the same eigenvalues as $G$ which must have eigenvalues with a positive real part if a steady state regime can be attained. Consequently under these assumptions $K$ must be a positive operator and we must have that the change in the diagonal terms of the diffusion tensor due to the presence of drift or inclusions is negative, i.e.

$$ D_{ii} - \langle \kappa_{ii} \rangle_0 = - \iint_{\Omega} d\mathbf{x} d\mathbf{y} \tilde{u}_i(\mathbf{x}) \tilde{u}_i(\mathbf{y}) K(\mathbf{x}, \mathbf{y}) < 0, \quad (95) $$

that is to say that diffusion in equilibrium systems is slowed down by the presence of the drift. This is physically obvious for say diffusion in periodic potentials where

the diffusing particle becomes trapped in local minima of the potential. This slowing down due to drift also occurs for diffusion in a medium of spatially varying diffusivity where

$$ H_{\mathbf{x}} = -\partial_{x_i} \kappa_{ij}(\mathbf{x}) \partial_{x_j}, \quad (96) $$

here $P_0 = 1/|\Omega|$ and $J_i = 0$. The drift field in this case is, in our notation, given by

$$ u_i(\mathbf{x}) = \partial_{x_j} \kappa_{ij}(\mathbf{x}), \quad (97) $$

and in terms of the pseudo-Greens' function $G$ we find

$$ D_{ij} = \frac{1}{|\Omega|} \int_{\Omega} d\mathbf{x} \, \kappa_{ij}(\mathbf{x}) \\ - \frac{1}{|\Omega|} \iint_{\Omega} d\mathbf{x} d\mathbf{y} [\partial_{x_k} \kappa_{ik}(\mathbf{x})][\partial_{x_l} \kappa_{il}(\mathbf{y})] G(\mathbf{x}|\mathbf{y}), \quad (98) $$

which recovers the formula derived in Ref. [21]. We see that the first term in the above is the spatial average of the local diffusivity tensor while the second terms gives a negative contribution when $i = j$. Finally we remark that the case of diffusion with constant diffusivity $\kappa_0$ in a periodic potential $\Phi$ gives $u_i = -\kappa_0 \beta \partial_{x_i} \Phi(\mathbf{x})$, where $\beta$ is the inverse temperature. Here the steady state distribution $P_0$ is given by the Gibbs-Boltzmann distribution on the unit cell $\Omega$

$$ P_0(\mathbf{x}) = \frac{\exp[-\beta\Phi(\mathbf{x})]}{Z_{\Omega}}, \quad (99) $$

where $Z_{\Omega}$ is the partition function normalizing the distribution, thus recovering the result of Ref. [21] for this particular case.

## VII. GENERAL RESULTS IN ONE DIMENSION

We consider a generic one dimensional system with periodicity denoted by $L$. Here we can give explicit formula for all terms in the static Kubo formula, although the computation is surprisingly involved. In what follows we will use notation similar to Refs. [26, 28] to facilitate comparison with their results. In one dimension the steady state current is constant in space and the steady state probability distribution is given by

$$ P_0(x) = JI_+(x), \quad (100) $$

where the term $I_+$ reads

$$ I_+(x) = \frac{\exp(\Gamma(x))}{\kappa(x)} \int_x^\infty dx' \exp(-\Gamma(x')), \quad (101) $$

with

$$ \Gamma(x) = \int_0^x dx' \frac{u(x)}{\kappa(x)}. \quad (102) $$
Due to the periodicity of $u$ and $\kappa$ the function $\Gamma$ obeys the relation

$$ \Gamma(x + L) = \Gamma(x) + \Gamma(L). \quad (103) $$

When $\Gamma(L) = 0$ the system clearly has a steady state equilibrium distribution with no current. In writing Eq. (101) we have assumed, without loss of generality, that $\Gamma(L) > 0$ so that the integral on the right hand side converges. The steady state current is then obtained from the condition of normalization of $P_0$ and is thus given by

$$ J = \frac{1}{\int_{0}^{L} dx \ I_{+}(x)}, \quad (104) $$

and consequently the effective drift is given by $V = JL$.

To compute the effective diffusion constant we use the representation for $D$ given in Eq. (74) in its one dimensional version, that is to say

$$ D = \int_{0}^{L} dx \ \kappa(x) P_{0}(x) + \int_{0}^{L} dx f^{*}(x) u^{*}(x). $$

We now define $f^*(x) = g(x)P_0(x)$ and write the diffusion constant as

$$ D = \int_{0}^{L} dx \, \kappa(x) P_{0}(x) + \int_{0}^{L} dx g(x) \left\{ J - \frac{d}{dx}[\kappa(x) P_{0}(x)] \right\}. $$

Defining $g(x) = -x + h(x)$ we find that $h$ obeys

$$ \kappa \frac{d^2 h}{dx^2} + u \frac{dh}{dx} = JL. \quad (105) $$

A first integration of this equation for $h$ yields

$$ \frac{dh(x)}{dx} = JLI_{-}(x), \quad (106) $$

where

$$ I_{-}(x) = \exp(-\Gamma(x)) \int_{-\infty}^{x} dx' \frac{\exp(\Gamma(x'))}{\kappa(x')}, \quad (107) $$

and where we have used the fact that $dg/dx = -1+dh/dx$ is periodic and thus $dh/dx$ must be periodic. Integrating again yields

$$ h(x) = J_0 L \int_0^x dx' I_-(x') + C \quad (108) $$

where $C$ is an integration constant that will be determined from the orthogonality condition $\int_0^L dx \, P_0(x)g(x) = 0$. It is not immediately obvious from Eq. (108) that $g(x) = -x + h$ is periodic, however it is easy to see that

$$ g(x + L) = g(x) - L + JL \int_{0}^{L} dx \, I_{-}(x). \quad (109) $$

Now by integration by parts one finds

$$ \int_{0}^{L} dx \, I_{-}(x) = \int_{0}^{L} dx \, I_{+}(x), \quad (110) $$

and using this and Eq. (104) in Eq. (109) yields $g(x) = g(x + L)$. Incidentally this shows that we may also write

$$ J = \frac{1}{\int_{0}^{L} dx \, I_{-}(x)}. \quad (111) $$

To proceed it is useful to define the two new functions

$$ A_{+}(x) = \int_{-\infty}^{x} dx' \frac{\exp(\Gamma(x'))}{\kappa(x')}, \quad (112) $$

$$ A_{-}(x) = \int_{x}^{\infty} dx' \exp(-\Gamma(x')). \quad (113) $$

Further more one can show that

$$ A_{+}(x + L) = \exp(\Gamma(L)) A_{+}(x), \quad (114) $$

$$ A_{-}(x + L) = \exp(-\Gamma(L)) A_{-}(x). \quad (115) $$

In terms of these two functions we can write $I_+$ and $I_-$ as

$$ I_{+}(x) = \frac{dA_{+}(x)}{dx} A_{-}(x), \quad (116) $$

$$ I_{-}(x) = - \frac{dA_{-}(x)}{dx} A_{+}(x). \quad (117) $$

Computing the constant $C$ in Eq. (108) using the orthogonality relation then yields

$$ D = \int_{0}^{L} dx \left[ \kappa(x) P_{0}(x) \frac{dh}{dx} + J [-x + h(x) + LxP_{0}(x) - Lh(x)P_{0}(x)] \right], \quad (118) $$

where we have carried out an integration by parts in the last term of Eq. (105). The function $h$ is explicitly given by

$$ h(x) = JL \int_{0}^{x} dx' I_{-}(x'), \quad (119) $$

and using the expression Eq. (117) for $I_-$ it can be written as

$$ h(x) = -JL \left[ A_{+}(x)A_{-}(x) - A_{+}(0)A_{-}(0) - \int_{0}^{x} dx' \frac{P_{0}(x')}{J} \right]. \quad (120) $$

We also, by integration by parts, note the identity

$$ \int_{0}^{L} dx \int_{0}^{x} dx' P_{0}(x') = L - \int_{0}^{L} dx \, xP_{0}(x). \quad (121) $$

Using these relations in Eq. (118), after some algebra, we obtain the compact expression

$$ D = \frac{L^2 \int_{0}^{L} dx \, \kappa(x) I_{\pm}(x)^2 I_{\mp}(x)}{\int_{0}^{L} dx \, I_{\pm}(x)^3}, \quad (122) $$
where ± indicates that the index can be taken to be + or −. One can verify that Eq. (122) agrees with the results of Reimann et al. [26, 28] for the case of diffusion in a tilted potential $\phi(x) = V(x) - fx$, where V is periodic. It also agrees with the formulas presented in Ref. [31] where the local mobility is also varying. The result given here shows how the diffusion constant of any one dimensional system with periodic diffusivity and drift can be obtained.

## VIII. CONCLUSIONS

The problem of computing effective transport coefficients occurs in many settings of fundamental physics and has many applications. Many results exist in the literature, for example in fluid mechanics using homogenization theory [10–20] as well as in statistical physics [21–25] where both exact and approximate results exist for problems such as diffusion in periodic or random potentials and diffusion in media with locally varying diffusivity. In these latter studies most results have been obtained for systems with equilibrium (currentless) steady states, given for example by the Gibbs Boltzmann distribution for diffusion in a periodic potential.

More recently results have been derived for diffusion in systems with non zero currents [26–32] In such non equilibrium systems interesting new physics such as massively increased diffusivity due to an external applied force has been discovered. This enhancement occurs, in the weak noise limit, at the external field where all minima disappear from the overall potential to be replaced by points of inflection [26]. The difference between different trajectories at this critical point is enhanced and dispersion in thus increased. Using the Kubo formula, given here, one could investigate the effect of an externally applied field on diffusion in higher dimensional periodic potentials. The formula given here cannot in general be evaluated analytically in higher than one dimension, however their numerical evaluation is straightforward and the transport coefficients can be accurately determined using standard packages to solve partial differential equations. Of course many of the systems considered here can be studied via numerical simulations where one integrates the corresponding SDE, however the simulation approach suffers from the need to estimate the errors due to statistical fluctuations and in the case of obstacles the correct imposition of the no-flux boundary condition requires careful treatment and is far from being obvious [40–42].

Many other problems can be tackled using the present approach, in Ref. [33] it was shown how the presence of an externally applied uniform field on a system with spatially varying diffusivity modifies the dispersion of a cloud of tracer particles. It was shown that the diffusion in the direction of the applied force could be hugely increased for systems in two or more dimensions. In addition it was shown that the components of the effective

diffusion tensor can exhibit counter intuitive monotonic behavior. The results thus imply that one can control the dispersion of a diffusing cloud with an external field, and such an effect may well have useful applications.

The formula given here are extremely general and should be valuable for the interpretation of experiments. For instance, potential fields in which colloidal particles diffuse are often generated by laser [43–46], in which case absorption of laser light can also lead to temperature gradients. The variation of the local temperature will clearly influence local transport properties through the temperature dependence of the surrounding fluid viscosity as well as due to thermodynamic forces associated with temperature gradients, i.e the Soret effect [5].

Finally there has been much recent study of fluctuation dissipation relations (FDR) in non-equilibrium systems [47–50]. A notable example of this type of FDR is the Stokes-Einstein relation between the differential response of the mean velocity $V_i$ and the diffusion tensor. In our preliminary study [33], we showed how the presence of a steady state current leads to corrections to the usual Stokes-Einstein formula. However the general analysis of the correction term needs further study in order to fully understand the non-equilibrium physics it encodes, as well as to make contact with the literature on the subject [47–50]. In the literature in question results are most often given in terms of an integral over a time dependent violation factor, the results given here could be useful to understand these results in terms of the static, time independent quantities, used here.

## Appendix A: Details on the numerical simulations

Here we briefly describe the algorithm used to simulate the motion of a Brownian walker in a 2D array of circular disks, in the presence of a force F, leading to the results presented in Fig.1(b). Let $\mathbf{X}$ be the position at time t of the walker, then we define $\mathbf{e}_\alpha$ the unit vector with direction $\mathbf{X} - \mathbf{X}_c$, with $\mathbf{X}_c$ the nearest disk center. $\mathbf{e}_\alpha$ makes an angle $\alpha$ with the x-direction. During the time step $\Delta t$, the motion in the direction $\mathbf{e}_\alpha$ is modified by the increment

$$dX_{\perp} = \sqrt{\Delta t} f_1(r_{\perp}/\sqrt{\Delta t}) + u_{\parallel} f_2(r_{\perp}/\sqrt{\Delta t}) + h \cos \alpha \Delta t \quad (A1)$$

Here $h = \beta FD$, $r_\perp$ is the distance to the nearest disk surface, and $f_1$ and $f_2$ are the functions calculated by Peters et al. in the case of reflecting boundaries [40], and $u_\parallel = \pm\sqrt{D\Delta t}$ with probabilities 1/2. In the direction parallel to the nearest surface obstacles, the increment is $dX_\parallel = -h \sin \alpha \Delta t + u$, where u is a stochastic variable with zero mean and variance $\langle u^2 \rangle = 2D\Delta t$. With this algorithm, when the particle is close to the boundaries, $r_\perp \sim \sqrt{\Delta t}$, the terms in $f_1$ and $f_2$ dominate the increment in the perpendicular direction in Eq. (A1) if $\Delta t$ is small enough. These terms account for the fact that,
since trajectories that cross the surface of the obstacles are forbidden, the averaging over authorized trajectories between $t$ and $t + \Delta t$ gives rise to an effective drift oriented towards the exterior of the obstacles, described by the term $f_1$, and also modifies the variance of these trajectories (term $f_2$), see Ref [40] for details. If the particle is far from the obstacles, $r_⊥ \gg \sqrt{\Delta t}$, then $f_1 \to 0$ and $f_2 \to \sqrt{2}$ and one recovers Brownian motion under force in the absence of obstacles. Although the algorithm of Ref [40] was proposed only for one-dimensional geometries, we expect that it is also valid in the present 2D case, since for sufficiently small $\Delta t$ the curvature of the surface should play no role. For the results of Fig. 1(b), the time step was $\Delta t = 10^{-5}L^2/D$ and was checked to be sufficiently small so that convergence was reached for the calculation of $D_{ij}$ which were estimated over 500,000 trajectories. The results of simulations match with the predictions of our (exact) Kubo formulas.

## Appendix B: Boundary condition (72) for the function $f_i$

In this appendix, we show that the formulation of the problem under the form of Eqs. (70,71,72) is correct by checking that Eq. (68) is the actual solution of these equations. We consider the quantity

$$M_i(\mathbf{x}) = \int_{\Omega} dy \{[H_y^\dagger G(\mathbf{x}|y)]f_i(\mathbf{y}) - G(\mathbf{x}|y)[H_y f_i(\mathbf{y})]\} \quad (B1)$$

where the adjoint operator reads

$$H_y^\dagger = -u_i(\mathbf{y})\partial_{y_i} - \kappa_{ij}(\mathbf{y})\partial_{y_i}\partial_{y_j}. \quad (B2)$$

The pseudo-Green function $G$ can be shown to satisfy the adjoint equation

$$H_y^\dagger G(\mathbf{x}|\mathbf{y}) = \delta(\mathbf{x} - \mathbf{y}) - P_0(\mathbf{x}), \quad (B3)$$

The following condition also holds

$$\int_{\Omega} dy P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) = 0. \quad (B4)$$

Using the above relations and Eqs. (69) and (71) we deduce that

$$M_i(\mathbf{x}) = f_i(\mathbf{x}) + \int_{\Omega} dy P_0(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) u_i^*(\mathbf{y}). \quad (B5)$$

On the other hand, using the explicit expressions of $H$ [Eq.(1)] and integrating by parts two times the term containing $H$ in Eq. (B1), we obtain

$$M_i(\mathbf{x}) = \int_{S_{\Omega}} dS_j(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \{-u_j(\mathbf{y})f_i(\mathbf{y}) + \partial_{y_k}[\kappa_{jk}(\mathbf{y})f_i(\mathbf{y})]\} - \int_{S_{\Omega}} dS_j(\mathbf{y}) f(\mathbf{y})\kappa_{jk}(\mathbf{y})\partial_{y_k}G(\mathbf{x}|\mathbf{y}). \quad (B6)$$

The last term of this equation vanishes, since the relation $dS_j\kappa_{kj}(\mathbf{y})P(\mathbf{x}, t|\mathbf{y}) = 0$ is known for reflecting boundaries [3]. Then, using the boundary conditions (72), we obtain

$$M_i(\mathbf{x}) = \int_{S_{\Omega}} dS_j(\mathbf{y}) G(\mathbf{x}|\mathbf{y}) \kappa_{jk}(\mathbf{x}), \quad (B7)$$

which, when compared to Eq. (B5), is exactly Eq. (68), thereby proving that the formulation of Eqs. (70,72,71) is correct.

[1] N. Van Kampen, *Stochastic Processes in Physics and Chemistry*, Third Edition (North-Holland, Amsterdam, 2007).

[2] B. Øksendal, *Stochastic differential equations* (Springer, New-York, 2003).

[3] C. Gardiner, *Handbook of stochastic methods for physics, chemistry and the natural sciences, second edition* (1985).

[4] H. Brenner, *Chem. Eng. Sci.* **16**, 242 (1961).

[5] A. Würger, *Rep. Progr. Phys.* **73**, 126601 (2010).

[6] F. Daumas, N. Destainville, C. Millot, A. Lopez, D. Dean, and L. Salomé, *Biophys. J.* **84**, 356 (2003).

[7] E. A. Reits and J. J. Neefjes, *Nat. Cell Biol.* **3**, E145 (2001).

[8] S. Condamin, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, *Nature* **450**, 77 (2007).

[9] G. Taylor, *Proc. R. Soc. Lon. A* **219**, 186 (1953).

[10] M. Rosenbluth, H. Berk, I. Doxas, and W. Horton, *Phys. Fluids* **30**, 2636 (1987).

[11] B. I. Shraiman, *Phys. Rev. A* **36**, 261 (1987).

[12] P. McCarty and W. Horsthemke, *Phys. Rev. A* **37**, 2112 (1988).

[13] A. J. Majda and P. R. Kramer, *Phys. Rep.* **314**, 237 (1999).

[14] H. Brenner, *Philos. Tr. R. Soc. A* **297**, 81 (1980).

[15] J. Rubinstein and R. Mauri, *SIAM J. Appl. Math.* **46**, 1018 (1986).

[16] M. Quintard and S. Whitaker, *Adv. Water. Ress.* **17**, 221 (1994).

[17] A. Alshare, P. Strykowski, and T. Simon, *Int. J. Heat Mass Transfer* **53**, 2294 (2010).

[18] H. P. A. Souto and C. Moyne, *Phys. Fluids* **9**, 2253 (1997).
[19] H. P. A. Souto and C. Moyne, Phys. Fluids 9, 2243 (1997).

[20] M. Quintard and S. Whitaker, Chem. Eng. Sci. 48, 2537 (1993).

[21] D. S. Dean, I. Drummond, and R. Horgan, J. Stat. Mech. 2007, P07013 (2007).

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

[23] R. Zwanzig, Proc. Natl. Acad. Sci. U S A 85, 2029 (1988).

[24] P. De Gennes, J. Stat. Phys. 12, 463 (1975).

[25] S. Lifson and J. L. Jackson, J. Chem. Phys. 36, 2410 (1962).

[26] P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J. Rubi, and A. Pérez-Madrid, Phys. Rev. Lett. 87, 010602 (2001).

[27] B. Lindner, M. Kostur, and L. Schimansky-Geier, Fluct. Noise Lett. 01, R25 (2001).

[28] P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J. Rubi, and A. Pérez-Madrid, Phys. Rev. E 65, 031104 (2002).

[29] P. Reimann and R. Eichhorn, Phys. Rev. Lett. 101, 180601 (2008).

[30] B. Lindner and L. Schimansky-Geier, Phys. Rev. Lett. 89, 230602 (2002).

[31] D. Reguera, G. Schmid, P. S. Burada, J. Rubi, P. Reimann, and P. Hänggi, Phys. Rev. Lett. 96, 130603 (2006).

[32] P. Burada, G. Schmid, P. Talkner, P. Hänggi, D. Reguera, and J. Rubi, BioSystems 93, 16 (2008).

[33] T. Guérin, T. and D. S. Dean, Phys. Rev. Lett. 115, 020601 (2015).

[34] G. Barton, Elements of Green's functions and propagation (Clarendon Press, Oxford, 1989).

[35] R. Carbonell and S. Whitaker, Chem. Eng. Sci. 38, 1795 (1983).

[36] I. Drummond and R. Horgan, J. Phys. A- Math. Gen. 20, 4661 (1987).

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

[38] D. S. Dean and T. Guérin, Phys. Rev. E 90, 062114 (2014).

[39] R. Stratonovich, Radiotekh Elektron. (Moscow) 3, 497 (1958).

[40] E. Peters and T. Barenbrug, Phys. Rev. E 66, 056701 (2002).

[41] T. Barenbrug, E. Peters, and J. Schieber, J. Chem. Phys. 117, 9202 (2002).

[42] G. Lamm and K. Schulten, J. Chem. Phys. 78, 2713 (1983).

[43] C. Dalle-Ferrier, M. Krüger, R. D. Hanes, S. Walta, M. C. Jenkins, and S. U. Egelhaaf, Soft Matt. 7, 2064 (2011).

[44] M. Evstigneev, O. Zvyagolskaya, S. Bleil, R. Eichhorn, C. Bechinger, and P. Reimann, Phys. Rev. E 77, 041107 (2008).

[45] F. Evers, R. Hanes, C. Zunke, R. Capellmann, J. Bewerunge, C. Dalle-Ferrier, M. Jenkins, I. Ladadwa, A. Heuer, R. Castañeda-Priego, et al., Eur. Phys. J. Special Topics 222, 2995 (2013).

[46] R. D. Hanes, C. Dalle-Ferrier, M. Schmiedeberg, M. C. Jenkins, and S. U. Egelhaaf, Soft Matter 8, 2714 (2012).

[47] T. Speck and U. Seifert, Europhys. Lett. 74, 391 (2006).

[48] V. Blickle, T. Speck, C. Lutz, U. Seifert, and C. Bechinger, Phys. Rev. Lett. 98, 210601 (2007).

[49] M. Baiesi, C. Maes, and B. Wynants, Proc. Roy. Soc. London A: Math. Phys. 467, 2792 (2011).

[50] C. Maes, S. Safaverdi, P. Visco, and F. Van Wijland, Phys. Rev. E 87, 022125 (2013).