1. Introduction
Over long time scales, the dispersion of particles and passive scalars in periodic incompressible flows is asymptotically a pure diffusive process, with an effective diffusivity tensor which can be computed, e.g. using the method of homogenisation (Majda & Kramer Reference Majda and Kramer1999; Pavliotis & Stuart Reference Pavliotis and Stuart2005). The flows that have attracted most attention in this context are shear flows and cellular flows because of the remarkably different dependence of their effective diffusivities on the Péclet number and the availability of closed-form results. The cellular flow – on which this paper focuses – is the two-dimensional flow with the streamfunction
\begin{equation} \psi(x,y)= U a \sin (x/a)\sin (y/a), \end{equation}
where 
$U$ is the maximum flow speed and with cell period 
$2{\rm \pi} a$. This flow consists of a doubly periodic array of cells containing four vortices, with the fluid rotating alternatively clockwise and anticlockwise in each quarter-cell. A classic result, due to Childress (Reference Childress1979), Rosenbluth et al. (Reference Rosenbluth, Berk, Doxas and Horton1987), Shraiman (Reference Shraiman1987) and Soward (Reference Soward1987), gives the asymptotic form of the (isotropic) effective diffusivity for non-inertial particles when advection dominates over diffusion. Using 
$a$ and 
$U/a$ as reference length and time, it reads
\begin{equation} \bar{D} \sim 2 \mathscr{S} Pe^{-1/2} \quad \text{for}\ Pe\gg 1, \end{equation}
where 
$Pe=Ua/D$ is the Péclet number and 
$D$ is the diffusion coefficient. The constant 
$\mathscr {S} \approx 0.5327\cdots$ was determined by Soward (Reference Soward1987) using Wiener–Hopf techniques. (See also the mathematical literature, e.g. Heinze (Reference Heinze2003) and Novikov, Papanicolaou & Ryzhik (Reference Novikov, Papanicolaou and Ryzhik2005), for rigorous bounds.)
In this paper, we investigate how the dispersion is affected by small-but-finite particle inertia. Numerical simulations (e.g. Pavliotis, Stuart & Band Reference Pavliotis, Stuart, Band, Niederreiter and Talay2006; Pavliotis, Stuart & Zygalakis Reference Pavliotis, Stuart and Zygalakis2009) show a strong enhancement of the effective diffusivity compared with the non-inertial case when the particles are denser than the fluid. This is because inertia expels heavy particles away from high vorticity regions, towards the cellular flow's separatrices, enhancing transport between cells and hence large-scale dispersion. Conversely, particles less dense than the fluid tend to accumulate in the (high vorticity) cell centres, leading to a much smaller effective diffusivity. Our aim is to quantify this by generalising (1.2) to inertial particles.
The dynamics of non-inertial particles is Brownian, governed by the non-dimensional equation
\begin{equation} \textrm{d}\boldsymbol{X}=\boldsymbol{u}\,\textrm{d} t+(2/Pe)^{1/2}\,\textrm{d}\boldsymbol{W}, \end{equation}
where 
$\boldsymbol {X}$ is the particle position, 
$\boldsymbol {W}$ a two-dimensional Wiener process and 
$\boldsymbol {u}=(\partial _y\psi ,-\partial _x\psi )$ is the fluid velocity. We model the dynamics of inertial particles by the Langevin equation based on the Maxey–Riley (Reference Maxey and Riley1983) model,
\begin{gather} \textrm{d}\boldsymbol{X}=\boldsymbol{V}\,\textrm{d} t, \end{gather}
\begin{gather} St\, \textrm{d}\boldsymbol{V}=(\boldsymbol{u}(\boldsymbol{X},t)-\boldsymbol{V})\,\textrm{d} t+St\beta \, \textit{D}_t\boldsymbol{u}\,\textrm{d} t+(2/Pe)^{1/2}\,\textrm{d} \boldsymbol{W}, \end{gather}
where 
$\boldsymbol {V}$ is the particle velocity, 
${D}_t=\partial _t+\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material derivative along the fluid velocity 
$\boldsymbol {u}$, 
$St=\tau U/a$ is the Stokes number, with 
$\tau =r_{p}^2/(9\beta \nu )$ the Stokes time scale, 
$\beta =3\rho _{f}/(2\rho _{p}+\rho _{f})\in [0,3]$ is the density parameter with 
$\rho _{f}$ and 
$\rho _{p}$ the mass density of the fluid and particle, 
$r_{p}$ is the particle radius and 
$\nu$ is the kinematic viscosity of the fluid. The momentum equation (1.4b) includes Stokes drag (first term on the right-hand side) and Auton's added-mass force (second term, see Auton, Hunt & Prud'Homme (Reference Auton, Hunt and Prud'Homme1988)) but neglects the Boussinesq–Basset force and Faxén correction as well as gravity.
Typical examples of single trajectories, obtained by solving (1.4) for light (
$\beta > 1$) and heavy (
$\beta < 1$) particles and (1.3) for non-inertial particles, are shown in figure 1(a–c). The early-time dispersion is illustrated by figure 1(d–f) which shows the location at a fixed time 
$t$ of an ensemble of particles initially concentrated within a single quarter-cell. The figure demonstrates the expected enhancement of dispersion for heavy particles and inhibition for light particles.
Single particle trajectory 
$\boldsymbol {X}(t)$ for 
$0\leq t\leq 500$ (a–c) and position at 
$t=30$ of 
$10^5$ particles initially distributed uniformly in a quarter-cell (d–f), with thin lines indicating the streamlines 
$\psi =0,{\pm }0.5$. The non-inertial results (b,e) are obtained by integrating (1.3) with 
$Pe=1000$; the results for inertial particles are obtained by integrating (1.4) with 
$Pe=1000$, 
$St=0.1$ and 
$\beta =2$ (a,d) or 
$\beta =0$ (c,f).
The Brownian dynamics (1.3), in which the velocity rather than the acceleration is white in time, is best thought of as a limit the Langevin dynamics (1.4), obtained when inertial effects are negligible compared with diffusion, that is, for 
$St \ll Pe^{-1}\ll 1$. We consider the more general, distinguished regime 
$St \ll 1$, 
$Pe \gg 1$ with 
$St \, Pe =O(1)$, where both inertial and diffusive effects enter the problem at the same order. We emphasise that analysing this distinguished regime makes it possible to treat the entire range of relative values of 
$St$ and 
$Pe$ at once, including the limiting cases 
$St \ll Pe^{-1}\ll 1$ and 
$Pe^{-1} \ll St \ll 1$. Our main result is an asymptotic formula for the effective diffusivity in this regime:
\begin{equation} \bar{D} \sim \frac{2\mathscr{S} Pe^{-1/2}}{Z(\alpha)} \quad \text{for} \ St, \, Pe^{-1} \ll 1,\ St \, Pe =O(1). \end{equation}Here,
\begin{equation} \alpha = St\,Pe \, (1-\beta) \end{equation}
and 
$Z(\alpha )$ is a function given explicitly in terms of elliptic integrals in (4.12) below and shown in figure 2. Comparing (1.5) with its non-inertial counterpart (1.2) shows that 
$Z(\alpha )$ completely captures the effect of inertia. Since it is monotonically decreasing and satisfies 
$Z(0)=1$, inertia is confirmed to decrease 
$\bar {D}$ for light particles (
$\alpha <0$) and to increase 
$\bar {D}$ for heavy particles (
$\alpha >0$).
Log–log plot of 
$Z(\alpha )$ defined by (4.12). Asymptotic expressions for 
$\alpha \to \pm \infty$ obtained in appendix C are shown as dashed and dot-dashed lines (
$\gamma _{E} = 0.577\cdots$ is the Euler–Mascheroni constant and 
$\varUpsilon =0.655\cdots$ is defined in (C 7)).
To derive (1.5), we first reduce the Langevin dynamics (1.4) to an effective Brownian dynamics which captures weak inertia (§ 2). We then apply the method of homogenisation, formulate its cell problem in § 3, and solve it asymptotically in § 4, recasting the results of Childress (Reference Childress1979), Rosenbluth et al. (Reference Rosenbluth, Berk, Doxas and Horton1987), Shraiman (Reference Shraiman1987) and Soward (Reference Soward1987) in the language of homogenisation along the way. We discuss the results, derive simplified versions of (1.5) valid when inertia dominates diffusion (
$St \, Pe \gg 1$), and conclude in § 5. The reader uninterested in the details of the computation may want to skip to § 5.
2. Effective Brownian dynamics
In this section, we derive an effective Brownian dynamics capturing inertial effects in the distinguished regime 
$St, \, Pe^{-1} \ll 1$ with 
$St \, Pe = O(1)$ and 
$\beta =O(1)$. The derivation is similar to that carried out in Pavliotis et al. (Reference Pavliotis, Stuart and Zygalakis2009) in the case 
$\beta = 0$. It is convenient to introduce a small parameter 
$\epsilon \ll 1$ such that 
$Pe^{-1}=\epsilon ^2$, let 
$St=\gamma \epsilon ^2$ with 
$\gamma =O(1)$, and define the rescaled relative velocity 
$\boldsymbol {P}=\gamma \epsilon (\boldsymbol {V}-\boldsymbol {u})$ to rewrite (1.4) as
\begin{gather} \textrm{d}\boldsymbol{X}=(\boldsymbol{u}+\boldsymbol{P}/(\gamma\epsilon))\,\textrm{d} t, \end{gather}
\begin{gather} \textrm{d}\boldsymbol{P}-(\boldsymbol{P}/(\gamma\epsilon^2)+\gamma \epsilon(1-\beta)\textit{D}_t\boldsymbol{u}-(\boldsymbol{P}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u})\,\textrm{d} t+\sqrt{2} \, \textrm{d} \boldsymbol{W}. \end{gather}We reduce the dynamics of (2.1) by considering the corresponding backward Kolmogorov equation, namely
\begin{equation} \partial_t g-\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}g=\left(\frac{1}{\gamma\epsilon}\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}} -\left(\frac{1}{\gamma \epsilon^2}\boldsymbol{p}+\gamma \epsilon(1-\beta)\textit{D}_t\boldsymbol{u}-\left(\,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}\right)\boldsymbol{u}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}+{\nabla}_{{p}}^2\right)g \end{equation}for functions
\begin{equation} g(\boldsymbol{x},\boldsymbol{p},t)=\mathbb{E}[\,f(\boldsymbol{X}(t),\boldsymbol{P}(t))|\boldsymbol{X}(0)=\boldsymbol{x},\boldsymbol{P}(0)=\boldsymbol{p}], \end{equation}
where 
$\mathbb {E}$ denotes the expectation over the Wiener process 
$\boldsymbol {W}$ and 
$f$ is an arbitrary function (see e.g. Evans Reference Evans2013). We now introduce the multiple time scales 
$t_n=\epsilon ^n t$ with 
$n=0,1,\ldots$ and the expansion
\begin{equation} g(\,\boldsymbol{p},\boldsymbol{x},t)=g_0(\boldsymbol{x},t_0,t_1,t_2,\ldots)+\epsilon g_1(\,\boldsymbol{p},\boldsymbol{x},t_0,t_1,t_2,\ldots)+\cdots. \end{equation}
Substituting (2.4) into (2.2) yields, up to order 
$\epsilon ^2$,
\begin{gather} \frac{\boldsymbol{p}}{\gamma }\boldsymbol{\cdot}(\boldsymbol{\nabla}_{\boldsymbol{x}}g_0-\boldsymbol{\nabla}_{\boldsymbol p}g_1)=0,\end{gather}
\begin{gather} \frac{\boldsymbol{p}}{\gamma}\boldsymbol{\cdot}(\boldsymbol{\nabla}_{\boldsymbol{x}}g_1-\boldsymbol{\nabla}_{\boldsymbol{p}} g_2)=(\partial_{t_0}-\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})g_0,\end{gather}
\begin{gather} \frac{\boldsymbol{p}}{\gamma }\boldsymbol{\cdot}(\boldsymbol{\nabla}_{\boldsymbol{x}}g_2-\boldsymbol{\nabla}_{\boldsymbol p}g_3)=(\partial_{t_0}-\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})g_1+\partial_{t_1}g_0-{\nabla}_{{p}}^2g_1-(\,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}g_1 ,\end{gather}
\begin{align} \frac{\boldsymbol{p}}{\gamma }\boldsymbol{\cdot}(\boldsymbol{\nabla}_{\boldsymbol{x}}g_3-\boldsymbol{\nabla}_{\boldsymbol p}g_4)&=(\partial_{t_0}-\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})g_2+\partial_{t_1}g_1 - {\nabla}_{{p}}^2g_2- (\,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}g_2 \nonumber\\ &\quad +\partial_{t_2}g_0+\gamma (1-\beta)\textit{D}_t\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}g_0. \end{align}
Using that 
$g_0$ is independent of 
$\boldsymbol {p}$, we can solve (2.5) successively. Solving (2.5a) then (2.5b) yields
\begin{gather} g_1=\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}g_0,\quad g_2=\frac{\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}}{2}g_1,\end{gather}
\begin{gather} \partial_{t_0}g_0 =\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}g_0.\end{gather}Introducing (2.6a) into (2.5c) and using (2.6b) gives
\begin{equation} \frac{1}{\gamma}\left(\frac{(\,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}})^3}{2}g_0-\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}}g_3\right)=\partial_{t_1}g_0, \end{equation}
which is satisfied by 
$g_3={(\boldsymbol {p}\boldsymbol {\cdot }\boldsymbol {\nabla }_{\boldsymbol {x}})^3}g_0/{6}$ and
\begin{equation} \partial_{t_1}g_0=0. \end{equation}
Finally, introducing (2.6a) into (2.5d), using (2.6b) and (2.8) and setting 
$\boldsymbol {p}=0$ yields
\begin{equation} \partial_{t_2}g_0=-\gamma(1-\beta)\textit{D}_t\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}g_0+ {\nabla}_{{x}}^2g_0. \end{equation} From (2.6b), (2.8) and (2.9), it is clear that advection by 
$\boldsymbol {u}$ is the dominant process while inertia and diffusion arise at the same, lower order. We reconstitute the dynamics of 
$g_0$ capturing all three effects by adding (2.6b), (2.8) and (2.9), setting 
$\partial _{t_0}+\epsilon \partial _{t_1}+\epsilon ^2\partial _{t_2}= \partial _t$ and using 
$\epsilon ^2 = Pe^{-1}$ and 
$\gamma \epsilon ^2=St$ to obtain the backward Kolmogorov equation
\begin{equation} \partial_t g_0=(\boldsymbol{u}-St(1-\beta)\textit{D}_t\boldsymbol{u})\boldsymbol{\cdot} \boldsymbol{\nabla} g_0+\frac{1}{Pe}\nabla^2 g_0 \end{equation}
for 
$g_0$. This corresponds to the effective Brownian (or overdamped Langevin) dynamics
\begin{equation} \textrm{d}\boldsymbol{X}=(\boldsymbol{u}-St(1-\beta)\textit{D}_t\boldsymbol{u})\,\textrm{d} t+(2/Pe)^{1/2}\,\textrm{d}\boldsymbol{W}, \end{equation}
which includes inertial correction through the term 
$St(1-\beta ){D}_t\boldsymbol {u}$. In the absence of diffusion (
$Pe\to \infty$), (2.11) recovers the so-called first-order Eulerian approximation (e.g. Ferry & Balachandar Reference Ferry and Balachandar2001) or equivalently the slow manifold dynamics discussed by Rubin, Jones & Maxey (Reference Rubin, Jones and Maxey1995) and Haller & Sapsis (Reference Haller and Sapsis2008).
According to (2.11), inertial particles behave as though advected by the effective flow
\begin{equation} \boldsymbol{u}_{e}=\boldsymbol{u}-St(1-\beta)\textit{D}_t\boldsymbol{u}. \end{equation}
For incompressible fluids, 
$\boldsymbol {u}$ is divergence free but 
$\boldsymbol {u}_{e}$ rarely is. It can be shown that
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{e}=-St(1-\beta)(\|\boldsymbol{S}\|^2-\|\boldsymbol{\varOmega}\|^2), \end{equation}
where 
$\boldsymbol {S}$ and 
$\boldsymbol {\varOmega }$ are the local strain rate and rotation rate tensors of 
$\boldsymbol {u}$. Thus particles denser than the fluid (
$\beta <1$) tend to accumulate (
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_{e}<0$) in low vorticity and high strain regions, while less-dense particles (
$\beta >1$) tend to accumulate in high vorticity and low strain regions. A detailed analysis of this clustering effect without diffusion is given by Sapsis & Haller (Reference Sapsis and Haller2010). Isodense particles (
$\beta =1$) are unaffected by inertia to the order of accuracy of (2.12). Non-zero higher-order terms resulting from the coupling between inertia and noise appear when the expansion of the backward Kolmogorov equation is carried out to 
$O(St^2)$; these involve derivatives of 
$g$ of order higher than 2, characterising a non-diffusive correction to the trajectories 
$\boldsymbol {X}(t)$. Note that differences between the behaviour of rigid isodense particles and the fluid they replace also arise from finite-size (Faxén) effects which we neglect.
3. Periodic homogenisation
In this section, we apply the methodology of homogenisation to the Brownian dynamics (2.11) to obtain an expression for the corresponding effective diffusivity. The derivation is standard (see e.g. Vergassola & Avellaneda Reference Vergassola and Avellaneda1997; Pavliotis & Stuart Reference Pavliotis and Stuart2005) and recorded here for completeness and to set up notation.
We rewrite the associated backward Kolmogorov equation (2.10) as
\begin{equation} \partial_tg=\boldsymbol{u}_{e}\boldsymbol{\cdot}\boldsymbol{\nabla} g+\frac{1}{Pe}\nabla^2g, \end{equation}
where 
$\boldsymbol {u}_{e}=\boldsymbol {u}-St(1-\beta )\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}$ is the effective velocity field (2.12) which is steady, periodic and divergent. We introduce the small parameter 
$\delta \ll 1$ along with the variables 
$\boldsymbol {X}=\delta \,\boldsymbol {x}$, 
$t_1=\delta \,t$ and 
$t_2=\delta ^2t$. We seek a solution of (3.1) in the form
\begin{equation} g(\boldsymbol{x},t)=g_0(\boldsymbol{X},t_1,t_2)+\delta g_1(\boldsymbol{x},\boldsymbol{X},t_1,t_2)+\delta^2g_2(\boldsymbol{x},\boldsymbol{X},t_1,t_2) + \cdots\,, \end{equation}
where the 
$g_i, \, i=1,2,\ldots ,$ are periodic functions of 
$\boldsymbol {x}$. We introduce the expansion (3.2) into (3.1) and collect terms at each order in 
$\delta$. At orders 
$\delta$ and 
$\delta ^2$, we find
\begin{gather} \partial_{t_1}g_0=\mathcal{L}_1g_0+\mathcal{L}_0g_1,\end{gather}
\begin{gather} \partial_{t_2}g_0=\mathcal{L}_0 g_2+\mathcal{L}_1 g_1+\mathcal{L}_2g_0, \end{gather}with
\begin{equation} \mathcal{L}_0=\boldsymbol{u}_{e}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}+\frac{1}{Pe} {\nabla}_{{x}}^2,\quad \mathcal{L}_1=\boldsymbol{u}_{e}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{X}}+\frac{2}{Pe}\boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}\quad \text{and}\quad \mathcal{L}_{2}=\frac{1}{Pe} {\nabla}_{{X}}^2. \end{equation}
Let 
$\phi ^{\dagger }(\boldsymbol {x})$ be the invariant measure associated with 
$\mathcal {L}_0$, that is, the periodic solution of
\begin{equation} \mathcal{L}_0^{\dagger}\phi^{\dagger}= - \boldsymbol{\nabla}_x \boldsymbol{\cdot} (\boldsymbol{u}_e \phi^\dagger) +\frac{1}{Pe} {\nabla}_{{x}}^2\phi^{\dagger}=0 \end{equation}
normalised such that 
$\langle \phi ^{\dagger }\rangle =1$, with 
$\langle \boldsymbol{\cdot} \rangle$ denoting spatial average over the periodic cell. Multiplying (3.3a) by 
$\phi ^{\dagger }$ and averaging yields
\begin{equation} \partial_{t_1}g_0=\boldsymbol{c}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{X}}g_0, \end{equation}where
\begin{equation} \boldsymbol{c}=\langle\phi^{\dagger}\boldsymbol{u}_{e}\rangle \end{equation}is the effective drift velocity vector.
Using (3.5) and (3.3a), 
$g_1$ is found to satisfy
\begin{equation} \mathcal{L}_0g_1+\left(\boldsymbol{u}_{e}-\boldsymbol{c}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{X}}g_0=0. \end{equation}
The solution takes the form 
$g_1=\boldsymbol {\varphi }\boldsymbol {\cdot }\boldsymbol {\nabla }_{\boldsymbol {X}}g_0$, where 
$\boldsymbol {\varphi }$ obeys the cell problem
\begin{equation} \mathcal{L}_0\boldsymbol{\varphi}+\boldsymbol{u}_{e}-\boldsymbol{c}=0 \end{equation}
with periodic boundary conditions. Substituting this solution into (3.3b), multiplying by 
$\phi ^{\dagger }$ and averaging yields the effective equation
\begin{equation} \partial_{t_2}g_0=\boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{\cdot}(\boldsymbol{\mathsf{D}}_{e}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{X}})g_0, \end{equation}where
\begin{equation} \boldsymbol{\mathsf{D}}_{e}=\left\langle\phi^{\dagger}\left(\frac{1}{Pe}\left(\mathbb{I}+2\boldsymbol{\nabla}_{\boldsymbol{x}}\otimes\boldsymbol{\varphi}\right)+\left(\boldsymbol{u}_{e}-\boldsymbol{c}\right)\otimes\boldsymbol{\varphi}\right)\right\rangle, \end{equation}
is the effective diffusivity tensor, with 
$\mathbb {I}$ denoting the identity tensor and 
$\otimes$ the tensorial product. Finally, we gather (3.5) and (3.9), set 
$\delta \partial _{t_1}+\delta ^2\partial _{t_2}=\partial _t$ and 
$\delta \boldsymbol {\nabla }_{\boldsymbol {X}} = \boldsymbol {\nabla }$ to obtain the effective backward Kolmogorov equation
\begin{equation} \partial_t g_0=\boldsymbol{c}\boldsymbol{\cdot}\boldsymbol{\nabla} g_0+\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\boldsymbol{\mathsf{D}}_{e}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)g_0. \end{equation}
This corresponds to a diffusive process, characterised by the effective drift velocity 
$\boldsymbol {c}$ and diffusivity tensor 
$\boldsymbol{\mathsf{D}}_{e}$, which approximates the dynamics (2.11) over long time scales.
Computing 
$\boldsymbol {c}$ and 
$\boldsymbol{\mathsf{D}}_{e}$ requires solving both (3.4) for 
$\phi ^{\dagger }$ and (3.8) for 
$\boldsymbol {\varphi }$. In the non-inertial case, 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}_{e}=0$ and the invariant measure 
$\phi ^{\dagger }$ is simply a constant. For inertial particles, 
$\phi ^{\dagger }$ is non-trivial and reflects the clustering caused by the divergence of 
$\boldsymbol {u}_{e}$.
We remark that we obtained (3.6) and (3.10) by first deriving the effective Brownian dynamics (2.11) then applying homogenisation. Alternatively, we could have first applied homogenisation to the Langevin dynamics (1.4), then exploited the asymptotic parameters to simplify the cell problems. We have checked that this alternative route yields the same cell problems; in other words, the limits 
$St \to 0$ and 
$\delta \to 0$ commute. Pavliotis & Stuart (Reference Pavliotis and Stuart2005) and Martins Afonso, Mazzino & Muratore-Ginanneschi (Reference Martins Afonso, Mazzino and Muratore-Ginanneschi2012) take 
$\delta \to 0$ first followed by 
$St \to 0$, but they do not consider the small-
$Pe$ asymptotics. We note that the cell problems for the Langevin dynamics are difficult to solve as they are defined within domains including unbounded velocity spaces. The numerical sampling of trajectories is also delicate for the Langevin dynamics when 
$St \ll 1$ since it requires exceedingly small time steps; we therefore make use of the effective Brownian dynamics for this purpose below.
In the next section, we solve the cell problems (3.4) and (3.8) asymptotically in the distinguished regime 
$St=O(Pe^{-1})\ll 1$ and compute the leading-order expression of the effective diffusivity (3.10) for the cellular flow.
4. Effective diffusivity
In non-dimensional variables, the effective velocity is
\begin{equation} \boldsymbol{u}_{e} = \boldsymbol{u} - \alpha \epsilon^2 \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \end{equation}with
\begin{equation} \boldsymbol{u} = (\partial_y \psi,-\partial_x \psi) \quad \text{and}\quad \psi(x,y)=\sin x\sin y, \end{equation}
where 
$\alpha = (1-\beta ) \gamma = (1-\beta ) St \, Pe = O(1)$. The inertial contribution turns out to have a gradient structure: we can write (4.1) as
\begin{equation} \boldsymbol{u}_{e} = \boldsymbol{u} + \alpha\epsilon^2\boldsymbol{\nabla}\varPhi \end{equation}with
\begin{equation} \varPhi(x,y)=-(\cos(2x)+\cos(2y) )/4. \end{equation} Figure 3 displays the 
${\rm \pi}$-periodic potential 
$\varPhi$. In effect, inertia adds a weak potential flow which, depending on the sign of 
$\alpha$, attracts particles to either the minima or the maxima of 
$\varPhi$ located at the corners or centres of the quarter-cells.
Potential 
$\varPhi$ defined in (4.2) (
$a$) and invariant measure 
$\phi ^{\dagger }$ for 
$\alpha =-10$, 
$0$ and 
$100$ computed from (4.8) (
$b$). Both fields are 
${\rm \pi}$-periodic and shown here in a quarter-cell. The thin dark lines represent the streamlines 
$|\psi |=0$, 
$0.25$ and 
$0.75$.
To compute the effective diffusivity (3.10), we first derive the leading-order approximation to 
$\phi ^{\dagger }$, then to 
$\boldsymbol {\varphi }$, and finally combine the resulting expressions.
4.1. Asymptotic computation of 
$\phi ^{\dagger }$
We rewrite (3.4) for 
$\phi ^{\dagger }$ as
\begin{equation} \epsilon^2 \boldsymbol{\nabla}\phi^{\dagger}-\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^{\dagger}+\alpha\epsilon^2\boldsymbol{\nabla}\boldsymbol{\cdot}(\phi^{\dagger}\,\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u})=0 \end{equation}
and expand 
$\phi ^{\dagger }$ in powers of 
$\epsilon ^2$ according to 
$\phi ^{\dagger }=\phi ^{\dagger }_0+\epsilon ^2\phi ^{\dagger }_1+\cdots$ to obtain
\begin{gather} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^{\dagger}_0= 0, \end{gather}
\begin{gather} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^{\dagger}_1=\nabla^2\phi^{\dagger}_{0}+\alpha\boldsymbol{\nabla}\boldsymbol{\cdot}(\phi^{\dagger}_0 \, \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}). \end{gather}
Equation (4.4a) implies that 
$\phi ^{\dagger }_0$ must be constant along streamlines: 
$\phi ^{\dagger }_0=\phi _0^{\dagger }(\psi )$. Equation (4.4b) can be solved for 
$\phi ^{\dagger }_1$ provided that a solvability condition, obtained by integrating (4.4b) along streamlines, is satisfied. We show in appendix A that this solvability condition is
\begin{equation} \frac{\textrm{d}}{\textrm{d}\psi}\left(a(\psi)\frac{\textrm{d}\phi^{\dagger}_0}{\textrm{d}\psi}+\alpha b(\psi)\phi^{\dagger}_0\right)=0, \end{equation}where
\begin{equation} a(\psi)=\oint_{\psi}\,\textrm{d} l\, |\boldsymbol{u}| \quad \text{and}\quad b(\psi)=\oint_{\psi}\,\textrm{d} l\, \boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}), \end{equation}
with 
$l$ denoting arclength and 
$\boldsymbol {n}=\boldsymbol {\nabla }\psi /|\boldsymbol {\nabla }\psi |$ the normal to the streamline. Note that 
$a(\psi )$ can be recognised as the circulation around each streamline, while 
$b(\psi )$ is the flux induced by inertial effects through the streamline and can be rewritten as
\begin{equation} b(\psi) = \oint_{\psi}\,\textrm{d} l\, k |\boldsymbol{u}|^2 = \iint_{D_{\psi}}\,\textrm{d} \boldsymbol{x}\, \det \boldsymbol{\nabla} \boldsymbol{\nabla} \psi , \end{equation}
where 
$k$ denotes the curvature of the streamline, 
$D_{\psi }$ the area enclosed by the streamline and 
$\boldsymbol {\nabla } \boldsymbol {\nabla } \psi$ the Hessian. We integrate (4.5) once, taking the integration constant to 
$0$ as required for a regular solution at the extrema 
$\psi = {\pm }1$. Integrating once more we find that the (bounded) normalised solution to (4.5) is
\begin{equation} \phi^{\dagger}_0\left(\psi\right)=\frac{1}{Z(\alpha)}\exp\left(-\alpha\int_0^{\psi}\,\textrm{d} \psi'\,\frac{b(\psi')}{a(\psi')}\right), \end{equation}
where 
$Z(\alpha )$ is a normalisation constant. After some algebra detailed in appendix A, we obtain
\begin{equation} Z(\alpha)=\frac{1}{{\rm \pi}^2}\int_0^1\,\textrm{d}\psi\,c(\psi)\exp\left(-\alpha\int_0^{\psi}\,\textrm{d}\psi'\frac{b(\psi')}{a(\psi')}\right), \end{equation}where
\begin{equation} c(\psi)=\oint_{\psi}\frac{\textrm{d} l}{|\boldsymbol{u}|} \end{equation}is the orbital time around the closed streamlines.
Using the form (1.1) for the streamfunction we can evaluate the integrals 
$a(\psi )$, 
$b(\psi )$ and 
$c(\psi )$ in terms of (complementary) complete elliptic integrals to find
\begin{gather} a(\psi)=8(E'(\psi)-\psi^2\,K'(\psi)), \end{gather}
\begin{gather} b(\psi)=8\psi(K'(\psi)-E'(\psi)), \end{gather}
\begin{gather} c(\psi)=4\,K'(\psi), \end{gather}
with 
$K'(\psi )=K(\sqrt {1-\psi ^2})$ and 
$E'(\psi )=E(\sqrt {1-\psi ^2})$ (see DLMF Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Miller, Saunders, Cohl and McClain2019). This yields the closed-form
\begin{equation} Z(\alpha)=\frac{4}{{\rm \pi}^2}\int_{0}^1\,\textrm{d}\psi\,K'(\psi)\exp\left(-\alpha\int_{0}^{\psi}\,\textrm{d}\psi'\,\frac{\psi'\left(K'(\psi')-E'(\psi')\right)}{E'(\psi')-\psi'^2K'(\psi')}\right). \end{equation} Figure 3 shows the function 
$\phi ^{\dagger }_0$ for three values of 
$\alpha$ in one quarter-cell. This is of interest because 
$\phi ^{\dagger }$ describes locally the spatial structure of the particle density which is given by 
$\phi ^{\dagger }(\boldsymbol {x})$ times a large-scale diffusive envelope. The function is even and monotonic with 
$|\psi |$, decreasing from the centre of quarter-cells for light particles (
$\alpha <0$) and increasing from the centre for heavy particles (
$\alpha >0$) reflecting the expected clustering induced by inertia.
4.2. Asymptotic computation of ${\varphi}$
It is clear from the symmetries of the cellular flow that 
$\boldsymbol {c}=0$ and that the effective diffusivity is isotropic. Therefore, we focus on a single component of 
$\boldsymbol {\varphi }$, say 
$\varphi =\boldsymbol {\varphi }\boldsymbol {\cdot }\boldsymbol {e}_x$, which satisfies
\begin{equation} \epsilon^2\nabla^2\varphi+(\boldsymbol{u}-\alpha\epsilon^2 \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{\nabla}\varphi+u-\alpha\epsilon^2 \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} u=0, \end{equation}
where 
$u=\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {e}_x$. Making use of the symmetries 
$(x,y,\varphi )\mapsto ({\rm \pi} -y,{\rm \pi} -x,-\varphi )$, 
$(x,y,\varphi )\mapsto (x,-y,\varphi )$ and 
$(x,y,\varphi )\mapsto (-x,y,-\varphi )$ we can focus on the quarter-cell 
$[0,{\rm \pi} ]^2$ and look for a solution satisfying the boundary conditions
\begin{equation} \varphi(0,y)=\varphi({\rm \pi},y)=0\quad \text{and}\quad\partial_y\varphi(x,0)=\partial_y\varphi(x,{\rm \pi})=0. \end{equation}The general solution of the problem is identical to the solution of this simpler problem up to an irrelevant constant.
We now introduce 
$\varTheta =2(\varphi +x)/{\rm \pi} -1$ satisfying
\begin{equation} \epsilon^2\nabla^2\varTheta+(\boldsymbol{u}-\alpha\epsilon^2 (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{\nabla}\varTheta=0, \end{equation}with the boundary conditions
\begin{equation} \varTheta(0,y)=-1,\quad \varTheta({\rm \pi},y)=1,\quad \partial_y\varTheta(x,0)=0\quad \text{and}\quad\partial_y\varTheta(x,{\rm \pi})=0. \end{equation}
We obtain an approximation for 
$\varTheta$ for 
$\epsilon \ll 1$ using matched asymptotics.
4.2.1. Interior solution
We first consider the solution in the quarter-cell interior. Introducing the expansion 
$\varTheta =\varTheta _0+\epsilon ^2 \varTheta _1+\cdots$, we obtain 
$\varTheta _0=\varTheta _0(\psi )$ at leading order and, at the next order,
\begin{equation} \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\varTheta_1=-\nabla^{2}\varTheta_0+\alpha ( \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} ) \boldsymbol{\cdot}\boldsymbol{\nabla}\varTheta_0. \end{equation}Integrating along a streamline yields the solvability condition
\begin{equation} \frac{\textrm{d}}{\textrm{d}\psi}\left(a(\psi)\frac{\textrm{d}\varTheta_0}{\textrm{d}\psi}\right)-\alpha b\left(\psi\right)\frac{\textrm{d}\varTheta_0}{\textrm{d}\psi}=0, \end{equation}
with 
$a(\psi )$ and 
$b(\psi )$ given in (4.11). The only bounded solution of (4.18) is a constant. This constant interior solution must be matched with a boundary layer solution around the separatrix 
$\psi =0$ so as to satisfy the boundary conditions (4.16a–d).
4.2.2. Boundary layer solution
Following Childress (Reference Childress1979), we introduce the variables
\begin{equation} \zeta=\epsilon^{-1}\psi \quad \text{and}\quad \sigma=-\int_0^{l}\,\textrm{d} l\,|\boldsymbol{u}|, \end{equation}
where 
$l$ is the arclength along a streamline. At the separatrix, 
$0<\sigma <8$ parameterises the boundary of the quarter-cell clockwise, starting from 
$0$ at 
$(0,0)$ and taking values 
$2$, 
$4$ and 
$6$ at successive corners. Introducing the expansion 
$\varTheta =\varTheta _0(\sigma ,\zeta )+O(\epsilon )$, we obtain at leading order, away from the corners,
\begin{equation} \partial_{\zeta}^2\varTheta_0-\partial_{\sigma}\varTheta_0=0, \end{equation}with the boundary conditions
\begin{gather} \varTheta_0(\sigma,0)=-1 \quad \text{for} \ \sigma\in[0,2], \end{gather}
\begin{gather} \varTheta_0(\sigma,0)=1 \quad \text{for} \ \sigma\in[4,6], \end{gather}
\begin{gather} \partial_{\zeta}\varTheta_0(\sigma,0)=0 \quad \text{for} \ \sigma\in[2,4]\cup[6,8]. \end{gather}
Equations (4.20) make up the so-called Childress problem, solved in closed form by Soward (Reference Soward1987). Using the symmetry of the boundary conditions (4.20b) and (4.20c), it can be shown that 
$\varTheta _0$ can only match the constant interior solution if this vanishes. Thus we conclude that, to leading order, 
$\varTheta$ is non-zero only within the boundary layer. As a result, as we now show, we can compute the effective diffusivity by combining (4.8) for 
$\phi ^{\dagger }$ with the solution of the Childress problem (4.20).
4.3. Computation of the effective diffusivity $\boldsymbol{\mathsf{D}}_{e}$
Using the symmetry of the problem, (3.8) and integration by parts, the effective diffusivity (3.10) can be recast into the isotropic form
\begin{equation} \boldsymbol{\mathsf{D}}_{e}=\bar{D}\,\mathbb{I} \end{equation}with
\begin{equation} \bar{D}=\epsilon^2\langle\phi^{\dagger}(2\partial_{x}\varphi+|\boldsymbol{\nabla}\varphi|^2+1)\rangle. \end{equation}
Since the four quarter-cells are equivalent, we can focus on 
$[0,{\rm \pi} ]^2$ where 
$\varphi ={\rm \pi} (1+\varTheta )/2-x$ to obtain
\begin{equation} \bar{D}=\frac{{\rm \pi}^2\epsilon^2}{4}\langle\phi^{\dagger}|\boldsymbol{\nabla}\varTheta|^2\rangle_{[0,{\rm \pi}]^2}, \end{equation}
where 
$\langle \cdot \rangle _{[0,{\rm \pi} ]^2}$ denotes the spatial average over the quarter-cell 
$[0,{\rm \pi} ]^2$. In the limit 
$\epsilon \ll 1$, 
$\varTheta$ becomes a boundary layer term localised in the 
$O(\epsilon )$ boundary layer around the separatrix 
$\psi =0$. Within this boundary layer, 
$\phi ^{\dagger }\sim Z(\alpha )^{-1}$ as (4.8) shows. As a result, to leading order in 
$\epsilon$, the effective diffusivity reduces to
\begin{equation} \bar{D} \sim \frac{\epsilon}{4Z(\alpha)}\int_{0}^8\,\textrm{d}\sigma\int_0^{\infty}\,\textrm{d}\zeta(\partial_{\zeta}\varTheta_0)^2, \end{equation}
using that 
$\textrm {d} x \textrm {d} y = \epsilon ^{-1}|\boldsymbol {\nabla } \psi |^{-2}\,\textrm {d} \sigma \, \textrm {d} \zeta$. In appendix B, we relate the integral in (4.23) to Soward's constant 
${\mathscr {S}}$ appearing in the asymptotic expression of the non-inertial effective diffusivity (1.2) to rewrite the effective diffusivity in the fully explicit form (1.5). Note that it is possible to bypass the computation in appendix B by observing that the integral in (4.23) is independent of 
$\alpha$ so that it can be determined by setting 
$\alpha =0$ in (4.23) and matching to the familiar non-inertial result (1.2), noting that 
$Z(0)=1$.
5. Discussion and conclusion
The key result of the paper is the asymptotic expression (1.5) for the effective diffusivity 
$\bar {D}$ in the distinguished regime 
$St = O(Pe^{-1}) \ll 1$. We test it against the direct numerical sampling of the dynamics of the particles. The Langevin equation (1.4) can be costly for 
$St \ll 1$ so we restrict its use to moderately small 
$St$; for smaller 
$St$, we instead integrate the effective Brownian dynamics (2.11). We use numerical schemes based on those of Pavliotis et al. (Reference Pavliotis, Stuart and Zygalakis2009) which split the flow over one time step between two advection steps (each carried out exactly using the variables 
$x+y$ and 
$x-y$) and a diffusion step. The exact area-preservation of the advection steps is essential to enable large time steps and long integration times. The results are summarised in figure 4. They demonstrate a very good agreement between direct estimations and asymptotic predictions of 
$\bar {D}$. The near coincidence of estimates obtained with the Langevin and effective Brownian dynamics also confirms the validity of the latter.
Effective diffusivity 
$\bar {D}$ as a function of 
$St$ for 
$Pe=1000$ (
$a$) and as a function of 
$Pe$ for 
$St=0.1$ (
$b$). The values estimated by direct numerical sampling of the Langevin dynamics (1.4) for 
$\beta =0$ (
$\star$), the effective Brownian dynamics (2.11) for 
$\beta =0$ (
$\triangledown$) and 
$\beta =2$ (
$\vartriangle$), and the non-inertial dynamics (1.3) (
$\circ$) are compared with the asymptotic prediction (1.5). The simplified asymptotic expressions (5.1) and (5.2) valid for 
$St \, Pe \gg 1$ are shown (dotted lines), as well as the refinement of (5.1) that includes the logarithmic correction in (C 6) (dot-dashed line).
The impact of inertia is entirely captured by the function 
$Z(\alpha )$ defined by (4.12) and shown in figure 2. We now analyse the behaviour of this function in detail. As noted in § 1, 
$Z(\alpha )$ is decreasing with 
$\alpha$, with 
$Z(\alpha )<1$ for 
$\alpha <0$ and 
$Z(\alpha )>1$ for 
$\alpha >0$. Consequently, the effective diffusivity of particles denser (respectively less dense) than the fluid is always larger (respectively smaller) than that of non-inertial particles. This can be attributed to the divergence (2.13) of the effective velocity which, when averaged along streamlines, leads to an accumulation (respectively depletion) of particles in the separatrix region (cf. figure 3), which controls the cell-to-cell and hence global transport.
It is interesting to consider the limiting behaviour of 
$Z(\alpha )$ and hence 
$\bar {D}$ as 
$\alpha = St \, Pe (1-\beta ) \to \pm \infty$, that is, when inertia dominates over diffusion. In the heavy-particle case 
$\alpha \to \infty$, the asymptotics of 
$Z(\alpha )$ derived in appendix C gives
\begin{equation} \bar{D} \sim {\mathscr{S}}{\rm \pi}^{3/2}\left(\frac{St(1-\beta)}{\log \left(Pe \, St (1-\beta) \right)}\right)^{1/2} \quad \text{for} \ Pe^{-1}\ll St \ll 1 \ \text{and} \ \beta < 1. \end{equation}
This corresponds to an effective diffusivity that depends only weakly on the Péclet number or, equivalently, on molecular diffusivity, and is instead controlled by inertia through the dependence on 
$St^{1/2}$. Note that the asymptotics (5.1) is rather poor because it ignores logarithmic corrections that are negligible only for exceedingly large 
$\alpha$. A more accurate formula can be obtained by using an improved approximation to 
$Z(\alpha )$ given in (C 6) and shown in figure 2. The predictions of (5.1) and its improvement are compared against simulation results and the full asymptotic approximation (1.5) in figure 4.
In the light-particle case 
$\alpha \to - \infty$, the asymptotics of 
$Z(\alpha )$ in appendix C gives
\begin{equation} \bar{D} \sim {\mathscr{S}}{\rm \pi} Pe^{1/2}St(\beta-1)\,\exp({-\varUpsilon \, St\, Pe(\beta-1)}) \quad \text{for} \ Pe^{-1}\ll St \ll 1 \ \text{and} \ \beta > 1, \end{equation}
where 
$\varUpsilon \approx 0.655$, which is also shown in figure 4. The effective diffusivity decreases exponentially with 
$Pe \, St$ corresponding to a dramatic inhibition of dispersion caused by inertia. Physically, particles are trapped by inertia near the quarter-cell centres and only escape by crossing the separatrix for rare realisations of the noise. The rate of these rare escapes and hence 
$\bar {D}$ could be estimated using small-noise large-deviation theory (e.g. Freidlin & Wentzell Reference Freidlin and Wentzell2012). Note that the key parameter 
$St \, Pe = U^2 \tau /D$ (which is independent from the flow scale 
$a$) can be rewritten as the Arrhenius-like number
\begin{equation} St \, Pe = {m U^2}/(k_B T), \end{equation}
where 
$m=m_{p}+m_{f}/2$ is the effective mass of the particles, 
$k_{B}$ is the Boltzmann constant and 
$T$ the temperature, on using the Einstein–Smoluchowski relation 
$D = \tau k_{B} T / m_{p}$. Thus (5.2) can be interpreted as a form of Arrhenius law, with the particle kinetic energy playing the role of activation energy.
We can assess the relative importance of inertia and diffusion by thermal noise for particles in a flow using (5.3). For instance, at room temperature, for particles with an effective density similar to that of water, we find that 
$St\,Pe \approx r_{p}^3 U^2 \times 10^{24}$, in SI units. Thus, inertia dominates diffusion (
$St\,Pe \gtrsim 1$) for particles of radius 
$r_{p}=10^{-6}\ \textrm {m}$ and 
$10^{-3}\ \textrm {m}$ when 
$U \gtrsim 10^{-3}\ \textrm {m}\,\textrm {s}^{-1}$ and 
$U \gtrsim 10^{-8}\ \textrm {m}\,\textrm {s}^{-1}$, respectively. Since typical flow velocities are likely to exceed these small values, the simplified asymptotic expressions (5.1) and (5.2) will be valuable. We emphasise that, though subdominant, diffusion plays a crucial role in setting the effective diffusivity of cellular flows, as the dependence in 
$Pe$ of (5.1) and (5.2) indicates. This is because diffusion is indispensible for particles to move between cells. We also note that, while (5.3) is restricted to diffusion by thermal noise, our results are also relevant to dispersion problems where diffusion models small-scale turbulent mixing, in which case the (turbulent) diffusivity is orders of magnitude larger than the molecular one and 
$St \, Pe$ is not necessarily large.
We remark that the 
$Pe^{-1/2}$ scaling for the effective diffusivity of non-inertial particles in (1.2) is a universal feature of periodic flows with closed streamlines (Heinze Reference Heinze2003; Novikov et al. Reference Novikov, Papanicolaou and Ryzhik2005). We expect that the conclusion that this is corrected to 
$Pe^{-1/2}/Z(\alpha )$ when inertia is taken into account, which we draw for the cellular flow (1.1), also generalises. Specifically, we expect that
\begin{equation} \bar{D}(Pe,\alpha) = {\bar{D}(Pe,0)}/Z(\alpha) \end{equation}
for all cellular flows. The form of 
$Z(\alpha )$ is specific to each flow but its qualitative dependence on 
$\alpha$ and its asymptotic scalings for 
$\alpha \to \pm \infty$ are likely to be as in the case of the canonical cellular flow (1.1). To illustrate this, we consider another cellular flow in appendix D and confirm (5.4). Flows with open streamlines behave very differently, however. In shear flows, numerical results from Pavliotis et al. (Reference Pavliotis, Stuart, Band, Niederreiter and Talay2006) suggests that inertia only has a negligible effect on the (Taylor) effective diffusivity. It would be of interest to examine the impact of inertia on more complex flows such as the cat's-eye flows of Childress & Soward (Reference Childress and Soward1989).
It would also be desirable to assess the effect of the Boussinesq–Basset force which we neglected. Although Manton (Reference Manton1974) argues that it is less significant than other inertial effects, more recent work by Daitche & Tél (Reference Daitche and Tél2011) and Langlois, Farazmand & Haller (Reference Langlois, Farazmand and Haller2015) suggests otherwise.
Acknowledgements
This work was supported by EPSRC Programme Grant EP/R045046/1: Probing Multiscale Complex Multiphase Flows with Positrons for Engineering and Biomedical Applications (PI: Professor. M. Barigou, University of Birmingham). We thank the two anonymous referees for useful comments.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the solvability condition (4.5)
We introduce the time-like coordinate 
$s$ such that
\begin{equation} \frac{\textrm{d}}{\textrm{d} s}=\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}. \end{equation}Integrating (4.4b) along streamlines gives
\begin{equation} \oint_{\psi}\,\textrm{d} s\,\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\nabla}\phi_0+\alpha\phi_0\, \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u})=0. \end{equation}
Following Haynes & Vanneste (Reference Haynes and Vanneste2014), we use the fact that the change of variables 
$(x,y)\mapsto (s,\psi )$ is area preserving, and the divergence theorem, to write
\begin{equation} \oint_{\psi}\,\textrm{d} s\, \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{f}=\frac{\textrm{d}}{\textrm{d}\psi}\left(\iint\, \textrm{d} x\,\textrm{d} y\, \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{f}\right)=\frac{\textrm{d}}{\textrm{d}\psi}\left(\oint_{\psi}\,\textrm{d} s\, \boldsymbol{\nabla} \psi\boldsymbol{\cdot} \boldsymbol{f}\right) \end{equation}
for an arbitrary vector field 
$\boldsymbol {f}$. This reduces (A 2) to
\begin{equation} \frac{\textrm{d}}{\textrm{d}\psi}\left(a(\psi)\frac{\textrm{d}\phi_0}{\textrm{d}\psi}+\alpha b(\psi)\phi_0\right)=0, \end{equation}where
\begin{equation} a(\psi)=\oint_{\psi}\,\textrm{d} s\, |\boldsymbol{\nabla}\psi|^2\quad \text{and} \quad b(\psi)=\oint_{\psi}\,\textrm{d} s\,\boldsymbol{\nabla}\psi\boldsymbol{\cdot}\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right). \end{equation}
Using that 
$\textrm {d} s=|\boldsymbol {\nabla }\psi | \, \textrm {d} l$ with 
$l$ the arclength, we obtain (4.6a,b). The derivation of (4.11) for 
$a(\psi )$ from (1.1) can be found in Haynes & Vanneste (Reference Haynes and Vanneste2014). To compute 
$b(\psi )$, we use the symmetry of the streamline, (A 1) and 
$\psi = \sin x\sin y$, to write
\begin{equation} b(\psi)=\psi\int_{\sin^{-1}\psi}^{{\rm \pi}/2}\frac{\cos^2 x \, \textrm{d} x}{\sqrt{\sin^{2}x-\psi^2}}. \end{equation}
The substitution 
$t^2 = (1-\sin ^2x)/(1-\psi ^2)$ then gives
\begin{equation} b(\psi)=\psi(1-\psi^2)\int_{0}^{1}\frac{t^2 \, \textrm{d} t}{\sqrt{1-t^2}\sqrt{1-(1-\psi^2)t^2}}. \end{equation}Using formula (19.2.6) in DLMF (Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Miller, Saunders, Cohl and McClain2019), this reduces to the expression in (4.11).
Appendix B. Computation of the integral in (4.23)
We show that the integral
\begin{equation} \mathcal{K}=\int_0^8\,\textrm{d}\sigma\int_0^{\infty}\,\textrm{d}\zeta\,(\partial_{\zeta}\varTheta_0)^2, \end{equation}
where 
$\varTheta _0$ solves (4.20), is related to
\begin{equation} \int_{0}^{\infty}\,\textrm{d}\zeta\,\varTheta_0(2,\zeta)=-2{\mathscr{S}} \end{equation}with
\begin{equation} {\mathscr{S}}=\sqrt{\frac{2}{{\rm \pi}}}\sum_{n=0}^{\infty}\frac{(-1)^n}{\sqrt{2n+1}}, \end{equation}
as calculated by Soward (Reference Soward1987). Integrating (B 1) by parts in 
$\zeta$ gives
\begin{equation} \mathcal{K}=-\int_{0}^8\,\textrm{d}\sigma\,\varTheta_0(\sigma,0)\partial_{\zeta}\varTheta_0(\sigma,0)-\int_{0}^{8}\,\textrm{d}\sigma\int_0^{\infty}\,\textrm{d}\zeta\,\varTheta_0\partial_{\zeta}^2\varTheta_0. \end{equation}
The second term can be shown to vanish using (4.20) and periodicity in 
$\sigma$. Using the boundary conditions (4.20b) and (4.20c) reduces the first term to
\begin{equation} \mathcal{K}=\int_0^2\,\textrm{d}\sigma\,(\partial_{\zeta}\varTheta_0(\sigma,0)-\partial_{\zeta}\varTheta_0(\sigma+4,0)). \end{equation}
Now, integrating (4.20) for 
$\zeta \in [0,\infty )$ gives 
$\int _{0}^{\infty }\,\textrm {d}\zeta \,\partial _{\sigma }\varTheta _0=-\partial _{\zeta }\varTheta _0(\sigma ,0)$ which can be introduced in (B 4) to obtain
\begin{equation} \mathcal{K}=\int_{0}^{\infty}\,\textrm{d}\zeta\,(\varTheta_{0}(0,\zeta)-\varTheta_{0}(2,\zeta)+\varTheta_{0}(6,\zeta)-\varTheta_{0}(4,\zeta)). \end{equation}
Finally, using that 
$\int _{0}^{\infty }\,\textrm {d}\zeta \,\varTheta _0(\sigma ,\zeta )$ is constant for 
$\sigma \in [2,4]$ and 
$[4,8]$, and the symmetry 
$\varTheta _0(\sigma +4,\zeta )=-\varTheta _0(\sigma ,\zeta )$, we obtain
\begin{equation} \mathcal{K}=-4\int_0^{\infty}\,\textrm{d}\zeta\,\varTheta_0(2,\sigma)=8{\mathscr{S}}. \end{equation}Appendix C. Asymptotic calculation of $Z(\alpha )$ for $\alpha \to \pm \infty$
C.1. Limit $\alpha \to \infty$
The outer integral in (4.12) is dominated by a neighbourhood of the minimum of the integral multiplying 
$\alpha$. This minimum can be verified to be the left endpoint 
$\psi =0$. Using asymptotic expansions for 
$K'(\psi )$ and 
$E'(\psi )$ as 
$\psi \to 0^+$, we find the approximation
\begin{equation} Z(\alpha) \sim \frac{4}{{\rm \pi}^2}\int_0^1\,\textrm{d}\psi\,\left(-\log\frac{\psi}{4}\right)\exp\left(\alpha\frac{\psi^2}{4}\left(1+2\log\frac{\psi}{4}\right)\right). \end{equation}
Introducing the new integration variable 
$x = -\alpha \psi ^2(1+2\log (\psi /4))/4$ yields, after cumbersome calculations (carried out using the symbolic algebra software Mathematica),
\begin{equation} Z(\alpha) \sim \frac{2}{{\rm \pi}^2\sqrt{\alpha}}\int_0^{\alpha c} \frac{(w(x)+1)\sqrt{w(x)} \mathrm{e}^{-x}}{w(x)-1}\,\textrm{d} x, \end{equation}
where 
$c=\log 2-1/4>0$ and 
$w(x)=-W_{-1}(-\mathrm {e}\,x/(4\alpha ))$, with 
$W_{-1}$ the 
$-1$ branch of the Lambert 
$W$ function, namely the inverse function of 
$x\mathrm {e}^{x}$ on 
$(-\infty ,-1]$ (see DLMF Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Miller, Saunders, Cohl and McClain2019). Using the asymptotic expansion
\begin{align} W_{-1}(-x) = \log x+\log(-\log x)+\frac{\log(-\log x)}{\log x}+o\left(\frac{\log(\log x)}{\log x}\right) \quad \text{as} \ x \to 0^+, \end{align}
we expand the integrand in (C 2) for 
$\alpha \to \infty$ then extend the integration interval to 
$[0,\infty )$ to find
\begin{equation} Z\left(\alpha\right) \sim \frac{2\sqrt{\log\alpha}}{{\rm \pi}^2\sqrt{\alpha}}\left(\left(1+\frac{3+2\log2+\log \log\alpha}{2\log\alpha}\right)I_{0}+\frac{1}{2\log\alpha}I_{1}\right), \end{equation}where
\begin{align} I_{0}=\int_{0}^{\infty}\frac{\mathrm{e}^{-x} \, \textrm{d} x}{\sqrt{x}}=\sqrt{{\rm \pi}} \quad \text{and}\quad I_{1}=-\int_0^{\infty}\frac{\mathrm{e}^{-x}\log x \, \textrm{d} x}{\sqrt{x}}=\sqrt{{\rm \pi}}(\gamma_{E}+2\log 2) \end{align}
with 
$\gamma _{E} = 0.577\cdots$, the Euler–Mascheroni constant. Combining (C 4) with (C 5a,b) finally yields
\begin{equation} Z(\alpha) \sim \frac{2\sqrt{\log\alpha}}{{\rm \pi}^{3/2}\sqrt{\alpha}}\left(1+\frac{3+\gamma_{E}+\log(16\log\alpha)}{2\log\alpha}\right) \quad \text{as} \ \alpha \to \infty. \end{equation}
Note that the corrections up to 
$(\log \alpha )^{-1}$ are necessary for numerical applications as they decay very slowly with 
$\alpha$ (see figure 2).
C.2. Limit $\alpha \to -\infty$
In this case, the outer integral in (4.12) is dominated by a neighbourhood of the maximum of the integral multiplying 
$\alpha$, located at the right endpoint 
$\psi =1$. Expanding the inner integral near 
$\psi =1$ we find
\begin{equation} Z(\alpha) \sim \frac{2}{{\rm \pi}|\alpha|}\mathrm{e}^{\varUpsilon|\alpha|} \quad \text{as} \ \alpha \to - \infty, \end{equation}where
\begin{equation} \varUpsilon=\int_{0}^{1}\frac{b(\psi)}{a(\psi)} \, \textrm{d}\psi = 0.655\cdots, \end{equation}
using that 
$K'(1)={\rm \pi} /2$ and 
$b(\psi )/a(\psi )\to 1$ as 
$\psi \to 1$.
Appendix D. Another cellular flow
In this appendix, we examine a cellular flow different from the canonical cellular flow (1.1) in order to verify (5.4). We consider the flow with streamfunction
\begin{equation} \psi(x,y)=\tfrac{3}{4}\sin x\sin y+\tfrac{1}{4}\sin (3x)\sin (3y). \end{equation}A period of this flow, displaying the streamfunction and corresponding velocity field, is shown in figure 5(a).
$(a)$ Streamfunction (D 1). 
$(b)$ Function 
$b(\psi )$ for the streamfunction (D 1) computed from (4.6b). The streamlines associated with the roots 
$\psi _1\approx 0.42$ and 
$\psi _2\approx 0.52$ are represented with solid lines in panel 
$(a)$. 
$(c)$ Ratio of inertial over non-inertial particle effective diffusivities 
$\bar {D}$ as a function of 
$\alpha$ for 
$Pe=10^3$ (circles), 
$10^4$ (triangles) and 
$10^6$ (stars). The effective diffusivities are estimated by direct numerical sampling of the effective Brownian dynamics (2.11). Results are compared against a numerical evaluation of 
$1/Z(\alpha )$ from (4.9).
In figure 5(c), we compare direct numerical estimates of the effective diffusivity ratio 
$\bar {D}(Pe,\alpha )/\bar {D}(Pe,0)$ with the value of 
$1/Z(\alpha )$ computed from (4.9), where the integrals along streamlines 
$a(\psi )$, 
$b(\psi )$ and 
$c(\psi )$ are evaluated from their definitions (4.6a,b) and (4.10). The good match confirms the validity of (5.4).
We also note here that for (D 1) the function 
$b(\psi )$ becomes negative for a range of streamlines inside the quarter-cells. Within the quarter-cells with anticlockwise rotation, these streamlines are characterised by 
$0<\psi _1\leq \psi \leq \psi _2<1$ (see figure 5b). This is significant because 
$b(\psi )$ represents a (suitably averaged) cross-streamline drift induced by the weak inertia, with 
$\alpha b(\psi )>0$ corresponding to drift towards the separatrix, and 
$\alpha b(\psi )<0$ towards the centres of the quarter-cells. Therefore, in the absence of diffusion particles heavier than fluid (
$\alpha >0$) initially on streamlines 
$0<\psi <\psi _1$ are attracted by the separatrix 
$\psi =0$, while particles on streamlines 
$\psi _1 < \psi < 1$ are attracted by the streamline 
$\psi =\psi _2$. In the presence of diffusion, the overall effect remains an expulsion of the particles from the cells towards the separatrices leading to enhanced dispersion because 
$Z(\alpha )<1$, but this is only established on long time scales due to the presence of the attracting streamlines 
$\psi =\psi _2$.
The change of sign of 
$b(\psi )$ raises the intriguing possibility that for certain flows inner streamlines can be sufficiently attractive that 
$Z(\alpha )>1$ for some 
$\alpha >0$, so that heavy particles concentrate in the interior of the cells and dispersion is inhibited by inertia. We have neither been able to find such flows, nor ruled out their existence.
