Emergence of long-range correlations and thermal spectra in forced turbulence

Abstract

Recent numerical studies have shown that forced, statistically isotropic turbulence develops a ‘thermal-equilibrium’ spectrum, $\mathcal {E}(k) \propto k^2$, at large scales. This behaviour presents a puzzle, as it appears to imply the growth of a non-zero Saffman integral, which would require the longitudinal velocity correlation function, $\chi (r)$, to satisfy $\chi (r\to \infty )\propto r^{-3}$. As is well known, the Saffman integral is an invariant of decaying turbulence, precisely because non-local interactions (i.e. interactions via exchange of pressure waves) are too weak to generate such correlations. Subject to certain restrictions on the nature of the forcing, we argue that the same should be true for forced turbulence. We show that long-range correlations and a $k^2$ spectrum arise as a result of the turbulent diffusion of linear momentum, and extend only up to a maximum scale that grows slowly with time. This picture has a number of interesting consequences. First, if the forcing generates eddies with significant linear momentum (as in so-called Saffman turbulence), a thermal spectrum is not reached – instead, a shallower spectrum develops. Secondly, the energy of turbulence that is forced for a while and then allowed to decay obeys Saffman's decay laws for a period that is much longer than the duration of the forcing stage.

1. Introduction

Probably the best-known result in the theory of turbulence is Kolmogorov's law for the spectral energy density in the inertial range, $\mathcal {E}(k)\propto k^{-5/3}$. This law follows from the conjecture of a constant flux of energy in $k$-space, from the large scales at which it is injected, to the small scales at which it is dissipated by molecular viscosity (Kolmogorov 1941b). However, a power-law spectrum can also be found at scales larger than the outer scale of the turbulence, if that scale is small compared with the system's size. Unlike the inertial-range spectrum, this small-$k$ spectral tail does not correspond to large eddies with size $k^{-1}$ – instead, it is controlled by statistical properties of the eddies at the outer scale (Davidson 2015). Provided 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 }$, decays sufficiently quickly with distance, a purely kinematic calculation shows that the energy spectrum of statistically isotropic and homogeneous turbulence satisfies

(1.1)\begin{equation} \mathcal{E}(k\to 0) = \frac{L k^2}{4 {\rm \pi}^2} + \frac{I k^4}{24 {\rm \pi}^2} + O(k^5), \end{equation}

where

(1.2)\begin{equation} L = \int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle, \end{equation}

and

(1.3)\begin{equation} I ={-} \int \mathrm{d}^3 \boldsymbol{r} r^2 \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle,\end{equation}

are known as the Saffman (the Saffman integral is sometimes known as the Saffman–Birkhoff integral, in recognition of the work by Birkhoff (1954) – for convenience, we shall use the more economical standard name in this work) and Loitsyansky integrals, respectively (Loitsyansky 1939; Saffman 1967). These integrals encode information about the distribution of linear and angular momentum in real space (Landau & Lifshitz 1959; Saffman 1967; Davidson 2009). Owing to the conservation of these momenta, it turns out that $L$ and $I$ are invariants of unforced, decaying turbulence (more precisely, $I$ is related to a weighted integral of angular momentum density, and is invariant only if correlations decay sufficiently rapidly with distance (Davidson 2009)), leading to a phenomenon often called the ‘permanence of large-scale structure’ – as turbulence decays, the small-$k$ part of the spectrum remains unchanged. This observation, together with the assumption of self-similarity, allows the decay of kinetic energy to be computed as a function of time (Saffman (1967), Batchelor & Proudman (1956); see Davidson (2015) for a review).

While these results are well established in the theory of decaying turbulence, the large-scale properties of forced turbulence, i.e. one into which energy is continually injected, are usually described in very different terms. In that context, the small-$k$ part of the energy spectrum has received particular attention in recent years, owing to an attractive analogy with statistical mechanics. It has been shown in numerical simulations that there is no net $k$-space energy flux to these scales (Dallas, Fauve & Alexakis 2015), as is to be expected on physical grounds. Accordingly, it has been argued that the largest scales of steady-state forced turbulence might constitute a subsystem in thermal equilibrium with the separate, non-equilibrium subsystem represented by the rest of the flow (Dallas et al. 2015; Cameron, Alexakis & Brachet 2017; Alexakis & Biferale 2018; Alexakis & Brachet 2019). This idea leads immediately to a prediction for the large-scale spectrum: energy should be equipartitioned between Fourier modes (if the large-scale Fourier modes of Navier–Stokes turbulence are taken to constitute a separate system to their smaller-scale forced and dissipating counterparts, then their thermal-equilibrium spectrum follows formally from the statistical mechanics of the truncated Euler equations, $\partial _t \boldsymbol {u} +\mathcal {P}_{K}[\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } p] = 0, \label {truncatedeuler}$ where $\boldsymbol {u}$ is the incompressible velocity field, $p$ is the pressure, determined by $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$, and $\mathcal {P}_K$ is a truncation operator that sets to zero all Fourier modes with $k>K$. This system satisfies a Liouville's theorem, and has an absolute equilibrium state that satisfies (1.4) in the absence of net kinetic helicity (Lee (1952), Orszag (1977) and Kraichnan (1973), see § 1 of Alexakis & Brachet (2019) for a review)), so, in three dimensions,

(1.4)\begin{equation} \mathcal{E}(k) \propto k^2. \end{equation}

In reality, the large-scale modes do not constitute an isolated system, but, if their nonlinear interaction with the turbulent scales is weak, it may be expected that they should develop a close-to-equilibrium state (Alexakis & Brachet 2019). Indeed, (1.4) is well supported by a number of numerical studies conducted in recent years (Dallas et al. 2015; Cameron et al. 2017; Alexakis & Biferale 2018; Alexakis & Brachet 2019). Furthermore, the validity of thermal-equilibrium spectra in more general types of turbulence has been demonstrated experimentally for capillary-wave turbulence by Michel, Pétrélis & Fauve (2017). An experiment to study the large scales of hydrodynamic turbulence is also in development by the same group.

However, like any statistical-mechanics argument, the reasoning outlined above does not elucidate the mechanism by which the equilibrium spectrum is attained. Furthermore, it is unclear what relation (1.4) has to the expansion of $\mathcal {E}(k)$ in terms of the Saffman and Loitsyansky integrals, (1.1). Until now, it has been assumed that the connection between forced turbulence and the concept of decaying ‘Saffman turbulence’ (i.e. that with $L\neq 0$) is superficial, despite both having the same large-scale spectral power law. This is because (i) analysis of the former is mostly concerned with the statistical steady state, obtained by taking a long-time average, while decaying turbulence is, by definition, transient; and (ii) large scales in the former may interact with the forcing, which is absent from decaying turbulence (Alexakis & Brachet 2019). Nonetheless, it should be noted that (1.1) is a purely kinematic result, and must, therefore, apply equally to the forced and decaying cases.

The central goal of the present work is to reconcile the kinematic and statistical-mechanical points of view. This problem turns out to be non-trivial, because of the invariance of the Saffman integral. As we shall show in § 2, this invariance is not restricted to decaying turbulence, but should also apply to forced turbulence, subject to certain reasonable conditions on the nature of the forcing. In particular, if the forcing is solenoidal and sufficiently local in real space (i.e. if its correlations decay sufficiently quickly), then non-local interactions via pressure waves are too weak to generate the long-range longitudinal velocity correlations, $\chi (r\to \infty )\propto r^{-3}$, required for a non-zero Saffman integral, as is the case in decaying turbulence (Batchelor & Proudman 1956; Saffman 1967; Davidson 2015). Thus, the naïve conclusion that the equilibrium spectrum (1.4) simply corresponds to $L\neq 0$ cannot be correct. For consistency with (1.1), therefore, it must always be the case that the equilibrium, $\propto k^2$, part of the spectrum terminates at some large cutoff scale, provided it is smaller than the system size. Above the cutoff scale, (1.1) demands that $\mathcal {E}(k\to 0)\propto k^4$.

In § 3, we shall argue that the physical mechanism by which the equilibrium part of the spectrum develops is the stochasticisation of the distribution of linear momentum, an inevitable consequence of interactions between eddies, even if each of them individually has zero net momentum when it forms. We shall show that this process leads naturally to a split-power-law spectrum at the large scales, with (1.4) satisfied up to a cutoff scale that grows with time, corresponding to the largest scale at which eddies have been able to stochasticise their momentum distribution. The requirement of momentum conservation in these interactions means that different eddies become correlated, which generates the long-range correlations, $\chi (r)\propto r^{-3}$, required for $\mathcal {E}(k)\propto k^2$, although only up to the cutoff scale, above which correlations decay rapidly.

In § 4, we propose a simple, though non-rigorous, model of this phenomenon, in which the large-scale momentum distribution of the flow evolves due to turbulent diffusion caused by flow-scale structures. Under this model, we find that the development of a $k^2$ spectrum is recovered for local, solenoidal forcing, with the cutoff scale growing like $t^{1/2}$. This prediction, along with a number of others, is borne out well in the numerical simulations that we present. Under the same model, we also consider forcing that is local in real space, but not solenoidal – arguably, a more generic situation. Making use of a theorem due to Saffman (1967), we show that such turbulence need not equilibrate at large scales, on account of the long-range real-space correlations present in the solenoidal part of the forcing. Instead, the turbulent diffusion of injected momentum leads to a shallower spectrum than (1.4).

Finally, in § 5, we investigate the implications of the equilibration phenomenon for decaying turbulence. We find inter alia that the energy $E$ of turbulence forced solenoidally without long-range correlations and then allowed to decay satisfies Saffman's law ${E\propto t^{-6/5}}$ (Saffman 1967) for a time period that is much larger than the period of forcing if the latter is large compared with the turnover time of the largest eddies.

Section 6 contains a short summary of our findings, followed by a discussion of their possible applications, implications and extensions in both hydrodynamical contexts and beyond – viz. in astrophysical MHD turbulence.

2. Long-range correlations and the invariance of Saffman's integral

Let us begin by reviewing an important kinematic result: turbulence with an energy spectrum satisfying $\mathcal {E}(k\to 0) \propto k^2$ necessarily has strong long-range correlations in real space (Batchelor & Proudman 1956; Saffman 1967; Davidson 2015).

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 (2015) 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 1959; Davidson 2015), 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 1938), 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 (1956) and Saffman (1967) 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. 2015; Cameron et al. 2017; Alexakis & Biferale 2018; Alexakis & Brachet 2019)? 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. 2015; Cameron et al. 2017; Alexakis & Brachet 2019). The algorithm that we employ to generate $\boldsymbol {f}$ is described in Appendix D of Hosking & Schekochihin (2021). 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 2015) 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 (2019), 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. The large-scale spectrum and linear momentum

Assuming the equivalence of volume and ensemble averages, the definition of the Saffman integral, (1.2), is equivalent to

(3.1)\begin{equation} L = \lim_{V\to \infty}\frac{1}{V}\left\langle \left[\int_V \mathrm{d}^3 \boldsymbol{r}\, \boldsymbol{u}\right]^2 \right\rangle \equiv \lim_{V\to \infty} \frac{\langle \boldsymbol{P}_V^2 \rangle}{V}.\end{equation}

The Saffman integral, therefore, is a measure of how much linear momentum $\boldsymbol {P}_V$ is contained in a large control volume $V$ (Saffman 1967; Davidson 2015). For instance, in Saffman turbulence, where each eddy in $V$ has non-zero, but random, linear momentum, $\langle \boldsymbol {P}_V^2 \rangle \propto V$ (accumulating as a random walk), so $L$ is finite. If, instead, each eddy has vanishing total momentum, as in Batchelor turbulence, then (3.1) is dominated by the contributions of eddies at the surface of $V$. In that case, $\langle \boldsymbol {P}_V^2 \rangle \propto V^{2/3}$, and so $L=0$.

This idea immediately provides a physical explanation for why long-range correlations (2.5) are required for a finite $L$. Consider an isolated turbulent eddy in an otherwise quiescent fluid. The linear momentum contained in a large sphere $V$ of radius $R$, centred on the eddy, is

(3.2)\begin{equation} \boldsymbol{P}_{\textrm{eddy}}=\int_V \mathrm{d}^3 \boldsymbol{r}\,\boldsymbol{u} = \int_{\partial V} \mathrm{d} \boldsymbol{S} \times \boldsymbol{A}, \end{equation}

where we have represented the solenoidal velocity field as $\boldsymbol {u} = \boldsymbol {\nabla } \times \boldsymbol {A}$. Clearly, $\boldsymbol {P}_{\textrm{eddy}}$ vanishes unless the average of $\boldsymbol {A}$ over $\partial V$ scales as $R^{-2}$ as $R\to \infty$, implying that the mean velocity on $\partial V$ scales as $R^{-3}$. One can imagine building a synthetic $L\neq 0$ turbulence by superimposing such eddies with random positions and orientations; this velocity field will necessarily have long-range correlations, owing to the long tails of the component eddies.

In fact, there is a deep connection between the linear-momentum content of the turbulence and the large-scale spectrum, that goes beyond the finiteness of the Saffman integral and the Taylor expansion (1.1). Davidson (2015) has shown that, in incompressible, homogenenous and isotropic turbulence, $\langle \boldsymbol {P}_V^2 \rangle$ is a functional of $\chi (r)$: if $V$ is a sphere of radius $R$

(3.3)\begin{equation} \langle \boldsymbol{P}_V^2 \rangle= 4{\rm \pi}^2 R^2 u^2 \int^{2R}_0 \mathrm{d} r r^3 \chi(r) \left[1-\left(\frac{r}{2R}\right)^2 \right]. \end{equation}

It is convenient to integrate this formula by parts, which gives

(3.4)\begin{equation} \langle \boldsymbol{P}_V^2 \rangle = 2{\rm \pi}^2 u^2 \int^{2R}_0 \mathrm{d} r r \int_0^r \mathrm{d} r' r'^3 \chi(r'),\end{equation}

boundary terms having vanished exactly. From (3.4), it is clear that $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ only if $\chi (r'\to \infty )\propto r'^{-3}$. If, instead, $\chi (r')$ decays quickly, viz., $\chi (r'\to \infty )=o(r'^{-4})$, the $r'$ integral in (3.4) is dominated by small $r'$, and hence $\langle \boldsymbol {P}_V^2 \rangle \propto R^2$, which is the scaling $\langle \boldsymbol {P}_V^2 \rangle \propto V^{2/3}$ obtained above. Equation (3.4) is also readily inverted, to yield

(3.5)\begin{equation} u^2 \chi(2R) = \frac{1}{128 {\rm \pi}^2} \frac{1}{R^3} \frac{\partial}{\partial R} \frac{1}{R} \frac{\partial \langle \boldsymbol{P}_V^2 \rangle}{\partial R}. \end{equation}

Therefore, full knowledge of $\langle \boldsymbol {P}_V^2 \rangle$ as a function of $R$ is equivalent to full knowledge of $\chi (r)$, and hence, via (2.1) and (2.3), of the energy spectrum. This observation suggests that one might seek the explanation of the growth of the thermal spectrum in figure 1 as a consequence of the evolution of $\langle \boldsymbol {P}_V^2 \rangle$.

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 (2015) 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 2015) 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 (2011). 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 (1967), Davidson (2011); 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 (2009, 2011) 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 2000; Lesieur, Métais & Comte 2005; Lesieur 2008).

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 1941a; Landau & Lifshitz 1959; Ishida, Davidson & Kaneda 2006; Davidson 2015). While this point is not universally accepted (e.g. $I$ is not conserved under the popular EDQNM closure; see Lesieur (2008) 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. 2006). 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 (1956), Davidson (2015); 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. 2006), 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. A solvable model of passive momentum diffusion

In § 3, we argued that the development of a thermal $k^2$ spectrum over a finite, but growing, large-scale band is a consequence of the quasi-randomisation of the linear-momentum distribution. In this section, we consider a model of this process in which the momentum density is a passive quantity, in which case its randomisation can be understood as a consequence of turbulent diffusion.

To motivate the model, let us consider the evolution of a velocity field $\boldsymbol {w}$ under the Navier–Stokes equations. Let $\boldsymbol {w} = \bar {\boldsymbol {w}} + \tilde {\boldsymbol {w}}$, where $\bar {\boldsymbol {w}}$ is the large-scale part of $\boldsymbol {w}$, formally defined as the result of applying a Fourier-space filter to $\boldsymbol {w}$ to isolate only those modes with $k< K$, for some $K$ much smaller than the characteristic wavenumber of the forcing, while $\tilde {\boldsymbol {w}}$ is the remaining smaller-scale part, consisting of modes with $k>K$. Then the evolution of $\boldsymbol {w}$ proceeds according to

(4.1)\begin{equation} \frac{\partial \boldsymbol{w}}{\partial t} + \mathcal{P}[\bar{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{w}} + \bar{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{w}} + \tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{w}} + \tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{w}}] = \nu \nabla^2 \boldsymbol{w} + \boldsymbol{f}, \end{equation}

where $\mathcal {P}$ is the Fourier-space projection operator that returns the solenoidal part of the field on which it operates: $[\mathcal {P}\boldsymbol {w}]_i \equiv (\delta _{ij}-k_ik_j/k^2)w_j$. Let us assume that, because the large-scale modes are energetically subdominant to the rest of the flow, advection by them is unimportant. Then we are left with

(4.2)\begin{equation} \frac{\partial \boldsymbol{w}}{\partial t} + \mathcal{P}[\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w}] = \nu \nabla^2 \boldsymbol{w} + \boldsymbol{f},\end{equation}

so the only important nonlinearity is advection by the small-scale part of $\boldsymbol {w}$. The small-scale part of (4.2) is

(4.3)\begin{equation} \frac{\partial \tilde{\boldsymbol{w}}}{\partial t} + \mathcal{P}[\widetilde{\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{w}}} + \widetilde{\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{w}}}] = \nu \nabla^2 \tilde{\boldsymbol{w}} + \tilde{\boldsymbol{f}}. \end{equation}

Again, owing to the energetic subdominance of $\bar {\boldsymbol {w}}$ (and its small gradients), let us assume that the term involving $\bar {\boldsymbol {w}}$ is negligible compared with the other term inside the square brackets, so

(4.4)\begin{equation} \frac{\partial \tilde{\boldsymbol{w}}}{\partial t} + \mathcal{P}[\widetilde{\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{w}}} ] = \nu \nabla^2 \tilde{\boldsymbol{w}} + \tilde{\boldsymbol{f}}.\end{equation}

Equation (4.4) shows that, under the approximations outlined so far, the evolution of $\tilde {\boldsymbol {w}}$ is entirely decoupled from that of $\bar {\boldsymbol {w}}$. Taking the large-scale part of (4.2), we find that $\bar {\boldsymbol {w}}$ satisfies

(4.5)\begin{equation} \frac{\partial \bar{\boldsymbol{w}}}{\partial t} + \mathcal{P}[\overline{\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{\boldsymbol{w}}} + \overline{\tilde{\boldsymbol{w}}\boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\boldsymbol{w}}} ] = \nu \nabla^2 \bar{\boldsymbol{w}} + \bar{\boldsymbol{f}}.\end{equation}

In a sense, therefore, $\bar {\boldsymbol {w}}$ is a passive field: although (4.5) shows that its evolution is affected by $\tilde {\boldsymbol {w}}$, $\tilde {\boldsymbol {w}}$ is not affected by $\bar {\boldsymbol {w}}$, according to (4.4).

Motivated by this property, we propose to replace $\tilde {\boldsymbol {w}}$ with an artificial field, $\boldsymbol {u}$, wherever the former appears as an advecting field in (4.4) and (4.5). This model can be summarised by

(4.6)\begin{equation} \frac{\partial \boldsymbol{w}}{\partial t} + \mathcal{P}[\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w}] = \nu \nabla^2 \boldsymbol{w} + \boldsymbol{f}.\end{equation}

This equation is sometimes called the ‘linear pressure model’ of the Navier–Stokes equation; a number of its properties have been studied by Benzi, Biferale & Toschi (2001), Adzhemyan et al. (2001a), Adzhemyan, Antonov & Runov (2001b), Antonov et al. (2003) and Arponen (2009). Physically, (4.6) describes ‘eddies’ of the field $\boldsymbol {w}$ interacting nonlinearly with eddies of the field $\boldsymbol {u}$, rather than other $\boldsymbol {w}$-eddies. The $\boldsymbol {w}$-eddies can receive momentum from their interaction with the $\boldsymbol {u}$-eddies, whereas the latter do not get anything back, as their motion is externally prescribed. Nonetheless, the receipt of momentum by $\boldsymbol {w}$-eddies still occurs in a way that locally conserves momentum, because $\int \mathrm {d}^3\boldsymbol {x} \boldsymbol {w}$ is an invariant of (4.6). Ultimately, then, the $\boldsymbol {w}$-eddies do have local interactions that satisfy net-momentum conservation, which is the key ingredient for the stochasticisation of their momentum distribution. If the field $\boldsymbol {u}$ is chosen so that its statistical properties are close to those of real turbulence, it may be hoped that the evolution of $\bar {\boldsymbol {w}}$ should mimic that of the large-scale part of a real velocity field (the same need not be true of $\tilde {\boldsymbol {w}}$, although see Benzi et al. 2001 for some similarities in small-scale properties).

In § 4.2, we shall present an analytic treatment of (4.6), taking $\boldsymbol {u}$ to be the so-called Kraichnan ensemble (Kraichnan 1965, 1994). First, however, we present numerical simulations to demonstrate the validity of the model (4.6).

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 (1968, 1994), 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 (1968) 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 2001 and Rincon 2019). 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. (2001), Adzhemyan et al. (2001a,b) and Arponen (2009). 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 (1967), 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 (1956). A proof (closely following the one presented by Saffman 1967) 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 (2019). 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 2019), 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 (2019).

To conclude this section, we note that a relation equivalent to (4.16) was derived by Lesieur (2008) under the EDQNM closure scheme by assuming (i) Markovian dynamics – i.e. neglecting finite-correlation-time effects and (ii) Kraichnan's (1987a,b) 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. Decay of initially forced turbulence

In this section, we consider how turbulence that has been forced for a long period decays after the forcing is removed. This situation is somewhat different to the one usually considered in theoretical treatments of decaying turbulence (see Davidson (2015) for a review and, e.g. Panickacheril John, Donzis & Sreenivasan (2022) for a recent numerical study), where it is common to consider an initial condition that was generated effectively instantaneously (physically, over a period $\lesssim$ the initial eddy turnover time). For an initial condition generated in this way, the decay is usually believed to be governed by the principle of the ‘permanence of the large-scale spectrum’: while the energy contained at the flow scale cascades to small scales and is dissipated, the large-scale power law is preserved, i.e. ${\mathcal {E}(k\ll l(t)^{-1})\simeq \mathrm {const.}}$ (see figure 7a). In other words, there is no thermalisation of the large scales and no quasi-randomisation of the momentum distribution. In this section, we shall argue that the same is true for previously forced, now decaying turbulence, even if $k_c^{-1}$ is initially much larger than the outer scale. This means that the growing integral scale eventually becomes comparable to $k_c^{-1}$, and, when it exceeds $k_c^{-1}$, the decay laws change – although, as we shall see, this can happen after a long time compared with the duration of the forced stage.

Figure 7. Schematic diagrams of the evolution of the energy spectrum of decaying isotropic turbulence initially forced for a period $t_0$, where (a) $t_0 \lesssim$ the initial eddy-turnover time and (b) $t_0 \gg$ the initial eddy-turnover time. The exponent $p$ is either $4$ or the asymptotic spectral exponent of the forcing as $k\to 0$, whichever is smaller (cf. (4.17)). Blue $\to$ red indicates earlier $\to$ later times.

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 1967; Davidson 2015) – 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 2015). 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 (1967), 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 (1941a). 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).

6. Conclusion

In this work, we have addressed the apparent disconnect between the ‘decaying-turbulence view’ of the large-scale structure of turbulence, i.e. the notion that it is determined kinematically by the values of certain invariants that describe statistical properties of the flow field, and the increasingly popular idea that the large-scale spectral tail might (in some cases) constitute an isolated subsystem in thermal equilibrium. If the latter were true, equipartition of energy between Fourier modes would imply $\mathcal {E}(k\to 0)\propto k^2$, which corresponds to a non-zero Saffman integral, $L$ (see (1.1)). However, as we found in § 2, $L$ is an invariant not only of decaying turbulence, but also of forced turbulence, provided that the forcing is solenoidal and sufficiently localised in real space. This invariance is a manifestation of the conservation of linear momentum: solenoidal, localised forcing can only generate eddies with vanishing total momentum, so the total momentum contained within a sufficiently large volume of the turbulence must always vanish (in the sense that only surface contributions matter to the total), in which case $L=0$ (see (3.1)).

Nonetheless, the total momentum contained within a finite volume of turbulence need not vanish, provided that one waits long enough for the momentum distribution to become stochastic on the relevant scale. As explained in § 3, this scale-dependent stochasticisation is the net result of the momentum-conserving interactions between eddies, and leads to a broken-power-law spectrum at large scales. The ‘knee’ wavenumber that separates the developing $k^2$ spectrum from the ‘Batchelor’ asymptotic $\mathcal {E}(k\to 0) \propto k^4$ is a decreasing function of time: numerically, we have found it to decrease like $k_c \propto t^{-1/2}$, as would be expected for a turbulently diffusing field (of which a good solvable model can be constructed using the Kraichnan flow: see § 4). In a finite system, such as a numerically simulated turbulence in a periodic box, the spectral knee eventually reaches the box scale, at which point the turbulence is essentially identical in character to a finite-system-size approximation of ‘Saffman turbulence’, i.e. turbulence with $L\neq 0$. In particular, $\mathcal {E}(k)\propto k^2$ at all resolved scales larger than the forcing scale. However, in an infinite domain, the spectral knee grows indefinitely, and increasingly slowly, as increasingly distant points become correlated. Interestingly, these conclusions may be modified somewhat if the forcing is non-solenoidal and thus momentum injecting, as explained in § 4.3 (see also below), and also if it is in a narrow spectral band, as explained in § 4.4.

We anticipate that there may be a number of applications of the ideas developed in this work to more complex variants of fluid turbulence – in particular, to naturally occurring astrophysical turbulence. One such astrophysical application might be to the evolution of primordial magnetic fields in the early universe. Typically, these magnetic fields are assumed to have a magnetic-energy spectrum $\propto k^4$ at large scales, because it is usually thought that long-range correlations in the magnetic field, of the sort required for a shallower spectrum, are excluded by causality constraints imposed by cosmological inflation models (Durrer & Caprini 2003; Brandenburg, Kahniashvili & Tevzadze 2015; Brandenburg & Kahniashvili 2017; Reppin & Banerjee 2017). In the statistically isotropic case, their decay then proceeds via reconnection of magnetic-field lines, conserving magnetic helicity, either in a ‘net’ (Hatori 1984; Biskamp & Müller 1999; Brandenburg & Kahniashvili 2017; Hosking & Schekochihin 2021) or ‘fluctuating’ sense (Hosking & Schekochihin 2021). However, intuitively, it seems likely that a process akin to the one described in this work may be able to induce a magnetic-energy spectrum $\propto k^2$ over a finite range of scales, if the magnetic energy is maintained by an effective ‘forcing’ from the velocity field (i.e. magnetic dynamo). In this case, the process of momentum diffusion would be replaced by ‘flux diffusion’, which presumably can occur due to magnetic reconnection. While a dedicated study would be necessary for a complete understanding of this effect, we note that numerical simulations of the magnetohydrodynamic (MHD) fluctuation dynamo that have scale separation between the box size and forcing scale do indeed appear to saturate with a magnetic-energy spectrum $\propto k^2$ at large scales (Maron & Blackman 2002; Brandenburg, Zhou & Sharma 2023). This would motivate consideration of a $k^2$ large-scale spectrum in the primordial magnetic field. In the process of decay of such a field, the ‘Saffman flux invariant’,

(6.1)\begin{equation} L_{\boldsymbol{B}} = \int \mathrm{d}^3 \boldsymbol{r} \langle \boldsymbol{B}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{B} (\boldsymbol{x} + \boldsymbol{r})\rangle,\end{equation}

which is the analogue of the Saffman integral (1.2) but for magnetic flux, should be conserved – the interested reader will find discussion of the effect of $L_{\boldsymbol {B}}$ on MHD decay laws in the Supplementary Information of Hosking & Schekochihin (2022). We also note that momentum and flux diffusion, in the sense described here, provides a mechanism for the transfer of some energy to large spatial scales. It might therefore be useful to consider this effect in the context of the fluctuation (non-helical) dynamo, whose saturated energy spectrum is an outstanding theoretical problem (see, e.g., Galishnikova, Kunz & Schekochihin (2022), Rincon (2019) and references therein).

Another application of the ideas presented here to MHD turbulence concerns the two-dimensional spectra of the latter in the presence of a strong mean field, in which case the idea of thermalisation turns out to have some traction in treating the scales perpendicular to the mean field that lie in the inertial range but are larger than the ‘critical-balance’ scale – an intrigued reader will find the details in Appendix B of Schekochihin (2022). Such two-dimensional spectral scalings turn out to be of some consequence also in the theory of phase-space turbulence in kinetic plasmas (Schekochihin et al. 2016).

Returning to hydrodynamics, one intriguing finding of this study is that turbulence does not develop a thermal large-scale spectrum when the forcing generates eddies with non-zero linear momentum directly. This can occur even with a spatially localised forcing, provided it is non-solenoidal. In that case, long-range correlations can be generated when the non-solenoidal part is removed, an idea that underpins the realisability of decaying Saffman turbulence (Saffman 1967). The measurement of a flatter-than-$k^2$ large-scale spectrum in a forced-turbulence experiment, such as the one under construction at ENS Paris, would therefore represent a neat direct demonstration of the physics underlying Saffman's theory of decaying turbulence, as would direct measurement of the Saffman decay laws (see § 5.1) in turbulence forced locally and solenoidally until a large-scale thermal spectrum developed, and then allowed to decay. From the statistical-mechanics perspective, the thermodynamical motivation for an equilibrium spectrum (see § 1) relies on weak interaction between the forcing scales and the much larger ones. This condition is satisfied when the forcing is solenoidal, because then the fluid response is essentially local in real space (see § 2.2), but it is violated when the forcing is non-solenoidal, because then long-range interactions between distant points via exchange of pressure waves become an important feature of the dynamics, even though the forcing itself might be local in real space.

In light of this observation, we suggest that an interesting topic for further study would be the large-scale structure of compressible turbulence, in which pressure (sound) waves propagate with finite velocity. In the large-Mach-number limit of highly supersonic motions, turbulent diffusion should be the dominant mechanism of momentum transport, so equilibration of the large scales should be possible. At finite Mach number, however, sound waves may correlate distant points before the turbulent diffusion of momentum can, owing to the fact that sound waves propagate ballistically, rather than diffusively. Whether this precludes thermalisation of the large scales in forced compressible turbulence is an intriguing question.