[

# [

[
###### Abstract

The dispersion of a passive scalar in a fluid through the combined action of advection and molecular diffusion is often described as a diffusive process, with an effective diffusivity that is enhanced compared to the molecular value. However, this description fails to capture the tails of the scalar concentration distribution in initial-value problems. To remedy this, we develop a large-deviation theory of scalar dispersion that provides an approximation to the scalar concentration valid at much larger distances away from the centre of mass, specifically distances that are rather than , where is the time from the scalar release.

The theory centres on the calculation of a rate function characterising the large-time form of the scalar concentration. This function is deduced from the solution of a one-parameter family of eigenvalue problems which we derive using two alternative approaches, one asymptotic, the other probabilistic. We emphasise the connection between the large-deviation theory and the homogenisation theory that is often used to compute effective diffusivities: a perturbative solution of the eigenvalue problems in the appropriate limit reduces at leading order to the cell problem of homogenisation theory.

We consider two classes of flows in some detail: shear flows and periodic flows with closed streamlines (cellular flows). In both cases, large deviation generalises classical results on effective diffusivity and captures new phenomena relevant to the tails of the scalar distribution. These include approximately finite dispersion speeds arising at large Péclet number (corresponding to small molecular diffusivity) and, for two-dimensional cellular flows, anisotropic dispersion. Explicit asymptotic results are obtained for shear flows in the limit of large . (A companion paper, Part II, is devoted to the large- asymptotic treatment of cellular flows.) The predictions of large-deviation theory are compared with Monte Carlo simulations that estimate the tails of concentration accurately using importance sampling.

Dispersion in the large-deviation regime I]Dispersion in the large-deviation regime. Part I: shear flows and periodic flows P. H. Haynes and J. Vanneste]P. H. Haynes and J. Vannestethanks: Email address for correspondence: J.Vanneste@ed.ac.uk

## 1 Introduction

Taylor (1953) identified the phenomenon of shear dispersion in which a passive scalar, e.g. a chemical pollutant, released in a pipe Poiseuille flow spreads along the pipe according to a diffusion law. The corresponding diffusivity, often termed effective diffusivity to distinguish it from molecular diffusivity, is inversely proportional to molecular diffusivity when the latter is small (see also Aris 1956; Young & Jones 1991). This effective diffusivity is associated with a random walk along the pipe that results from the random sampling of the Poiseuille flow by molecular Brownian motion across the pipe. The diffusive description of this random walk, and the corresponding Gaussian profile of the scalar concentration, of course only apply on time scales that are much longer than the Lagrangian correlation time scale.

Shear dispersion is a striking example of a broad class of phenomena in which the interaction between fluid motion and Brownian motion leads to a strong enhancement of dispersion and to effective diffusivities that are orders of magnitude larger than molecular diffusivity. The importance of these phenomena in applications, in particular industrial, biological and environmental applications, is obvious. This has motivated studies of effective diffusivity in many different flows (see Majda & Kramer 1999, for a review). These include spatially periodic flows which can be analysed using the method of homogenisation. This method, which exploits the separation between the (small) scale of the flow and the (large) scale of the scalar field that emerges in the long-time limit, has proved highly valuable: it applies to more complicated flows, including time-dependent and random flows, and provides a unifying framework for methods used earlier. Shear dispersion, in particular, can be regarded as a special case of homogenisation applied to periodic flows, where cells repeat in the along pipe direction and the flow in each cell is simple Poiseuille flow.

In the large literature on shear dispersion, efforts have been made to overcome the restriction to large times that underlies the diffusive approximation, and improved asymptotic estimates that capture some of the early-time behaviour have been obtained (see Young & Jones 1991 for a review and Camassa et al. 2010 for more recent results). For periodic flows, because the effective diffusivity is more difficult to compute, the focus has mainly remained on the derivation of asymptotic estimates and bounds, in particular in the limit of small molecular diffusivity (e.g. Majda & Kramer 1999; Novikov et al. 2005).

Here we consider a different aspect. The characterisation of dispersion in the long-time limit by an effective diffusivity and hence by a Gaussian scalar distribution holds only close to the centre of mass of the distribution: the results of homogenisation are in essence a manifestation of the central-limit theorem and apply only to particles displaced from the mean by distances. Our aim is to go beyond this and describe the concentration far from the mean. To achieve this, we derive large-deviation estimates for the concentration, that is, we derive the rate function in an approximation of the form for the scalar concentration at position and time .

Large-deviation theory extends the central-limit theorem and applies to numerous probabilistic problems (e.g. Dembo & Zeitouni 1998; den Hollander 2000). When applied to the stochastic differential equations governing the motion of fluid particles advected and diffused in a fluid flow, it naturally yields an improved approximation to the scalar concentration (interpreted as a particle-position probability function, cf. Jansons & Rogers 1995). This approximation is valid for distances from the mean that are rather than and therefore captures the tails of the distribution. These are typically non-Gaussian and not adequately represented by the diffusive approximation. This is illustrated in Figure 1 by the example of dispersion in a plane Couette flow, one of the shear flows considered in detail in this paper. The top panel shows the profile along the flow of the cross-stream averaged concentration at four successive times in the case of small molecular diffusivity. The figure compares the averaged concentration obtained numerically using a Monte Carlo simulation (symbols) with the Gaussian, diffusive approximation (dashed lines) and the large-deviation approximation derived in §§2–3 (solid lines). The units of and have been chosen so that the maximum flow velocity and (Taylor) effective diffusivity are both . The inadequacy of the diffusive approximation in describing the tails of the concentration and the superiority of the large-deviation approximation are apparent in the top panel for the earliest profile . They are obvious for all the profiles in the bottom panel which displays the results using logarithmic scale for . This emphasises the tails of to reveal how the diffusive prediction overestimates dispersion and to demonstrate the effectiveness of the large-deviation approximation. We note that while large deviation formally applies for , it appears here remarkably accurate for moderate . (The discrepancies between large-deviation and Monte Carlo results for are mainly attributable to the limitations of the straightforward Monte Carlo method used here and are much reduced with the more sophisticated methods discussed in §3.)

As the Couette-flow example illustrates, large-deviation theory provides estimates of the low scalar concentrations in the tails, where the diffusive approximation fails. This makes it relevant to a range of applications in which low concentrations matter. Examples include the prediction of the first time at which the concentration of a pollutant released in the environment exceeds a low safety threshold, and the quantification of the impact of stirring on chemical reactions in a fluid. In such examples, there is a strong sensitivity of the response (physiological or chemical) to low scalar concentrations that makes the logarithm of the concentration, and hence the rate function , highly relevant quantities. This broad observation can be made precise for the certain classes of chemical reactions. For F-KPP reactions (e.g. Xin 2009), the combination of diffusion and reaction leads to the formation of concentration fronts that propagate at a speed that turns out to be controlled by the large-deviation statistics of the dispersion and given explicitly in terms of the rate function (Gärtner & Freidlin 1979; see also Freidlin 1985, Ch. 7, Xin 2009, Ch. 2, and Tzella & Vanneste 2014a).

The present paper starts in §2 with a relatively general treatment of the large-deviation theory of dispersion which applies to time-independent periodic flows and to shear flows. The key result is a family of eigenvalue problems parameterised by a variable . The principal eigenvalue, , is the Legendre transform of the rate function . These eigenvalue problems can be thought of as generalised cell problems in that they resemble and extend the cell problem that appears when homegenization is used to compute effective diffusivities. In §§2.12.2 we present two alternative derivations of the the eigenvalue problems: the first is a direct asymptotic method that treats the large-deviation form of the concentration as an ansatz (see Kuske & Keller 1997); the second follows the standard probabilistic approach based on the Ellis–Gärtner theorem and considers the cumulant generating function of the particle position (e.g. Ellis 1995; Dembo & Zeitouni 1998; den Hollander 2000; Touchette 2009). We then discuss the relation between large deviation and homogenisation (§2.3). Homogenisation, and the corresponding diffusive approximation, are shown to be recovered when the eigenvalue problems yielding are solved perturbatively for small up to errors. Carrying out the perturbation expansion to higher orders provides a systematic way of improving on the diffusive approximation; in the case of shear dispersion, this recovers earlier results (Mercer & Roberts 1990; Young & Jones 1991).

The rest of the paper is devoted to dispersion in specific shear and periodic flows. We compute the functions and for the classical Couette and Poiseuille flows in §3 by solving the relevant one-dimensional eigenvalue problem numerically. We also obtain asymptotic results for the concentration at small and large distances from the centre of mass. While the first limit recovers the well-known expression for the effective diffusivity of shear flows, the second captures the finite propagation speed that exists when diffusion along the pipe is neglected. This provides a transparent example of the limitations of the diffusive approximation. Section 4 is devoted to a standard example of periodic flow, the two-dimensional cellular flow with streamfunction . The numerical solution of the corresponding eigenvalue problems for specific values of the Péclet number (measuring the relative strength of advection and diffusion) reveals interesting features of the dispersion, such as anisotropy, that are not captured in the diffusive approximation. Using a regular perturbation expansion, we derive explicit results in the limit of small . We examine the opposite, large-Péclet-number limit in a companion paper (Haynes & Vanneste 2014, hereafter Part II). We conclude the paper with a Discussion in §5.

Throughout the present paper and Part II, we verify the predictions of large-deviation theory against direct Monte Carlo simulations of particle dispersion. This is not without challenges since this requires estimating the tails of distributions which are associated with rare events and are, by definition, difficult to sample. We have therefore used importance sampling and implemented two methods that are applicable broadly. These are described in Appendix B. Two other Appendices are devoted to technical details of certain asymptotic limits.

## 2 Formulation

We start with the advection–diffusion equation for the concentration of a passive scalar. Using a characteristic spatial scale as reference length and the corresponding diffusive time scale , where is the molecular diffusivity, as a reference time, this equation can be written in the non-dimensional form

 ∂tC+Peu⋅∇C=∇2C, (2.0)

where is the Péclet number. Here is the typical magnitude of the velocity field, which is assumed to be time independent, , and divergence free, .

Equation (2) can be considered as the Fokker–Planck equation associated with the stochastic differential equation (SDE) which governs the position of fluid particles,

 dX=Peu(X)dt+√2dW, (2.0)

where denotes a Brownian motion. In this interpretation and with , the initial condition for the concentration is and the concentration at later times can then be thought of as the transition probability for a particle to move from at to at . We focus on this initial condition and use the notation when the dependence on needs to be made explicit.

In this paper we consider two somewhat different flow configurations. The first, relevant to Taylor dispersion, corresponds to parallel shear flows, with unidirectional and varying in the cross-flow direction only, and a domain that is bounded in this direction. The concentration then satisfies a no-flux condition at the boundary. The second configuration corresponds to a periodic in an unbounded domain. In both cases, our interest is in the dispersion in the unbounded directions of the domain. The shear-flow configuration can essentially be regarded as a particular case of the more general periodic-flow configuration, with the domain extending over only one period in the streamwise direction and no-flux boundary conditions replacing periodicity conditions. Because of this, we consider the two configurations together when developing the general large-deviation approach in the rest of this section. Any ambiguity that may arise as a result will be clarified in §3 and §4 when we apply the approach separately to shear flows and to two-dimensional periodic flows and obtain explicit results. Mixed configurations, in which the flow is periodic in certain directions and bounded in others, could also be treated with no essential changes.

### 2.1 Large-deviation approximation

We are interested in the form of for . Under the assumption that , the solution to (2) can be sought as the expansion

 C(x,t|x0)=t−d/2e−tg(ξ)(ϕ0(x,ξ)+t−1ϕ1(x,ξ)+⋯),where  ξ=(x−x0)/t, (2.0)

where is the number of spatial dimensions. This can be considered to be a WKB expansion with as large parameter. The leading-order approximation

 C(x,t|x0)∼t−d/2ϕ(x,ξ)e−tg(ξ), (2.0)

has the characteristic large-deviation form in which is the Cramér or rate function (e.g. Dembo & Zeitouni 1998; Touchette 2009, and references therein). The conservation of total mass – the spatial integral of – imposes that

 minξg(ξ)=0 (2.0)

and explains the presence of the prefactor in (2.1), as an application of Laplace’s method shows. Note that we concentrate on this leading-order approximation throughout and hence omit the subscript from .

Introducing the expansion (2.1) into (2) and retaining only the leading order terms gives

 (ξ⋅∇ξg−g)ϕ=∇2ϕ−(Peu+2∇ξg)⋅∇ϕ+(Peu⋅∇ξg+|∇ξg|2)ϕ. (2.0)

Letting

 q=∇ξgandf(q)=q⋅ξ−g, (2.0)

this equation reduces to

 ∇2ϕ−(Peu+2q)⋅∇ϕ+(Peu⋅q+|q|2)ϕ=f(q)ϕ, (2.0)

where can be regarded as a parameter. This can be rewritten compactly as

 eq⋅x(∇2−Peu⋅∇)(e−q⋅xϕ)=f(q)ϕ, (2.0)

in which the form of the operator on the left-hand side makes transparent the connection to the advection–diffusion operator . The function satisfies no-flux boundary conditions when impermeable boundaries are present or periodic boundary conditions in the case of unbounded domains with periodic .

Equation (2.1) is central to this paper. Together with its associated boundary conditions, it gives a family of eigenvalue problems for parameterised by , with as the eigenvalue. Solving these eigenvalue problems (numerically in general) provides as the principal eigenvalue, that is, the eigenvalue with largest real part. The rate function is then recovered by noting from (2.1) that and are related by a Legendre transform

 f(q)=supξ(q⋅ξ−g(ξ))andg(ξ)=supq(ξ⋅q−f(q)). (2.0)

The fact that the critical points of are suprema and the convexity of can be deduced from the probabilistic interpretation of discussed below.Note that the second equality assumes that is differentiable (e.g. Touchette 2009). It follows that

 ξ=∇qf, (2.0)

which gives a one-to-one map between the parameter and the physical variable . The eigenfunction of (2.1) associated with can therefore be equivalently thought of as a function of , as in (2.1), or of , as in (2.1). Note that the maximum principle can be used to show that is real and that is sign definite (e.g. Berestycki et al. 1994). This is consistent with the asymptotics (2.1) and the observation that the concentration is positive for all time if it is initially positive.

To summarise, solving the eigenvalue problem (2.1) for arbitrary and performing a Legendre transform of the principal eigenvalue yields the large- approximation (2.1) of the concentration. This approximation is valid for and thus, as discussed below, extends the standard diffusive approximation which requires . The eigenvalue problem (2.1) can be thought of as a generalised cell problem since, as discussed in § 2.3, it generalises the cell problem of homogenisation theory. Bensoussan et al. (1989, §4.3.1) derive this eigenvalue problem as part of a Floquet–Bloch theory for linear equations with periodic coefficients and term it ‘shifted cell problem’ (see also Papanicolaou 1995, §3.6, and §4 below).

### 2.2 Probabilistic derivation

An alternative view of the problem considers the moment generating function

 w(q,x,t)=Eeq⋅X,with  X(0)=x (2.0)

for the position of the fluid particles satisfying (2). Here denotes the expectation over the Brownian process in (2). The generating function obeys the backward Kolmogorov equation

 ∂tw=Peu⋅∇w+∇2w,with  w(q,x,0)=eq⋅x (2.0)

(e.g. Øksendal 1998; Gardiner 2004). A solution can be sought in the form

 w(q,x,t)=eq⋅x+f(q)tϕ†(q,x), (2.0)

where the function remains to be determined but will shortly be identified with that in (2.1).

Introducing (2.2) into (2.2) leads to

 ∇2ϕ†+(Peu+2q)⋅∇ϕ†+(Peu⋅q+|q|2)ϕ†=f(q)ϕ†, (2.0)

with no-flux or periodic boundary conditions. This corresponds to a family of eigenvalue problems, again parameterised by , which are the adjoints of those in (2.1), and hence have the same eigenvalues and in particular the same principal eigenvalue , justifying the notation in (2.2). This eigenvalue controls for . As a result, it can alternatively be defined by

 f(q)=limt→∞1tlogEeq⋅X(t) (2.0)

and interpreted as the limit as of the cumulant generating function scaled by . This function is convex by definition.

The relationship between the large- asymptotics of encoded in and that of can be made obvious. Noting from the definition (2.2) that is the Legendre transform with respect to of with the variable dual to , we apply Laplace’s method to obtain

 w(q,x,t)=∫eq⋅x′C(x′,t|x)dx′≍∫et(q⋅(ξ+x/t)−g(ξ))dξ≍eq⋅x+tsupξ(q⋅ξ−g(ξ)),

where denotes the asymptotic equivalence of the logarithms as and we use (2.1) to write .

From (2.2) we obtain the first part of (2.1). Under the assumption of differentiability of , which ensures that is convex, the second part follows, allowing the computation of the rate function. The argument used in this subsection, which relies on Laplace’s method to establish a connection between rate function and scaled cumulant generating function , is an instance of the Gärtner–Ellis theorem, a fundamental result of large-deviation theory which extends Cramér’s treatment of the sum of independent random numbers (see, e.g., Ellis 1995; Dembo & Zeitouni 1998; Touchette 2009). Rigorous results for a problem very similar to that defined above can be found in Freidlin (1985, Ch. 7). It may be worth contrasting the large-time () large deviations discussed in this paper, with the small-noise () large deviations developed by Freidlin & Wentzell (see Freidlin & Wentzell 2012): while for small noise a single (maximum-likelihood or instanton) trajectory controls the rate function , this is not generally the case for large time. As we discuss in the case of shear flows in §3, it is only for and sufficiently large that can be expressed in terms of single trajectory and that the two forms of large deviations intersect.

Some properties of and are useful to infer properties of the dispersion directly from without the need to carry out the Legendre transform explicity. As noted, and are convex. Therefore, from (2.1), increasing correspond to increasing , and can be thought of as a proxy for the more physical variable . It is clear from (2.2) that ; correspondingly,

 ∇qf(0)=ξ∗, (2.0)

defines which, by (2.1), minimizes . Eq. (2.1) then indicates that the maximum of and its centre of mass are located at . Qualitatively the Legendre transform implies that a slow growth of away from its minimum corresponds to a rapid growth of and vice versa. In particular, linear asymptotes for , say as in the one-dimensional case, correspond to vertical asymptotes for , as . This implies that vanishes for , reflecting a finite maximum transport speed for the scalar. Exactly linear asymptotes do not arise for because the eigenvalue problem (2.1) for has the simple solution which corresponds to a purely diffusive behaviour. However, for large , there can be a range of values of for which is approximately linear and a finite transport speed controls scalar dispersion.

### 2.3 Relation with homogenisation and its extensions

Much of the literature on scalar dispersion focuses on the computation of an effective diffusivity governing the dispersion for and . In this approximation, (2) reduces to the diffusion equation

 ∂tC+Pe⟨u⟩⋅∇C=∇⋅(k⋅∇C), (2.0)

where is the spatial average of , and is an effective diffusivity tensor. Alternatively, and can be obtained from the particle statistics using

 limt→∞1tEX=Pe⟨u⟩andlimt→∞12tE(X−Pe⟨u⟩t)⊗(X−Pe⟨u⟩t)=k. (2.0)

The form of has been derived for a variety of flows using several essentially equivalent methods, starting with Taylor’s (1953) work on shear flows. In the last 20 years, homogenisation, as reviewed in Majda & Kramer (1999) and Pavliotis & Stuart (2007), has become the systematic method of choice.

The diffusive approximation (2.0) can be recovered from the more general large deviation results: since the assumption implies that and hence that , we can expand according to

 f(q)=ξ∗⋅q+12q⋅Hf⋅q+O(|q|3), (2.0)

where is the Hessian of evaluated at . Taking the Legendre transform gives

 g(ξ)∼12(ξ−ξ∗)⋅H−1f⋅(ξ−ξ∗). (2.0)

In this approximation the concentration is

 C(x,t|x0)≍e−(x−ξ∗t)⋅H−1f⋅(x−ξ∗t)/(2t) (2.0)

corresponding to the solution of (2.0) with

 Pe⟨u⟩=ξ∗andk=Hf/2. (2.0)

This result also follows from (2.0) noting that the mean and covariances that appear on the left-hand sides are given by the first and second derivatives with respect to of the cumulant generating function evaluated .

Since the diffusive approximation is recovered from the large-deviation results by an expansion for small , it can be expected that the method of homogenisation is equivalent to the perturbative solution of the eigenvalue problem (2.1) or (2.2). This is plainly the case. Consider the periodic-flow configuration and assume that for simplicity. Expanding

 ϕ=1+|q|ϕ1+|q|2ϕ2+⋯andf=|q|α1+|q|2α2+⋯, (2.0)

and introducing this into (2.1) yields at ,

 ∇2ϕ1−Peu⋅∇ϕ1+Peu⋅^q=α1,

where is a unit vector. Averaging this equation gives that . The solution is then written as

 ϕ1=−^q⋅χ

in terms of the periodic, zero-average solution of the so-called cell problem

 ∇2χ−Peu⋅∇χ=Peu. (2.0)

(see Majda & Kramer 1999, §2.1). At order , the eigenvalue problem reduces to

 ∇2ϕ2−Peu⋅∇ϕ2−2^q⋅∇ϕ1+Pe(u⋅^q)ϕ1=α2.

Averaging gives

 α2=1+Pe⟨(u⋅^q)ϕ1⟩=1+^qi⟨∇χi⋅∇χj⟩^qj,

where the second equalities follows after some manipulations using (2.0) (see Majda & Kramer 1999, p. 251 for details). This corresponds to an effective diffusivity with components

 kij=12(Hf)ij=δij+⟨∇χi⋅∇χj⟩,

which is the standard homogenisation result. An analogous computation detailed in Appendix A shows how the homogenisation results for shear flows are recovered from the large-deviation calculation.

The perturbative solution of the eigenvalue problem (2.1) offers a route for the systematic improvement of the diffusive approximation. Such improvements, which have been derived for shear flows by Chatwin (1970, 1972), Mercer & Roberts (1990) and others (see Young & Jones 1991, for a review), extend the diffusion equation (2.0) to include higher-order spatial derivatives and increase the accuracy of the approximation for . They lead to effective equations of the form

 ∂tC+Pe⟨u⟩∇⋅C=kij∂ijC+k(3)ijk∂ijkC+k(4)ijkl∂ijlkC+⋯, (2.0)

where summation over repeated indices is understood and we have introduced higher-order effective tensors , etc. The behaviour of the large-deviation function as encodes all these tensors. This can be deduced from the large-deviation form (2.1) of the concentration which implies that and . Combining these formally leads to the effective equation

 ∂tC=f(−∇)C. (2.0)

Comparison with (2.0) shows that the various effective tensors that appear are given as derivatives of at . Hence they can be computed by continuing the perturbative solution of the eigenvalue problem (2.1) to higher orders in . This is demonstrated to for shear flows in Appendix A.

Another kind of improvement captures finite-time effects, specifically the fact that the mean and variance of the particle position have corrections to their linear growth which depend on initial conditions. These corrections have been computed for some shear flows (Aris 1956; Mercer & Roberts 1990; Young & Jones 1991) and termed ‘initial displacement’ and ‘variance deficit’. Although we do not consider them further in what follows, it can noted that Eq. (2.2) for the moment generating function is exact. Its solution for finite time can be expressed as a series of the form , where and denote the complete set of eigenvalues and eigenfunctions of (2.2). The constants can be determined from the initial condition of the concentration. It is clear, then, that the first 2 terms in the Taylor expansion of , where the mode corresponds to the eigenvalue , determine the initial displacement and variance deficit; the other eigenvalues contribute to exponentially small corrections.

In the rest of the paper, we apply the results of this section to several specific shear and periodic flows. We start with the case of shear flows for which the eigenvalue problems (2.1) and (2.2) simplify considerably.

## 3 Shear flows

Consider the advection by a parallel shear flow in two dimensions, in a channel of width corresponding to for the dimensionless coordinate . Without loss of generality (exploiting a suitable Galilean transformation as necessary) the velocity can be assumed to satisfy

 ⟨u⟩=12∫1−1u(y)dy=0. (3.0)

Because it is the longitudinal dispersion that is of interest, we modify (2.1) and take the large-deviation form of the concentration to be

 C(x,t)∼t−1/2ϕ(y,ξ)e−tg(ξ),% where  ξ=Pe−1x/t, (3.0)

assuming . Similarly, we write the moment generating function as

 (3.0)

Note that and depend only on the longitudinal variables and and that can be taken -independent because of the -independence of the flow. The factors are introduced in (3.0)–(3.0) for convenience: they lead to a Legendre pair of functions and that are independent of in the limit , at least for . The eigenvalue problem (2.1) then reduces to the Schrödinger form

 d2ϕdy2+(qu(y)+Pe−2q2)ϕ=f(q)ϕ. (3.0)

This one-dimensional eigenvalue problem is completed by the no-flux boundary conditions

 dϕdy(−1)=dϕdy(1)=0. (3.0)

Note that the operator in (3.0) is self adjoint and hence the same equation arises for the eigenvalue problem (2.2) for associated with the moment generating function. Note also that (3.0) can be derived more directly using the Feynman–Kac formula. To see this, write (2) explicitly as

 dX=Peu(Y)dt+√2dW1,dY=√2dW2, (3.0)

and note that . The generating function (3.0) then becomes

 w(q,x,t)=Eeq(Pe−1(x+√2W1)+∫t0u(y+√2W2)dt′)=ePe−1qx+Pe−2q2tEeq∫t0u(y+√2W2)dt′.

Using the Feynman–Kac formula (e.g. Øksendal 1998), is seen to satisfy

 ∂tw=∂yyw+(qu(y)+Pe−2q2)w

and hence, for , to depend on as with the principal eigenvalue in (3.0).

Alternatively, (3.0) is obtained when seeking normal-mode solutions of the form to the advection–diffusion equation (2) provided that the identification and is made. The large-deviation form of is then recovered by applying the steepest-descent method to the normal-mode expansion of . The large-deviation approach makes it clear that the saddle point in the plane is on the imaginary axis with a purely imaginary associated frequency .

Below we solve (3.0)–(3.0) numerically for some classical shear flows. Several general remarks can already be made. First, the term proportional to in (3.0) is associated with longitudinal (molecular) diffusion. For , it can be neglected for , leading to the simpler eigenvalue problem

 d2ϕdy2+qu(y)ϕ=f(q)ϕ (3.0)

which makes clear that and hence are independent of in the limit with . The large-deviation form of can be written in terms of dimensional variables and as

 C(x∗,t∗)≍e−a−2κt∗g(x∗/(Ut∗)), (3.0)

and its range of validity as and . In what follows, we mostly concentrate on the limit and solve (3.0) rather than (3.0): the effect of the neglected longitudinal diffusion on is straightforward, since it simply adds , but the corresponding change in is somewhat more complicated. It is nonetheless a simple matter to estimate the size of for which the neglect of longitudinal diffusivity ceases to be a good approximation.

Second, the perturbative solution of eigenvalue problem (3.0) for , provides an effective diffusivity as sketched in §2.3. In terms of , the dimensional effective diffusivity is expressed from (3.0) as

 k∗=a2U22κf′′(0), (3.0)

and is inversely proportional to the molecular diffusivity in the limit . The perturbative calculation carried out in Appendix A gives

 12f′′(0)=⟨(∫y−1u(y′)dy′)2⟩. (3.0)

and recovers the explicit form of as obtained using homogenisation (e.g. Majda & Kramer 1999; Camassa et al. 2010). The first of the corrections to the diffusive approximation of Mercer & Roberts (1990) and Young & Jones (1991) is also computed in Appendix A.

Third, the asymptotics of (3.0) indicates that tends to as , where denote the maximum and minimum velocities in the channel. This can be seen by noting that is the lowest eigenvalue of a Schrödinger operator which, in the semiclassical limit , is given by the minimum of the potential (e.g. Simon 1983). The implication, as discussed in §2.2, is that as . Physically, this corresponds to the fact that fluid particles have longitudinal velocities in the range ; changes in the concentration therefore propagate at finite speeds and the concentration is compactly supported for . This is only an approximation of course: when longitudinal molecular diffusion is taken into account, there is no limit on the propagation speed. It is readily seen that the term becomes comparable to in for and that the rate function is approximately the diffusive for near () or larger (smaller). This form of can also be shown to arise from an application of the Freidlin & Wentzell (2012) small-noise large-deviation theory and is controlled by a single maximum-likelihood trajectory. (This applies only when is sufficiently large: the dimensional expression (3.0) makes this clear, with an argument of the exponential that scales like whereas the small-noise large deviation necessarily leads to a scaling, corresponding to a factor with our non-dimensionalisation.)

Finally, we note that the eigenfunctions , where the dependence is inferred from the -dependence using , have a simple interpretation. For the amount of scalar at for can be approximated as

 ∫∞ξtC(x,y,t)dx≍ϕ(ξ,y)e−tg(ξ), (3.0)

since, by the convexity of , the integral is dominated by the contribution of the endpoint . Therefore gives the scalar distribution across the shear flow of particles with average speed greater than . Similarly, for , gives the distribution of particles with speed less than .

### 3.1 Couette flow

We now examine classical shear flows, starting with the plane Couette flow

 u(y)=y. (3.0)

The dispersion in this flow is illustrated in Figure 1. The figure shows how the diffusive and large-deviation approximations provide a good approximation in the core of the scalar distribution and how only large deviation captures the tails. Figure 1 does not resolve the tails of with sufficient detail to assess the validity of the large-deviation approximation fully, however. In what follows, we test systematically the large-deviation prediction for , defined as

 f(q)=limt→∞1tlogEePe−1qX(t) (3.0)

with our shear-flow scaling, by comparing the value obtained by solving the eigenvalue problem (3.0) for a range of with careful Monte Carlo estimates. The eigenvalue problem is solved using a finite-difference scheme. (An exact solution can be written in terms of Airy functions, but it is not particularly illuminating). The Monte Carlo estimates approximate the right-hand side of (3.0) as an average over a large number of solutions of (3.0). However, a straightforward implementation does not provide a reliable estimate for except for small values of . This is because for moderate to large is controlled by rare realisations which are not sampled satisfactorily. To remedy this, it is essential to use an importance-sampling technique which concentrates the computational effort on these realisations. For the results reported in this paper, we have implemented a version of Grassberger’s (1997) pruning-and-cloning technique which we describe in Appendix B.1.

Results for the plane Couette flow are displayed in the leftmost panels of Figure 2. The top panel shows the eigenvalue and Monte Carlo approximations of along with asymptotic approximations valid for small and large . The small- approximation for is found from (3.0) as

 f(q)∼215q2as  q→0. (3.0)

The large- approximation is obtained by noting that for , the solution to (3.0) is localised in boundary layers near . Concentrating on , we introduce and into (3.0). To leading order, this gives

 d2ϕdY2−Yϕ=μϕ, (3.0)

with solution decaying as . Imposing the boundary condition at gives the equation for . Hence we have

 f(q)∼|q|−1.019|q|2/3as  |q|→∞, (3.0)

using symmetry to deal with .

The top left panel of Figure 2 confirms the validity of the eigenvalue calculation and of the asymptotic estimates. In the case of the estimates, a constant is added to (3.0) to ensure a good match; with this correction, the asymptotic formula appears to be accurate for as small as , say. The dispersive approximation corresponding to the parabola (3.0) overestimates for all , indicating that this approximation overestimates the speed of dispersion or equivalently the magnitude of the tails of the distribution.

The rate function is shown in the second row of Figure 2. The solid curve is obtained by Legendre transforming the function computed by numerical solution of the eigenvalue problem. This is compared with direct Monte Carlo estimates. Again, it is crucial to use importance sampling to obtain a reliable estimate of for not small. We have chosen to integrate a modified dynamics in which particles, instead of simply diffusing in the -direction, also experience of drift towards the wall at (or ). A better sampling is obtained because the wall regions control for large ; the method is described in Appendix B.2. The Figure also shows the asymptotic approximations for deduced from (3.0) and (3.0) by Legendre transform and given by

 g(ξ)∼158ξ2  as  ξ→0andg(ξ)∼4⋅1.019327(1∓ξ)2  as  ξ→±1. (3.0)

The match between the values of derived from the eigenvalue problem and those obtained by Monte Carlo sampling provides a direct check on the validity of the large-deviation theory. The discrepancy between the exact and its diffusive approximation confirms that diffusion overestimates the dispersion speed, as inferred already from the plot of . The finite support of the concentration distribution for , arising from the neglect of longitudinal molecular diffusion, is also hinted at by the large slopes of for . The large- approximation to (with term fixed by inspection) is seen to be accurate for and could be combined with the small approximation to provide a satisfactory uniform approximation.

The third panel on the left of Figure 2 shows the map between that arises as part of the Legendre transform. This identifies the location which control the corresponding exponential moment for large . Finally, the fourth panel shows profiles of the eigenfunctions of (3.0) for several values of . According to (3.0), these give the structure of the concentration profile for larger than . Thus, for instance, the eigenfunction for approximately corresponds to (see third panel). As and hence increase (or decrease) the profile becomes more and more localised in the region of maximum (or minimum) velocity, that is, near (). The eigenfunctions for finite are to be contrasted with the standard (homogenisation) results on Taylor dispersion which correspond to eigenfunctions that are small, perturbations to the uniform eigenfunction .

### 3.2 Plane Poiseuille flow

We next examine the plane Poiseuille flow

 u(y)=1/3−y2. (3.0)

The small- approximation in this case is readily found from (3.0) to be

 f(q)∼8945q2as  q→0. (3.0)

For , the solution is localised around the maximum of the velocity at . For the required boundary-layer analysis, we let and and obtain

 d2ϕdY2−Y2ϕ=μϕ. (3.0)

The solution corresponding to the largest eigenvalue is the Gaussian , leading to and

 f(q)∼q/3−q1/2as  q→∞. (3.0)

For , the asymptotic treatment is similar to that of the Couette flow: we let and and find that and hence . This gives the approximation

 f(q)∼−2q/3−1.617q2/3as  q→−∞. (3.0)

The corresponding rate function is derived by Legendre transform, yielding the asymptotic behaviours

 g(ξ)∼94532ξ2  as  ξ→0, (3.0)
 g(ξ)∼14(1/3−ξ)  as  ξ→1/3,  and  g(ξ)∼4⋅1.617327(2/3+ξ)2  as  ξ→−2/3. (3.0)

The numerical and asymptotic results obtained for the plane Poiseuille flow are displayed in the second column of Figure 2. As for the Couette flow, the diffusive approximation (3.0) and (3.0) is seen to overestimate the speed of dispersion, leading to an overestimate of and an underestimate of . The concentration distribution for the Poiseuille flow is skewed, with increasing faster for than corresponding to smaller concentrations for than for . The eigenfunctions shown in the bottom panel illustrate how for large (small ) and hence for large (small ) are controlled by motion near the centre (periphery) of the flow. This culminates in the limits () as the boundary-layer form of the eigenfunctions derived above indicates.

### 3.3 Pipe Poiseuille flow

We conclude this section by considering the Poiseuille flow in a pipe, with velocity

 u(r)=1/2−r2, (3.0)

where . This flow is three-dimensional, with particles diffusing across the flow in both the - and -directions. While the eigenfunctions for axisymmetric flows can in principle depend on and independently, the principal eigenvalue determining is obtained for axisymmetric : . Correspondingly, the eigenvalue problem (3.0) of plane shear flows is replaced by

 1rddr(rdϕdr)+qu(r)ϕ=f(q)ϕ (3.0)

with boundary conditions at .

The small-, diffusive approximation for general axisymmetric shear flows is quoted in Appendix A as (A 0). For the Poiseuille flow, this gives

 f(q)∼1192q2as  q→0. (3.0)

For , an approximation to is derived from (3.0) using a boundary-layer approach: we let and to find the leading-order equation

 1RddR(RdϕdR)−R2ϕ=μϕ, (3.0)

with solution , corresponding to . Therefore,

 f(q)∼q/2−2q1/2as  q→∞. (3.0)

The analysis for is almost identical to that carried out for the plane Poiseuille flow and leads to

 f(q)∼−q/2−1.617q2/3as  q→−∞. (3.0)

Computing the Legendre transform of (3.0), (3.0) and (3.0) yields the corresponding asymptotics results for the rate function, namely

 g(ξ)∼48ξ2  as  ξ→0, (3.0)
 g(ξ)∼1(1/2−ξ)  as  ξ→1/2,  and  g(ξ)∼4⋅1.617327(1/2+ξ)2  as  ξ→−1/2. (3.0)

Note that (3.0) recover’s Taylor’s original result (Taylor 1953).

The numerical and asymptotic results for the pipe Poiseuille flow are shown in the rightmost panels of Figure 2. The diffusive approximation is seen to mostly overestimate the dispersion speed, although it turns out to be remarkably accurate for . Close inspection shows in fact that there is a range of values of for which diffusion underestimates somewhat the concentration, in contrast to the other cases considered. Note that the skewness for the pipe Poiseuille flow is opposite to that of the plane Poiseuille flow, with larger concentrations predicted for than .

## 4 Periodic flows

We now turn to two-dimensional periodic flows. The formalism of § 2 applies directly: is obtained by solving the eigenvalue problem (2.1) with periodic boundary conditions for . Eq. (2.1) can also be obtained in an alternative manner: because the advection–diffusion equation (2) has periodic coefficients, its solutions can be sought in the Floquet–Bloch form , which leads to (2.1) with and (Bensoussan et al. 1989; Papanicolaou 1995). This approach gives a representation of the concentration as an integral over whose large- asymptotics, derived using the steepest-descent method, reduces to the large-deviation form (2.1).

We focus our attention on the cellular flow

 u(x,y)=(−∂yψ,∂xψ)withψ=−sinxsiny. (4.0)

This flow, with period in both the - and -direction, consists of a regular array of cells in which the fluid is rotating alternatively clockwise and counterclockwise along closed streamlines; see Figure 3. It has received a great deal of attention, most of it devoted to the properties of the effective diffusivity that can be computed by homogenisation, especially in the limit of large Péclet number; see Majda & Kramer (1999, §2) for a review, and Novikov et al. (2005) and Gorb et al. (2011) for more recent references.

To illustrate the dispersion of a passive scalar in this flow, we show in Figures 45 the concentration field obtained by solving numerically the advection–diffusion equation (2) for and . Only the first quadrant is shown since the field has a four-fold symmetry. For , molecular diffusion plays a major part across the domain, leading to a smooth evolution, with only some modulations in the form of diagonal bands in the central sector of the quadrant and of cells located near the coordinate axes. For , advection dominates, resulting in an apparent finite propagation speed and the obvious mark of the flow structure on the scalar field. The importance of the separatrices, around which boundary layers of high concentrations are established, is clear. As the distance from the origin increases, there is gradual change in the scalar distribution within the cells, from almost uniform near the origin to essentially zero at large distance. This feature is discussed briefly below and fully elucidated in Part II.

Let us now turn to the predictions of large-deviation theory. We have developed a code for the numerical solution of the eigenvalue problem (2.1) for (4.0). This relies on a straightforward finite-difference discretisation and on the matlab routine ‘eigs’ for the solution of the resulting matrix eigenvalue problem. The convergence of the algorithm requires a good first guess for the eigenvalue; since we are interested in obtaining for a range of , the code performs an iteration over and , using at each step the previous value of as its first guess. Since satisfies the obvious symmetries , we concentrate on the first quadrant of the -plane. The symmetry can also be exploited.

The left panel of Figure 6 shows the numerical approximation to obtained using this code for . It is compared with the result of a Monte Carlo estimate which relies on the importance-sampling algorithm described in Appendix B.1. In addition to confirming the validity of the large-deviation approximation and of the numerical implementation, the figure illustrates general qualitative features of . For small , is approximately isotropic, consistent with the result of homogenisation theory which predicts a diagonal effective diffusivity tensor. For of order-one or larger, however, is anisotropic, taking smaller values along the axes and than along the diagonal . Physically, this implies that dispersion is slower along the axis than along the diagonal. This can be understood by considering the streamline geometry: continued advection along one of the axes requires particles to also meander in the perpendicular direction, resulting in a decrease in average speed by a factor ; by contrast, advection along the diagonal happens in staircase-like fashion which decreases the speed by a factor . That motion along the diagonal is faster is also apparent in the rate function obtained by Legendre transform and shown on the right panel of Figure 6: when is not small, the contours of , which directly correspond to concentration contours, are anisotropic with the larger scalar concentrations along the diagonal.

A direct Monte Carlo estimate of — as opposed to the indirect estimate deduced from Legendre transforming the Monte Carlo approximation to — proves difficult to compute reliably. Figure 7 illustrates this: even for a large number of realisations of , the direct Monte Carlo approach only provides a valid approximation for , in range where remains roughly isotropic. Attempts at implementing importance sampling in a manner analogous to that used for shear flows and described in Appendix B.2 did not lead to significant improvements in the estimation of in this direct manner. A conclusion, therefore, is that a more efficient Monte Carlo approximation to is achieved by sampling and taking a Legendre transform. Of course, for this problem the most efficient method for obtaining and remains the numerical solution of the eigenvalue problem (2.1).

It is interesting to examine the eigenfunctions associated with the eigenvalue for given since these provide the structure of the scalar field at position (with convex so that can be interpreted as a proxy for ). Figure 8 shows the eigenfunctions obtained by numerical solution of the eigenvalue problem for three values of . For small and hence small , is essentially constant over the whole cell, with only small modulations. This is consistent with the perturbative treatment of the eigenvalue problems for and , amounting to homogenisation, which indicates that . As and increase, the modulations, in the form of diagonal stripes, increase in amplitude so that, for large , is close to zero in wide stripes. The form of the eigenfunctions depends on the angle of , of course, and for or for instance, corresponding to dispersion along the and axes, they have a have a cellular rather than banded structure (not shown). The structure of the eigenfunctions is consistent with the concentration field shown in Figure 4. To see this, recall that the concentration depends on both and on the rate function ; across a single cell, the latter varies slowly and can be approximated by a Taylor expansion, leading to the spatial dependence , since . For large , the dominant effect is the exponential decay of the concentration in the direction of , with the form of introducing the banded modulations observed in Figure 4.

Some insight into the large-deviation behaviour of cellular flows can be gained by considering the regime corresponding to weak advection. The effective diffusivity in this limit was computed by Moffatt (1983, §7) and Sagues & Horsthemke (1986) who obtained (in our notation) the approximation . The generalisation to the large-deviation regime is straightforward and described in Appendix C. It leads to the asymptotic approximation

 f(q)=q21+q22+Pe28q21+q22+q41+6q21q22+q421+2(q21+q22)+(q21−q22)2+O(Pe3) (4.0)

whose small- limit is consistent with the effective diffusivity just quoted. This approximation is tested against numerical results in Figure 9 which shows the correction to purely diffusive behaviour for . The figure confirms the validity of (4.0); it also shows that dispersion is fastest along the diagonal, as noted for . The correction to behaves in fact very differently for than it does for : whereas is i