Turbulent particle pair diffusion: A theory based on local and non-local diffusional processes

A re-appraisal of the Richardson’s 1926 dataset [Richardson, L. F. Proc. Roy. Soc. Lond. A 100, 709–737, (1926)] displays an unequivocal non-local scaling for the pair diffusion coefficient, K∼σl1.564, quite different to the previously assumed locality scaling law ∼σl4/3, where σl is the pair separation. Consequently, the foundations of turbulent pair diffusion theory are re-examined here and it is shown that pair diffusion is governed by both local and non-local diffusional processess inside the inertial subrange. In the context of generalised energy spectra, E(k) ∼ k−p for 1 < p ≤ 3, the new theory predicts two non-Richardson regimes depending on the size of the inertial subrange: (1) in the limit of asymptotically infinite subrange, non-local scaling laws is obtained, K∼σlγ, with γ intermediate between the purely local and the purely non-local scalings, i.e. (1 + p)/2 < γ ≤ 2; and (2) for finite (short) inertial subrange, local scaling laws are obtained, K∼σl(1+p)/2. The theory features a novel mathematical approach expressing the pair diffusion coefficient through a Fourier integral decomposition.


Introduction
Turbulent transport and mixing play an essential role in many natural and industrial processes [1], in cloud formation [2], in chemically reacting flows [3] and combustion systems [4,5], and atmospheric and oceanographic turbulence determines the spread of pollutants and biological agents in geophysical flows [6][7][8][9]. Turbulent concentration fluctuations often play a critical role in such systems, and this is related to the separation of nearby fluid particles.
Turbulent particle pair diffusion (or relative diffusion) was introduced by L. F. Richardson [10], who laid the foundations for a theory of how ensembles of pairs of fluid particles (tracers) initially close together move apart due to the effects of atmospheric winds and turbulence. Pair diffusion is usually classified into three regimes depending on the pair separation distance l relative to the Kolmogorov length scale η of the turbulence: (i) the dissipation subrange where l ( η; (ii) the inertial subrange when η ( l ( L, where L is some outer length scale of the turbulence (such as the integral length scale or the Taylor microscale); (iii) at much larger times a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 when l ) L, their motions become independent and the pair diffusion collapses to twice the one-particle Taylor diffusion, hl 2 i ! 2hx 2 i $ t [11].
Batchelor introduced a fourth regime, the so-called ballistic regime [12]. He noted that turbulence is correlated in time and space, so for very short times after release the motion should be well correlated with its initial conditions, thus the relative velocity is approximately constant for a short time which leads to hl 2 i $ l 2 0 þ v 2 0 ðt À t 0 Þ 2 , where v 0 = v(t = t 0 ) is the rms pair velocity at initial time t = t 0 . In fact, this is true in any spatio-temporally correlated velocity field and not just in turbulence.
In this work, we focus our interest on inertial subrange scaling laws, and we will ignore the other regimes which are fairly well understood at the present time.
Richardson argued that as particle pairs separate in a field of turbulence, the rate at which they move apart is affected by different scales of motions, and the energies in the scales of motion (eddies) of the same scale as the pair separation are the most effective in the diffusional process-this is the so-called locality hypothesis. Most theories of turbulent pair diffusion since then have assumed locality.
Richardson was also motivated by a desire to bring molecular and turbulent pair diffusional processes into a unified picture through the use of a non-Fickian diffusion equation with scale dependent diffusion coefficient, K(s). Assuming homogeneous isotropic turbulence, Richardson posed the problem in 3D in terms of the probability density function (pdf) of the pair separation, q = q(s, t), subject to the normalization, R 1 0 4ps 2 qðs; tÞds ¼ 1, he suggested the following diffusion equation to describe q, Here, s is the sample space variable for the random pair separation distance. K(s) is also called the pair diffusion coefficient, or the pair diffusivity.
From observational data of turbulent pair diffusion coefficients collected from different sources, Richardson assumed an approximate constant power law fit to the data, K(l) $ l 4/3 . This is equivalent to [12,13] often referred to as the Richardson-Obukov t 3 -regime. l(t) is the pair separation, at time t, and the angled brackets is the ensemble average over particle pairs. It is no longer believed that it is possible to unify molecular and turbulent diffusional processes because their physics are fundamentally different; Brownian motion characterizes molecular diffusion, while convective gusts of winds that increase the pair separation in surges characterizes turbulent diffusion [14][15][16]; but the idea of a scale dependent turbulent diffusivity has survived.
Richardson's assumed 4/3-scaling law is consistent with Kolmogorov turbulence K41 [17] in the following manner. The K41 energy spectrum in the inertial subrange is, E(l) $ ε 2/3 l 5/3 , from which it follows that the pair diffusion coefficient depends only upon l and ε (the rate of kinetic energy dissipation per unit mass). The typical pair relative velocity in eddy scales of order l is vðlÞ $ The admittance of a solution to Eq (1) from a point source implies that the initial separation of paricle pairs is effectively zero, and therefore the Kolmogorov scale must also asymptote to zero. Richardson's theory is thus true strictly only in the asymptotic limit of infinite Reynolds number, Re ! 1 (which is equivalent to infinite inertial subrange, R k ! 1), and also l 0 ! 0, as t ! 0.
With the 4/3-scaling for K, an explicit self-similar solution of Eq (1) exists for diffusion from a point source with boundary conditions q(0, t) = q(1, t) = 0, qðr; tÞ ¼ 429 70 ffiffiffiffiffiffiffi ffi 143 2 Turbulence is both scale dependent and time correlated (non-Markovian), and it is not clear whether the pdf in Eq (3), which describes a local and Markovian process, can accurately represent the turbulent pair diffusion process. However, attempts have been made to derive alternative non-Markovian models for pair diffusion, such as Levy Flight type models, [18][19][20][21][22][23], and this remains a topic of active research in the field.
However, a time dependent diffussion coefficient is hard to justify physically if we assume steady state equilibrium turbulence, because it implies unrealistic values for K as time progresses: for b > 0 the pair diffusivity increases in time without limit; for b < 0 the pair diffusivity approaches zero in time and the separation process effectively stops. Both cases seem unlikely, and in the ensuing we will restrict the discussion to steady state equilibrium turbulence and consider only the case, b = 0.
It is worth remarking, however, that time-dependent diffusivities like Eq (4) may be appropriate in the context of non-equilibrium turbulence. Furthermore, Hentchel & Procaccia [26] discuss time-dependent diffusivity it the context of clouds of particles; in this case, we are not dealing with point-sources and the description of the pair diffusion by a diffusion-type equation like Eq (1) is questionable.
Some of the questions that we address in this work are as follows. What evidence is there for the locality hypothesis? Can we develop a theory for turbulent pair diffusion without introducing the assumption of locality from the beginning? What scalings laws for the pair diffusion coefficients does this yield, firstly in the limit of infinite inertial subranges (infinite Reynolds number), and secondly for finite inertial subranges (moderate Reynolds numbers)?
Here, we construct a new theoretical framework through a novel method of analysis namely the Fourier decomposition of the relative pair velocity. We derive the scaling laws for the pair diffusion coefficient without making the assumption of locality from the outset.
In [27], we investigate this new theory through numerical simulations and where all predictions of the theory presented here have been verified.
The rest of the paper is organised as follows. In Section 2 we re-examine the body of evidence available on turbulent pair diffusion especially from large scale turbulence containing large inertial subranges. In Section 3 we focuss upon a reappraisal of Richardson's original 1926 dataset. In Section 4, a new theory that does not a priori make the assumption of locality is constructed through a novel mathematical method from which new scaling laws for the pair diffusivity is derived. In Section 5 some precitions from the new theory is derived. We discuss the implications of the theory in Section 6.

What is the evidence for locality?
The general consensus among scientists in the field at the current time is that the collection of observational data, experimental data, and Direct Numerical Simulation, suggests a convergence towards a Richard-Obukov locality scaling. However, the relatively low Reynolds numbers in the experiments and DNS, and the high error levels in collecting the data means that this is by no means a forgone conclusion. As noted by Salazar and Collins [25], "‥ there has not been an experiment that has unequivocally confirmed R-O scaling over a broad-enough range of time and with sufficient accuracy." It is not known what size of the inertial subrange is required to observe unequivocally the pair diffusion scaling, but it is widely assumed to be many orders of magnitude. Only geophysical turbulence, such as in the atmosphere and in the oceans, can produce such extended inertial subranges. We define the size of the inertial subrange R k to be, where the inertial subrange part of the energy spectrum E(k) is defined in the wavenumber range k 1 k k η . Observations of approximate R-O t 3 -regimes in geophysical flows, have been reported by Tartarski [28], Wilkins [14], Sullivan [29], and Morel & Larchaveque [30]. More recent observations include Julian et al. [31] in the atmosphere, and Lacasce & Ohlmann [32], and Ollitrault & Gabillet [33] in the oceans. But the high error levels and other assumptions, such as two-dimensionality, made in these observations means that a firm conclusions cannot be drawn from them regarding pair diffusion theories.
Cases in point are the recent experiments in the atmosphere [34] and in the Nordic sea [35]. In [34], the authors revisited The EOLE experiment in 1973 to study turbulent processes in the lower stratosphere circulation Relative dispersion of balloon pairs was studied by calculating the finite-scale Lyapunov exponent, an exit-time-based technique. The improved analysis supports a k −5/3 behavior in the mesoscale range 100-1000 km. However, they were unable to deduce the origin of this spectrum-whether it concerns 2D inverse energy cascade, gravity wave breaking with direct energy cascade, or shear (zonal) dispersion in a diffusive (meridional) field.
In [35], the authors examine the relative dispersion of surface drifters deployed in the POLEWARD experiment in the Nordic Seas during 2007-2008. The authors found some evidence for Richardson pair diffusion but could not rule out that this may be due an inverse energy cascade as this is a quasi-2D system. The deformation circle when shear effects become important is small in the Nordic sea, so the Richardson diffusion was observed only in a short range over one decade of scales, 10-100 km, in separation distance.
Direct Numerical Simulations (DNS) is inconclusive at the current time because it does not produce a big enough inertial subrange in order to test pair diffusion laws convincingly. For example, Ishihara et al. [36] perform a DNS with 4096 3 , at Taylor scale based Reynolds numbers R λ % 1200 showing an approximate inertial subrange energy spectrum over a very short range of just 40. Other DNS of particle pair studies at low Reynolds numbers are, Yeung [37] at R λ = 90, Boffeta & Sokolov [38] at R λ = 200, Ishihara & Kaneda [39] at R λ = 283, Yeung& Borgas [40] at R λ = 230, Sawford et al. [41] at R λ = 650. Scatamacchia et al. [42] at R λ = 300. See also Bitane et al. [43], and Biferale et al. [44]. The maximum separation of time scales between the integral time scale and the Kolmogorov times scale observed to date in DNS is about a factor of, 100. This is still too small to fully test inertial subrange pair diffusion scalings. For 2D turbulence, see [45,46].
Particle Tracking Velocimetry (PTV) laboratory experiments Mass et al. [47], and Malik et al. [48] have been providing pair diffusion statistics at low to moderate Reynolds numbers. Like DNS, these Reynolds numbers are too small to reliably test pair diffusion laws. Virant & Dracos [16] obtain pair diffusion measurements from PTV. More recently, Berg et al. [49] obtained measurements in a water tank at R λ = 172, and Bourgoin et al. [50] and Ouellete et al. [51] tracked hundreds of particles at high temporal resolution at R λ = 815. Although higher resolution tracking experiments using high-energy physics methods have been performed for single particle trajectories [52], they have not yet been applied to particle pair studies.
Apart from experiments and DNS, many models of turbulent relative particle diffusion have been proposed from classical stochastic random flight type models [53,54], to more recent multifractal models of turbulence [55]. However, none of these models address the problem of local and non-local diffusional processes directly, and as such it is not possible to compare them with the new model developed here. We will therefore not refer to such models any further in this study.
Datasets, both numerical and experimental, on relative dispersion are being made available on a number of databases such as [56] and [57].
In summary, at the current time the scaling laws for inertial subrange pair diffusion remain inconclusive.

A reappraisal of the 1926 dataset
In 1926 Richrdson reported in the Proceedings of the Royal Society of London data on turbulent pair diffusivities collected from different sources [10], which is reproduced here with a brief description in Table 1. He plotted the turbulent diffusivity against the pair separation in log-log scale, shown as the red and black filled circles in Fig 1. Motivated by an attempt to unify pair diffusional processes across all possible scales in the limit of infinite Reynolds number, he made two important assumptions: firstly that the pair diffusion can be modeled by a diffusion-type equation, Eq (1), and secondly that the diffusion coefficient is scale dependent and can be modelled as a unique power law across all scales. From the collected data, he assumed the scaling K $ l 4/3 as a reasonable fit (Fig 1, dotted blue line). This is not the least squares line of best fit to the data, it is just an approximation to the data. The actual line of best fit is, K $ l 1.248 (Fig 1, solid red line). Later, Richardson and Stommel [58,59] commenting on diffusion of floats in the sea noted that their new measurements were roughly consistent with the 1926 data. At the end of their paper they wrote, 'Note added in proof.-After this manuscript was submitted the writers have read two unpublished manuscripts by C. L. von Weisaecker and W. Heisenberg in which the problem of turbulence for large Reynolds number is treated deductively with the result that they arrive at the 4/3-law. The agreement between von Weisaecker and Heisenberg's deduction and our quite independent induction is a confirmation of both', see [60].
But they also observed that, 'any power between $ l 1.3 and $ l 1.5 would be a tolerable fit to the data'. Indeed, the 1926 data are from very different sources varying in length and time scales by orders of magnitude. Apart from this there are other physical processes whose impact on pair diffusion cannot be readily assessed, such as the effects of bouyancy, and the quasi two-dimensionality of the atmosphere. Similar concerns can be made regarding the 1948 and 1949 data. No attempt was made to assess the magnitude of these effects, so it is not surprising that such a wide error band in the power law scaling was given.
There is a more fundamental problem with the data. Data-point number N1 in Table 1 ( Fig  1, red filled circle) is not from turbulence measurements at all, rather it is from studies of the molecular diffusion of oxygen into nitrogen and whose length scale is stated to be of the order of 10 −2 cm. At such a small scale it cannot be turbulent, and certainly cannot contain an inertial subrange. This data-point is therefore disregarded as an outlier in the current investigation which strictly demands the existence of turbulence with an inertial subrange.
The remaining six data-points (Fig 1, black filled circles) are sound, coming from geophysical turbulence settings and certainly containing extended inertial subranges. The line of best fit to this new improved dataset (N2 to N7) displays an unequivocal non-local scaling, K $ l 1.564 (Fig 1, solid black line). The coefficient of determination, R 2 = 0.97, thus all the data-points lie close to line of best fit, Fig 1. The correct scaling must be close this; but even taking a reasonably wide margin of error, say between $ l 1.5 and $ l 1.6 , this is outside the range noted by Richardson & Stommel.
The 1926 dataset is used as a guide to inspire a new idea-locality in the case of Richardson. Richardson and Stommel were aware of the wide margin of error in the data, but this did not prevent them from making a reasoned guess on the scaling power based on this data. The situation is similar here, except that the dataset has been corrected and now shows a clear shift towards non-locality, indicating that non-local diffusional processes cannot be ignored a priori in a general theory of turbulent pair diffusion, which is where we turn to in the next section.

A new theory
Like Richardson the focus here is on the diffusion coefficient K. Although the mean square separation hl 2 i is often the focus of diffusion studies, it is related directly to the diffussion coefficent by the exact relation, K = 0.5dhl 2 i/dt. If a general scaling, K $ s g l , is assumed then it yields hl 2 i $ t χ . For steady-state equilibrium turbulence χ is given by the relation, As such hl 2 i does not provide any additional information. We will therefore focus our attention mainly on K in the ensuing work, and refer to hl 2 i only where necessary. Our interest is scaling laws inside the inertial subrange, so we will ignore discussion of the viscous diffusion range and the long time Taylor diffusion. Furthermore, as our focus is on scalings rather than exact quantities we will supress unimportant constants wherever possible.
We will consider turbulence with generalized energy spectra of the form, E(k) $ k −p . The term 'wavenumber' will be used to denote both the vector, k, as well as its magnitude, k = |k|, and the context will make it clear. Generalised spectra of this type are routinely used in pair diffusion studies; it helps in understanding the balance of diffusional processes as the turbulence spectra changes. In the current work, generalised spectra will play an important role in developing the new theory, and also for testing the predictions of the theory [27]. The pair diffusion coefficient K(l, p) is now also a function of p. Eq (7) is true for generalised spectra as well, with γ = γ(p), and χ = χ(p).

The statement of the problem
The problem is to determine the pair diffusivity, K = h l Á v i, of an ensemble of pairs of fluid particles in a field of homogeneous turbulence with an energy density spectrum, E(k) containing an generalised inertial subrange, E(k) $ k −p , k 1 k k η , for 1 < p 3, and such that E(k) ! 0 as k ! 0. The particles in a pair are located at x 1 (t) and x 2 (t) at time t, the pair displacement vector is l(t) = x 2 (t) − x 1 (t), and the pair separation is lðtÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The initial separation at some earlier time, t 0 , is denoted by l 0 = |x 2 (t 0 ) − x 1 (t 0 )|. The turbulent velocity field is, u(x, t), and the particle velocities at time t are, respectively, u 1 (t) = u(x 1 (t), t) and u 2 (t) = u (x 2 (t), t), and the pair relative velocity is v(t) = u 2 (t) − u 1 (t).
We assume point source release, which in practical terms means that the initial pair separation must be close to the Kolmogorov length scale, l 0 % η. The particles will diffuse apart and eventually decorrelate with the initial conditions-they will 'forget' their initial conditions, (l 0 , v 0 ), as Batchelor put it-after some travel time, t l 0 , when the pair is inside the inertial subrange. During this travel time, the pair will display ballistic motion with essentially constant velocity equal to the initial relative velocity v = v 0 ; The transition from the ballistic regime to the explosive inertial subrange regime occurs on a time scale of the order of the eddy turnover time scale, t l 0 % t 0 ðl 0 Þ $ ε 1=3 l 1=3 0 . If l 0 % η then the travel time is approximately equal to the Kolmogorov time microscale, t l 0 % ε 1=3 Z 1=3 , which is very short. At much longer times we can ignore the ballistic regime because t ) t l 0 and hl 2 i ) l 0 .
Without loss of generality, it will also be assumed that, t 0 = 0.

The mathematical framework
The main task here is to uncover what scales of turbulent motions contribute to the pair separation process when the pair separation l(t) is inside the inertial subrange. For this purpose, we decompose the pair diffusion coefficient in to its component scales using Fourier transforms. This is be obtained as follows.
The pair relative velocity is, v(x) = d l / dt, and K is defined as the ensemble average of the scalar product of v with l, For homogeneous, isotropic, incompressible, reflectional and statistically stationary turbulence, the Fourier expression for the velocity field u is [69], where A(k) is the Fourier transform of the flow field, k is the associated wavenumber. The relative velocity v across a finite displacement l, is Using Eq (9) this is gives, Taking the scalar product of v with l, and then the ensemble average hÁi over particle pairs yields an expression for hl Ávi. But the left hand side is a Lagrangian quantity, while the right hand side is an Eulerian quantity. We assume that the Lagrangian ensemble scales like the Eulerian ensemble; such a closure is often made in diffusion studies. Thomson & Devenish [70] for example make this assumption implicitly in their analysis of Lagrangian diffusion models.
We thus obtain a scaling for the pair diffusivity, Because of homogeneity, the ensemble average removes the factor exp (i k Á x 1 ) without altering the scaling upon l. This gives, Let k l = 1/l be the pair separation wavenumber; we follow the usual convention and replace l with the scaling, l $ σ l , throughout this work, so that where s 2 l ¼ hl 2 i. Note that k l (t) changes with time.

Turbulent transport processes
There are several transport processes that must to be considered. The largest eddies carry the smaller scales of motion, as depicted in Richardson's poem 'Big whorls have little whorls'. This is known as the sweeping effect, which means that in a frame of reference moving with these large scale convective motions the small scale transport process should proceed as if there were no large scale sweeping motions at all in a statistical sense. Thus pair diffusion inside the inertial subrange, which is a small scale process, should be the same in both frames of reference. This allows us to simplify the problem statement by working in the swept frame of reference by eliminating the large scale part of the energy spectrum by setting E(k) = 0 for k < k 1 . In this frame of reference, E(k) $ k −p , k 1 k k η , and E(k) = 0 outside of this range of wavenumbers.
In the swept frame of reference, consider a particle pair whose separation distance is l(t), and the associated wavenumber is k l in Eq (14). In principle, all scales of motion in the inertial subrange could affect the pair diffusion process.
The eddies at wavenumbers that are much bigger than k l (the smallest scales of motion) probably act like molecular diffusion motion since the time scales are very small and the motions are somewhat randomised due to the high frequencies at these scales.
The eddies at wavenumbers close to k l are considered to be in the local wavenumber range to k l .
The large scale motions at wavenumbers that are much smaller than k l will be less correlated with the motions at scale k l . The largest of these scales become the sweeping scales which are statistically decoupled from the local wavenumbers other than carrying them en block. As the wavenumber spectrum is continuous, it follows that between the convecting scales and the wavenumbers local to k l , there could be a range of wavenumbers that are weakly correlated with the motions near k l while still being inside the inertial subrange, and this we call the nonlocal wavenumber range to k l .
In this picture, the hypothesis of locality corresponds to the assumption that the non-local scales of motion do not contribute to the inertial pair diffusion process.
Here, we do not make such an assumption. The physical assumption adopted here about the nature of the diffusional processes that are occurring in the system is that there exist three broadly independent diffusional processes within the inertial subrange that potentially contribute to the pair diffusion process as a whole, each process acting from its own range of wavenumbers relative to the inverse pair separation wavenumber k l .
The three physical processes operate, respectively, from (i) the wavenumbers that are larger than k l , (ii) the wavenumbers that are local to k l , and (iii) the wavenumbers that are non-local to k l . The associated frequencies are, oðkÞ / ffiffiffiffiffiffiffiffiffiffiffiffiffi k 3 EðkÞ p , according to the usual assumption that the frequencies scale with the inverse turnover time of the eddy at wavenumber k. The local eddy frequency is, o l / ffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 3 l Eðk l Þ p , and the local eddy turnover time is, T l $ 1/ω l . The integral in Eq (13) is then the sum of three integrals over different wavenumber ranges which are defined as follows: s: the small scales such that k ) k l , and |k Á l| ) 1, and the associated frequencies are much larger than ω l , ω(k) ) ω l .
l: the local scales such that k % k l , and |k Á l| % 1, and the associated frequencies are of the same order as ω l , ω(k) % ω l .
nl: the non-local scales such that k ( k l , and |k Á l| ( 1, and the associated frequencies are much smaller than ω l , ω(k) ( ω l . Eq (13) then becomes, which we rephrase as, We might approximate the integrand in Eq (15) by expanding the exponential term such that, exp (i k Á l) − 1 % i k Á l. However, such an expansion is accurate only for low wavenumbers k ( k l where, |k Á l| ( 1. For local wavenumbers k % k l where |k Á l| % 1, this expansion is only approximately true. For high wavenumbers k ) k l where |k Á l| ) 1, this expansion is not accurate.

The physics of the small scales of motion
A simplification can be made with respect to the small scales of motion whose contribution to the diffusion process is, There is no need to evaluate this integral directly because the net ensemble effect can be assessed on physical grounds alone. K s is an integral over high wavenumbers and represents the contribution from scales of turbulent motion which are much smaller than the pair separation, i.e., from k ) k l . The energy contained in these scales is very small if the energy spectrum decreases as k increases, such as an inverse power law of the type E(k) $ k −p , with p > 1.
Furthermore, these small scales are associated with the unsteadiness of high frequencies, ω ) ω l . Statistically, these high frequency motions induce random and rapid changes in the direction and the magnitude of the pair displacement vector.
Overall, the net statistical impact of these high frequency, low energy, and random turbulent velocity fluctuations on the pair diffusion process is expected to be extremely small, i.e., K s ( Max(K l , K nl ). Thus, we assume that K s % 0 and K s will be ignored from now on.

The physics of the local and non-local scales of motion
With the effect of the small scale contributions eliminated, the simplified expression for the pair diffusivity is, The expansion of the exponential in the integrand to leading order is accurate only in the non-local range because |k Á l| ( 1. In the local range where |k Á l| % 1, such an expansion is only approximately true because the local eddies are moderately unsteady with frequencies that are of the same order of magnitude as, ω l . The effect on the local diffusion process is assumed to be likewise moderate. We will assume that the ensemble effect of the unsteadiness from the local wavenumbers is to reduce the magnitude of K l , but without altering its overall scaling behaviour. Then, the expansion, exp (i k Á l) − 1 % i k Á l, can be used in Eq (18) but with the magnitude of K l reduced by some factor, F l ≲ 1-this constant is smaller than unity, but not too small. Then (18) becomes, is not expected to be a universal constant because it will depend upon various parameters, like p, and R l = k l /k 1 (which is the size of the inertial subrange relative to the particle pair separation), and also implicitly upon the size of an eddy in wavenumber space, C (to be defined later). After absorbing constants, the integrands in Eq (19) become, where α is the angle between l and A, and β is the angle between l and k. For isotropic random fields averaging (20) over all directions, again, does not affect the scaling behaviour. α and β are not uniformly distributed in all direction, it is well known that l aligns preferentially in the positive strain directions, and this ensures that the ensemble average above is non-zero.
Retaining the angled brackets hÁi to include averaging over all directions, Eq (19) with Eq (20) simplifies to, where a = |A|, and dA(k) is the element of surface area at radius k in wavenumber space. If the closure, hl 2 aki $ hl 2 ihaki, is assumed, then upon integrating over the surface area this integral becomes Note that haki 6 ¼ 0 even though the vectors k and A are orthogonal because a and k are magnitudes.
R ha 2 idA(k) is the energy density per unit wavenumber averaged over all directions [69] and scales like $ E(k)/k. If the closure, R hakidAðkÞ $ k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R ha 2 idAðkÞ p is assumed then Eq (22) becomes, As a further check, this can also be derived as follows. The velocity variance from the scales k to k + dk is E(k)dk, and the variance of velocity gradient is k 2 E(k)dk. The particle pair velocity variance is, $ hl 2 ik 2 E(k)dk. The time scale of eddies of wavenumber k is 1= . So the incremental contribution to the diffusivity from these scales is the pair velocity variance times the time scale, dK $ hl 2 i ffiffiffiffiffiffiffiffiffiffiffi kEðkÞ p dk, which leads to Eq (23). To make further progress the actual form of the turbulence spectrum must be specified. In this work, the focus is on turbulence which contains an inertial subrange. For pair diffusion statistics, as mentioned earlier, this is implemented by working in the swept frame of reference by setting the spectrum in the large energy scales to zero, E(k) = 0 for k < k 1 , and assuming an inverse power-law energy spectrum in the inertial subrange, where C k is a constant. A large length scale L is necessary for dimensional consistency. L scales with some length scale that is characteristic of the large energy scales, such as the integral length scale, or the Taylor length scale. With the spectrum in Eq (24), and with s 2 l ¼ hl 2 i, Eq (23) becomes, This is the most general expression for K that can be derived from the present analysis without any a priori assumption regarding locality.

Validation: The locality limit
To check the effectiveness of the mathematical approach adopted here in deriving Eq (25), the locality limit from this expression must first be validated against Richardson's locality hypothesis for which there is a known theory.
The assumption of locality means that the non-local (first) term on the right hand side in Eq (25) is ignored. To evaluate the remaining local integral we assume a cut-off wavenumber k Ã l that separates the local and non-local ranges, and such that k 1 ( k Ã l < k l . Thus, k 1 k < k Ã l defines the non-local kernel (range) to k l in wavenumber space, and k Ã l k k l defines the local kernel to k l .
Using the scaling hl 2 i $ 1=k 2 l the local integral in (25) yields, Let the size of the locality kernel with respect to k l be defined by, where C is finite and greater than unity. The size of the inertial subrange with respect to k l is, Note that R l (t) changes with time because k l (t) changes with time. R l , is related to a local Reynolds number through, Re l $ R 4=3 l . Thus, all dependencies on R l could be replaced by dependencies on Re l , e.g. C = C(p, Re l ). We will use R l in the current analysis.

The influence of non-local scales
A priori there is no reason to neglect the non-local term, which is the first term on the right hand side in Eq (25), This is equivalent to strained relative motion where each scale of turbulence in the nonlocal wavenumber range, k 1 < k < k Ã l , will set up a straining field in the neighbourhood of k l which will alter the rate of increase of the pair separation. Previous theories have assumed that such non-local effects are negligible. However, there are three factors that suggest that this may be an oversimplification.
Firstly, the non-local wavenumbers, k 1 < k < k Ã l , possess much greater energies than at the local separation wavenumber k l , and this will increase their relative influence in the pair diffusion process.
Secondly, the time scale of the non-local scales, T nl (k), are much larger than the local turnover time scale, T nl (k) ) T l $ 1/ω l , so the straining fields set up by non-local wavenumbers will persist for longer times than T l , which will enhance their effectiveness.
Thirdly, although an individual non-local wavenumber contribution is weak, the integral is over a large part of the energy spectrum. Again, this will enhance the effectiveness of the nonlocal scales in the pair diffusion process.
Taken together, there is a fair chance that the total non-local contributions could be significant at least in some parameter range.
Changing variables in the integrand in Eq (31) to s = k/k 1 , yields Inside the integral the non-local wavenumbers in range of integration ½k 1 ; k Ã l do not scale with k l so the integral is a definite integral producing a non-dimensional number, S nl = S nl (p, R l , C). S nl is not expected to be a universal constant. This gives, If the upper end of the inertial subrange is assumed to scale with the large scales, k 1 $ 1/L, then this simplifies further to, γ nl (p) is the non-locality scaling, and it is always equal to 2, independent of p. K nl is thus always strain dominated being proportional to s 2 l .

A general expression for the pair diffusivity
The overall expression for the turbulent pair diffusion coefficient is therefore, or simply, where, Eq (36) with Eqs (37) and (38) means that turbulent pair diffusion is governed by a local diffusional process, and a non-local diffusional process which is produced by straining fields set up by the large scales.
But, what is the relative balance of their contributions? What does Eq (36) imply for the overall scaling for K?

The scaling laws for K
To address these questions, we need to first establish some general properties of the function C in Eq (27).
C ¼ k l =k Ã l is a fundament physical quantity. It represents the range of wavenumbers close to k l over which motions are well correlated with motions at scale k l . It is a type of non-dimensionalised correlation scale for particle pairs. We call it the locality kernel, or we may call it the locality correlation length. In essence, this means that motions with scales within this kernel are deemed to scale locally with k l , but motions with scales outside this kernel do not scale with k l .
C is finite, 1 < C < 1, and it depends on R l , and on the energy spectrum E(k) $ k −p . E(k) becomes less steep as p ! 1, so there is relatively more energy in the smaller scales when p approaches 1; this indicates that C must increase in this limit, C ) 1. On the other hand, as p ! 3, then C ! 1 because all scales of motion act non-locally so that k Ã l ! k l . The balance of the local and non-local diffusion is defined as the ratio, M K = K nl /K l . M K is a function of, p, R l , and C, and from Eqs (29) to (34) we obtain, Usually, we expect that C ) 1, which simplifies the above to, Furthermore, in subranges that are bigger than the locality kernel, R l /C ) 1, this simplifies further to, We can now obtain some properties of M k , and hence predict some scalings for the pair diffusion coefficient.

Scaling for Kolmogorov turbulence p = 5/3
Before we look at the general case, 1 < p 3, it is instructive to look first at the Kolmogorov turbulence case p = 5/3. Then in the limit of infinite subrange R l /C ! 1, Eq (41) is, This is a finite limit and we can therefore expect the balance between local and non-local processes to befinite. In other words, neither one process nor the other will be completely dominant.
But in the limit of short finite subranges such that R l /C ! 1, i.e. R k /C > R l /C % 1, Eq (40) gives, Thus, for reasonably short inertial subranges, the non-local processes are negligible and we expect locality to be completely dominant, at least approximately.
For an ad hoc value of F l = 0.25, and C = 100, we obtain M k % 0.15 in the limit of R l /C ! 1. This is illustrated in Fig 2, which shows plots of M k against log(R l ), for p = 5/3, C = 100, and for three choices of F l = 0.1, 0.25, 0.5. These are only estimates, but it illustrates the finite balance between the local and non-local diffusional processes.
If the inertial subrange is too short, R l < C, then the above analysis breaks down and we do not observe any kind of inertial range scaling.
In summary, for asymptotically infinite inertial subrange (Reynolds number) we obtain a finite balance between local and non-local diffusional processes; but in the limit of short inertial subrange locality is restored because the non-local processes are negligible.
Do these limiting cases also exist for general spectra?. We investigate this in the next section.

Infinite inertial subrange R k ! 1 (Re ! 1)
First, we derive the scaling laws for the pair diffusion coefficient in the limit of infinite inertial subrange (infinite Reynolds number). Henceforth we assume that, F l < 1, and R l > C > 1, and R k > R l , unless otherwise stated.
From Eq (40), as p ! 1, then M K ! (1/C − 1/R l )/F l . For large inertial subranges, R k /C > R l /C ) 1, we obtain M K ! 1/F l C ( 1, and therefore locality dominates in this limit as expected. From Eq (39), it is easily shown that as p ! 3, then M K ! 1, and therefore non-locality always dominates in this limit, again as expected. (Exact non-locality corresponds to C = 1).
In the intermediate range, 1 < p < 3, although the balance M K cannot be obtained quantitatively without the explicit form of F l and C as a functions of p and R l , it is instructive to investigate M K with some test functions for C to uncover some general trends in M K .
For this purpose, we simplify further by assuming an ad hoc value of F l = 0.25. We choose smooth exponential type test functions, C(p) = 1 + 100(exp(−A(p − 1)) − exp(−2A)), such that C ) 1 at p = 1, and C = 1 at p = 3. A > 0 is a given constant. M K was calculated using C(p) in (39). https://doi.org/10.1371/journal.pone.0202940.g003 Turbulent particle pair diffusion indicating that both local and non-local diffusional processes could play significant roles in large inertial subranges.
We can draw the following conclusions for large inertial subranges where R k /C > R l /C ) 1. Firstly, as p ! 1 then M K ( 1, and therefore K nl ( K l , yielding the locality limit, Fourthly, The influence of the non-local process increases at any p in the rangel 1 < p < 3 as the size of the inertial subrange increases. From these considerations, there is a fair chance that in the critical range of spectra close to Kolmogorov p = 5/3 both the local and non-local processes could exert significant influence in the diffusional process. We now apply an extension of Richardson's second hypothesis, that at any p the diffusion coefficient is described by a single power, say γ(p), across all scales in the limit of infinite inertial subrange. This is plausible because self-similarity is exact in this limit-that is you cannot tell the difference in the scaling at different scales, so you must observe the same power at all scales. https://doi.org/10.1371/journal.pone.0202940.g011

Turbulent particle pair diffusion
Since M K is a smooth and continuous function of p in the range 1 < p 3, then K(l, p) must also be a smooth and continuous function of p in this range and it must display a smooth transition between its asymptotic limits, Kðl; 1Þ $ s 1 l and Kðl; 3Þ $ s 2 l , as p passes smoothly from 1 to 3.
Is it possible that either one of the local or non-local processes dominates throughout the inertial subrange for any given p? That would imply a discontinuous jump between locality and non-locality scalings at some value of p in order to satisfy the asymptotic limiting cases in (44) and (45), but this would violate the continuity in K as a function of p.
Eq (46) is the main prediction of the new theory in the limit R k /C ! 1, with γ l (p) and γ nl (p) given by (37) and (38). This is equivalent to the mean square separation scaling, hl 2 i $ t χ(p) , with χ(p) given by Eq (7). Eq (7) is a non-linear relation between γ and χ, so small changes in γ could produce large changes in χ.
For Kolmogorov turbulence, E $ k −5/3 , the new theory predicts that, γ > 4/3, and χ > 3. For turbulence with intermittency μ I > 0, such that E $ k −(5/3+μ I ) , the scaling is again greater than from a purely local theory, It is interesting to note that even under the classical locality assumption, in real turbulence with intermittency p = 1.72 we should obtain the scaling γ l = 1.36, and χ l = 3.125; thus, the https://doi.org/10.1371/journal.pone.0202940.g012 Turbulent particle pair diffusion classical RO-t 3 regime does not actually exist! This is important for another reason because it gives us a way of testing how close current DNS is to at the locality limit-see Discussion, Section 6.
We define M γ (p) to be the ratio of the scaling power γ(p) to and local scaling powers γ l (p), M γ (p) is equal to 1 at both p = 1 and p = 3, and since M γ > 1 in the range 1 < p < 3, then there must be a maximum in M γ at some p = p m for an energy spectrum E $ k −pm , where 1 < p m < 3.
Richardson's 1926 dataset, Fig 1, from real geophysical turbulence (i.e. including intermittency) suggests a scaling of, g m I % 1:564, i.e. K m I $ s 1:564 l . The theory developed here in Eq (46) is more consistent with this data than previous theories. However, not all large-scale measurements agree with this, as discussed in Section 2.
The preditions for asymptotically infinite inertial subrange R k ! 1 is summarised in the sketch in Fig 13 which shows the predicted log-log plots of K against σ l as p increases from 1 to 3.
Finally, we remark that corrections from the ends of the inertial subrange may modify some of the regimes predicted above. Ultra-violet corrections from the high wavenumber close to k η , and infra-red corrections from the low wavenumbers close to k 1 may penetrate some way in to the inertial subrange, but the degree of penetration is unknown. Nevertheless, for very large subranges, we still expect to observe unambigous inertial range scalings over an extended inner part of the subrange.

Exposing the local process-Short inertial subrange, R k ( 1 (Re ( 1)
Eq (36) means that turbulent pair diffusion is governed by two broadly independent local and non-local processes. This excites a compelling question; are there circumstances in which we can 'turn off' either one of the processes thus exposing the other process explicitly?
The non-local process, which is the first term in Eq (36), exists in the wavenumber range ½k 1 ; k Ã l ; but suppose that the inertial subrange itself is so small that it is close to the size of the locality kernel? Then the non-local range of scales would be absent to good approximation, and we might then observe a purely local diffusional process, as illustrated in the sketch in Fig 14. However, due to the small size of the inertial subrange, the infra-red and ultra-violet corrections from the ends of the inertial subrange may have a relatively greater influence than in very large inertial subranges. So such a regime may only be approximately local in nature-i.e. a quasi-local regime, As we progressively increase the size of the inertial subrange we would expect to see a smooth transition from the locality regime at moderate inertial subrange to the non-locality regime at very large (effectively infinite) subrange, as illustrated in the sketch in Fig 15. Such a quasi-local regime if it exists is non-Richardson in character because locality was hypothesised for strictly infinite inertial subranges by Richardson.

5.4
Exposing the non-local process-Very small initial separation l 0 ( η We can 'turn off' the local process in Eq (36) by simply removing the 'local' part of the specturm-this is equivalent to taking a very small initial separation l 0 ( η. Then there is a spectral gap between the k 0 = 1/l 0 and k η ( k 0 where E(k) = 0. If this gap is large enough then the spectrum between k 1 k k η will in efffect be non-local to the pair separation process so long as σ l (t) ( η. Strictly speaking, this is not inertial range scaling any more because the separation is outside the inertial subrange, but because this system can be used to test the fundamental premise of the new theory, we will consider it here. According to the theory described in Eq (36), with the local range removed we should now observe purely strain dominated pair separation which has the diffusion coefficient scaling, for all p. In fact, this scaling is independent of the size of the inertial subrange and also of the form of the spectrum, so long as σ l ( η. Eq (50) is equivalent to exponential growth in time, σ l (t) $ l 0 exp(St), where S depends upon the form of E(k). But as σ l approaches η from below, i.e. σ l /η ! 1, the inertial subrange scaling will begin to take effect and we expect the scaling in (50) to change over.

Discussion
Richardson's hypothesis of locality was based on the 1926 data-set of pair diffusion coefficients. However, the reappraised 1926 data presented here shows an unequivocal non-local scaling for the turbulent pair diffusivity, K $ s 1:564 l , Fig 1. Consequently, the foundations of turbulent pair diffusion theory have been re-examined here in an effort to resolve one of the most important and enduring problems in turbulence.
A novel mathematical approach has been developed by expressing the pair diffusion coeffiicient through a Fourier integral decomposition. a priori assumptions regarding locality have not been made, and this has led to an expression for K as the sum of local and non-local contributions in Eq (36).
The main contribution of this investigation is to propose a new local-non-local theory of turbulent particle pair diffusion based upon the principle that the turbulent pair diffusion process in statistically stationary homogeneous turbulence is governed by both local and nonlocal diffusional processes. The theory preditcs the existance of two non-Richardson pair diffusion regimes in turbulence with generalized energy spectra, E(k) $ k −p .
For asymptotically infinite inertial subrange the pair diffusion coefficient scales like, Kðl; pÞ $ s gðpÞ l , with (1 + p)/2 < γ(p) 2, in the range 1 < p 3, Eq (46), which is intermediate between the purely local and purely non-local scaling power laws. The reappraised 1926 geophysical data, Fig 1, provides support for the new theory in this limit.
For short inertial subrange quasi-local regimes are predicted in which the pair diffusion coefficient scales approximately locally like, Kðl; pÞ $ s ð1þpÞ=2 l , in the range 1 < p 3. The theory also predicts that the existence of the local and non-local diffusional processes can be demonstrated in principle by islolating each one; first by taking short inertial subranges which would expose the local process, and then by taking very small initial separations which would expose the non-local process.
In summary, the new theory predicts the following two non-Richardson pair diffusion regimes: local : KðlÞ ¼ s g l ðpÞ l ; g l ðpÞ ¼ ð1 þ pÞ=2; R k =C % 1 ðRe ( 1Þ ð51Þ non À local : KðlÞ ¼ s gðpÞ l ; g l ðpÞ < gðpÞ 2; R k =C ! 1 ðRe ! 1Þ ð52Þ The new theory could explain why DNS studies in pair diffusion have failed to observe non-local regimes of pair diffusion-because DNS would have to generate inertial subranges orders of magnitude bigger than currently possible. It was shown in Section 5.2 that in the locality limit, we should obtain the scalings KðlÞ $ s 1:36 l and hl 2 i $ t 3.125 if we assume an intermittent spectrum of p = 1.72. However, no DNS study has so far reported a scaling greater than hl 2 i $ t 3 , which indicates that the size of the inertial subrange produced in DNS has not reached the minimum needed to observe locality scaling at the current time. It is important to note that the theory presented here is based upon fundamental physical principles, and as such it is not dependent on any specific numerical method of simulation. However, the theory cannot predict quantitatively the power laws γ(p) themselves, or the precise size of the inertial subrange that will yield either of the two non-Richardson regimes. For these quantities, and in order to investigate other predictions of the new theory, ideally we would need experiments or DNS.
Neither laboratory experiments nor DNS can produce extended inertial subranges at the current time, which is an essential requirement in the present theory. However, some Lagrangian diffusion models can generate large inertial subranges. In a companaion paper, [27], such a numerical simulation method is used to investigate the new theory where all predictions of the new theory have been verified.