{'fraction_non_alphanumeric': 0.04501420067696378, 'fraction_numerical': 0.014725907481616932, 'mean_word_length': 4.151002004008016, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalar Decay in Chaotic Mixing

I review the local theory of mixing, which focuses on infinitesimal blobs of scalar being advected and stretched by a random velocity field. An advantage of this theory is that it provides elegant analytical results. A disadvantage is that it is highly idealised. Nevertheless, it provides insight into the mechanism of chaotic mixing and the effect of random fluctuations on the rate of decay of the concentration field of a passive scalar.


Introduction
The equation that is in the spotlight is the advection-diffusion equation for the time-evolution of a distribution of concentration θ(x, t), being advected by a velocity field v(x, t), and diffused with diffusivity κ. The concentration θ is called a scalar (as opposed to a vector). We will restrict our attention to incompressible velocity fields, for which ∇ · v = 0. For our purposes, we shall leave the exact nature of θ nebulous: it could be a temperature, the concentration of salt, dye, chemicals, isotopes, or even plankton. The only assumption for now is that this scalar is passive, which means that its value does not affect the velocity field v. Clearly, this is not strictly true of some scalars like temperature, because a varying buoyancy influences the flow, but is often a good approximation nonetheless. The advection-diffusion equation is linear, but contrary to popular belief that does not mean it is simple! Because the velocity (which is regarded here as a given vector field) is a function of space and time, the advection term (the second term in (1)) can cause complicated behaviour in θ. Broadly speaking, the advection term tends to create sharp gradients of θ, whilst the diffusion term (the term on the right-hand side of (1)) tends to wipe out gradients. The evolution of the concentration field is thus given by a delicate balance of advection and diffusion.
The advection term in (1) is also known as the stirring term, and the interplay of advection and diffusion is often called stirring and mixing. As we shall see, the two terms have very different rôles, but both are needed to achieve an efficient mixing.
To elicit some broad features of mixing, we will start by deriving some properties of the advection-diffusion equation. First, it conserves the total quantity of θ. If we use angle brackets to denote the average of θ over the fixed domain of interest V , i.e.
then we find diretly from (1) that Because the velocity field is incompressible, we have v · ∇θ = ∇ · (θ v), and also ∇ 2 θ = ∇ · (∇θ). Thus, we can use the divergence theorem to write (2) as where S is the surface bounding V , and dS is the element of area, andn outward-pointing normal to the surface. For a closed flow, two possibilities are now open to us: (i) the domain V is periodic; or (ii) v and ∇θ are both tangent to the surface S. In the first case, the terms on the right-hand side of (3) vanish because boundary terms always vanish with periodic boundary conditions (a bit tautological, but true!). In the second case, both v·n and ∇θ· n vanish. Either way, ∂ t θ = 0 (4) so that the mean value of θ is constant. Since V is constant, this also implies that the total amount of θ is conserved. The second set of boundary conditions we used implies that there is no fluid flow or flux of θ through the boundary of the volume. It is thus natural that the total θ is conserved! For periodic boundary conditions, whatever leaves the volume re-enters on the other side, so it also makes sense that θ is conserved. Because of (4), and because we can always add a constant to θ without changing its evolution (only derivatives of θ appear in (1)), we will always choose θ = 0 (5) without loss of generality. In words: the mean of our scalar vanishes initially, so by (4) it must vanish for all times. Now let's look at another average of θ: rather than averaging θ itself, which has yielded an important but boring result, we average its square. The variance is defined by where the second term on the right vanishes by (5). To obtain an equation for the time-evolution of the variance, we multiply (1) by θ and integrate, We rearrange on the left and integrate by parts on the right, to find Now there are some boundary terms that vanish under the same assumptions as before, and we get ∂ t Var = −2κ |∇θ| 2 .
Notice that, once again, the velocity field has dropped out of this averaged equation. However, now the effect of diffusion remains. Moreover, it is clear that the term on the right-hand side of (7) is negative-definite (or zero): this means that the variance always decreases (or is constant). The only way it can stop decreasing is if ∇θ vanishes everywhere, that is, θ is constant in space. But because we have assumed θ = 0, this means that θ = 0 everywhere. In that case, we have no choice but to declare the system to be perfectly mixed : there are no variations in θ at all anymore. Equation (7) tells us that variance tends to zero, which means that the system inexorably tends to the perfectly mixed state, without necessarily ever reaching it. Variance is thus a useful measure of mixing: the smaller the variance, the better the mixing.
There is a problem with all this: equation (7) no longer involves the velocity field. But if variance is to give us a measure of mixing, shouldn't its time-evolution involve the velocity field? Is this telling us that stirring has no effect on mixing? Of course not, as any coffee-drinker will testify, whether she likes it with milk or sugar: stirring has a huge impact on mixing! So what's the catch?
The catch is that (7) is not a closed equation for the variance: the righthand side involves |∇θ| 2 , which is not the same as θ 2 . The extra gradient makes all the difference. As we will see, under the right circumstances the stirring velocity field creates very large gradients in the concentration field, which makes variance decrease much faster than it would if diffusivity were acting alone. In fact, when κ is very small, in the best stirring flows the gradients of θ scale as κ −1/2 , so that the right-hand side of (7) becomes independent of the diffusivity. This, in a nutshell, is the essence of enhanced mixing.
Several important questions can now be raised: -How fast is the approach to the perfectly-mixed state? -How does this depend on κ? -What does the concentration field look like for long times? What is its spectrum? -How does the probability distribution of θ evolve? -Which stirring fields give efficient mixing?
The answers to these questions are quite complicated, and not fully known. In the following sections we will attempt to give some hints of the answers and give some references to the literature. This is not meant to be a comprehensive review article, so entire swaths of the literature are missing. We focus mainly on local or Lagrangian theories, which involve deterministic and stochastic approaches for quantifying stretching using a local idealisation of the flow. The essential feature here is that the advection-diffusion equation is solved along fluid trajectories. These theories trace their origins to Batchelor [1], who treated constant matrices with slow time dependence, and Kraichnan, who introduced fast (delta-correlated) time dependence [2,3]. Zeldovich et al. [4] approached the problem from the random-matrix theory angle in the magnetic dynamo context. More recently, techniques from large-deviation theory [5][6][7][8] and path integration [9][10][11][12][13] have allowed an essentially complete solution of the problem. It is this work that will be reviewed here, as it applies to the decay of the passive scalar (and not the PDF of concentration or its power spectrum). We will favour expediency over mathematical rigour, and try to give a flavour of what these local theories are about without describing them in detail.
The story will proceed from here as follows: in Section 2 advection of a blob by a linear velocity field is considered, with diffusion included. This problem has an exact solution, but it can be made simpler in the limit of small diffusivity. Solutions are examined for a straining flow in two and three dimensions (Sections 2.2 and 2.4), as well as a shear flow in two dimensions (Section 2.3). Randomness is added in Section 3: the strain associated with the velocity field is assumed to vary, and the consequences of this for a single blob (Section 3.1) and a large number of blobs (Section 3.2) are explored. Practical implementation is discussed in Section 4, and a simple model for a micromixer is analysed in Section 4.1. Finally, the limitations of the theory presented herein (Section 4.2).

Advection and Diffusion in a Linear Velocity Field
We will start by considering what happens to a passive scalar advected by a linear velocity field. The overriding advantage of this configuration is that it can be solved analytically, but that is not its only pleasant feature. Like most good toy models, it serves as a nice prototype for what happens in more complicated flows. It also serves as a building block for what may be called the local theory of mixing (Section 3).
The perfect setting to consider a linear flow is in the limit of large Schmidt number. The Schmidt number is a dimensionless quantity defined as where ν is the kinematic viscosity of the fluid and κ is the diffusivity of the scalar. The Schmidt number may be thought of as the ratio of the diffusion time for the scalar to that for momentum in the fluid. Alternatively, it can be regarded as the ratio of the (squared) length of the smallest feature in the velocity field to that in the scalar field. This last interpretation is due to the fact that if θ varies in space more quickly than √ κ, then its gradient is large and diffusion wipes out the variation. The same applies to variations in the velocity field with respect to ν. Hence, for large Schmidt number the scalar field has much faster variations than the velocity field. This means that it is possible to focus on a region of the domain large enough for the scalar concentration to vary appreciably, but small enough that the velocity field appears linear. Because there are many cases for which Sc is quite large, this motivates the use of a linear velocity field. In fact, large Sc number is the natural setting for chaotic advection. It is also the regime that was studied by Batchelor and leads to the celebrated Batchelor spectrum [1]. The limit of small Sc is the domain of homogenization theory and of turbulent diffusivity models. We shall not discuss such things here.

Solution of the Problem
We choose a linear velocity field of the form where σ is a traceless matrix because ∇·v must vanish. Inserting this into (1), we want to solve the initial value problem Here the coordinate x is really a deviation from a reference fluid trajectory.
(In Appendix A we derive (8) from (1) by transforming to a comoving frame and assuming the velocity field is smooth.) We will follow closely the solution of Zeldovich et al. [4], who solved this by the method of "partial solutions." Consider a solution of the form where k 0 is some initial wavevector. We will see if we can make this into a solution by a judicious choice ofθ(k 0 , t) and k(t). The time derivative of (9) is Putting these together into (8) and cancelling out the exponential gives This must hold for all x, and neitherθ nor k depend on x, so we equate powers of x. This gives the two evolution equations We can write the solution to (10a) in terms of the fundamental solution T(t, 0) as and Id is the identity matrix. The advantage of doing this is that we can use the same fundamental solution for all initial conditions. We will usually write T t to mean T(t, 0). Note that because Tr σ = 0, we have This is a standard result that is proved in Appendix B for completeness. If σ is not a function of time, then the fundamental solution is simply a matrix exponential, but in general the form of T t is more complicated. Now that we know the time-dependence of k, we can express the solution to (10) as We can think of T t as transforming a Lagrangian wavevector k 0 to its Eulerian counterpart k. Thus (13b) expresses the fact thatθ decays diffusively at a rate determined by the cumulative norm of the wavenumber k experienced during its evolution. The full solution to (8) is now given by superposition of the partial solutions, whereθ 0 (k 0 ) is the Fourier transform of the initial condition θ 0 (x). 1 Assuming the spatial mean of θ vanishes, the variance (6) is which from (13b) becomes We thus have a full solution of the advection-diffusion equation for the case of a linear velocity field and found the time-evolution of the variance. But what can be gleaned from it? We shall look at some special cases in the following section.

Straining Flow in 2D
We now take an even more idealised approach: consider the case where the velocity gradient matrix σ is constant. Furthermore, let us restrict ourselves to two-dimensional flows. After a coordinate change, the traceless matrix σ can only take two possible forms, Case (2a) is a purely straining flow that stretches exponentially in one direction, and contracts in the other. Case (2b) is a linear shear flow in the x 1 direction. We assume without loss of generality that λ > 0 and U ′ > 0. The form σ (2b) is known as the Jordan canonical form, and can only occur for degenerate eigenvalues. Since by incompressibility the sum of these identical eigenvalues must vanish, they must both vanish. The corresponding fundamental matrices T t = exp(−σ t) are These are easy to compute: in the first instance one merely exponentiates the diagonal elements, in the second the exponential power series terminates after two terms, because the square of σ (2b) is zero. Let us consider Case (2a), a flow with constant stretching (the case considered by Batchelor [1]). The action of the fundamental matrix on k 0 for Case (2a) isθ for the Fourier transform in d dimensions.
The wavevector k(t) = T (2a) t · k 0 grows exponentially in time, which means that the length scale is becoming very small. This only occurs in the direction x 2 , which is sensible because that direction corresponds to a contracting flow. Picture a curtain being closed: the bunching up of the fabric into tight folds is analogous to the contraction. (Of course, it is difficult to close a curtain exponentially quickly forever!) The component of the wavevector in the x 1 direction decreases in magnitude, which corresponds to the opening of a curtain.
Let's see what happens to one Fourier mode. By inserting (19) in (13b), we haveθ The time-integral can be done explicitly, and we find For moderately long times (t λ −1 ), we can surely neglect e −2λt compared to 1, and 1 compared to e 2λt , Actually, this assumption of moderately long time is easily justified physically. If κk 2 /λ ≪ 1, where k is the largest initial wavenumber (that is, the smallest initial scale), then the argument of the exponential in (20) is small, unless e 2λt Pe (21) where the Péclet number is Pe = λ κ k 2 .
Thus the assumption that e 2λt is large is a consequence of Pe being large, since otherwise the exponential in (20) is near unity and can be ignoredvariance is approximately constant. We can turn (21) into a requirement on the time, It is clear from (23) that λ −1 sets the time scale for the argument of the exponential in (20) to become important. The Péclet number influences this time scale only weakly (logarithmically). This is probably the most important physical fact about chaotic mixing: Small diffusivity has only a logarithmic effect. Thus vigorous stirring always has a chance to overcome a small diffusivity, no matter how small: we need just stir a bit longer. Note that the variance is given by which is approximately constant for t ≪ λ −1 log Pe 1/2 . This does not mean that the concentration field · k 0 is a function of time from (18). This time dependence becomes important for t λ −1 .
The Péclet number may be thought of as the ratio of the advection time of the flow to the diffusion time for the scalar. It is usually written as where U is a typical velocity and L a typical length scale. Our velocity estimate in (22) is λ/k, and our length scale is k, which are both natural for the problem at hand. Just like large Sc, large Pe is the natural setting for chaotic advection. In fact, if Pe is small then diffusion is faster than advection, and stirring is not really required! Large Pe means that diffusion by itself is not very effective, so that stirring is required. We shall always assume that Pe is large.
We return to (20): the striking thing about that equation is its prediction for the rate of decay of the concentration field. Roughly speaking, (20) predicts for λt ≫ 1. Eq. (26) is the exponential of an exponential-a superexponential decay. This is extremely fast decay. In fact, unnaturally so: it is hard to imagine a physically sensible system that could mix this quickly. Something more has to be at work here. If we examine (20) closely, we see that the culprit is the term which grows exponentially fast. This term has its origin in the Laplacian in the advection-diffusion equation (1): the contracting direction of the flow (the x 2 direction) leads to an exponential increase in the wavenumber via the curtain-closing mechanism. This is exactly the mechanism for enhanced mixing we advertised on p. 3: very large gradients of concentration are being created, exponentially fast. This mechanism is just acting too quickly for our taste! So what's the problem? We are doing the wrong thing to obtain our estimate (26). This estimate tells us how fast a typical wavevector decays, and it says that this occurs very quickly. What we really want to know is what modes survive superexponential decay the longest, and at what rate they decay. Clearly the concentration in most wavenumbers gets annihilated almost instantly, once the condition (23) is satisfied. But a small number remains: those are the modes with wavevector closely aligned to the x 1 (stretching) direction, or equivalently that have a very small projection on the x 2 (contracting) direction. To overcome the exponential growth in (27), we require that is at any given time we need consider only wavenumbers satisfying (28), since the concentration in all the others has long since been wiped out by diffusion. The consequence is that the k 02 integral in (24) is dominated by these surviving modes. To see this, we blow up the k 02 integration by making the coordinate change k 02 = k 02 e λt in (24), The decay factor e −λt has appeared in front. For small diffusivity, we can neglect the k 0 2 1 term in the exponential (it just smooths out the initial concentration field a little). 2 We can then take the inverse Fourier transform ofθ 0 (k 01 , k 02 e −λt ), and insert this into (29). We then interchange the order of integration: the k 01 integral gives a δ-function, and the k 02 integral gives a Gaussian. The final result is where is a normalised Gaussian distribution with standard deviation ℓ, and we defined the length If the initial concentration decays for large |x 2 | (as when we have a single blob of dye), then (30) can be simplified to So the x 1 dependence in (32) is given by the "stretched" initial distribution, averaged over x 2 . The important thing to notice is that This is a much more reasonable estimate for the decay of concentration than (26)! The concentration thus decays exponentially at a rate given by the rate-of-strain (or stretching rate) in our flow. The exponential decay is entirely due to the narrowing of the domain for eligible (i.e., nondecayed) modes. This "domain of eligibility" is also known as the cone or the cone of safety [4,14]. (In two dimensions it is more properly called a wedge.) The concentration associated with wavevectors that fit within this cone is temporarily shielded from being diffusively wiped out, but as the aperture of the cone is shrinking exponentially more and more modes leave the safety of the cone as time progresses. Notice that (33) is independent of κ. This brings us to the second most important physical fact about chaotic mixing: The asymptotic decay rate of the concentration field tends to be independent of diffusivity. But note that a nonzero diffusivity is crucial in forcing the alignment (28). The only effect of the diffusivity is to lengthen the wait before exponential decay sets in, as given by the estimate (23). But this effect is only logarithmic in the diffusivity.
We can also try to think of (32) in real rather than Fourier space. Consider an initial distribution of concentration. Our straining flow will stretch this distribution in the x 1 direction, and contract it in the x 2 direction. Gradients in x 2 will thus become very large, so that eventually diffusion will limit further contraction in the x 2 direction and the distribution will stabilise with width κ/λ (see Fig. 1). This is what the Gaussian prefactor in (32) is telling us: the asymptotic distribution has "forgotten" its initial shape in x 2 . We say that the contracting direction has been stabilised.

Shear Dispersion in 2D
So far we have only considered case (2a) in (16). For case (2b), we have from (17) T Inserting (34) in (14), we have We can then explicitly do the time integral in the exponential, The enhancement to diffusion in this case is reflected in the cubic power of time in the exponential. This is not as strong as the exponential enhancement of case (2a), but is nevertheless very significant. This phenomenon is known as shear dispersion or Taylor dispersion. The mechanism is often called the venetian blind effect. Assuming the initial distribution θ 0 depends only on k 01 , then lines of constant concentration which are initially vertical are tilted by the shear flow, in a manner reminiscent of venetian blinds. The distance between the lines of constant concentration decreases with time as (U ′ t) −1 , which gives an effective enhancement to diffusion. The time required to overcome a weak diffusivity is thus If we use k 0 −1 1 as a length scale and U ′ as a time scale, we can define a Péclet number Pe := U ′ /(k 01 2 κ) and rewrite (36) as which should be compared to (23), the corresponding expression for the case (2a). Here there is a power law dependence on the Péclet number, rather than logarithmic, so we may have to wait a long time for diffusion to become important. This makes the linear velocity field approximation more likely to break down. Let us consider the time-asymptotic limit U ′ t ≫ 1: we might then be tempted to neglect everything but the U ′ t k 01 term in the argument of the exponential in (35). However, this would be a mistake. To see more clearly what happens, define the dimensionless time τ := U ′ t and the length The first two terms in the exponential are just what is expected of regular diffusion in the absence of flow. The next term is the enhancement to diffusion along the x 1 direction: it will force the modes k 01 ∼ τ −3/2 to be dominant, since everything else will be damped away. Similarly, the last term forces k 02 ∼ τ −1/2 . Assuming these scalings, the only term that can be neglected for τ ≫ 1 is the very first one, k 0 2 1 τ . We make the change of variable k 01 = τ 3/2 k 01 , k 02 = τ 1/2 k 02 in (38), If we approximateθ 0 ( k 01 τ −3/2 , k 02 τ −1/2 ) ≃θ 0 (0, 0), we can do the integrals in (39) and find For moderate values of The width in the x 2 direction of an initial distribution thus increases as χτ 1/2 = √ κ t. This is independent of U ′ and is exactly the same as expected from pure diffusion. The width in the x 1 direction in (40) increases as χτ 3/2 = U ′ t √ κ t (see Fig. 2).

Three Dimensions
In three dimensions, there are three basic forms for the matrix σ: with corresponding fundamental matrices We can assume without loss of generality that λ 1 ≥ 0, λ 1 ≥ λ 2 , and U ′ > 0, but the sign of λ 2 and λ is arbitrary; however, we must have The case of greatest interest to us is (3a). The relevant k(t), corresponding to (18), is which is used in (14) to give What happens next depends on the sign of λ 2 : the question is whether e −2λ2t grows or decays for t ≫ |λ 2 | −1 . If λ 2 > 0, then we have whilst for λ 2 < 0 Both approximations are valid when t ≫ max(λ −1 1 , |λ 2 | −1 ). For λ 2 = 0 the situation is similar to the two-dimensional case (2a): valid when t ≫ λ −1 1 . The rest of the calculation is very similar to the two-dimensional case (2a), in going from (24) to (32). In both (46) and (47) the x 3 direction is stabilised, that is we need to blow up the k 03 integral to remove the time dependence from the exponential, and find that the integral is dominated by k 03 ≃ 0. The x 2 direction is also stabilised in (47), so we can set k 02 ≃ 0. We thus find for λ 2 ≥ 0, and for λ 2 < 0, where ℓ i := κ/|λ i |. Contracting directions have their spatial dependence given by a time-independent Gaussian, with an overall exponential decay; stretching directions do just that: they stretch the initial distribution, with no diffusive effect. Solutions of the form (48) are called pancakes, and those of the form (49) are called ropes or tubes.
There is another way of thinking about the asymptotic forms (32), (48), and (49) [15]: contracting directions are stabilised near some constant width ℓ j , and expanding directions lead to exponential growth of the width of an initial distribution along the direction. Thus, the volume of the initial distribution grows exponentially at a rate given by the sum of λ i 's associated with stretching directions, but the total amount of θ remains fixed (the mean is conserved). Hence, the concentration at a point should decay inversely proportional to the volume, which is exactly what (32), (48), and (49) predict.

Random Strain Models
In Section 2 we analysed the deformation of a patch of concentration field (a 'blob') in a linear velocity field. Though this is interesting in itself, it is a far cry from reality. We will now inch slightly closer to the real world by giving a random time dependence to our velocity field.

A Single Blob
Consider a single blob in a two-dimensional linear velocity field of the type we treated in Section 2.2 (case (2a)). Now assume the orientation and stretching rate λ of the straining flow change randomly every time τ . This situation is depicted schematically in Figure 3. We assume that the time τ is much larger than a typical stretching rate λ, so that there is sufficient time for the blob to be deformed into its asymptotic form (32) at each period, which predicts that at each period the concentration field will decrease by a factor exp(−λ (i) τ ), where λ (i) is the stretching rate at the ith period. The concentration field after n periods will thus be proportional to the product of decay factors, PSfrag replacements θ We may rewrite this as where t = nτ , and is the 'running' mean value of the stretching rate at the nth period. As we let n become large, how de we expect the concentration field to decay? We might expect that it would decay at the mean valueλ of the stretching rates λ (i) . This is not the case: the running mean (51) does not converge to the meanλ. Rather, by the central limit theorem its expected value isλ, but its fluctuations around that value are proportional to 1/ √ t. These fluctuations have an impact on the decay rate of θ.
The ensemble of variables λ (i) is known as a realisation. Now let us imagine performing our blob experiment several times, and averaging the resulting concentration fields: this is known as an ensemble average over realisations. Ensemble-averaging smooths out fluctuations present in each given realisation. We may then replace the running mean Λ n by a sample-space variable Λ, together with its probability distribution P (Λ, t). The mean (expected value) of the αth power of the concentration field is then proportional to The overbar denotes the expected value. The factor e −αΛt gives the amplitude of θ α given that the mean stretching rate at time t is Λ, and P (Λ, t) measures the probability of that value of Λ occuring at time t.
The form of the probability distribution function (PDF) P (Λ, t) is given by the central limit theorem: that is, a Gaussian distribution (31) with meanΛ and standard deviation β/t. The quantity β is a measure of the nonuniformity (or fluctuations) of stretching in the flow, both spatially and temporally. The decrease of the standard deviation β/t with time reflects the convergence of the average stretching rate to the Lyapunov exponentΛ. I will say more about β in Section 4.1, when we look at a practical example. Actually, the central limit theorem only applies to values of Λ that do not deviate too much from the mean. The theorem understimates the probability of rare events; a more general form of the PDF of Λ comes from large deviation theory [16,17], (A derivation of (54) is given in Appendix C.) The function S(x) is known as the rate function, the entropy function, or the Cramér function, depending on the context (that is, which literature one is reading). It is a time-independent convex function with a minimum value of 0 at 0: S(0) = S ′ (0) = 0. If Λ is near the mean, we have which recovers the Gaussian result (53) with β = 1/S ′′ (0). Both (53) and (54) are only valid for large t (which in our case means t ≫ τ , or equivalently n ≫ 1). We can now evaluate the integral (52) with the PDF (54), where we have omitted the nonexponential prefactors, and defined Since t is large, the integral is dominated by the minimum value of H(Λ): this is the perfect setting for the well-known saddle-point approximation. The minimum occurs at Λ sp where H ′ (Λ sp ) = α + S ′ (Λ sp −Λ) = 0, and is unique because S is convex and has a unique minimum. The decay rate is then given by There's a caveat to this: for α large enough the saddle point Λ sp is negative. This is not possible: the stretching rates are defined to be nonnegative (the integral (56) involves only nonnegative Λ). Hence, the best we can do is to choose Λ sp = 0-the integral (56) is dominated by realisations with no stretching. Thus, in that case γ α = H(0), or We re-emphasise: for small enough α, the saddle point is positive and the decay rate is given by (57). Beyond that, we must choose zero as the saddle point and the decay rate is given by (58). To find the critical value α crit where we pass from (57) to (58), observe that this happens as the saddle point nears zero. Thus, we may solve our saddle point equation H ′ (Λ sp ) = 0 by Taylor expansion, But the saddle point will not be small unless the first terms cancel in (59), that is α crit = −S ′ (−Λ). We may thus recapitulate the result for the decay rate, Clearly γ α is continous, and it can be easily shown that dγ α /dα is also continuous.
As an illustration, we use the Gaussian approximation (55) for the Cramér function, with β = 1/S ′′ (0). The critical α is α crit = −S ′ (−Λ) =Λ/β. The saddle point is positive for α <Λ/β, so from (60) we get This is plotted in Figure 4. Notice that the solid curve (for a random flow) lies below the dashed line (for a nonrandom flow). This is a general result: if f (x) is a convex function and x a random variable, Jensen's inequality says that Now, e −αtΛ is a convex funtion of Λ, so we have which means that the rate of decay satisfies γ α ≤ αΛ, which is exactly what is seen in Figure 4. Thus, fluctuations in Λ inevitably lead to a slower decay rate γ α .
Stronger fluctuations also means that the decay rate γ α saturates more quickly with α. Clearly, in the absence of fluctuations we recover the nonrandom result:Λ/β is infinite and only the α <Λ/β case is needed in Eq. (61). If there are lots of fluctuations,Λ/β is small, and there is a greater probability of obtianing a realisation with no stretching. For large enough fluctuations this exponentially-decreasing probability dominates, and we obtain the second case in (60).

Many Blobs
In Section 3.1 we considered the evolution of the concentration of a single blob of concentration in a random straining field. Now we turn our attention to a large number of blobs, homogeneously and isotropically distributed, with random concentrations. We assume that the mean concentration over all the blobs is zero. A simplified view of this initial situation is depicted in Figure 5(a), with shades of gray indicating different concentrations. If we now apply a uniform straining flow of the type (2a) (see Section 2.2), the blobs are all stretched horizontally (the x 1 direction) and contracted in the vertical (x 2 ) direction, as shown in Figure 5(b). They are pressed together in the x 2 direction until diffusion becomes important ( Figure 5(c)). The effect of diffusion is to homogenise the concentration field until it reaches a value which is the average of the concentration of the individual blobs. This is depicted by the long gray blob in Figure 5(d), which will itself keep contracting until it reaches the diffusive length ℓ.
Of course, here the initial concentration field θ 0 represents the concentration of all the homogeneously-distributed blobs together, so it does not decay at infinity: we must thus use Eq. (30) rather than (32). The summing (and hence averaging) over blobs is manifest in Eq. (30), which contains an PSfrag replacements integral over the initial distribution θ 0 in the x 2 direction, windowed by a Gaussian.
In practice, this implies that the expected value of the concentration at a point x on the gray filament is given by where θ (i) 0 is the initial concentration of the ith blob, and · blobs denotes the expected value of the sum over the overlapping blobs at point x (not the same as spatial integration · ). We assume that N ≫ 1 blobs have overlapped. Equation (63) gives the concentration at a point, summed over N overlapping blobs. Of course, Eq. (63) converges to zero for large N , because the blobs average out. Not so for the fluctuations at that point: by the central limit theorem, we have since the initial blobs have identical distributions. The blob-summed fluctuation amplitude θ 2 blobs is thus proportional to the number N of overlapping blobs. But the number of overlapping blobs is proportional to e Λt : as time increases more and more blobs converge to a given x in the contracting direction and overlap diffusively (this can be seen in Eq. (30): the width of the windowing region grows as e λt ). Assuming the variance of θ (i) 0 is finite, we conclude from (64) that Compare this to (51) for the single-blob case: the overlap between blobs has led to an extra square root. Thus, the ensemble averages θ 2 (x, t) α blobs for the overlapping blobs are computed exactly as in Section 3.1, resulting in (60). Because of the assumption of homogeneity, the point-average is the same as the average over the whole domain (see Section 4 for more on this), and we have 3 with γ α given by (60). (In (66) the angle brackets denote spatial averaging, not spatial integration, because the total variance is infinite in this case.)

Three Dimensions
In three dimensions, we will only treat Case (3a) (a purely straining flow) of Section 2.4. For λ 2 < 0, where the asymptotic concentration is given by (49) (ropes), the situation is basically identical to the 2D case of Sections 3.1-3.2: the statistics of the stretching direction λ 1 determine γ α from (60). The contracting directions x 2 and x 3 are stabilised by diffusion. For λ 2 ≥ 0, the asymptotic concentration is given by (48) (pancakes). We have two fluctuating quantities to worry about (λ 1 and λ 2 ). But since the decay rate in (48) only depends on λ 3 , we can instead focus on its fluctuations only. For a single blob, the average (52) is then replaced by where of course Λ 3 is the average of λ 3 . This PDF achieves a distribution of the large-deviation form (54). The analysis thus follows exactly as in Sections 3.1 and 3.2, except the Cramér function for |Λ 3 | must be used. 4

Practical Considerations
One may rightly wonder if the blobs in a random uniform straining flow depicted in Section 3 bear any resemblance to reality. The single-blob scenario doesn't, but the many-blobs scenario has a fighting chance, as we will try to justify here. There are two important considerations: where does the ensemble-averaging come from, and what are the stretching rates given by?
The decay rate (60) depends crucially on ensemble-averaging: with that averaging the decay rate fluctuates wildly for a given realisation. At the end of Section 3.2 we assumed that homogeneity allowed us to generalise from the average at a point to the average over the whole domain. But the average over the whole domain can actually do a lot more for us: it can provide the ensemble of blobs that we need for averaging! Thus, we can forget about speaking of realisations as if we were running many parallel experiments, and instead speak of the moments of the concentration field as given by an average over randomly-distributed blobs. The decay rate will then be naturally smoothed-out over blobs experiencing different stretching histories. The saturation of the decay rate with α in Eq. (60) is due to θ 2α being dominated by the fraction of blobs that have experienced no stretching.
What about the stretching rates λ? Luckily, it is not them but their timeaverage Λ that matters. If we imagine following a blob as it moves through the flow, we can see that this time-averaged stretching rate is nothing but the finite-time Lyapunov exponent associated with this blob and its particular initial condition. A given blob will be constantly reoriented as it moves along in the flow, so its finite-time Lyapunov exponent is not just the average of the stretching rates (in fact, it must be strictly less than this average). But in a chaotic system we are guaranteed that, on average, these reorientations do not lead to a vanishing (infinite-time) Lyapunov exponent. This is guaranteed by the celebrated Oseledec multiplicative theorem for random matrices [18]. 5 We may thus use for P (Λ, t) the distribution of finite-time Lyapunov exponents, which is well-known to have the large-deviation form (54) [19].
The result of these considerations is the local theory of passive scalar decay. It is called local because of the reliance of such a local concept as the finite-time Lyapunov exponents, which come from a linearisation near fluid element trajectories. In Section 4.1 we discuss a specific example. We postpone a discussion of the validity of the local theory until Section 4.2, but for now we point out that it is known to be exact at least in some simple model flows [15,20].
The derivation presented in this Section was based on the work of Balkovsky and Fouxon [15], who used a slightly more rigorous approach. Son [21] also obtained the decay rate (60) using path-integral methods. Earlier, Antonsen et al. [8] derived the decay rate for the second moment θ 2 in terms of the Cramér function, using a different (and not quite equivalent) approach, though they did not allow for the second case in (60).

An Example: Flow in a Microchannel
We illustrate how to compute the decay rates γ α with a practical problem. Specifically, we will use a three-dimensional model of a microchannel. The system is shown in Figure 6. It consists of a narrow channel, roughly 100 µm wide and slightly shallower. These types of channel are widely used in microfluidics applications ("lab-on-a-chip"), and often one wants to achieve good mixing in the lateral cross-section of the channel. This is difficult, since the Reynolds number of the flow varies between 0.1 and 100-far from turbulent. Clever techniques have to be used to induce chaotic motion of the fluid particle trajectories in order to enhance mixing. Stroock et al. [22] used patterned grooves at the bottom of the channel to induce vortical motions, and found that the mixing efficiency was dramatically increased. Here we use a variation on this where the bottom is pattern with an electro-osmotic coating, which induces fluid motion near the wall [23]. The effect of the electro-osmotic coating is well-approximated by a moving wall boundary condition. The pattern is chosen in a so-called herringbone pattern to maximise the mixing efficiency (though not in a staggered herringbone, which is even better but is more difficult to model). Rather than solving the full equations numerically, we adopt here an analytical model based on Stokes flow in a shallow layer [24]. The longitudinal (x) direction is taken to be periodic. The flow is steady, but because it is three-dimensional it can still exhibit chaos. Figure 7 shows two Poincaré sections for the flow. These are taken at two constant x planes, one at x = 0 and the other at the mipoint of the x-periodic pattern. The two colours represent two trajectories that have periodically punctured those planes many times over. It is clear from the Figures that the flow contains large chaotic regions, as well as smaller regular regions (known as islands). We focus here on the chaotic regions. Now that we have established (or at least strongly suspect) the existence of chaotic regions, we can compute the distribution of finite-time Lyapunov exponents. There are many ways of doing this: because we are not interested in exetremely long times, the most direct route may be used. We have an analytical form for the velocity field, so the velocity gradient matrix is easily computed. This allows us to linearise about trajectories in the standard manner [19,25]. Each trajectory will thus have a finite-time Lyapunov exponent associated with it, which shows the tendency of infinitesimally close trajectories to diverge exponentially. This is then repeated over many different trajectories within the same chaotic region, and a histogram is made of the finite-time Lyapunov exponents. This histogram changes with time, as shown in Figure 8. For these relatively early times, it changes dramatically The evolution of the mean and standard deviation of the distribution is shown in Figure 9. The mean is converging to a constantΛ ≃ 0.116, and the standard deviation is decreasing as β/t, with β ≃ 0.168. (BothΛ and β are fitted values, and arise from the complicated nonlinear nature of fluid particle trajectories-they cannot be predicted or calculated from first principles except in the simplest cases.) These facts taken together are strongly indicative that the distribution is converging to a Gaussian of the form (53). This is easily confirmed by plotting the PDFs at different times and rescaling the horizontal axis by √ t, as shown in Figure 10. Note that this case exhibits a particularly nice Gaussian form, which is not necessarily the norm for all chaotic flows.
Using the values forΛ and β we just obtained, we can calculate the decay rates γ α with the Gaussian approximation. The ratioΛ/β is 0.69, so the change in character in (61) occurs at α ≥ 0.69. Since from (66) the decay of θ 2 α is given by γ α , this means that moments of order 2α ≥ 1.38 will decay at the same rate. This includes the variance θ 2 , so we have from (66) The mixing time is thus γ −1 1 ≃ 25 seconds. This is about a factor of four improvement over the purely diffusive time for, say, DNA molecules (κ ≃   10 −10 m 2 s −1 ). This is not spectacular, but can be greatly increased by staggering the herringbone pattern. The mixing time assuming the decay proceeds at the rate of the mean Lyapunov exponentΛ is roughly 9 seconds, so that the fluctuations multiply this by a factor of three! Of course, we do not know if this is actually a good estimate for the mixing time, since we haven't directly solved the advection-diffusion equation numerically: this is prohibitive in a three-dimensional domain for such a small diffusivity. This is one of the advantages of the local theory: it is usually less expensive to compute the distribution of finite-time Lyapunov exponents than it is to solve the advection-diffusion equation equation directly. We will say more on the validity of the local theory in Section 4.2.

Limitations of the Local Theory
So is this local theory of mixing correct? Well, certainly not always, even in the Batchelor regime. There are many assumptions underlying the model, some of them difficult to verify. (Do blobs really undergo a series of stretching events as described here? Do correlations between these events matter?) My feeling is that sometimes it will, but most of the time it won't. More experiments and numerical simulations are needed to get to the bottom of this. For a detailed discussion of possible problems with the local theory, see Fereday and Haynes [26]. They make a good case that the theory must break down for long times: the blobs discussed here meet the boundaries of the fluid domain and must begin to fold. The folding forces them to interact with themselves in a correlated fashion. We enter the regime of the strange eigenmode [27], which has received a lot of attention lately [14,26,[28][29][30][31][32][33][34][35]. Maybe we'll hear more about that in ten years. . . . 6 dinates x to coordinates r measured from a reference fluid trajectory x 0 (t). The coordinates r are not quite material (Lagrangian) coordinates, since we follow the trajectory of only one fluid element.
We thus let and write the concentration field as θ(x, t) =θ(r, t).
The time derivative of θ can be written where ∂/∂t| x denotes a derivative with x held constant. Now from (68) ∂r Spatial derivatives are unchanged by (68): ∇ x θ = ∇ rθ . Hence, inserting (69) into (1), we find We Taylor expand the velocity field in (70) to get where we neglected terms of order |r| 2 . This is only valid if the velocity field changes little over the region we consider (i.e., if it is smooth enough), which is true for large Schmidt number. Equation (71) is the same as (8), and tells us how to find σ(t).

B Volume Preservation
It is useful to know how to prove that a divergence-free vector field will lead to volume-preservation of a small blob of integrated trajectories. Mathematically, we want to go from (11) to (12). We start from the definition of the determinant of a d-dimensional matrix, where we have dropped the 't' subscript from T for this section, and ǫ is the fully-antisymmetric Levi-Civita symbol. Taking a time derivative of (72), since the d terms obtained after using the product rule for derivatives are the same. We use the equation of motion (11) for T, We rearrange the sums on the right, and recognise that, up to a factor of det T, the term in the parentheses is the cofactor representation of the inverse of T: We conclude that if Tr σ = 0 then the determinant of T is constant. Since its initial condition is unity (from Eq. (11)), then it must remain so for all time.

C Large Deviation Theory
In this Appendix we will justify the large-deviation form of the PDF, Eq. (54), assuming little prior knowledge of probability theory. First, we define the characteristic function e −s(k) of a random variable x by e −s(k) = p(x) e −i kx dx , that is, the characteristic function is simply the Fourier transform of the PDF of x. We have s(0) = 0,x = −i s ′ (0), and x 2 −x 2 = s ′′ (0). Now define the random variable X to be the mean of n variables, where the x i are independent and identically distributed with PDF p(x i ) = p(x). How do we find the PDF P (X n ) of X n , in the limit where n is large? First observe that (from here on we drop the subscript on X n ) P (X) = p(x 1 ) · · · p(x n ) δ x 1 + · · · + x n n − X dx 1 · · · dx n since the joint PDF p(x 1 , . . . , x n ) = p(x 1 ) · · · p(x n ) by independence of the x i . The characteristic function e −S(k) for P (X) is then e −S(k) = P (X) e −ikX dX = p(x 1 ) · · · p(x n ) δ x 1 + · · · + x n n − X e −ikX dx 1 · · · dx n dX We do the X integral, and then observe that we get a product of n identical x i integrals, each of which is equal to e −s(k/n) : e −S(k) = p(x 1 ) · · · p(x n ) e −i(k/n)(x1+···+xn) dx 1 · · · dx n = e −ns(k/n) .
Thus the characteristic function for X is the nth power of the characteristic function for x. We can invert the Fourier transform to find the PDF P (X): where K = k/n. We let H(K, X) = s(K) − i KX ; For large values of n the integral in (74) is dominated by the stationary points K sp (X) of H(K, X) (saddle-point approximation): K sp (X) such that ∂H ∂K (K sp , X) = s ′ (K sp ) − i X = 0.
(I will often leave out the X dependence of K sp (X) to shorten the expressions.) In that case we can approximate the integrand in (74) using which allows us to do the integral explicitly: P (X) = n 2π s ′′ (K sp (X)) e −nH(Ksp(X),X) where K sp (X) is given by (75). As a final step, let us calculate the mean of X using this PDF: X = XP (X) dX = X n 2π s ′′ (K sp (X)) e −nH(Ksp(X),X) dX .
Again, for large n we can use the saddle-point method to evaluate this integral. The important observation is that the saddle-point X 0 of H sp (K(X), X) satisfies The ∂H/∂K term vanishes because it is evaluated at K sp ; hence, K sp (X 0 ) = 0, which implies H(K sp (X 0 ), X 0 ) = 0. Inserting this into the integral (77), we find X = X 0 : the mean of X and the minimum of H coincide. This means that it makes sense to define S(X − X) := H(K sp (X), X), with S(0) = 0 and S ′ (0) = 0, which is the sought-after Cramér function. Note also that S ′′ (X − X) = 1/s ′′ (K sp (X)), and that for large n the non-exponential coefficient in (76) can thus be approximated by evaluating it at the saddle-point K sp (X) = 0, with s ′′ (0) = 1/S ′′ (0). The final form of our large-deviation result is thus which is the same as Eq. (54). As a simple example (treated in every textbook, see for example [17]), consider a random variable x with PDF where x + > x − are constants-this is a binomial distribution (or Bernoulli distribution in this case). If we take the mean X of n such variables, what is the PDF of X for large n? First, we compute the characteristic function for x, We take the logarithm to obtain s(k) and find K sp (X) by solving the saddlepoint equation (75), where ∆ := x + − x − and we restrict x − ≤ X ≤ x + . Inserting this into H(K sp (X), X), we find from (78) It is easy to verify that, since X = (1 − ε)x + + ε x − , we have S(0) = S ′ (0) = 0 and X 2 − X 2 = 1/S ′′ (0) = ε(1 − ε) ∆ 2 .
The binomial distribution (80) is a useful model of stretching of an infinitesimal line segment by a uniform incompressible straining flow in two dimensions, assuming the straining axes of the flow change direction randomly at regular intervals τ . If we set x ± = ±λτ = ±β, where λ is the strain rate, then X is the averaged logarithm of the length ℓ of the segment, i.e. ℓ = e nX . Thus, the mth power of the length of the segment will on average grow as ℓ m = e mnX = e −S(i mn) = e −ns(i m) = (1 − ε) e βm + ε e −βm n .
We know that ℓ −2 must be constant in a 2D incompressible flow [13], so that the term in braces in (82) must be unity. We use this to solve for ε, ε = (1 + e 2β ) −1 , which then allows us to use (82) to write the growth rate χ m of line segments as The Lyapunov exponent, which is given by dχ m /dm at m = 0, has a value of λ tanh β for this flow: it is less than for a uniform straining flow because of the time taken for the segment to realign with the new straining axis when its direction changes.