Dispersion in rectangular networks: effective diffusivity and large-deviation rate function

The dispersion of a diffusive scalar in a fluid flowing through a network has many applications including to biological flows, porous media, water supply and urban pollution. Motivated by this, we develop a large-deviation theory that predicts the evolution of the concentration of a scalar released in a rectangular network in the limit of large time $t \gg 1$. This theory provides an approximation for the concentration that remains valid for large distances from the centre of mass, specifically for distances up to $O(t)$ and thus much beyond the $O(t^{1/2})$ range where a standard Gaussian approximation holds. A byproduct of the approach is a closed-form expression for the effective diffusivity tensor that governs this Gaussian approximation. Monte Carlo simulations of Brownian particles confirm the large-deviation results and demonstrate their effectiveness in describing the scalar distribution when $t$ is only moderately large.

In this Letter, we investigate the dispersion of a diffusive scalar released in a fluid flowing through a rectangular network (see Fig. 1). A vivid example of application -and a motivation for our work -is the spreading of a pollutant released suddenly in the streets of a city with a regular grid plan such as Manhattan. The primary question concerns the form taken by the scalar concentration C(x, t) long after release, when the disparity between the (large) scale of the scalar patch and the (small) scale of the network makes the problem challenging. The question arises in numerous applications across science and engineering besides urban pollution [1,2]: vascular and respiratory flows [3], microfluidic devices [4,5], porous media [6][7][8] and water distribution [9] for example. Its answer sheds light on the subtle interplay between advection, diffusion and geometry that controls dispersion in networks.
As is typical for advection-diffusion problems, C(x, t) for t 1 can be approximated by a Gaussian, parameterised by an effective diffusivity tensor [10,11]. This approximation applies only to the core of the scalar distribution, specifically to distances O(t 1/2 ) away from the centre of mass: the network geometry leads to non-Gaussian behaviour in the tails of the distribution. These tails are important in applications where low concentrations are critical, e.g., for highly toxic chemicals or in the presence of amplifying chemical reactions. To capture both the Gaussian core and the tails, we develop a large-deviation theory [12,13] that leads to a general approximation, of the form C(x, t) ∝ exp(−tg(x/t)) and holds for distances up to O(t) away from the centre of mass [14]. Here g is a rate function which we compute (by solving a transcendental equation) and approximate explicitly in asymptotic limits. Its quadratic approximation gives a closed-form expression for the effective diffusivity controlling the core of the scalar distribution. Monte Carlo simulations confirm the effective-diffusivity and large-deviation results and demonstrate the benefits of the latter, particularly for moderately large t when the non-Gaussian behaviour is most conspicuous.
Model.-We consider the rectangular network in Fig.  1 composed of one-dimensional edges of length L and βL in the xand y-directions along which fluid flows with uniform velocity U and V . This simple model has proved remarkably effective in describing pollution spreading through dense city centres [1,2]. More broadly, it provides an excellent prototype for geometry-induced non-Gaussianity and its description by large deviations.
Taking L as reference length and L 2 /κ as reference time, the one-dimensional advection-diffusion equations for the scalar concentration C read in edges oriented along x and y [15]. The non-dimensional parameters U and V are Péclet numbers measuring the strength of advection relative to diffusion. These equations are supplemented by boundary conditions applied at the vertices separated by distances 1 in x and β in y. The boundary conditions express (i) continuity of C, where the subscripts denote the limiting value to the west, east, etc. of the vertex, and (ii) vanishing of the net concentration flux which, on using (2), simplifies into Eqs. (19)-(3) form a closed system which can be solved numerically to predict the evolution of C for arbitrary initial conditions (e.g. using Laplace transforms [16][17][18]).
Here we consider a scalar initially released at a vertex taken to be the origin so that C(x, y, 0) = δ(x)δ(y). Large deviations.-Analytic progress is possible using the theory of large deviations [12,13]. This describes the concentration in the long-time limit t 1 as [14] The rate (or Cramér) function g(ξ) provides a continuous approximation for the most rapid changes in C and is the main object of interest. The function φ is supported on the network and has periods 1 and β in x and y. The factor t −1 is imposed by normalisation. Introducing (18) into ( To write these we have defined q = (q x , q y ) T = ∇ ξ g and f (q) = ξ · q − g(ξ), (7) which implies that f and g are Legendre transforms of one another, with q and ξ the dual independent variables. Eqs. (23)- (24) are supplemented by the boundary conditions inferred from (2) and (3): continuity of φ and Together, (23)-(22) form a family of eigenvalue problems parameterized by q, with f (q) as the eigenvalue. We solve Eqs. (23)- (24) to obtain explicit expressions for the eigenfunction φ in the 4 edges incident to the vertex (0, 0) using periodicity [19]. Introducing the solution into (22) gives where α U = f (q) + U 2 /4 and similarly for α V . This transcendental equation for f (q) is our central result. It can be solved numerically for a range of q to obtain f (q); the rate function g(ξ) is deduced by Legendre transform. We start our analysis by considering the behaviour of g(ξ) near its minimum. This provides a closed-form expression for the effective diffusivity of the network.
Effective diffusivity.-The Gaussian, diffusive approximation is deduced from (18) by Taylor expanding g(ξ) around its minimum ξ * , identified as the velocity of the centre of mass of the scalar. It can be shown [19] that ξ * = ∇ q f (0), and that the effective diffusivity tensor Introducing the Taylor expansion of f around q = 0 into (9) and solving gives, after lengthy manipulations, and the components . Note that effective diffusivities are more commonly derived using homogenization [20][21][22][23] or the method of moments [24,25]: solving their cell problem amounts to a perturbative solution of (23)-(22) [14]. The explicit expressions (11)- (14) illustrate the complex interplay between advection and diffusion that determines dispersion in networks. They are visualised for a range of U and V and two values of β as ellipses of constant x T K −1 x (corresponding to constant concentration) in Fig. 2. For U, V 1 (small Péclet number), the asymptotic formula h(x) = 1/12 + O(x 2 ) as x → 0 provides the approximation K 11 ∼ 1/(1 + β) + γU 2 , These grow linearly in U and V which dimensionally corresponds to components that are independent of the molecular diffusivity κ. This is characteristic of a regime termed geometric [26] or mechanical [8,27] dispersion. The tensor K is singular to leading order in U and V : effective diffusion is strong in the direction (−U, V ) but weak in the perpendicular direction (V, U ) (see Fig. 2) For U 1 and V 1, i.e., strong flow in the x-direction, This corresponds to a mostly longitudinal diffusivity with a κ −1 scaling characteristic of Taylor dispersion [10]. Note that the geometric and Taylor regimes can be understood in terms of a random-walk model with correlation time determined by advection in the first case and molecular diffusion in the second [17].
Rate function.
-Effective diffusivity provides a partial description of dispersion: the rate function g obtained from (25) is much more informative. This is demonstrated in Fig. 3 which shows typical examples of g obtained numerically for two values of (U, V ) (for β = 1) and its quadratic approximation corresponding to the Gaussian (30). This approximation is excellent in the vicinity of ξ * , with circular ( Fig. 3(a)) and elliptical contours ( Fig. 3(b)). Beyond the vicinity of ξ * , it is inadequate, failing for instance to capture the anisotropy of g and hence of C for U = V = 0, or underestimating g in large portions of the ξ-plane (hence overestimating C by an exponentially large factor) for U = V = 0. Note however that for U = V and β = 1, f is exactly quadratic along the line q x = q y (see (25)); hence, g coincides with its quadratic approximation for ξ x = ξ y , as evident in Fig. 3.
The limitations of the quadratic (Gaussian) approximation are best demonstrated by considering the large-ξ behaviour of g(ξ) or, equivalently, the large-q behaviour of f (q). In this regime, and with the distinguished scaling U, V = O(|q|), (25) reduces to Either term on the right-hand side is exponentially large, precluding the solution of (15) unless This gives a leading-order approximation to f which, remarkably, is independent of β. The Legendre transform of (16) is cumbersome for arbitrary U and V , but physical insight is gained by considering limiting cases. For U = V = 0, (16) leads to g(ξ) ∼ (|ξ x | + |ξ y |) 2 /4, Selected contours (with values 0.1, 1, 5 and 10) compare g (black) with its quadratic, Gaussian approximation (white) (30). This approximation is clearly valid near the minimum ξ * of g. in accordance with the diamond-shaped contours of g for large ξ in Fig. 3(a). This implies a concentration C ∼ exp(−(|x| + |y|) 2 /(4t)), which can be interpreted as a generalised form of diffusion with the Euclidian distance replaced by the L 1 (or Manhattan) distance. When U q x and V q y dominate in (16), the linear dependence of f on q implies that g → ∞ as ξ x → U and as ξ y → V , reflecting the finite propagation speed of the scalar when molecular diffusion is neglected against advection.  (25) becomes which implies a concentration independent of molecular diffusivity, generalising the notion of geometric or mechanical dispersion to the large-deviation regime. Monte Carlo simulations.-We now test our predictions against Monte Carlo simulations of Brownian particles. The concentration, derived as the probability density function (PDF) of their positions X(t), is compared with the large-deviation estimate in Fig. 4. The PDF is obtained from an ensemble of N = 10 6 particles by integrating the stochastic differential equations associated with (19), with the additional microscopic rule that particles entering a vertex exit through a random edge [28]. Although formally valid for t 1, the large-deviation approximation (18) is remarkably accurate for the moderate values of t = 1 and 5 considered. Its relevance is clear at t = 1, when comparing Figs. 3 and 4. As time progresses, the Gaussian approximation (30) becomes sufficient to describe the bulk of the scalar patch which assumes a characteristic elliptical form (cf. Fig. 3).
A detailed assessment of the large-deviation approximation requires a careful numerical evaluation of the rate function g. This is achieved by estimating its Legendre transform as the scaled cumulant generating function [29][30][31] f (q) = lim t→∞ t −1 log E e q·X(t) , where E is the expectation over the Brownian motion. To reduce sampling error to an acceptable level, we have adopted the pruning-cloning technique described in [14] based on [32]. Fig. 5 shows an excellent agreement between largedeviation predictions and numerical results (obtained for N = 10 3 particles at t = 5) and illustrates the restricted range of validity of the Gaussian approximation.
Conclusion.-We characterise the dispersive properties of a rectangular network by a rate function g deduced from (25). This describes the scalar concentration over a broad range of distances |x − ξ * t| = O(t) which proves particularly pertinent for moderately long times. In the narrower range |x − ξ * t| = O(t 1/2 ), it recovers the Gaussian, diffusive approximation and provides a convenient route to derive the effective diffusivity. Several conclusions can be drawn from the results: (i) in the absence of advection, the dispersion switches from a standard, L 2 diffusion with diffusivity κ/2 near the point of release, to an L 1 diffusion with diffusivity κ at large distances; (ii) correspondingly, the Gaussian approximation misrepresents the shape of the scalar patch and underestimates its area (by a factor π/4); (iii) advection leads to a complex, anisotropic behaviour, even in the Gaussian regime, with an enhancement of dispersion in the direction (−U, V ) corresponding to a constant advective travel time x/U + y/V ; (iv) strong advection (large Péclet number) leads to a geometric-dispersion regime, in which the rate function, and hence the effective diffusivity, are independent of the molecular diffusivity κ; (v) advection aligned with one of the axes of the network is anomalous in this respect, with an effective diffusivity that instead scales like κ −1 as in Taylor dispersion; (vi) the Gaussian approximation can under-and overpredict the scalar concentration for |x| = O(t), depending on x, by a factor that is exponentially large in t.
We emphasise that our large-deviation approach generalises straightforwardly to other periodic networks. It can capture anomalous diffusion [26] (when g is not quadratic near its minimum) and be further extended to fractal and random networks [16,17,33,34]. Our results are also directly applicable to reactive fronts: a Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) reaction, which adds Da C(1 − C) to (19) (with the Damköhler number Da as non-dimensional reaction rate), leads to the emergence of a concentration front. Its long-time speed of propagation v is determined by g through the condition g(v) = Da [13,35].

SUPPLEMENTAL NOTE
In this supplemental note, we provide details of the derivation of the eigenvalue problem determining the rate function g. We also deduce the corresponding effective diffusivity.

Eigenvalue problem
As discussed in the Letter, in the large-deviation regime, the concentration C of a dispersing scalar takes the asymptotic form Substituting (18) into the advection-diffusion equations and equating powers of t −1 yields, at leading order, Eqs.
which further simplifies into (Eq. (8) in the Letter) once we use continuity of φ and ∇ ξ g. We let q = ∇ ξ g(ξ) and f (q) = ξ · q − g(ξ) (Eq. (7) in the Letter). Treating q as a parameter, we obtain the eigenvalue problem with f (q) as the eigenvalue (Eqs. (5)-(6) in the Letter). The focus is on the principal eigenvalue (that with maximum real part) because it corresponds to the slowest decaying solution of (18). The Krein-Rutman theorem implies that this eigenvalue is unique, real and isolated, with a positive associated eigenfunction φ > 0. Moreover, f (q) ≥ 0 and is convex so that f (q) and g(ξ) are related by a Legendre transform We now solve (23)-(24) explicitly. Consider the intersection at (x, y) = (0, 0) and denote by φ E (x) the eigenfunction in the street to the east of it, and by A the value of φ at the intersection (0, 0) and hence, by periodicity, at all intersections. Solving (23) with the boundary conditions φ E (0) = φ E (1) = A, we find that where α U = f (q) + U 2 /4. The solution to the west of (0, 0) is found by substituting x → x + 1 in this expression to obtain φ W (x) = A sinh α U −e (qx+U/2)(x+1) sinh(α U x) + e (qx+U/2)x sinh(α U (x + 1)) .
We can now apply the boundary condition (22) by evaluating the derivatives of the solution. After some simplifications, this leads to which is Eq. (9) in the Letter.

Effective diffusivity
The rate function g has a single minimum, ξ * say, around which it can be expanded according to where ∇ ξ ∇ ξ g(ξ * ) is the Hessian of g (matrix of second derivatives) evaluated at ξ * . It follows from the Legendre transform that q = ∇ ξ g = 0 corresponds to the minimum of g, hence ξ * = ∇ q f (0). Taking the gradient of the relation q = ∇ ξ g with respect to ξ and evaluating at ξ * gives ∇ ξ ∇ ξ g(ξ * ) = ∇ ξ q(ξ * ).
On the other hand, the gradient with respect to ξ of ξ = ∇ q f and the chain rule give where I is the identify matrix. Evaluating at ξ = ξ * and, correspondingly, q = 0 leads to the standard relation between the Hessians of g and f , ∇ ξ ∇ ξ g(ξ * ) = (∇ q ∇ q f (0)) −1 .
Introducing this into (26) and using in (18) yields the Gaussian approximation C(x, t) ∼ t −1 e −(x−ξ * t) T K −1 (x−ξ * t)/(4t) , for the concentration, with the effective diffusivity tensor K = ∇ q ∇ q f (0)/2. This is Eq. (10) of the Letter. In order to deduce explicit expressions for ξ * = ∇ q f (0) and K from (25), we introduce the expansion into (25), expand in powers series of q and solve at order O(|q|) for ξ * and O(|q| 2 ) for K. This is best carried out using a symbolic-algebra package.