2.1. A $k^2$ spectrum requires strong long-range correlations

The energy spectrum is the Fourier transform of the two-point velocity correlation function, $\langle \boldsymbol {u} (\boldsymbol {x}) \boldsymbol {\cdot } \boldsymbol {u} (\boldsymbol {x} + \boldsymbol {r}) \rangle \equiv \langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$, where angle brackets indicate an ensemble average. For statistically homogeneous and isotropic turbulence, $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ is a function of $r= |\boldsymbol {r}|$ only, and then the energy spectrum is

(2.1)\begin{equation} \mathcal{E}(k) = \frac{k^2}{4{\rm \pi}^2} \int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle \, {\rm e}^{-{\rm i}\boldsymbol{k} \,\boldsymbol{\cdot}\, \boldsymbol{r}} = \frac{1}{\rm \pi} \int^\infty_0 \mathrm{d} r \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle kr \sin (kr).\end{equation}

If correlations between points separated by distances much greater than the energy-containing scale of the turbulence, $l$, decay sufficiently quickly, then (2.1) may be Taylor expanded for small $k$. Namely, if $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle = o(r^{-5})$ as $r \to \infty$, then (1.1) holds.

From (1.1), it would appear that the ‘thermal’ $k^2$ spectrum corresponds to $L\neq 0$. However, this conclusion is problematic, because $L$ is an invariant. This fact is well known in the context of decaying turbulence, for which the conservation of $L$ implies a meaningful distinction between turbulence with finite $L$, called ‘Saffman turbulence’, and that with $L=0$, called ‘Batchelor turbulence’. These two canonical types of turbulence have a number of differences, chief among them their laws for the decay of energy with time (see Davidson (Reference Davidson2015) for a review). As we shall show in § 2.2, conservation of $L$ should also be expected in forced turbulence, provided that long-range correlations in the forcing function are sufficiently weak to prohibit injection of $L$. As a result, if $L=0$ at $t=0$, $L=0$ at all subsequent times.

The relevance of correlations in the forcing function is that sufficiently strong long-range correlations in the velocity field are required for $L$ to be non-zero. Statistical isotropy and homogeneity, together with incompressibility, restrict the allowed form of the two-point velocity correlation tensor $\langle u_i u_j'\rangle$ to

(2.2)\begin{equation} \langle u_i u_j'\rangle = \frac{u^2}{2r}[(r^2 \chi)' \delta_{ij} - \chi'(r) r_i r_j], \end{equation}

where $\chi (r) = \langle u_r (\boldsymbol {x}) u_r (\boldsymbol {x}+\boldsymbol {r}) \rangle /u^2$ is the longitudinal correlation function (see, e.g. Landau & Lifshitz Reference Landau and Lifshitz1959; Davidson Reference Davidson2015), and we follow the convention ${u^2 \equiv \langle u_x^2 \rangle = \langle |\boldsymbol {u}|^2\rangle }$/3. Equation (2.2) implies

(2.3)\begin{equation} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle = \frac{1}{r^2}\frac{\partial}{\partial r} (r^3 u^2 \chi).\end{equation}

Integrating (2.3) over all space, we find that the Saffman integral, (1.2), is

(2.4)\begin{equation} L = 4{\rm \pi} u^2 \lim_{r\to \infty} r^3 \chi (r). \end{equation}

Thus, $L$ is finite if and only if

(2.5)\begin{equation} \chi(r\to\infty) = O(r^{{-}3}).\end{equation}

Note that, somewhat counter-intuitively, (2.5) need not mean that $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ decays slowly with $r$, as may be shown by substituting (2.5) in (2.3). As a consequence, the long-range correlations implied by (2.5) do not necessarily invalidate the expansion (1.1), which required $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle = o(r^{-5})$. An extreme example is a white-noise velocity field

(2.6)\begin{equation} \langle \boldsymbol{u} (\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x} +\boldsymbol{r}) \rangle \propto \delta^3 (\boldsymbol{r}) \implies \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle (r) \propto \frac{\delta (r)}{r^2}. \end{equation}

It follows immediately from (1.2) and (2.6) that $L\neq 0$ for such a field (in this case, the $k^2$ spectrum extends to all scales). However, we see from (2.3) that $\chi (r)=1/r^3$, and so, from (2.2), $\langle u_i u_j'\rangle = 3u^2 r_i r_j/2r^5$ for $i\neq j$. This means that even a white-noise velocity field, if incompressible, must have long-range correlations hidden in the off-diagonal components of its spectral tensor.

2.2. Non-local fluid processes are insufficient to generate long-range correlations

Intuitively, no local (in real space) forcing mechanism can set up correlations between infinitely separated points, at least in the absence of non-local fluid processes. Of course, this need not be an obstacle to the development of a non-zero Saffman integral, and hence a thermal-equilibrium $k^2$ spectrum, because incompressible turbulence is subject to non-local interactions: physically, incompressibility is enforced via the action of pressure waves, which propagate at infinite speed through the fluid. Let us estimate the strength of correlations that can develop as a consequence of the pressure-mediated interactions.

Taking the divergence of the Navier–Stokes equation

(2.7)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla} p + \nu \nabla^2 \boldsymbol{u} + \boldsymbol{f},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0,\end{equation}

and assuming that $\boldsymbol {f}$ is solenoidal (i.e. that $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {f} = 0$ – we shall return to the case of $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {f} \neq 0$ in § 4.3. The reader may wonder why this distinction is necessary. After all, only the solenoidal part of $\boldsymbol {f}$ is dynamically significant; the compressive part is negated by the pressure in an incompressible fluid. The problem is that when $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {f} \neq 0$, the solenoidal part of $\boldsymbol {f}$ is not necessarily local in real space, even if $\boldsymbol {f}$ is. Remarkably, we shall find in § 4.3 that when $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {f}\neq 0$, locally forced turbulence does not generically tend to equilibrate towards $\mathcal {E}(k) \propto k^2$ at large scales), we have

(2.8)\begin{equation} \nabla^2 p ={-} \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}) , \end{equation}

so the pressure is always exactly what is required to negate the non-solenoidal part of the inertial force. Inverting (2.8), we find that the far-field pressure generated by an eddy localised at $\boldsymbol {x}=0$ in an otherwise quiescent fluid is

(2.9)\begin{align} p(\boldsymbol{x}) & = \frac{1}{4{\rm \pi}} \int \frac{\mathrm{d}^3 \boldsymbol{x}'}{|\boldsymbol{x}' -\boldsymbol{x}|} \frac{\partial}{\partial x_i'}\frac{\partial}{ \partial x_j'}u_i(\boldsymbol{x}') u_j(\boldsymbol{x}') \end{align}

(2.10)\begin{align} & =\frac{1}{4{\rm \pi}} \frac{\partial}{\partial x_i}\frac{\partial}{ \partial x_j}\frac{1}{x}\int \mathrm{d}^3 \boldsymbol{x}' u_i(\boldsymbol{x}')u_j(\boldsymbol{x}') + O(x^{{-}4}) \end{align}

(2.11)\begin{align} & =O(x^{{-}3}), \end{align}

where we have Taylor expanded the Green's function $|\boldsymbol {x}'-\boldsymbol {x}|^{-1}$ in (2.11) in small $|\boldsymbol {x}'|$. Thus, a localised eddy generates a pressure field that extends to arbitrarily large distances, falling off as $x^{-3}$, with the corresponding $\boldsymbol {\nabla } p$ force decaying as $x^{-4}$.

Informally, we can imagine constructing a homogeneous and isotropic turbulence as a random assembly of many such eddies. Each would exert a force on distant ones that scales with their separation $r$ as $r^{-4}$. The strength of correlations that would develop due to these pressure-mediated interactions may be estimated using the von Kármán–Howarth equation (von Kármán & Howarth Reference von Kármán and Howarth1938), which follows from (2.7) under the assumptions of statistical isotropy and homogeneity,

(2.12)\begin{equation} \frac{\partial}{\partial t} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle = \frac{1}{r^2} \frac{\partial}{\partial r} \frac{1}{r} \frac{\partial}{\partial r} (r^4 u^3 K) + 2 \nu \nabla^2 \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle + 2 \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle, \end{equation}

where ${K(r) = \langle u_r(\boldsymbol {x}) u_r(\boldsymbol {x}) u_r(\boldsymbol {x}+\boldsymbol {r})\rangle / u^{3/2}}$ is the longitudinal triple-correlation function. Pressure does not appear in (2.12) – this is a consequence of statistical homogeneity, which demands that $\partial \langle u_i p'\rangle /\partial r_j = 0$. Instead, pressure enters implicitly via the coupling of (2.12) to higher-order correlators, i.e. the term containing $K(r)$. The analogue of (2.12) for triple correlations is

(2.13)\begin{align} &\frac{\partial }{\partial t} \langle u_i u_j u'_k\rangle= \frac{\partial}{\partial r_l}\langle u_i u_j u'_k u_l\rangle - \frac{\partial}{\partial r_l}\langle u_i u_j u'_k u'_l\rangle - \left\langle u_i u_j \frac{\partial p'}{\partial x'_k}\right\rangle - \left\langle u'_k\left( u_i \frac{\partial p}{\partial x_j}+u_j\frac{\partial p}{\partial x_i}\right)\right\rangle \nonumber\\ &\quad + \mathrm{viscous\ terms} +\langle u_i u_j f'_k \rangle+\langle (u_i f_j+f_i u_j) u'_k \rangle, \end{align}

where the terms involving $p$ do not vanish. In particular, for our ensemble of randomly distributed eddies, the correlator $\langle u_i u_j \partial _k' p'\rangle$ is $O(r^{-4})$ as $r\to \infty$, because the part of $p'=p(\boldsymbol {x}+\boldsymbol {r})$ that is correlated with the velocity field at position $\boldsymbol {x}$ decays like $r^{-3}$. This suggests that

(2.14)\begin{equation} K(r\to \infty) = O(r^{{-}4}). \end{equation}

The above argument for the scaling (2.14) is informal: there is an obvious inconsistency in evaluating $\langle u_i u_j \partial _k' p'\rangle$ by assuming that different patches of the turbulence are uncorrelated with the conclusion that correlations $K(r\to \infty ) = O(r^{-4})$ must develop. However, the argument can be formalised – we prove in Appendix A that (2.14) holds for real turbulence provided spatial correlations in the forcing decay sufficiently quickly (viz., exponentially) under suitable assumptions (this proof is a generalisation to forced turbulence of arguments presented by Batchelor & Proudman (Reference Batchelor and Proudman1956) and Saffman (Reference Saffman1967) for decaying turbulence). Importantly, (2.14) is too weak a correlation to permit $L\neq 0$: integrating (2.12) over all $\boldsymbol {r}$, we find

(2.15)\begin{equation} \frac{\mathrm{d} L}{\mathrm{d} t} = 4{\rm \pi} \lim_{r\to\infty}\left[\frac{1}{r} \frac{\partial}{\partial r} (r^4 u^3 K)\right] +2 \int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle, \end{equation}

where the term involving $K$ vanishes for $K(r\to \infty ) = O(r^{-4})$.

2.3. Correlations generated directly by the forcing

The argument in § 2.2 indicates that non-local interactions between fluid elements are too weak to allow the Saffman integral to change with time. Another concern is that correlations in the forcing function itself might decay sufficiently slowly with distance to permit development of $L\neq 0$; this effect is encoded in the second term on the right-hand side of (2.15). We can estimate how slowly these correlations need to decay for $\mathrm {d} L/\mathrm {d} t$ to be non-zero by examining the formal solution for $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {f}' \rangle$ that is obtained by integrating the Navier–Stokes equation in time:

(2.16)\begin{align} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle &= \int^t_0 \mathrm{d} s \left[ \frac{\partial}{\partial r_j}\langle u_i(s) u_j(s) f'_i(t) \rangle - \frac{\partial}{\partial r_i} \langle p(s) f_i'(t) \rangle \right. \nonumber\\ & \quad + \left.\vphantom{\frac{\partial}{\partial r_j}} \nu \nabla^2 \langle \boldsymbol{u}(s) \boldsymbol{\cdot} \boldsymbol{f}'(t)\rangle + \langle \boldsymbol{f}(s) \boldsymbol{\cdot} \boldsymbol{f}'(t) \rangle\right]. \end{align}

Of the four terms in the right-hand side of (2.16), the first and third give rise to surface terms in (2.15), which vanish provided the relevant correlators fall off faster than $r^{-2}$ and $r^{-1}$, respectively (a proof that they do, under the assumption that correlations in the forcing decay exponentially with distance, is presented in Appendix A). The second term is identically zero by the solenoidality of $\boldsymbol {f}$. The fourth term is a two-point, two-time correlation function of $\boldsymbol {f}$, which, because $\boldsymbol {f}$ is a solenoidal, statistically isotropic vector field, satisfies (cf. (2.3))

(2.17)\begin{equation} \int^t_0 \mathrm{d} s \langle \boldsymbol{f}(s) \boldsymbol{\cdot} \boldsymbol{f}'(t)\rangle = \frac{1}{r^2}\frac{\partial}{\partial r} [r^3 H(t,r)], \end{equation}

where $H(t,r)$ is the time-integrated longitudinal correlation function of $\boldsymbol {f}$. From (2.15), we find that the contribution of this term to the rate of change of the Saffman integral vanishes if

(2.18)\begin{equation} H(t,r)=o(r^{{-}3}),\end{equation}

which is unsurprising, given (2.5).

The arguments presented in §§ 2.2 and 2.3 indicate that, if the forcing function is solenoidal and sufficiently localised, then correlations between infinitely separated points that are strong enough to change the Saffman integral cannot arise in finite time, even accounting for the non-local nature of the pressure force (the reader used to thinking of forcing whose properties are specified in spectral, rather than real, space, might wonder whether the condition of ‘sufficient localisation’ is satisfied for the common choice of forcing in a finite spectral band. In Appendix B, we show that a finite-band forcing with a smooth spectrum has very weak correlations at the largest scales (it decays faster with $r$ than any power law), as is intuitive, given the absence of energy in large-scale modes. If the spectrum of $\boldsymbol {f}$ is not smooth, but instead has sharp discontinuities at the edges of the band, correlations are induced in $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ that oscillate in $r$ at the wavenumbers of the edges, and decay in amplitude as $r^{-3}$. While these correlations may, in principle, propagate into all other correlators, we show in Appendix B that any oscillatory component of $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ always has a negligible effect on $\mathcal {E}(k)$ at small $k$, so these oscillatory correlations are of little dynamical significance). Because the Saffman integral was zero at $t=0$, when $\boldsymbol {u} = 0$, it remains zero at all times, and therefore it might appear that the system is forbidden from developing a $k^2$ spectrum at $k \to 0$.

2.4. Long-range correlations as a cumulative effect of short-range interactions

How, then, does one explain the numerical evidence for a thermal-equilibrium $k^2$ spectrum in forced turbulence (Dallas et al. Reference Dallas, Fauve and Alexakis2015; Cameron et al. Reference Cameron, Alexakis and Brachet2017; Alexakis & Biferale Reference Alexakis and Biferale2018; Alexakis & Brachet Reference Alexakis and Brachet2019)? The answer is that the $k^2$ spectrum is established not by non-local processes (in real space), but by the cumulative effect of local ones. Then, while infinitely separated points can never be strongly correlated enough to induce a $k^2$ spectrum, points separated by a large but finite distance can be (as long as one is prepared to wait long enough), leading to a $k^2$ spectrum that spans a finite, time-dependent range of wavenumbers.

Let us now demonstrate that this is indeed the typical behaviour of forced turbulence, by means of a numerical simulation. We take the forcing function $\boldsymbol {f}$ to be a solenoidal, Gaussian random field (as is a common choice for forced-turbulence studies), and to be delta correlated in time, so the spectrum of energy injection is

(2.19)\begin{equation} F(k) \equiv \frac{k^2}{2{\rm \pi}^2} \int_0^t \mathrm{d} s\int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{f}(t) \boldsymbol{\cdot} \boldsymbol{f}'(s)\rangle \, {\rm e}^{-{\rm i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}} \propto k^4 \exp({-}k^2/k_p^2), \end{equation}

where the peak wavenumber $k_p=80$. Because the power injected into the $k=0$ mode is always zero, the average of the velocity (momentum) over the periodic box is zero for all $t$. The large-scale $k^4$ tail of $F(k)$ is consistent with the generic spectral tail of an isotropic field with short spatial correlations (an expansion of $F(k\to 0)$ analogous to (1.1) yields $F(k\to 0)\propto k^4$ if $\langle \boldsymbol {f} \boldsymbol {\cdot } \boldsymbol {f}'\rangle$ decays rapidly with $r$. A faster decay of $F(k\to 0)$ would require the equivalent of the Loitsyansky integral (1.3) for $\boldsymbol {f}$, $I_{\boldsymbol {f}}\equiv -\int ^t_0\mathrm {d} s \int \mathrm {d}^3 \boldsymbol {r} r^2 \langle \boldsymbol {f}(t) \boldsymbol {\cdot } \boldsymbol {f}'(s)\rangle$, to be zero, which is an artificial situation), although this choice is not essential to observe the development of a $k^2$ band – other studies have used finite-band forcing (Dallas et al. Reference Dallas, Fauve and Alexakis2015; Cameron et al. Reference Cameron, Alexakis and Brachet2017; Alexakis & Brachet Reference Alexakis and Brachet2019). The algorithm that we employ to generate $\boldsymbol {f}$ is described in Appendix D of Hosking & Schekochihin (Reference Hosking and Schekochihin2021). With this choice, we solve the Navier–Stokes equations (2.7) in a periodic domain of size $2{\rm \pi}$ using the pseudo-spectral code Snoopy (Lesur Reference Lesur2015) with $512^3$ resolution. We employ de-aliasing according to the $2/3$-rule, and use eighth-order hyper-dissipation, i.e. $\nu \nabla ^2$ is replaced by $\nu _8 \nabla ^8$ in (2.7), where $\nu _8 = 10^{-16}$. The use of hyper-dissipation ensures that the effect of viscosity on the development of the large-scale structure is negligible. The simulation time step ${\rm \Delta} t$ is chosen automatically by the code so as to be sufficiently small to maintain the stability of the simulation according to the Courant–Friedrichs–Lewy (CFL) criterion. The spectral scheme allows the viscous term to be solved exactly at each time step.

The results of this simulation are shown in figure 1. We observe that $\boldsymbol {u}$ gradually develops a $k^2$ spectrum at large scales, with a spectral knee at a time-dependent wavenumber $k_c(t)$ separating the $\propto k^4$ and $\propto k^2$ parts, as anticipated. By fitting the numerical spectrum to a trial function of the form $k^2[1-\exp (-k^2/k_c(t)^2)]$, we find that $k_c(t)\propto t^{-1/2}$ (see inset to figure 1). At large enough times, the $k^2$ part of the spectrum extends all the way to the box size, which is the steady state (close to the box scale, i.e. say, at $k<4$, we observe some deviation from the $\propto k^2$ scaling at late times, which presumably is due to the absence of statistical isotropy at these scales). The ability of the system to reach a steady state hinges on the finite size of the simulation box – in an infinite system, $\mathcal {E}(k)\propto k^2$ would only ever be satisfied in an ever-broadening but finite band of wavenumbers.

Figure 1. Saturation of the large scales in simulated Navier–Stokes turbulence forced by a delta-correlated, Gaussian random field with weak long-range spatial correlations (so that $F(k\to 0\propto k^4)$, as explained in the text). Displayed spectra are logarithmically spaced in time, with blue $\to$ red indicating earlier $\to$ later times. The inset shows the evolution of the knee wavenumber, $k_c(t)$, that separates the $\propto k^4$ and $\propto k^2$ parts of the spectrum. In the chosen units, the energy-injection rate is $0.7$, and the root-mean-square (r.m.s.) velocity is $\simeq 0.5$.

We note that the turbulence in this simulation is not fully developed – we sacrifice the resolution of the $k^{-5/3}$ inertial range to facilitate resolving many forced structures, so that we may invoke statistical isotropy and homogeneity in our analysis, and also to allow a wide-band forcing function, so as to eliminate spurious effects that occur when forcing is restricted to a narrow spectral band (see § 4.4). We do not expect that the development of the $k^2$ band is a consequence of our failure to resolve the inertial range, as the small-$k$ spectral asymptotic is determined by the statistical properties of outer-scale structures (via (1.1)). This view is supported by the study of Alexakis & Brachet (Reference Alexakis and Brachet2019), which presents simulations of a turbulence that appears closer to being fully developed than ours (achieved by forcing in a narrow spectral band) but still develops the thermal spectrum. We do not expect that the use or order of hyperdiffusion affects the process of thermalisation, for the same reason.

To summarise our progress so far, we have seen that the law of conservation of the Saffman integral, ported from the theory of decaying turbulence, also holds for forced turbulence, and that this law prohibits the thermal equilibration of arbitrarily large scales in finite time. Nonetheless, thermal equilibration up to a large but finite scale is not prohibited, and indeed this is the behaviour that turbulence whose forcing is spatially decorrelated tends to adopt (as is shown by figure 1). However, we still lack a physical mechanism for the development of the thermal spectrum. In the next section, and the one that follows it, we shall argue that this mechanism is turbulent diffusion of linear momentum.

3.1. Broken-power-law spectra and their momentum content

The above discussion suggests that we might interpret the growth of a $k^2$ spectrum over a finite range of wavenumbers as indicating the development of random fluctuations in momentum that satisfy $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ over the corresponding range of scales. More precisely, these fluctuations would be quasi-random, in that the momenta of the eddies contained within $V$ would cancel more precisely when $R$ was large enough, so that $\langle \boldsymbol {P}_V^2 \rangle$ would be dominated by eddies at the surface of $V$, so that $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$. Then, ${\mathcal {E}(k\to 0)\propto k^4}$. A schematic representation of the distribution of momentum of this ‘quasi-Saffman turbulence’, similar to those presented by Davidson (Reference Davidson2015) for the true Saffman and Batchelor turbulence, is shown in figure 2.

Figure 2. Schematic of a ‘quasi-random’ distribution of linear momentum, i.e. one that would result in a broken-power-law spectrum, as in figure 1. Sufficiently large patches of turbulence have vanishing total momentum – a number of such patches (identified in a non-unique manner) are shown in different colours in the figure. For a control volume $V$ that is larger than the outer scale of the turbulence but smaller than the characteristic scale of the net-zero-momentum patches (e.g. the smaller circle in the figure), $\langle \boldsymbol {P}^2\rangle \propto R^3$ because the eddies contained by $V$ (represented by individual blobs) have uncorrelated, random momenta (represented by arrows). On the other hand, $\langle \boldsymbol {P}^2\rangle \propto R^2$ for $V$ much larger than the zero-net-momentum patches, because then only patches at the surface of $V$ contribute to the sum – in the figure, the central orange and yellow patches do not contribute to the total momentum contained within the larger circle.

Let us now check that these intuitive expectations hold up mathematically, i.e. that broken-power-law spectra do correspond to broken power laws in the dependence of $\langle \boldsymbol {P}_V^2\rangle$ on $R$. From (2.3) and

(3.6)\begin{equation} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle = 2 \int^{\infty}_0 \mathrm{d} k \mathcal{E}(k) \frac{\sin (kr)}{kr},\end{equation}

which is the inverse of (2.1), it follows that

(3.7)\begin{equation} u^2 \chi(r) = 2 \int^{\infty}_0\mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3}.\end{equation}

In Appendix C, we present a formal asymptotic expansion of (3.7), assuming that ${\mathcal {E}(k)\propto k^a}$ for ${k_1 \leq k \leq k_2}$ with ${k_2 \gg k_1}$ (we remain agnostic about $\mathcal {E}(k)$ outside of this range). This choice for ${\mathcal {E}(k)}$ models the broken-power-law spectrum shown in figure 1. We show that, for ${1/k_2 \ll r\ll 1/k_1}$

(3.8) \begin{equation} u^2 \chi(r) \simeq \begin{cases} \displaystyle \mathrm{constant} \sim \int^{k_1}_{0} \mathrm{d} k \mathcal{E}(k) & \text{if }a<{-}1; \\ \displaystyle \frac{1}{3} \frac{\ln(k_1 r)}{\ln (k_2/k_1)} \int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) & \text{if } a={-}1; \\ \displaystyle \mathrm{undetermined,} \lesssim (k_2 r)^{{-}1-a} \int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) & \text{if }a = 4,6,8\ldots;\\ \displaystyle -\varGamma(a-2)(a^2-1)\sin\left(\frac{a{\rm \pi}}{2}\right) (k_2 r)^{{-}1-a}\int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) & \text{otherwise.} \end{cases}\end{equation}

Let us explain qualitatively each case in turn.

If $a<-1$, $\chi (r)\sim \mathrm {const.}$, which is intuitive: the energy contained in the band $\{k_1, k_2\}$ is dominated by the largest structures, while we are looking at correlations on scales much smaller than them ($r\ll k_1^{-1}$).

If $a=-1$, then every scale in the band $\{k_1,k_2\}$ contributes equally to the energy contained within it – this energy diverges in the limit $k_2/k_1\to \infty$, explaining the factor of $\ln (k_2/k_1)$ in the denominator of (3.8). It turns out that the $r$ dependence of $\chi (r)$ is logarithmic in this case.

If $a = 4,6,8\ldots$, then although $\chi (r)$ must decay faster than $r^{-1-a}$ in the range ${1/k_2 \ll r\ll 1/k_1}$, its behaviour is not uniquely determined by our assumption of a power-law scaling for $\mathcal {E}(k)$. This was to be expected, because even-power spectra are precisely the ones generated in the expansion (1.1), and no specific strength of long-range correlations in the velocity field is required for the coefficients of $k^a$ with $a = 4,6,8\ldots$ in this expansion to be non-zero (unlike for $a=2$).

Finally, for all other cases, including that of $a=2$, $\chi (r)$ decays like $r^{-1-a}$ for ${1/k_2 \ll r\ll 1/k_1}$. In particular, note that (2.5) may be recovered from (3.8) for $a=2$, as ${\lim _{a\to 2}\varGamma (a-2)\sin (a{\rm \pi} /2)=-{\rm \pi} /2}$.

Our motivation in deriving (3.8) was to obtain the dependence of $\langle \boldsymbol {P}^2_V\rangle$ on $R$ that is associated with a finite-extent large-scale power law in $\mathcal {E}(k)$. Let us consider scales $1/k_2 \ll R \ll 1/k_1$, where $k_2$ is now identified with the outer scale of the turbulence, i.e. $k_2\sim 1/l$, and $k_1$ is identified with the scale of the spectral knee $k_c$ in figure (1). Then, from (3.4),

(3.9)\begin{equation} \langle \boldsymbol{P}_V^2 \rangle = 2{\rm \pi}^2 u^2 \int^{2R}_0 \mathrm{d} r r \left[\int_0^{X/k_2} \,\mathrm{d} r' r'^3 \chi(r')+ \int_{X/k_2}^r \,\mathrm{d} r' r'^3 \chi(r')\right],\end{equation}

where $X$ is chosen so that $1 \ll X \ll k_2/k_1$, in which case (3.8) is applicable in the second integral appearing inside the brackets in (3.9). This integral dominates over the first one for $r\gg X/k_2$ as long as $r^3 \chi (r)\geq O(1/r)$, which, according to (3.8), it will do if the spectrum follows a local power law with exponent $a \leq 3$. Otherwise, the first integral, which is independent of $r$, dominates. Thus, we have

(3.10)\begin{equation} \langle \boldsymbol{P}_V^2 \rangle \propto \begin{cases} R^2 & \text{if }a>3, \\ R^2 \ln R & \text{if } a=3, \\ R^{5-a} & \text{if } -1< a<3, \\ R^{6}\ln R & \text{if } a={-}1, \\ R^{6} & \text{if } a<{-}1. \\ \end{cases}\end{equation}

We note that the classical scalings (see Davidson Reference Davidson2015) are readily recoverable from (3.10): the intuitive ‘surface-term-dominated’ $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$ is recovered for steep spectral slopes, ${a>3}$, corresponding to weak long-range correlations, while the Saffman scaling ${\langle \boldsymbol {P}_V^2 \rangle \propto R^3}$ is recovered for $a=2$. The scaling $\langle \boldsymbol {P}_V^2 \rangle \propto R^6$ for $a<-1$ is also an intuitive one: such a spectrum is energetically dominated by structures with characteristic scale much larger than $R$, therefore control volumes $V\propto R^3$ will contain a total amount of momentum that is proportional to $V$. Although these results are familiar, (3.10) has the important new feature that it does not require the spectral power law to extend all the way to $k=0$ – it is sufficient for $\mathcal {E}(k)\propto k^a$ only for $k_1 \leq k\leq k_2$, as long as we restrict attention to volumes with $1/k_2 \ll R\ll 1/k_1$.

As an aside, we remark that, besides the generalisation of previous results to a finite-band power law in $\mathcal {E}(k)$, the other qualitatively new feature in (3.10) is that we have allowed for non-integer $a$. In this respect, (3.8) and (3.10) can be viewed as an extension of the results for $\chi (r\to \infty )\propto r^{-m}$ for integer $m$ derived by Davidson (Reference Davidson2011). While non-integer large-scale spectral power laws may be of limited applicability to real turbulence (though they can, of course, be manufactured numerically), they nonetheless have pedagogical value, particularly for $3< a<4$, when $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$, implying that arguments for the invariance of the large-scale spectrum in decaying turbulence that are based on momentum conservation (Saffman (Reference Saffman1967), Davidson (Reference Davidson2011); also see § 5) do not apply. If it is true that the large-scale asymptotic of the energy spectrum is indeed invariant in decaying turbulence with $3< a<4$, then this must be a result of the conservation of some other quantity. The arguments presented in Davidson (Reference Davidson2009, Reference Davidson2011) suggest that angular-momentum conservation, if it holds, would result in the invariance of this asymptotic; direct numerical simulations of turbulence with $3< a<4$ might therefore shed some light on the unsolved problem of angular-momentum conservation in turbulence in open domains. Interestingly, large-scale spectra with $3< a<4$ are not invariant under the Eddy Damped Quasi-Normal Markovian (EDQNM) closure model, whereas those with $a<3$ are (Eyink & Thomson Reference Eyink and Thomson2000; Lesieur, Métais & Comte Reference Lesieur, Métais and Comte2005; Lesieur Reference Lesieur2008).

3.2. The development of ‘quasi-random’ momentum fluctuations

Having confirmed that broken-power-law spectra, of the form shown in figure 1, do correspond to $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ over a finite range of scales, we now turn to the question of the physical mechanism that is responsible for the development of such a scaling.

Intuitively, $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ can arise as a simple consequence of momentum transport by the flow. Consider a localised fluid motion that develops at $t=0$ as a result of the forcing. While the total linear momentum associated with this structure will be zero, its momentum density will be transported under the action of the flow (i.e. the sum of the eddy's own motion and that of the rest of the flow), and therefore will become distributed over an ever-increasing volume as time advances. When this volume becomes large compared with the control volume $V$ for which we are interested in computing the total square momentum, this structure will contribute to $\boldsymbol {P}_V$ as a ‘volume term’, rather than as a surface one. The occurrence of this process at all points in space will then lead to a ‘quasi-random’ momentum distribution, of the form depicted in figure 2.

In figure 3, we present a simple toy model to illustrate this idea. In this model, turbulent eddies initialised by the forcing at $t=0$ are represented by pairs of particles. Each particle in the pair is initialised with the same random position in two-dimensional space, although they have opposite momenta – see panel (a). This means that at $t=0$, $\langle \boldsymbol {P}_V^2 \rangle \propto R$, because only ‘eddies’ at the boundary of $V$ contribute, much as in a real forced turbulence that has just reached saturation at the forcing scale (note that, naturally, the surface-dominated and volume-dominated scalings are different in two dimensions). Subsequently, the particles move ballistically, i.e. without interacting, all at the same speed, $u$, but in random directions. At later times, the distribution of their momenta becomes quasi-random, in the sense described above. For $R \ll ut\equiv R_c(t)$, $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$, because the control volume $V$ will only contain one particle from each pair. For $R\gtrsim R_c(t)$, $\langle \boldsymbol {P}_V^2 \rangle \propto R$, as the volume will contain both particles, whose contributions will cancel, unless they straddle the boundary.

Figure 3. A toy model to illustrate quasi-randomisation of eddy momentum. Eddies are represented by pairs of particles that are initialised with equal and opposite momenta, but at the same position in space, shown as red and blue arrows in panel (a). Panel (b) shows the state of the system at $t = l_{\textrm{box}}/5u$. Panel (c) shows the evolution of $\langle \boldsymbol {P}_V^2 \rangle$ with time (blue = early, red = late), demonstrating the development of the ‘stochastic’ momentum scaling, $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$, as explained in the text. Panel (d) shows that the position $R_c(t)$ of the ‘knee’ in the scaling behaviour of $\langle \boldsymbol {P}_V^2 \rangle$, between $\propto R$ and $\propto R^2$, grows linearly with time, $R_c = ut$.

While ballistic streaming is likely a poor model of the real motions of turbulent fluid structures, and while this toy model also neglects the effect of continuous forcing and energy dissipation, it captures the essential idea: transport of linear momentum means that highly ordered states where the total momentum of each flow structure vanishes cannot be maintained. Instead, it seems inevitable that the distribution of momentum will become quasi-random, i.e. that $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ (in three dimensions) up to a finite $R=R_c(t)$, which will grow with time.

Intuitively, ballistic streaming represents the upper bound permitted by causality on the rate at which the distribution of momentum can become stochasticised in the absence of significant non-local interactions. In Appendix D, we show how this causal bound, $R_c \sim ut$, can be recovered from the Navier–Stokes equation directly. In real turbulence, however, momentum is transported chaotically, rather than ballistically, so we should expect that $R_c \ll ut$. Owing to the long nonlinear time scale associated with interactions between structures on the largest scales, it is reasonable to suppose that ‘momentum density’ is transported passively by the turbulent diffusivity of the flow, at least insofar as the large scales are concerned. In that case, we expect a diffusive scaling $R_c \sim 1/k_c \propto t^{1/2}$, rather than the ballistic one, $R_c \propto t$. As we shall see in § 4, the diffusive scaling is indeed in excellent agreement with direct numerical simulations.

3.3. $R_c \propto t^{1/2}$ due to linear growth of the Loitsyanksy integral

Before we explore this topic further, however, let us pause to consider a tempting, if dangerous, argument that would appear to guarantee the scaling $R_c \sim 1/k_c \propto t^{1/2}$ without any further assumptions.

Let us accept, on the basis of the intuitive momentum-stochasticisation argument of § 3.2, that a $k^2$ spectrum will develop in a limited range of $k$, as depicted in figure 1.

It follows from integrating (2.12) over $r$ that the Loitsyansky integral (1.3) grows according to

(3.11)\begin{equation} \frac{\mathrm{d} I}{\mathrm{d} t} = 8{\rm \pi} u^3 \lim_{r\to\infty} [ r^4 K(r) ] - 12 \nu L - 2 \int \mathrm{d}^3 \boldsymbol{r} r^2 \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle. \end{equation}

Because the system always has a vanishing Saffman integral, the second term on the right-hand side of (3.11) is zero. For simplicity, let us assume that the forcing has a short correlation time, in which case (3.11) becomes, after substitution of (2.16) and (2.17),

(3.12)\begin{equation} \frac{\mathrm{d} I}{\mathrm{d} t} = 8{\rm \pi} u^3 \lim_{r\to\infty} [ r^4 K(r) ] + 8{\rm \pi} \int \mathrm{d} r r^4 H(r).\end{equation}

It is often conjectured that $I$ is conserved by an isotropic turbulence decaying from an initial state with a $k^4$ spectrum (Kolmogorov Reference Kolmogorov1941a; Landau & Lifshitz Reference Landau and Lifshitz1959; Ishida, Davidson & Kaneda Reference Ishida, Davidson and Kaneda2006; Davidson Reference Davidson2015). While this point is not universally accepted (e.g. $I$ is not conserved under the popular EDQNM closure; see Lesieur (Reference Lesieur2008) and references therein), the evidence from direct numerical simulations appears to support the conservation of $I$, at least after an initial transient period (Ishida et al. Reference Ishida, Davidson and Kaneda2006). As may be seen from (3.12), the invariance of $I$ in the absence of forcing requires that $K(r\to \infty )=o(r^{-4})$, i.e. long-range triple correlations must decay faster than the $K (r\to \infty ) = O(r^{-4})$ that follows from considering long-range pressure-mediated interactions (Batchelor & Proudman (Reference Batchelor and Proudman1956), Davidson (Reference Davidson2015); see our § 2.2). Supposing that a state with $K(r\to \infty )=o(r^{-4})$ can also arise in forced turbulence, (3.12) implies linear growth of $I$, whence, by (1.1),

(3.13)\begin{equation} \mathcal{E} (k \to 0) \propto k^4 t. \end{equation}

From our expectation that the system will saturate with a spectrum $\mathcal {E} (k) \propto k^2$, the wavenumber $k_c$ that has just saturated at time $t$ satisfies

(3.14)\begin{equation} k_c^4 t \propto k_c^2 \implies k_c \propto t^{{-}1/2},\end{equation}

which is precisely the diffusive scaling suggested above.

However, this argument should be treated with caution, because it seems unlikely that $K(r\to \infty )=o(r^{-4})$ could be realised in forced turbulence. Numerical evidence suggests that this condition is only satisfied in decaying turbulence after an initial transient period (Ishida et al. Reference Ishida, Davidson and Kaneda2006), during which the system loses memory of the initial conditions. Prior to this, growth of $I$ is observed, which requires $K(r\to \infty )=O(r^{-4})$. Forced turbulence, of course, is essentially always in this ‘transient’ regime, as the system never loses memory of the statistical properties of the forcing. Indeed, it is clear that $K(r\to \infty )=O(r^{-4})$ from the fact that $\mathcal {E}(k)\propto k^4$ does develop at the largest scales in turbulence forced with $I_{\boldsymbol {f}}=0$, which is the case, e.g. for forcing in a finite spectral band (see § 4.1).

On the other hand, the diffusive scaling (3.14) may still be obtained if ${\lim _{r\to \infty } r^4 K(r)}$ is constant in time. Admittedly, it is not a priori clear that this should be the case, because the value of this limit can change as a result of the long-range, pressure-mediated interactions between eddies, whose statistical properties do, after all, change with time as a result of the stochasticisation of linear momentum. Nonetheless, in the next section, we show that a passive model of the large-scale dynamics reproduces the linear growth (3.13) of $I$, indicating that ${\lim _{r\to \infty } r^4 K(r)}=\mathrm {const.}$ may be a reasonable approximation in real forced turbulence.

4.1. Assessing the passive-velocity model in simulated turbulence

In figure 4, we present results of simulations with both ‘active’ and ‘passive’ velocity fields. Specifically, we plot the evolution of the energy spectra of the fields $\boldsymbol {v}$ and $\boldsymbol {w}$, denoted $\mathcal {E}_{\boldsymbol {v}}(k)$ and $\mathcal {E}_{\boldsymbol {w}}(k)$, respectively, where $\boldsymbol {v}$ is determined by the forced Navier–Stokes equation,

(4.7)\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} ={-} \boldsymbol{\nabla} p_{\boldsymbol{v}} + \nu \nabla^2 \boldsymbol{v} + \boldsymbol{f}_{\boldsymbol{v}},\end{equation}

while $\boldsymbol {w}$ is governed by the passive-velocity equation (4.6), viz.

(4.8)\begin{equation} \frac{\partial \boldsymbol{w}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w} ={-} \boldsymbol{\nabla} p_{\boldsymbol{w}} + \nu \nabla^2 \boldsymbol{w} + \boldsymbol{f}_{\boldsymbol{w}},\end{equation}

where $\boldsymbol {u} = \mathcal {P}_K \boldsymbol {v}$, $\mathcal {P}_K$ is the truncation operator that removes all Fourier modes with $k< K=40$, and $\boldsymbol {f}_{\boldsymbol {v}}$ and $\boldsymbol {f}_{\boldsymbol {w}}$ are forcing functions that are delta correlated in time and inject an equal amount of energy into each Fourier mode in the band $40< k<80$ at every timestep (the box size is $2{\rm \pi}$). The other details of the simulations are as described in § 2.3.

Figure 4. Development of the ‘thermal’ $k^2$ spectrum by a Navier–Stokes velocity field, $\boldsymbol {v}$, described by (4.7), and a ‘passive velocity field’, $\boldsymbol {w}$, described by (4.8). Panel (a) shows the case where $\boldsymbol {w}$ and $\boldsymbol {v}$ are forced by the same function, $\boldsymbol {f}_{\boldsymbol {v}}=\boldsymbol {f}_{\boldsymbol {w}}$, while panel (b) shows the case where $\boldsymbol {w}$ and $\boldsymbol {v}$ are forced independently. Spectra of $\boldsymbol {w}$, $\mathcal {E}_{\boldsymbol {w}}(k)$, are plotted with dashed black lines, while spectra of $\boldsymbol {v}$, $\mathcal {E}_{\boldsymbol {v}}(k)$, are plotted with solid coloured lines: blue $\to$ red indicates earlier $\to$ later times. Panel (c) shows the evolution of the knee wavenumber $k_c(t)$ between the $k^4$ and $k^2$ parts of the spectrum. In the chosen units, the energy injection rate into each of $\boldsymbol {v}$ and $\boldsymbol {w}$ is $2.5$, and the r.m.s. values of all velocity fields are $\simeq 1.0$.

Panel (a) of figure 4 shows the results of a simulation where $\boldsymbol {f}_{\boldsymbol {v}}=\boldsymbol {f}_{\boldsymbol {w}}$. The only difference between $\boldsymbol {v}$ and $\boldsymbol {w}$ in this case is that $\boldsymbol {w}$ evolves without being advected by the modes in the large-scale tail of the spectrum of $\boldsymbol {v}$. We see that both $\boldsymbol {v}$ and $\boldsymbol {w}$ develop $k^2$ spectra at large scales gradually, with a spectral knee separating the $\propto k^4$ and $\propto k^2$ parts, as we anticipated in figure 1. The spectra $\mathcal {E}_{\boldsymbol {w}}(k)$ and $\mathcal {E}_{\boldsymbol {v}}(k)$ are almost the same at all scales and at all times. This finding demonstrates that the role of the large-scale structure of the turbulence in advecting itself and the small-scale flow is of negligible importance to the development of the thermal spectrum.

In panel (b), we plot the same spectra for a simulation where $\boldsymbol {f}_{\boldsymbol {v}}$ and $\boldsymbol {f}_{\boldsymbol {w}}$ are independent random variables. In this case, the small-scale field that advects $\boldsymbol {w}$ resembles the small-scale part of $\boldsymbol {w}$ only in a statistical sense. Nonetheless, $\boldsymbol {w}$ develops a $k^2$ band at roughly the same rate as $\boldsymbol {v}$. We view this finding as numerical justification of the passive-velocity model.

Finally, panel (c) shows the wavenumbers of the spectral knees in $\mathcal {E}_{\boldsymbol {v}}(k)$ and $\mathcal {E}_{\boldsymbol {w}}(k)$ as functions of time, computed by fitting a trial function of the form $k^2[1-\exp ({-k^2/k_c^2})]$ to the large-scale tail of the spectra. In all three cases, $k_c \propto t^{1/2}$, which is the diffusive scaling anticipated in § 3.2. In the next section, we shall derive this result (including the functional form of the knee) analytically from the passive-velocity equation (4.6).

4.2. Advection by a Kraichnan flow

In this section, we shall compute $\mathcal {E}_{\boldsymbol {w}}(k, t)$ from (4.6) analytically. The price we pay to do so is the need to make modelling assumptions about the advecting velocity field $\boldsymbol {u}$ – specifically, we take both $\boldsymbol {u}$ and $\boldsymbol {f}_{\boldsymbol {w}}$ to be delta correlated in time. Under this assumption, the $\boldsymbol {k}$-space correlation function of $\boldsymbol {u}$ is

(4.9)\begin{equation} \langle u_i(t, \boldsymbol{k}) u_j (t', \boldsymbol{k}')\rangle = (2{\rm \pi})^3 \delta (t-t') \delta (\boldsymbol{k} + \boldsymbol{k}') \kappa_{ij}(\boldsymbol{k}), \end{equation}

where the appearance of $\delta (\boldsymbol {k} + \boldsymbol {k}')$ in this expression is a consequence of statistical homogeneity. A similar expression is adopted for $\boldsymbol {f}_{\boldsymbol {w}}$ (see Appendix E). Together, incompressibility and isotropy further imply that

(4.10)\begin{equation} \kappa_{ij}(\boldsymbol{k}) = \kappa (k) \mathcal{P}_{ij}(\boldsymbol{k}), \end{equation}

where $\mathcal {P}_{ij}(\boldsymbol {k}) = \delta _{ij} - k_i k_j/k^2$ is the usual $\boldsymbol {k}$-space projection operator. A synthetic velocity field satisfying (4.9) is often called the (incompressible) Kraichnan ensemble, after Kraichnan (Reference Kraichnan1968, Reference Kraichnan1994), who proposed it as a model for studying the behaviour of a passive scalar advected by a turbulent flow. The same model was used independently by Kazantsev (Reference Kazantsev1968) to study the growth of magnetic fields via the turbulent dynamo effect. In both of these applications, the model gave rise to a lively analytical following (see reviews by Falkovich, Gawȩdzki & Vergassola Reference Falkovich, Gaweȩdzki and Vergassola2001 and Rincon Reference Rincon2019). The inertial-range statistics of the passive-velocity equation (4.6) with $\boldsymbol {u}$ the Kraichnan ensemble have also been studied in detail by Benzi et al. (Reference Benzi, Biferale and Toschi2001), Adzhemyan et al. (Reference Adzhemyan, Antonov, Mazzino, Muratore-Ginanneschi and Runov2001a,Reference Adzhemyan, Antonov and Runovb) and Arponen (Reference Arponen2009). The short-correlation-time approximation is a natural one for our problem, because the time scale on which large-scale structures diffuse is much longer than the correlation time of the outer-scale turbulence.

Finally, we assume further that

(4.11)\begin{equation} \kappa (k) = \kappa_0 \delta (k-k_f),\end{equation}

i.e. that the advecting field has a single wavenumber, $k_f$. While this assumption is not strictly required to produce a closed set of equations, it nonetheless greatly simplifies the calculation, and is not particularly limiting considering the simplifications already adopted. It should be noted that (4.11) does not restrict the applicability of the model to turbulence that is forced at a single scale, as $\boldsymbol {f}_{\boldsymbol {w}}$ can still be multi-scale. We also note that there is little to be gained by choosing $\boldsymbol {u}$ to have large-scale structure, because in any situation in which the large-scale structure of the advecting flow is important, the passive model of momentum diffusion will not be appropriate anyway.

We show in Appendix E that under these assumptions, the spectrum $\mathcal {E}_{\boldsymbol {w}}(k)$ of the passive-velocity field $\boldsymbol {w}$ satisfies the following mode-coupling equation

(4.12)\begin{equation} \partial_t \mathcal{E}_{\boldsymbol{w}}(k) + 2 [\nu + \nu_T(k)] k^2 \mathcal{E}_{\boldsymbol{w}}(k) = \frac{\kappa_0 k_f}{ (2{\rm \pi})^2} k \int^{k+k_f}_{|k-k_f|} \frac{\mathrm{d} k' }{k'} K(k',k) \mathcal{E}_{\boldsymbol{w}} ( k') + F_{\boldsymbol{w}}(k),\end{equation}

where the spectrum of energy injection is

(4.13)\begin{equation} F_{\boldsymbol{w}}(k) = \frac{k^2}{2{\rm \pi}^2} \int_0^t \mathrm{d} s\int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{f}_{\boldsymbol{w}}(t) \boldsymbol{\cdot} \boldsymbol{f}'_{\boldsymbol{w}}(s)\rangle \, {\rm e}^{-{\rm i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}}. \end{equation}

The turbulent viscosity $\nu _T(k)$ and the kernel $K(k', k)$ that appears in the mode-coupling integral are both unwieldy functions whose precise forms are given in Appendix E. However, because our interest is in wavenumbers $k \ll k_f$, we only need the small-$k$ part of (4.12), and thus only the small-$k$ limits of $\nu _T(k)$ and $K(k', k)$. These are

(4.14a,b)\begin{equation} \lim_{k\to 0}\nu_T(k)=\frac{\kappa_0 k_f^2}{10{\rm \pi}^2}, \quad \lim_{k,q\to 0}\left[\frac{k_f k}{k_f + q}K(k'=k_f+q,k)\right]=\frac{k^4 - q^4 }{2 k}.\end{equation}

Substituting (4.14a,b) into (4.12) yields, for $k \ll k_f$,

(4.15)\begin{equation} \partial_t \mathcal{E}_{\boldsymbol{w}}(k) + \beta k_f^2 k^2 \mathcal{E}_{\boldsymbol{w}}(k) = \frac{5}{8} \frac{\beta}{k} \int^{k}_{{-}k} \mathrm{d} q (k^4 - q^4) \mathcal{E}_{\boldsymbol{w}} ( k_f+q) + F_{\boldsymbol{w}}(k),\end{equation}

where we have defined $\beta = \kappa _0 / 5{\rm \pi} ^2$ and assumed that the turbulent viscosity dominates over the molecular one. Finally, if $k$ is small compared with the wavenumber scale on which $\mathcal {E}_{\boldsymbol {w}}(k)$ varies in the vicinity of $k_f$, then we may take $\mathcal {E}_{\boldsymbol {w}}(k_f+q)\simeq \mathcal {E}_{\boldsymbol {w}}(k_f)$ in (4.15) (we shall consider the effect of relaxing this assumption in § 4.4). Then (4.15) becomes

(4.16)\begin{equation} \partial_t \mathcal{E}_{\boldsymbol{w}}(k) + \beta k_f^2 k^2 \mathcal{E}_{\boldsymbol{w}}(k) = \beta k^4 \mathcal{E}_{\boldsymbol{w}} (k_f) + Ck^b, \end{equation}

where we have replaced $F_{\boldsymbol {w}}(k)$ by its small-$k$ asymptotic form, taken to be a power law with exponent $b$ (note that the case of finite-band forcing may be recovered by setting $C=0$ in what follows).

Equation (4.16) is coupled to the forcing-scale modes via the appearance of $\mathcal {E}_{\boldsymbol {w}}(k_f)$. Therefore, in order to calculate the growth of $\mathcal {E}_{\boldsymbol {w}}(k\ll k_f)$ from an initial state with $\mathcal {E}_{\boldsymbol {w}}(k)=0$, we should, strictly speaking, compute the evolution of $\mathcal {E}_{\boldsymbol {w}}(k_f)$ from (4.12) and substitute the result into (4.16). However, we expect the spectrum to saturate much more quickly at the forcing scale than at $k\to 0$, so we may, with negligible error, take $\mathcal {E}_{\boldsymbol {w}}(k_f)$ to be equal to its saturated value at all times. Then, solving (4.16) subject to $\mathcal {E}_{\boldsymbol {w}}(t=0,k)=0$ gives

(4.17)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k) = (1-{\rm e}^{-\beta k^2 k_f^2 t}) \frac{C k^b + \beta \mathcal{E}_{\boldsymbol{w}}(k_f)k^4}{\beta k^2 k_f^2}.\end{equation}

As anticipated, (4.17) exhibits a split power law. For any value of $b$, the critical wavenumber demarcating the two regimes is

(4.18)\begin{equation} k_c \sim \frac{1}{\sqrt{\beta k_f^2 t}} \sim \frac{1}{l} \sqrt{\frac{t_{\textrm{nl}}}{t}},\end{equation}

where $l\sim k_f^{-1}$ and $t_{\textrm{nl}}\sim 1/\beta k_f^4$ is the characteristic nonlinear advection time at the injection scale. This is the diffusive scaling for the spectral knee anticipated at the end of § 3.2.

For $k\ll k_c$, (4.17) reduces to

(4.19)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k) = [C k^b + \beta \mathcal{E}_{\boldsymbol{w}}(k_f)k^4]t.\end{equation}

Therefore, at small enough $k$ (or early enough times), $\mathcal {E}_{\boldsymbol {w}}(k)$ has a $k^b$ power law if $b\leq 4$, or a $k^4$ power law if $b>4$. In the case of solenoidal forcing that is local in real space, which has been our focus so far, $b=4$, so $\mathcal {E}_{\boldsymbol {w}}(k\to 0)\propto k^4$, consistent with the numerical results presented in figures 1 and 4. The development of a $k^4$ spectrum in the case of $b>4$ or $C=0$ reflects the fact that turbulence with zero Loitsyansky integral is unsustainable: even if the forcing has $I_{\boldsymbol {f}}=0$, the flow will develop $I\neq 0$ on a dynamical time scale, owing to interactions between eddies (cf. (3.12)). We note that the linear dependence of the right-hand side of (4.19) on $t$ indicates that our passive model of momentum diffusion corresponds to real turbulence with ${\lim _{r\to \infty } r^4 K(r)}$ constant in time (see § 3.3).

In the opposite limit, $k\gg k_c$, (4.17) becomes

(4.20)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k) = \frac{C k^b + \beta \mathcal{E}_{\boldsymbol{w}}(k_f)k^4}{\beta k^2 k_f^2},\end{equation}

so $\mathcal {E}_{\boldsymbol {w}}(k)\propto k^2$ if $b \geq 4$. Thus, we recover the expected thermal spectrum for $b\geq 4$, i.e. for real-space correlations in the forcing function that satisfy ${H_{\boldsymbol {w}}(r\to \infty )\leq O(r^{-5})}$ (where $H_{\boldsymbol {w}}$ is the analogue of $H$ for $\boldsymbol {f}_{\boldsymbol {w}}$) (this statement follows from a calculation directly analogous to the one that showed that $\mathcal {E}(k\to 0)\propto k^a \iff \chi (r\to \infty )\leq O (r^{5})$ for $a\geq 4$; see (3.8)). As explained in § 3.1, this corresponds to the development of a quasi-random momentum distribution, ${\langle \boldsymbol {P}_V^2\rangle \propto R^3}$. If long-range correlations in the forcing are stronger, i.e. $b<4$, then the thermal spectrum is not realised: instead, a shallower $k^{b-2}$ spectrum develops.

According to (3.10), (4.19) and (4.20) correspond to

(4.21) \begin{equation} \langle \boldsymbol{P}^2_V \rangle \propto \begin{cases} R^{7-b} & \text{if }l\ll R\ll R_c, \\ R^{5-b} & \text{if } R \gg R_c\ \&\ b<3, \\ R^{2} & \text{if }R \gg R_c\ \&\ 3< b<4, \\ \end{cases} \end{equation}

where $R_c \equiv k_c^{-1}$. The $R\gg R_c$ scalings are easily interpreted: for large volumes for which there has not been enough time for momentum diffusion to act, the turbulence inherits the momentum scaling dictated by the forcing. The $\langle \boldsymbol {P}^2_V \rangle \propto R^{7-b}$ scaling in the range $l\ll R\ll R_c$ is less intuitive. Interestingly, different values of $b$ in the range $3< b<4$ tend to saturate with different power laws, even though all $\mathcal {E}_{\boldsymbol {w}}\propto k^a$ spectra with $3< a<4$ have $\langle \boldsymbol {P}^2_V \rangle \propto R^2$, as (3.10) shows. To understand the origins of these scalings, it is convenient to consider the momentum-diffusion process as the net result of a series of instances of a decaying passive vector field. While the characteristic scale of the diffusing momentum grows like $t^{1/2}$ in all cases, the energy of the diffusing field decays at a rate that depends on the exponent $b$, and therefore the contribution of each instance of forcing to $\langle \boldsymbol {P}^2_V \rangle \propto R^2$ depends on $b$, even when $3< b<4$. In Appendix F, we show how the scalings (4.21) can be derived directly by thinking about the diffusion of momentum in such terms.

4.3. Local non-solenoidal forcing

While forcing with $b\leq 4$ is easy to implement in numerical simulations, where complete control of the forcing spectrum is possible, such forcing is artificial in the sense that $b<4$ corresponds to long-range correlations that decay with distance in a very particular way (see (3.8)). An important exception to this statement is the case of $b=2$. An expansion of ${F_{\boldsymbol {w}}(k\to 0)}$ analogous to (1.1) yields

(4.22)\begin{equation} F_{\boldsymbol{w}}(k\to0) = \frac{L_{\boldsymbol{f}_{\boldsymbol{w}}}k^2}{2{\rm \pi}^2}+\cdots , \end{equation}

so $b=2$ corresponds to a finite value of

(4.23)\begin{equation} L_{\boldsymbol{f}_{\boldsymbol{w}}} \equiv \int_0^t \mathrm{d} s \int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{f}_{\boldsymbol{w}}(t) \boldsymbol{\cdot} \boldsymbol{f}_{\boldsymbol{w}}'(s) \rangle,\end{equation}

which is the analogue of the Saffman integral for $\boldsymbol {f}_{\boldsymbol {w}}$. If the forcing is solenoidal, then strong long-range correlations in $\boldsymbol {f}_{\boldsymbol {w}}$ are required for $L_{\boldsymbol {f}_{\boldsymbol {w}}}$ to be finite, because then (2.17) yields

(4.24)\begin{equation} L_{\boldsymbol{f}_{\boldsymbol{w}}} = 4{\rm \pi} \lim_{r\to\infty} r^3 H_{\boldsymbol{w}}(r), \end{equation}

so $L_{\boldsymbol {f}_{\boldsymbol {w}}}\neq 0$ requires $H_{\boldsymbol {w}}(r\to \infty )=O(r^{-3})$. However, these long-range correlations need not be present if the forcing is non-solenoidal, as it turns out that they are generated naturally when the non-solenoidal part of $\boldsymbol {f}$ is removed by the action of the projection operator $\mathcal {P}$ – this effect is illustrated in figure 5. Physically, the correlations arise because non-solenoidal forcing generates pressure gradients that decay slowly with distance from the point at which an impulse is applied (these gradients are established instantaneously in an incompressible fluid). This result is due to Saffman (Reference Saffman1967), who used it to argue that naturally occurring decaying turbulence need not have ${\mathcal {E}(k\to 0)\propto k^4}$, as had been supposed by Batchelor & Proudman (Reference Batchelor and Proudman1956). A proof (closely following the one presented by Saffman Reference Saffman1967) and some further comments are given in Appendix G.

Figure 5. The effect of Fourier-space projection of a non-solenoidal forcing. Panel (a) shows a uniformly directed two-dimensional impulse that decays exponentially with distance from the origin. Panel (b) shows the result of removing the non-solenoidal part of this impulse by application of the Fourier-space operator $\mathcal {P}_{ij}= \delta _{ij}-k_ik_j/k^2$; the impulse now falls off much more slowly with distance from the origin.

In summary, there are two values of $b$ that are relevant to turbulence that is forced locally in real space. If the forcing is solenoidal (or non-solenoidal but with $L_{\boldsymbol {f}_{\boldsymbol {w}}}=0$), $b=4$, in which case our passive diffusion model predicts saturation with a thermal $k^2$ spectrum. If the forcing is non-solenoidal and has $L_{\boldsymbol {f}}\neq 0$, then $b=2$, and (4.17) predicts saturation with a flat spectrum: from (4.17), with $C= L_{\boldsymbol {f}_{\boldsymbol {w}}}/2{\rm \pi} ^2$,

(4.25)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k) \simeq (1-{\rm e}^{-\beta k^2 k_f^2 t})\frac{L_{\boldsymbol{f}_{\boldsymbol{w}}}}{2{\rm \pi}^2 \beta k_f^2}. \end{equation}

As we show in figure 6, these predictions hold up reasonably well in our numerical simulations of Navier–Stokes turbulence, as do the corresponding scalings for linear momentum. For the $b=4$ case, panel (a) shows that a finite-band $k^2$ spectrum develops over a wavenumber interval that widens over time, in close agreement with (4.18) and (4.20) (panel (a) is the same as figure 1, and is presented again here to facilitate comparison). We note that, at late times, the $k^2$ spectrum persists all the way to the box scale (to good approximation). Such scales are not strictly within the domain of validity of our theory (which employs isotropic and homogeneous statistics, only valid at scales much smaller than the box size). Nonetheless, such behaviour should be expected under the statistical-mechanical interpretation of the $k^2$ spectrum (see discussion around (1.4)), to which our theory is complementary. In principle, it would be desirable to have much larger simulations with increased separation between box and forcing scales in order to rule out any box-scale effects on the development of the $k^2$ spectrum; however, such simulations are unaffordable at present because of the stiff scaling of the simulation cost with size (accounting for the additional time for a larger simulation to reach saturation, which scales as box size squared by (4.18), the cost is proportional to the fifth power of the size).

Figure 6. Saturation of the large scales in Navier–Stokes turbulence with a delta-correlated Gaussian random forcing. Panel (a) (the same as figure 1) shows the evolution of the energy spectrum, while (b) shows the evolution of the mean square momentum $\langle \boldsymbol {P}^2_V \rangle$, here computed for cubic subvolumes of the box with side length $2R$. Simulations shown in panels (a) and (b) had the forcing spectrum $F(k)\propto k^4 \exp (-k^2/k_p^2)$. Panels (c) and (d) show the same quantities for $F(k)\propto k^2 \exp (-k^2/2k_p^2)$. In both cases, the peak of $F(k)$ is at $k_p = 80$. In panel (c), a numerical fit of the data to $(Ak^2/k_p^2+B)\exp (-Ck^2)$, as explained in the text, is plotted as a dotted line. Insets to panels (a) and (c) show the evolution of the spectral knee $k_c(t)$. In the chosen units, the energy injection rate is $0.7$ in both cases, and the r.m.s. velocities are $\simeq 0.5$. The plotted curves are logarithmically spaced in time, with blue $\to$ red indicating earlier $\to$ later times. Details of the numerical setup are described in §§ 4.1 and 2.4.

Panel (b) shows the corresponding development of $\langle \boldsymbol {P}^2_V \rangle \propto R^3$, although the split-power-law structure in $\langle \boldsymbol {P}^2_V \rangle$, between $\propto R^2$ and $\propto R^3$, is somewhat less pronounced than in the spectrum.

For $b=2$, panel (c) shows that the saturated spectrum is somewhat steeper than $k^0$, although still shallower than $k^2$. Likewise, the decrease in $k_c$ with time is somewhat faster than (4.18) predicts, as is shown by the inset to panel (c). It is plausible that these effects are a consequence of the scale separation between the forcing scale and the scale of the simulation box being insufficient to observe the true $k^0$ large-scale asymptotic: when $b=2$, unlike when $b=4$, (4.20) only reduces to the asymptotic behaviour $\mathcal {E}(k)\propto k^{b-2}$ when

(4.26)\begin{equation} \frac{k}{k_f} \ll \sqrt{\frac{L_{\boldsymbol{f}_{\boldsymbol{w}}}}{2{\rm \pi}^2 \beta k_f^2 \mathcal{E}_{\boldsymbol{w}}(k_f)}}.\end{equation}

If the right-hand side of (4.26) happened to be a moderately small number, then the $k^0$ and $k^2$ terms in (4.20) would be comparable over a range of $k\lesssim k_f$, giving the appearance of a steeper spectrum than $k^0$. To illustrate this possibility, we show in panel (c) of figure 6 a numerical fit of the function ${(A k^2/k_p^2 + B)\exp (-C k^2)}$, where $k_p=80$ is the peak forcing wavenumber, and $A$, $B$, $C$ are fitting parameters, to the final data curve. The result reproduces the data well with $A/B \simeq 20$. Alternatively, the discrepancy with the prediction of the passive-velocity model might be a result of the neglected effect of advection by large-scale modes, which do possess a significant proportion of the total energy for a close-to-flat spectrum. As a result, closure schemes such as the popular EDQNM model may be better at capturing this effect than the passive vector model – see the discussion at the end of § 4.4. Panel (d) shows that $\langle \boldsymbol {P}^2_V \rangle$ follows a scaling reasonably close to $R^{5}$ at late times, which is consistent with a $k^{0}$ spectrum at large scales, according to (3.10). We note that, for a forcing with $b=2$, running the simulation to later times than those shown in figure 6 tends to produce a build-up of energy in the largest-scale modes, making the large-scale asymptotic difficult to measure; a similar effect was found in some of the simulations of Alexakis & Brachet (Reference Alexakis and Brachet2019). Here, we ended our runs before this effect became significant.

4.4. Narrow-band forcing

Finally, we discuss the case of forcing in a narrow spectral band. As noted above, the prediction of the passive model for forcing in a finite band may be obtained by setting $C=0$ in (4.17). This is valid for $k$ much smaller than any other characteristic wavenumber associated with $F_{\boldsymbol {w}}(k)$, including the inverse characteristic width of the forcing spectrum; this assumption entered when we used

(4.27)\begin{equation} k \ll \left[ \frac{1}{\mathcal{E}_{\boldsymbol{w}}(k)}\frac{\mathrm{d} \mathcal{E}_{\boldsymbol{w}}(k)}{\mathrm{d} k} \right]^{{-}1} _{k=k_f},\end{equation}

in order to justify setting $\mathcal {E}_{\boldsymbol {w}}(k_f + q)\simeq \mathcal {E}_{\boldsymbol {w}}(k_f)$ in the $q$ integral in (4.15). Of course, (4.27) is always justified for $k\to 0$. However, if the forcing is concentrated in a narrow band of wavenumbers of width ${\rm \Delta} k \ll k_f$, which is somewhat artificial compared with the more obviously physically realisable ${\rm \Delta} k \sim k_f$, but a common choice for numerical simulations (see, e.g. Alexakis & Brachet Reference Alexakis and Brachet2019), then it is possible that the energy contained in the forced band will greatly exceed the energy contained by the nearby unforced modes. In that case, there will be an extended range of $k$ for which

(4.28)\begin{equation} \left[ \frac{1}{\mathcal{E}_{\boldsymbol{w}}(k)}\frac{\mathrm{d} \mathcal{E}_{\boldsymbol{w}}(k)}{\mathrm{d} k} \right]^{{-}1}_{k=k_f} \sim {\rm \Delta} k \ll k \ll k_f. \end{equation}

For $k$ in this range, (4.17) does not apply. Taking

(4.29)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k_f+q) \simeq \begin{cases} \dfrac{C'}{{\rm \Delta} k}, & \text{if } |q-k_f| < {\rm \Delta} k/2, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

where $C'$ is a constant, modifies the injection term on the right-hand side of (4.16) to be $\propto k^3$ rather than $\propto k^4$, viz.,

(4.30)\begin{equation} \partial_t \mathcal{E}_{\boldsymbol{w}}(k) + \beta k_f^2 k^2 \mathcal{E}_{\boldsymbol{w}}(k) = \tfrac{5}{8}\beta C' k^3.\end{equation}

The saturated spectrum for ${\rm \Delta} k \ll k \ll k_f$ is then readily obtained from (4.30) with $\partial _t\to 0$

(4.31)\begin{equation} \mathcal{E}_{\boldsymbol{w}}(k) \simeq \frac{5}{8} \frac{C' k}{k_f^2},\end{equation}

i.e. $\mathcal {E}_{\boldsymbol {w}}(k)\propto k$, not $k^2$, in this range.

How should we interpret this behaviour? These scalings are reminiscent of two-dimensional turbulence – in two dimensions, the expansion (1.1) of $\mathcal {E}_{\boldsymbol {w}}(k\to 0)$ becomes

(4.32)\begin{equation} \mathcal{E}(k) = \frac{L_{{2D}}k}{4{\rm \pi}}+\frac{I_{{2D}}k^3}{16{\rm \pi}}+\cdots, \end{equation}

where $L_{{2D}}$ and $I_{{2D}}$ are the two-dimensional analogues of the Saffman and Loitsyansky integrals. In the absence of the inverse cascade, therefore, the two-dimensional forced-turbulence spectrum would consist of a growing $k^3$ part at the largest scales, changing to a growing $k^1$ band at $k_c \propto t^{1/2}$, precisely as we have found for a narrow spectral band in three dimensions. While the latter system is not two-dimensional in real space, the same scalings are obtained because the Fourier modes that dominate the mode-coupling integral in (4.16) are confined to the $k=k_f$ surface in Fourier space.

In practice, saturation precisely according to (4.31) is unlikely, because the energy in forced modes usually does not greatly exceed the total energy at the flow scale at saturation. Nonetheless, there can still be some deviation from a precise $k^2$ spectrum, as reported by Alexakis & Brachet (Reference Alexakis and Brachet2019).

To conclude this section, we note that a relation equivalent to (4.16) was derived by Lesieur (Reference Lesieur2008) under the EDQNM closure scheme by assuming (i) Markovian dynamics – i.e. neglecting finite-correlation-time effects and (ii) Kraichnan's (Reference Kraichnan1987a,Reference Kraichnanb) distant-interaction algorithm, in which certain types of non-local (in Fourier space) interactions are discarded. We therefore might have derived the corollaries of (4.16) that were discussed in §§ 4.2, 4.3 and 4.4 as consequences of the EDQNM closure scheme (or other similar models). We also note that it might be possible to recover the development of a steeper-than-$k^0$ large-scale spectrum for $k^2$ forcing under the EDQNM scheme (without applying (ii)). We are grateful to an anonymous referee for pointing this out to us. We have not pursued that line of development in the present work, preferring to use the passive vector model introduced in § 4 as a simple model that clearly illustrates the role of turbulent diffusion of the momentum distribution of initially localised structures as the key process by which the $k^2$ spectrum forms.

5.1. Conservation of momentum in decaying turbulence

The large-scale part of the energy spectrum of instantaneously generated turbulence typically follows an unbroken power law, i.e. $\mathcal {E}(k\ll l^{-1})\propto k^a$. The ‘classical’ exponents $a=2$ and $a=4$ can each arise from initial impulses that do not have long-range spatial correlations (cf. (1.1)), although it is also possible to consider other values for $a$. For $a\leq 3$, the principle of the permanence of the large-scale spectrum that we outlined above is a consequence of the conservation of linear momentum (Saffman Reference Saffman1967; Davidson Reference Davidson2015) – let us briefly review how this works. According to (3.10), $\langle \boldsymbol {P}_V^2\rangle / R^2 \to \infty$ as $R\to \infty$ when $a\leq 3$, indicating that eddies throughout the volume $V$, not just those at its surface, contribute to $\langle \boldsymbol {P}_V^2\rangle$. On the other hand, $\mathrm {d} \langle \boldsymbol {P}_V^2\rangle /\mathrm {d} t = O(R^2)$, because $\langle \boldsymbol {P}_V^2\rangle$ only changes as a result of random fluxes through the surface of $V$ (formally, this requires that ${K(r\to \infty )=o(r^{-3})}$; see (D5) and Davidson Reference Davidson2015). Therefore,

(5.1)\begin{equation} \lim_{V\to\infty}\frac{\mathrm{d} \log \langle\boldsymbol{P}_V^2\rangle}{\mathrm{d} t} = 0,\end{equation}

so the large-$R$ scaling of $\langle \boldsymbol {P}_V^2\rangle$ vs $R$ is preserved and hence so is the large-scale spectral power law, which imposes the scaling

(5.2)\begin{equation} u^2 l^{1+a}\sim \mathrm{const.}\end{equation}

Assuming that the decay is self-similar and occurs on the turnover time scale $t_{nl}\sim l/u$ of the largest eddies, so that $\mathrm {d} u^2/\mathrm {d} t \propto -u^3/l$, it follows from (5.2) that

(5.3a,b)\begin{equation} u^2(t) \sim u^2(t_0)\left(\frac{t-t_0}{t_{{nl,}0}}\right)^{{-}2(1+a)/(3+a)}, \quad l(t)\sim l(t_0) \left(\frac{t-t_0}{t_{{nl,}0}}\right)^{2/(3+a)},\end{equation}

where $t_0$ is the time that the forcing ceases and decay commences and $t_{{nl,}0}\equiv t_{nl}(t_0)\ll t-t_0$. In the classical case of $a=2$ considered by Saffman (Reference Saffman1967), which can result from a non-solenoidal initial forcing without long-range correlations (see § 4.3), (5.3a,b) gives $u^2 \propto (t-t_0)^{-6/5}$, $l\propto (t-t_0)^{2/5}$. In the other classical case of $a=4$, (5.3a,b) indicates that $u^2 \propto (t-t_0)^{-10/7}$, $l\propto (t-t_0)^{2/7}$, which are the decay laws predicted by Kolmogorov (Reference Kolmogorov1941a). Note that the $a=4$ laws do not follow from the conservation of $\langle \boldsymbol {P}_V^2\rangle$, which constrains the decay only when $a \leq 3$. Instead, they are conventionally justified from the invariance of angular momentum via the Loitsyansky integral (1.3), although this idea has been challenged – for details, we refer the reader to the footnote on page 12 and references therein.

As we saw in § 3, the broken large-scale spectrum developed by forced turbulence is a consequence of changes in the local scaling of $\langle \boldsymbol {P}_V^2\rangle$ vs. $R$. This means that any local-power-law energy spectrum must be preserved during decay of previously forced turbulence if the local exponent is $\leq 3$. Therefore, the decaying turbulence must initially follow (5.3a,b) with $a$ set by the local power-law exponent on the infra-red (small-$k$) side of the wavenumber $k\sim l(t_0)^{-1}$ (see figure 7b). This decay will continue until such time when $l(t)^{-1}\sim k_c$.

What happens to $k_c$ during this period? If the $k\to 0$ asymptotic of the spectrum is shallower than $k^3$, it will be preserved by momentum conservation, so $k_c$ must be constant. However, if the $k\to 0$ asymptotic of the spectrum is steeper than $k^3$ – and it is $\propto k^4$ when the forcing is solenoidal and local in real space – then it is a priori unclear what the evolution of $k_c$ might be. In particular, it does not appear possible to argue for preservation of a $\propto k^4$ asymptotic from the conservation of the Loitsyansky integral (1.3) – as we found in § 4, the latter can grow during momentum diffusion, which plausibly could happen during the transient stage of decay, even if it is ruled out for decay with an unbroken power law. We shall therefore employ our passive model of momentum diffusion, (4.8), to determine the evolution of $k_c$ in decaying turbulence.

5.2. Passive evolution of the large scales in decaying turbulence

Let us consider how the small-$k$ part of $\mathcal {E}(k)$ evolves during the period of decay when $l(t) \ll k_c^{-1}$, assuming that the largest scales of the velocity field are advected passively by the decaying integral-scale flow, in the sense of (4.8) with $\boldsymbol {f}_{\boldsymbol {w}}=0$. We again model the energy-containing scales with the Kraichnan ensemble (4.9), taking $k_f \to l^{-1}$ and $\kappa _0 \sim u^3 l$ in (4.11), which follows from assuming that the energy-containing scales evolve in a self-similar manner. Under these choices, (4.16) becomes

(5.4)\begin{equation} \partial_t \mathcal{E}(k) + ul k^2 \mathcal{E} = ul^{3-a} A k^4, \end{equation}

where we have used (5.2) to obtain $\mathcal {E}(l^{-1}) \sim A l^{-a}$, where $A$ is a constant. We note that a relation similar to (5.2) may also be derived from the EDQNM closure approximation under certain assumptions – see the discussion at the end of § 4.4. Let us define $s = \int _{t_0}^t \mathrm {d} t' u(t')l(t')$, where $t_0$ is the time at which forcing ceases. Because $l\sim u(t-t_0)$ by (5.3a,b), $s\sim l^2$ for $t \gg t_{{nl,0}}$ independently of $a$. Changing variables from $t$ to $s$, we find that (5.4) can be rewritten as

(5.5)\begin{equation} \partial_s \mathcal{E}(k) + k^2 \mathcal{E}(k) = l^{2-a} A k^4.\end{equation}

The characteristic scale in $s$ that is associated with the diffusion term on the left-hand side of (5.5), $k^2 \mathcal {E}(k)$, is constant, $\sim 1/k^2$. However, this is also the value of $s$ at which the growing outer scale $l$ reaches $\sim 1/k$, at which point (5.5) ceases to be a valid description of $\mathcal {E}(k)$ at the particular $k$ under consideration. We conclude that the diffusion term can be neglected when the turbulence is decaying. It is straightforward to show that the solution for $t\gg t_{{nl},0}$ of (5.4) without the diffusion term $ul k^2 \mathcal {E}$ is

(5.6)\begin{equation} \mathcal{E}(k, t)\simeq \mathcal{E}(k, t_0) + \mathrm{const.} \times (kl)^4 \mathcal{E}(l^{{-}1}, t) \end{equation}

(when $a = 4$, there is an additional factor of log t in the second term in (5.6)). In words, $\mathcal {E}(k, t)$ is given by the larger of the initial condition or the Batchelor ($\propto k^4$) spectrum that would correspond to the instantaneous outer scale of the turbulence.

Thus, the energies of large-scale modes are preserved as the turbulence decays, provided that these exceed the energies that would correspond to a Batchelor spectrum. For unbroken spectra, this is simply a recovery of the permanence of the large-scale spectrum that was derived from the conservation of $\langle \boldsymbol {P}_V^2\rangle$ in § 5.1. However, for broken spectra (i.e. for $k_c \ll l^{-1}$), (5.6) also shows that there is no evolution of $k_c^{-1}$ before the integral scale reaches it (see figure 7b).

5.3. Decay laws for initially forced turbulence

We now report the analogues of the classical decay laws described above but for turbulence that was forced for a long period by a body force without long-range correlations in real space. If the forcing was solenoidal, so that the spectrum at $t=t_0$ is $\propto k^4$ for $k \ll k_c$ and $\propto k^2$ for $k_c \ll k \ll l(t_0)^{-1}$, then the decay is according to the Saffman laws, viz., (5.3a,b) with $a=2$, until $l\sim k_c^{-1}\sim l(t_0)(t_0/t_{{nl,0}})^{1/2}$, which occurs at $t = t_c \sim t_0 + t_0 (t_0/t_{{nl,0}})^{1/4} \gg t_0$. After this, the decay follows the Batchelor laws, viz., (5.3a,b) with $a=4$ and $t_0$ replaced by $t_c$. An interesting feature of these results is that memory of the forcing is retained for a long time – the Saffman laws are followed for a period that is asymptotically longer (when $t_0/t_{{nl,0}} \gg 1$) than the duration of the forcing stage.

The other interesting case is one with non-solenoidal forcing, for which the spectrum at $t=t_0$ is generically $\propto k^2$ for $k \ll k_c$ and, as per (4.17), $\propto k^0$ for $k_c \ll k \ll l(t_0)^{-1}$. (In § 4.3, we already saw numerical evidence that the $k^0$ law may not actually be developed by real turbulence. This is likely due to the importance of advection by large-scale modes, which is neglected in our model – see discussion in § 4.3. However, a power law reasonably close to $k^0$ does appear to be reached (see figure 6), so we assume the $k^0$ scaling here, for simplicity and lack of a better theory). In this case, the decay follows (5.3a,b) with $a=0$, viz., $u^2\propto (t-t_0)^{-2/3}$, $l\propto (t-t_0)^{2/3}$, until $l\sim k_c^{-1}$ at $t = t_c \sim t_0 + t_0 (t_0/ t_{{nl,0}})^{-1/4} \ll t_0$. Thus, unlike solenoidal forcing, non-solenoidal forcing tends to generate turbulence that has a short memory – the $a=0$ laws are followed for a period that is $\ll t_0$ when $t_0/t_{{nl,0}} \gg 1$. After $t=t_c$, the decay follows the Saffman laws, viz., (5.3a,b) with $a=2$ and $t_0$ replaced by $t_c$.

These results indicate that both types of forcing encourage turbulence to decay in the Saffman regime – for solenoidal forcing, this is because the transient period of Saffman-like decay is asymptotically long compared with the forcing period, while for non-solenoidal forcing, this is because the transient period of non-Saffman decay is asymptotically short compared with the same. A numerical study to test these predictions would be extremely valuable, although it would also incur significant numerical cost – it would be necessary to resolve two scale separations, first between the box size and $k_c(t_0)^{-1}$, and second between $k_c(t_0)^{-1}$ and $l(t_0)$, and to run the forced part of the simulation for long enough for the latter scale separation to be reached while furthermore also ensuring that the Reynolds number at the integral scale were always large enough for the Saffman decay laws ((5.3a,b) with $a=2$) to be valid. We therefore defer such a numerical study to future work (or invite the reader to undertake it).