Hostname: page-component-848d4c4894-nmvwc Total loading time: 0 Render date: 2024-07-06T00:30:59.362Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  16 October 2023

D.N. Hosking*
Affiliation:
Oxford Astrophysics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544, USA Merton College, Merton Street, Oxford OX1 4JD, UK Gonville & Caius College, Trinity Street, Cambridge CB2 1TA, UK
A.A. Schekochihin
Affiliation:
Merton College, Merton Street, Oxford OX1 4JD, UK Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Parks Road, Oxford OX1 3PU, UK
*
Email address for correspondence: dhosking@princeton.edu

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.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press

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 Reference Kolmogorov1941b). 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 Reference Davidson2015). 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 (Reference Birkhoff1954) – for convenience, we shall use the more economical standard name in this work) and Loitsyansky integrals, respectively (Loitsyansky Reference Loitsyansky1939; Saffman Reference Saffman1967). These integrals encode information about the distribution of linear and angular momentum in real space (Landau & Lifshitz Reference Landau and Lifshitz1959; Saffman Reference Saffman1967; Davidson Reference Davidson2009). 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 Reference Davidson2009)), 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 (Reference Saffman1967), Batchelor & Proudman (Reference Batchelor and Proudman1956); see Davidson (Reference Davidson2015) 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 Reference Dallas, Fauve and Alexakis2015), 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. Reference Dallas, Fauve and Alexakis2015; Cameron, Alexakis & Brachet Reference Cameron, Alexakis and Brachet2017; Alexakis & Biferale Reference Alexakis and Biferale2018; Alexakis & Brachet Reference Alexakis and Brachet2019). 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 (Reference Lee1952), Orszag (Reference Orszag1977) and Kraichnan (Reference Kraichnan1973), see § 1 of Alexakis & Brachet (Reference Alexakis and Brachet2019) 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 Reference Alexakis and Brachet2019). Indeed, (1.4) is well supported by a number of numerical studies conducted in recent years (Dallas et al. Reference Dallas, Fauve and Alexakis2015; Cameron et al. Reference Cameron, Alexakis and Brachet2017; Alexakis & Biferale Reference Alexakis and Biferale2018; Alexakis & Brachet Reference Alexakis and Brachet2019). 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 (Reference Michel, Pétrélis and Fauve2017). 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 Reference Alexakis and Brachet2019). 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 Reference Batchelor and Proudman1956; Saffman Reference Saffman1967; Davidson Reference Davidson2015). 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 (Reference Saffman1967), 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 Reference Saffman1967) 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 Reference Batchelor and Proudman1956; Saffman Reference Saffman1967; Davidson Reference Davidson2015).

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

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

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

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

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

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

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

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

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

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

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

Thus, $L$ is finite if and only if

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

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

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

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

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

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

Taking the divergence of the Navier–Stokes equation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3. Correlations generated directly by the forcing

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

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

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

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

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

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

which is unsurprising, given (2.5).

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

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

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

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

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

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

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

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

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

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

3. 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 Reference Saffman1967; Davidson Reference Davidson2015). 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 (Reference Davidson2015) 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 (Reference Davidson2015) for the true Saffman and Batchelor turbulence, is shown in figure 2.

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

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

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

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

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

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

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

Let us explain qualitatively each case in turn.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

which is precisely the diffusive scaling suggested above.

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

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

4. 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 (Reference Benzi, Biferale and Toschi2001), Adzhemyan et al. (Reference Adzhemyan, Antonov, Mazzino, Muratore-Ginanneschi and Runov2001a), Adzhemyan, Antonov & Runov (Reference Adzhemyan, Antonov and Runov2001b), Antonov et al. (Reference Antonov, Hnatich, Honkonen and Jurčišin2003) and Arponen (Reference Arponen2009). 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. Reference Benzi, Biferale and Toschi2001 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 Reference Kraichnan1965, Reference Kraichnan1994). 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 (Reference Kraichnan1968, Reference Kraichnan1994), who proposed it as a model for studying the behaviour of a passive scalar advected by a turbulent flow. The same model was used independently by Kazantsev (Reference Kazantsev1968) to study the growth of magnetic fields via the turbulent dynamo effect. In both of these applications, the model gave rise to a lively analytical following (see reviews by Falkovich, Gawȩdzki & Vergassola Reference Falkovich, Gaweȩdzki and Vergassola2001 and Rincon Reference Rincon2019). The inertial-range statistics of the passive-velocity equation (4.6) with $\boldsymbol {u}$ the Kraichnan ensemble have also been studied in detail by Benzi et al. (Reference Benzi, Biferale and Toschi2001), Adzhemyan et al. (Reference Adzhemyan, Antonov, Mazzino, Muratore-Ginanneschi and Runov2001a,Reference Adzhemyan, Antonov and Runovb) and Arponen (Reference Arponen2009). The short-correlation-time approximation is a natural one for our problem, because the time scale on which large-scale structures diffuse is much longer than the correlation time of the outer-scale turbulence.

Finally, we assume further that

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

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

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

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

where the spectrum of energy injection is

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3. Local non-solenoidal forcing

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.4. Narrow-band forcing

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5. 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 (Reference Davidson2015) for a review and, e.g. Panickacheril John, Donzis & Sreenivasan (Reference Panickacheril John, Donzis and Sreenivasan2022) 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 Reference Saffman1967; Davidson Reference Davidson2015) – let us briefly review how this works. According to (3.10), $\langle \boldsymbol {P}_V^2\rangle / R^2 \to \infty$ as $R\to \infty$ when $a\leq 3$, indicating that eddies throughout the volume $V$, not just those at its surface, contribute to $\langle \boldsymbol {P}_V^2\rangle$. On the other hand, $\mathrm {d} \langle \boldsymbol {P}_V^2\rangle /\mathrm {d} t = O(R^2)$, because $\langle \boldsymbol {P}_V^2\rangle$ only changes as a result of random fluxes through the surface of $V$ (formally, this requires that ${K(r\to \infty )=o(r^{-3})}$; see (D5) and Davidson Reference Davidson2015). Therefore,

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

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

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

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

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

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

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

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

5.2. Passive evolution of the large scales in decaying turbulence

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

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

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

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

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

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

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

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

5.3. Decay laws for initially forced turbulence

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

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

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

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 Reference Durrer and Caprini2003; Brandenburg, Kahniashvili & Tevzadze Reference Brandenburg, Kahniashvili and Tevzadze2015; Brandenburg & Kahniashvili Reference Brandenburg and Kahniashvili2017; Reppin & Banerjee Reference Reppin and Banerjee2017). In the statistically isotropic case, their decay then proceeds via reconnection of magnetic-field lines, conserving magnetic helicity, either in a ‘net’ (Hatori Reference Hatori1984; Biskamp & Müller Reference Biskamp and Müller1999; Brandenburg & Kahniashvili Reference Brandenburg and Kahniashvili2017; Hosking & Schekochihin Reference Hosking and Schekochihin2021) or ‘fluctuating’ sense (Hosking & Schekochihin Reference Hosking and Schekochihin2021). 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 Reference Maron and Blackman2002; Brandenburg, Zhou & Sharma Reference Brandenburg, Zhou and Sharma2023). 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 (Reference Hosking and Schekochihin2022). 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 (Reference Galishnikova, Kunz and Schekochihin2022), Rincon (Reference Rincon2019) 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 (Reference Schekochihin2022). 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. Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016).

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 Reference Saffman1967). 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.

Acknowledgements

This work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).

Funding

D.N.H. was supported by a UK STFC studentship. The work of A.A.S. was supported in part by the UK EPSRC grant EP/R034737/1.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Invariance of the Saffman integral in forced turbulence

In this appendix, we formalise the arguments presented in § 2.2 for the invariance of the Saffman integral in forced turbulence. Our key assumption is that the value of the forcing function at any given point in space and time is statistically independent of the value it has at all other finite times at arbitrarily distant points in space. More precisely, we assume

A-I: All cumulants of $f(\boldsymbol {x},t)$,

(A1)\begin{equation} \langle f_i(\boldsymbol{x}, t) f_j(\boldsymbol{x}', t')f_k(\boldsymbol{x}'', t'') \ldots \rangle_c, \end{equation}

decay sufficiently quickly with distance for their integral moments (with respect to $\boldsymbol {x}'-\boldsymbol {x}$, $\boldsymbol {x}''-\boldsymbol {x}$, etc.) to converge.

For $n>2$, the cumulant $\langle \cdots \rangle _c$ of the $n$th-order correlation function $\langle \cdots \rangle$ is the difference between $\langle \cdots \rangle$ and the value in terms of second-order correlators that it would have if the underlying statistics were Gaussian. For example, $\langle f_i f_j'f_k'' f_l''' \rangle _c \equiv \langle f_i f_j'f_k'' f_l''' \rangle - \langle f_if_j'\rangle \langle f_k''f_l'''\rangle - \langle f_if_k''\rangle \langle f_j'f_l'''\rangle - \langle f_if_l'''\rangle \langle f_j'f_k''\rangle$, where $f_i \equiv f_i(\boldsymbol {x}, t)$, $f_i' \equiv f_i(\boldsymbol {x}', t')$, etc. For $n\leq 2$, $\langle \cdots \rangle _c = \langle \cdots \rangle$.

A-I is a natural generalisation to forced turbulence of the assumption adopted by Batchelor & Proudman (Reference Batchelor and Proudman1956) for decaying turbulence, namely that integral moments of cumulants of $\boldsymbol {u}$ at the initial time converge (Our results reduce to theirs for $\boldsymbol {f}(\boldsymbol {x}, t)\to \boldsymbol {f}(\boldsymbol {x})\delta (t)$) (also see Saffman (Reference Saffman1967) for a similar analysis of the decaying-turbulence problem in Fourier space). Using methods analogous to theirs, we shall use A-I to determine from the forced Navier–Stokes equation all of the initial ($t=0$) time derivatives of the velocity correlation functions in the large-separation limit. We further follow Batchelor & Proudman (Reference Batchelor and Proudman1956) in assuming that

A-II: The large-separation asymptotics of correlators at $t>0$ can be written as convergent Taylor series in $t$, with derivatives evaluated at $t=0$.

The time derivative in such a Taylor series that decays slowest with the separation $r$ gives the large-$r$ asymptotic of the correlator at $t>0$.

We note that A-II is neither trivial nor unquestionable. Although he uses a similar assumption in his theory of decaying turbulence, Saffman (Reference Saffman1967) points out that there exist initial conditions for which the resulting velocity field is not an analytic function of time, which suggests that correlation functions need not be either. Furthermore, correlation functions are not, in general, analytic functions of time in diffusing systems. For example, with $u$ a statistically homogeneous and isotropic solution of the one-dimensional diffusion equation ${\partial u/\partial t = D\partial ^2 u/\partial x^2}$, the correlation function ${C(r,t)=\langle u(x,t)u(x+r,t)\rangle }$ also solves a diffusion equation: $\partial C/\partial t = 2D\partial ^2 C/\partial r^2$. Therefore, if $C(r,t=0)$ has compact support $r< R$, then $C(r>R,t)$ is non-analytic in $t$ – we see from the diffusion equation for $C$ that all of its $t$ derivatives are zero at $t=0$ for $r>R$, but $C(r>R,t)$ is not zero at all times: the support of $C(r,t)$ is only compact at $t=0$. We are grateful to an anonymous referee for pointing this example out to us. It should be noted, however, that such examples do not necessarily contradict A-II, as the latter is an assumption only about the behaviour of the large-separation asymptotics of correlators, rather than about the correlators (or solutions) in general. For example, even in the case of diffusion of a quantity with initially compact support, A-II does give the correct asymptotic behaviour at large separation: at arbitrary finite $t$, the large-$r$ limit of $C(r,t)$ is super-exponentially small. Batchelor & Proudman (Reference Batchelor and Proudman1956) similarly note that A-II is an assumption about asymptotics of correlation functions of solutions to the Navier–Stokes equation, which, being averaged properties, are likely to be better behaved than the solutions themselves. For these reasons, we shall take A-II as acceptable in what follows, but we do flag to the reader that it is far from trivial, and remains an assumption, rather than a theorem.

Equation (2.15) contains two terms that can cause growth of the Saffman integral $L$: one that involves the triple correlator $K(r)$ and the other involving $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {f}' \rangle$. From the isotropy and solenoidality of $\boldsymbol {u}$ and $\boldsymbol {f}$, we have that (cf. (2.2))

(A2)\begin{equation} \varPsi_{ij}\equiv \langle u_i f_j'\rangle = \frac{1}{2r}\left[\frac{\partial}{\partial r}(r^2 \psi) \delta_{ij} - \frac{\partial \psi}{\partial r} r_i r_j\right] \implies \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle = \frac{1}{r^2}\frac{\partial}{\partial r} (r^3\psi). \end{equation}

Therefore, (2.15) may be rewritten as

(A3)\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] +8{\rm \pi} \lim_{r\to\infty}(r^3 \psi),\end{equation}

which demands the examination of the large-$r$ asymptotics of $K(r)$ and $\psi (r)$. We consider each in turn, with the goal of showing that both of the limits appearing in (A3) are zero.

A.1. Asymptotic of $K(r)$

Consider the triple correlator $S_{ijk}\equiv \langle u_i u_j u'_k\rangle$. We wish to estimate all of its time derivatives at $t=0$. Let us assume that $\boldsymbol {u} = 0$ at $t=0$. It follows immediately that $S_{ijk}= \partial _t S_{ijk} = \partial _t^2 S_{ijk}=0$ at $t=0$, while we have from the Navier–Stokes equation (2.7) that $\partial _t^3 S_{ijk}=\langle f_i f_j f_k'\rangle$, which, by A-I, decays rapidly in space (i.e. faster than any power law). Similarly,

(A4)\begin{align} \frac{\partial^4S_{ijk}}{\partial t^4} &= \left\langle \left(\nu \nabla^2 f_i + \frac{\partial f_i}{\partial t}\right)f_j f_k' \right\rangle \nonumber\\ &\quad + \left\langle f_i\left(\nu \nabla^2 f_j + \frac{\partial f_j}{\partial t}\right) f_k'\right\rangle + \left\langle f_i f_j\left(\nu \nabla'^2 f'_k + \frac{\partial f'_k}{\partial t}\right)\right\rangle, \end{align}

which also must decay rapidly with increasing $r$.

In contrast, the fifth derivative in time of $S_{ijk}$ is not guaranteed by A-I to decay rapidly with $r$ at $t=0$, because it contains terms of the sort $\langle f_i f_j \partial _t^2 \partial _k' p'\rangle$. The pressure field $p'$ is determined non-locally by (2.9) (note that pressure did not appear in the first four time derivatives of $S_{ijk}$ because, according to (2.9), both $\boldsymbol {\nabla } p$ and its first derivative in time are zero at $t=0$). Differentiating (2.9) twice with respect to time, we have

(A5)\begin{equation} \frac{\partial^2 p(\boldsymbol{x})}{\partial t^2} = \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'} \frac{\partial u_i'}{\partial t} \frac{\partial u_j'}{\partial t}=\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'} f_i' f_j' , \end{equation}

where the second equality is valid only at $t=0$. It follows that

(A6)\begin{equation} \frac{\partial^5 S_{ijk}}{\partial t^5} \sim \left\langle f_i f_j \frac{\partial}{\partial x'_k}\frac{\partial^2 p'}{\partial t^2} \right\rangle = \frac{1}{4{\rm \pi}} \frac{\partial}{\partial r_k} \int \frac{\mathrm{d}^3 \boldsymbol{x}''}{|\boldsymbol{x}'' -\boldsymbol{x}'|} \frac{\partial}{\partial x_l''}\frac{\partial}{ \partial x_m''} \langle f_l'' f_m'' f_i f_j \rangle.\end{equation}

Here, and in what follows, we use ‘$\sim$’ to relate expressions that have a common asymptotic scaling with ${r=|\boldsymbol {x}' - \boldsymbol {x}|}$ at large $r$. To determine the large-$r$ asymptotic of (A6), we make use of the fact that $\langle f_l'' f_m'' f_i f_j \rangle - \langle f_l'' f_m''\rangle \langle f_i f_j \rangle =\langle f_l'' f_m'' f_i f_j \rangle _c$ + $\langle f_l''f_i \rangle \langle f_m'' f_j \rangle$ + $\langle f_l'' f_j \rangle \langle f_m'' f_i \rangle$, and, therefore, has convergent integral moments, by A-I. This means that we can make use of the large-$r$ expansion

(A7)\begin{equation} \frac{1}{|\boldsymbol{x}''-\boldsymbol{x}|} = \frac{1}{|\boldsymbol{r} + \boldsymbol{s}|}= \frac{1}{r}+s_i \frac{\partial}{\partial r_i}\frac{1}{r}+s_i s_j \frac{\partial}{\partial r_i}\frac{\partial}{\partial r_j}\frac{1}{r} + O\left(\frac{1}{r^4}\right),\end{equation}

where $\boldsymbol {s} \equiv \boldsymbol {x}''-\boldsymbol {x}'$, to write (A6) as

(A8)\begin{equation} \frac{\partial^5 S_{ijk}}{\partial t^5} \sim \frac{\partial}{\partial r_k}\frac{\partial}{\partial r_l}\frac{\partial}{ \partial r_m} \frac{1}{r} \int \mathrm{d}^3 \boldsymbol{s} (\langle f_l'' f_m'' f_i f_j \rangle - \langle f_l'' f_m''\rangle \langle f_i f_j \rangle) + O\left(\frac{1}{r^5}\right).\end{equation}

The expansion is justified because A-I guarantees convergence of the higher-order (in $1/r$) terms.

Equation (A8) shows that the large-$r$ asymptotic of $\partial _t^5 S_{ijk}$ decays like $1/r^4$. Thus, we expect the large-$r$ asymptotic of $S_{ijk}$ to decay like $1/r^4$ for $t>0$, unless even slower-decaying terms are present in higher time derivatives of $S_{ijk}$. Let us check that they are not. Using Leibnitz’ theorem, the derivatives in question can be written as

(A9)\begin{align} \frac{\partial^n S_{ijk}}{\partial t^n} &= \sum_{m_1+m_2+m_3 = n} \frac{n!}{m_1!\,m_2!\,m_3!} \left\langle \frac{\partial^{m_1} u_i}{\partial t^{m_1}} \frac{\partial^{m_2} u_j}{\partial t^{m_2}} \frac{\partial^{m_3} u_k'}{\partial t^{m_3}}\right\rangle \nonumber\\ & \sim \sum_{m_1+m_2+m_3=n}\frac{n!}{m_1!\,m_2!\,m_3!}\left\langle \frac{\partial^{m_1-1} f_i}{\partial t^{m_1-1}} \frac{\partial^{m_2-1} f_j}{\partial t^{m_2-1}} \frac{\partial^{m_3-1}}{\partial t^{m_3-1}}\frac{\partial p'}{\partial x'_k}\right\rangle. \end{align}

In writing the final expression in (A9), i.e. in identifying that terms of this form are the ones that decay slowest with $r$ at large $r$, we used the facts that (i) correlators involving multiple instances of pressure will contain more powers of $1/r$ after the expansion (A7) is taken than those with a single instance, (ii) correlators of $f'_k$ (or its time derivatives) with a product of $f_i$ and $\partial _j \partial _t^2 p$ (or their time derivatives) are vanishingly small (this is because, after (A5) is used to eliminate pressure, the correlation functions to be integrated in $\boldsymbol {x}''$ contain both $f$ and $f'$. They are thus vanishingly small in the $r\to \infty$ limit, for any value of $\boldsymbol {x}''$) and (iii) the terms in $\partial _t^n S_{ijk}$ that can be formed by eliminating time derivatives of $\boldsymbol {u}$ in favour of spatial ones via $\partial _t u_i = \partial _j(u_i u_j)+\cdots$ (thus producing a higher-order correlator in the sense of number of correlated fields) may be seen straightforwardly to be the same or even higher order in the expansion in $1/r$. Substituting for pressure using (2.9) and for $|\boldsymbol {x}''-\boldsymbol {x}'|^{-1}$ using (A7), we have

(A10) \begin{align} \frac{\partial^n S_{ijk}}{\partial t^n} &\sim \sum_{m_1+m_2+m_3 = n} \frac{n!}{m_1!m_2!m_3!} \sum_{m_4+m_5=m_3-3}\frac{(m_3-3)!}{m_4!m_5!}\nonumber\\ &\quad\times\frac{\partial}{\partial r_k} \int \frac{\mathrm{d}^3 \boldsymbol{x}''}{|\boldsymbol{x}'' -\boldsymbol{x}'|} \frac{\partial}{\partial x_l''}\frac{\partial}{ \partial x_m''}\left\langle \frac{\partial^{m_1-1} f_i}{\partial t^{m_1-1}} \frac{\partial^{m_2-1} f_j}{\partial t^{m_2-1}} \frac{\partial^{m_4}f_l''}{\partial t^{m_4}} \frac{\partial^{m_5}f_m''}{\partial t^{m_5}}\right\rangle\ \nonumber\\ & \sim \sum_{\{m_{\alpha}\}} \frac{\partial}{\partial r_k} \frac{\partial}{\partial r_l}\frac{\partial}{ \partial r_m} \frac{1}{r} \int \mathrm{d}^3 \boldsymbol{s} \left( \left\langle \frac{\partial^{m_1-1} f_i}{\partial t^{m_1-1}} \frac{\partial^{m_2-1} f_j}{\partial t^{m_2-1}} \frac{\partial^{m_4}f_l''}{\partial t^{m_4}} \frac{\partial^{m_5}f_m''}{\partial t^{m_5}}\right\rangle \right.\nonumber\\ &\quad\left.\vphantom{\sum_{\{m_{\alpha}\}} \frac{\partial}{\partial r_k} \frac{\partial}{\partial r_l}\frac{\partial}{ \partial r_m}\left(\frac{1}{r}\right)}- \left\langle \frac{\partial^{m_1-1} f_i}{\partial t^{m_1-1}} \frac{\partial^{m_2-1} f_j}{\partial t^{m_2-1}} \right\rangle \left\langle\frac{\partial^{m_4}f_l''}{\partial t^{m_4}} \frac{\partial^{m_5}f_m''}{\partial t^{m_5}}\right\rangle \right) + O\left(\frac{1}{r^5}\right)\nonumber\\ &=O\left(\frac{1}{r^4}\right), \end{align}

where, after the second equality, we have for brevity suppressed the combinatoric factors and the explicit forms of the sums over $\{m_{\alpha }\}$ that appear in the first line. In moving from the final expression in (A9) to the first equality of (A10), we used the fact that the slowest-decaying correlators that are obtained by eliminating $\partial u''_i/ \partial t$ using the Navier–Stokes equation are the ones containing $f''_i$, rather than the inertial terms or the pressure gradient. This is because the latter correlators are higher order in $1/r$, for the same reasons as explained in points (ii) and (iii) above. We conclude from (A10) that $S_{ijk}=O(r^{-4})$ for $t>0$, from which it follows that ${u^3 K(r) = S_{xxx} = O(r^{-4})}$. The term containing $K(r)$ in (A3) is therefore zero.

A.2. Asymptotic of $\psi (r)$

To establish the large-$r$ scaling of $\psi$, we evaluate time derivatives at $t=0$ of $\langle u_i f_j'\rangle$ in the same manner as we have done for the triple correlator $S_{ijk}$. Long-range pressure-induced correlations appear first in the third derivative of $\varPsi _{ij}$ with respect to time, in terms such as

(A11)\begin{align} \frac{\partial^3 \varPsi_{ij}}{\partial t^3}&\sim \left\langle f_i \frac{\partial}{\partial x'_j}\frac{\partial^2 p'}{\partial t^2} \right\rangle = \frac{1}{4{\rm \pi}} \frac{\partial}{\partial r_j} \int \frac{\mathrm{d}^3 \boldsymbol{x}''}{|\boldsymbol{x}'' -\boldsymbol{x}'|} \frac{\partial}{\partial x_l''}\frac{\partial}{ \partial x_m''} \langle f_l'' f_m'' f_i \rangle \nonumber\\ &\sim \frac{\partial}{\partial r_j}\frac{\partial}{\partial r_l}\frac{\partial}{ \partial r_m} \frac{1}{r} \int \mathrm{d}^3 \boldsymbol{s} \langle f_l'' f_m'' f_i \rangle + O\left(\frac{1}{r^5}\right). \end{align}

The final integral in (A11) vanishes because $\partial \langle f_l'' f_m'' f_i \rangle /\partial r_i=0$. Therefore, $\partial _t^3 \varPsi _{ij} \sim O(r^{-5})$.

As we did in the previous section, we now check that there are no slower-decaying higher-order time derivatives of $\varPsi _{ij}$. The $n$th derivative of $\varPsi _{ij}$ with respect to time satisfies

(A12)\begin{equation} \frac{\partial^n\varPsi_{ij}}{\partial t^n} = \sum^{n}_{m=0}\frac{n!}{m!(n-m)!} \left\langle \frac{\partial^{m} f_i}{\partial t^{m}} \frac{\partial^{n-m} u_j'}{\partial t^{n-m}}\right\rangle \sim \sum^{n}_{m=0}\frac{n!}{m!(n-m)!} \left\langle \frac{\partial^{m} f_i}{\partial t^{m}} \frac{\partial^{n-m}}{\partial t^{n-m}}\frac{\partial p'}{\partial x'}\right\rangle. \end{equation}

The final expression follows because eliminating time derivatives of $u_i$ in favour of spatial derivatives by means of the Navier–Stokes equation introduces terms that are higher-order in $1/r$. An expansion analogous to (A11) then shows that $\partial _t^n\varPsi _{ij}= O(r^{-5})$, whence $\varPsi _{ij}= O(r^{-5})$. Thus, $\psi =O(r^{-5})$, and hence the term containing $\psi$ in (A3) is zero.

We conclude that $L$ is constant under the assumptions A-I and A-II stated at the start of this appendix.

Appendix B. Long-range correlations induced by finite-band forcing

In § 2, we approached the problem of the conservation of the Saffman integral (i.e. the coefficient of $k^2$ in the expansion (1.1) of $\mathcal {E}(k)$ as $k\to 0$) in terms of the strength of real-space correlations possessed by the forcing function. This is natural for turbulence occurring in nature or in a laboratory, where the spatial profile of the forcing may be measured, or prescribed. However, in numerical studies, forcing is commonly implemented in spectral space – in particular, it is common to force in a finite band of wavenumbers. Here, we examine the strength of correlations in real space induced by such a forcing.

As we have seen, the strongest long-range correlations in $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {f}' \rangle$ are typically induced by the last term on the right-hand side of (2.16), so that, according to (2.17),

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

The function $H(t,r)$ is related to the spectral function $F(k)$ by a Fourier transform: dropping explicit dependence on time,

(B2)\begin{equation} \frac{1}{r^2}\frac{\partial}{\partial r} [r^3 H(r)] = 2 \int^{\infty}_0 \mathrm{d} k F(k) \frac{\sin (kr)}{kr}, \end{equation}

whence it follows that

(B3)\begin{equation} H(r) = 2 \int^{\infty}_0\mathrm{d} k F(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3}.\end{equation}

Let us now suppose that $F(k)$ is non-zero only for $k_1 \leq k \leq k_2$. Under the assumption that $F$ is smooth, it is possible to show that $H(r)$ decays arbitrarily quickly (i.e. more quickly than any power law) at large $r$. This is intuitive, given the absence of any non-zero large-scale Fourier modes of $\boldsymbol {f}$. A proof is as follows.

Integrating (B3) by parts, we have

(B4)\begin{equation} H(r) = \frac{2}{r^3} \int_0^{\infty}\mathrm{d} k\frac{1}{k} \left(\frac{\partial}{\partial k}\frac{F}{k}\right)\sin(kr) = \frac{2}{r^3} \,\mathrm{Im} \left[ \int_0^{\infty}\mathrm{d} k \frac{1}{k}\left(\frac{\partial}{\partial k}\frac{F}{k}\right)\,{\rm e}^{{\rm i}kr}\right], \end{equation}

where the boundary terms vanish because $F$ and all its derivatives vanish at $0$ and $\infty$. If we integrate by parts a further $n$ times, boundary terms continue to vanish, leaving

(B5)\begin{equation} H(r) = \frac{2}{r^{3+n}} \,\mathrm{Im} \left[i^{n} \int_0^{\infty}\mathrm{d} k \left(\frac{\partial^n}{\partial k^n}\frac{1}{k}\frac{\partial}{\partial k}\frac{F}{k}\right)\,{\rm e}^{{\rm i}kr}\right], \end{equation}

so that

(B6)\begin{equation} | H(r) | = \frac{2}{r^{3+n}} \left| \int_0^{\infty}\mathrm{d} k \left(\frac{\partial^n}{\partial k^n}\frac{1}{k}\frac{\partial}{\partial k}\frac{F}{k}\right)\,{\rm e}^{{\rm i}kr} \right| < \frac{2}{r^{3+n}} \int_0^{\infty}\mathrm{d} k \left| \frac{\partial^n}{\partial k^n}\frac{1}{k}\frac{\partial}{\partial k}\frac{F}{k}\right| = \frac{C_n}{r^{3+n}}, \end{equation}

where $C_n$ is a finite constant. Thus, as $r\to \infty$, $| H(r) | \leq \mathrm {const.} \times r^{-m}$, for any $m$. QED.

However, $F(k)$ is not typically chosen to be smooth in numerical studies. Instead, $F(k)$ is often discontinuous at $k_1$ and $k_2$, with ${F(k< k_1)=F(k>k_2)=0}$. In this case, (B3) becomes, after integrating by parts,

(B7)\begin{equation} H(r) = \frac{2}{r^3}\left[\frac{F(k_1)\sin(k_1 r)}{k_1^2} -\frac{F(k_2)\sin(k_2r)}{k_2^2}\right] + \frac{2}{r^3} \int_{k_1}^{k_2}\mathrm{d} k\frac{1}{k} \left(\frac{\partial}{\partial k}\frac{F}{k}\right)\sin(kr).\end{equation}

The boundary term that has arisen in (B7) is dominant over the integral, in the limit $r \gg 1/k_1$, $1/k_2$, by the Riemann–Lebesgue lemma. Subdominant terms in the expansion may be computed by continuing to integrate by parts: at every order, the resulting terms will be of the form (oscillating part at frequency $k_1$ or $k_2$$\times r^{-n}$. Thus, the effect of the discontinuity in $F(k)$ at $k_1$ and $k_2$ is to induce correlations that, at large distances, oscillate about zero with component wavenumbers $k_1$ and $k_2$, and an amplitude decreasing as a power law in $r$.

While a forcing of this sort will not inject energy into large-scale modes directly (i.e. via (B1)), oscillatory behaviour will inevitably propagate into other correlators via the von Kármán–Howarth equation (2.12) and its higher-order analogues. It may then be the case that the oscillatory correlations induced in $\chi (r)$ decay more slowly with distance than the monotonically decaying ones implied by (B7). However, this should not affect the properties of the small-$k$ part of the energy spectrum, as it turns out that oscillatory behaviour in $\chi (r\to \infty )$, with the amplitude of oscillations decaying as a power law, always has negligible effect on $\mathcal {E}(k\to 0)$, independently of the power-law exponent.

To see why, let us suppose that

(B8)\begin{equation} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{u}' \rangle(r\geq r_0)= A\left[\frac{\sin(k_f r)}{r^n} + \frac{1}{r_0^n}\left(\frac{r_0}{r}\right)^m \right], \end{equation}

i.e. that the behaviour of $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ at $r$ larger than some $r_0$ is the sum of an oscillation with frequency $k_f\gg 1/r_0$, with an amplitude that decays as $r^{-n}$, and a monotonically decaying part that decays as $r^{-m}$, where $m$ and $n$ are both larger than $1$ for finite $\mathcal {E}(k)$ as $k \to 0$. According to (2.1), the contribution of the oscillating part of $\langle \boldsymbol {u} \,\boldsymbol {\cdot } \,\boldsymbol {u}' \rangle$ to the energy spectrum at small $k$ is

(B9)\begin{equation} \frac{Ak}{\rm \pi}\int_{r_0}^{\infty}\mathrm{d} r \frac{1}{r^{n-1}} \sin(k_f r) \sin (kr) = \frac{Ak^2}{{\rm \pi} k_f r_0^{n-2}} \cos(k_f r_0)\left[ 1 + O\left(\frac{1}{k_f r_0}\right) \right],\end{equation}

where we have integrated $\sin (k_f r)$ by parts and applied the Riemann–Lebesgue lemma. Meanwhile, the contribution of the monotonically decaying part is

(B10)\begin{equation} \frac{A k r_0^{m-n}}{\rm \pi}\int_{r_0}^{\infty}\mathrm{d} r \frac{1}{r^{m-1}} \sin (kr) = \frac{A F_m }{{\rm \pi} k r_0^{n}} (k r_0)^{\min[m,3]} [ 1 + O(k r_0)], \end{equation}

where $F_m$ is a positive number that depends on $m$. The equality in (B10) is non-trivial: it is proven by direct application of asymptotic methods to the integral, with a case distinction for $m<3$ and $m \geq 3$, or else via evaluation of the integral directly (resulting in a hypergeometric function), and taking the limit of small $k r_0$.

The right-hand side of (B10) is always large compared with that of (B9); for $m > 3$, which is likely the only case of physical significance, the difference is $k_f r_0\gg 1$. This confirms our assertion that decaying oscillatory behaviour in $\chi (r\to \infty )$ has negligible effect on $\mathcal {E}(k\to 0)$: the latter is always dominated by any part of $\chi (r\to \infty )$ that decays monotonically. In conclusion, therefore, while forcing with a discontinuous spectrum does induce long-range correlations (as it must, owing to the sharp features induced in the spectrum of $\boldsymbol {u}$), these are likely of little significance to the dynamics of the large scales.

Appendix C. Asymptotic form of $\chi (r)$ for turbulence with a power-law energy spectrum

In this appendix, we present a derivation of (3.8), which gives the longitudinal correlation function, $\chi (r)$, corresponding to turbulence with spectrum

(C1)\begin{equation} \mathcal{E}(k) = C k^a, \end{equation}

at scales $k_1 \leq k \leq k_2$, for values of $r$ satisfying $k_2^{-1}\ll r \ll k_1^{-1}$. The constant $C$ may be expressed in terms of the total energy contained at scales $k_1 \leq k \leq k_2$: since

(C2)\begin{equation} \int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) \simeq \begin{cases} C k_1^{a+1}/|1+a| & \text{if }a<{-}1, \\ C \ln{(k_2/k_1)} & \text{if }a={-}1, \\ C k_2^{a+1}/(1+a) & \text{if }a>{-}1, \end{cases} \end{equation}

we have

(C3)\begin{equation} C = \int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) \times \begin{cases} |1+a| k_1^{{-}1-a} & \text{if }a<{-}1, \\ 1/\ln{(k_2/k_1)} & \text{if }a={-}1, \\ (1+a) k_2^{{-}1-a} & \text{if }a>{-}1. \end{cases}\end{equation}

Let us start from (3.7),

(C4)\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}

and consider first the part of the integral with $k< k_1$, which we denote $J_1$. A straightforward Taylor expansion in $kr \ll 1$ yields

(C5)\begin{equation} J_1 \equiv 2 \int^{k_1}_0 \mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3} = \frac{2}{3} \int_0^{k_1} \mathrm{d} k \mathcal{E}(k) + O[(k_1 r)^3]. \end{equation}

Then (C4) becomes

(C6)\begin{equation} u^2 \chi(r) = \frac{2}{3} \int_0^{k_1} \mathrm{d} k \mathcal{E}(k) + J_2 + J_3 + O[(k_1 r)^3],\end{equation}

where $J_2$ and $J_3$ are integrals that correspond to the parts of (C4) with $k_1< k< k_2$ and $k>k_2$, respectively, so that

(C7)\begin{equation} J_2 + J_3 = 2 \int^{\infty}_{k_1} \mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3}.\end{equation}

Depending on the particular value of $a$, it can be convenient to integrate (C7) by parts before choosing precise definitions for $J_2$ and $J_3$. This allows the boundary term at $k = k_2$ to be redistributed between $J_2$ and $J_3$. For some values of $a$ in (C1), it will then be possible to show that $J_2$ is the dominant contribution to (C6), allowing the dominant part of $\chi (r)$ to be computed despite our incomplete knowledge of $\mathcal {E}(k)$. We proceed by considering different ranges of $a$ in turn, defining and computing $J_2$ and $J_3$ for each case. The final asymptotic expressions for $\chi (r)$ will be assembled using (C6) at the end of this appendix.

C.1. Case of $a<2$

C.1.1. Calculation of $J_2$

In this case, it is unnecessary to integrate by parts. We define

(C8) \begin{equation} J_2 \equiv 2 \int^{k_2}_{k_1} \mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(k r)^3} = 2 C r^{{-}a-1}\int^{k_2 r}_{k_1 r} \mathrm{d}\kern 0.06em x x^{a} \frac{\sin x-x\cos x}{x^3}.\end{equation}

For $-1< a<2$, (C8) is convergent for $k_1 r \to 0$, $k_2 r \to \infty$: to leading order in $k_1 r$ and $1/ k_2 r$

(C9)\begin{equation} J_2 ={-} 2 C r^{{-}a -1}(a-1)\varGamma(a-2)\sin \frac{a {\rm \pi}}{2},\quad -1< a<2.\end{equation}

Therefore, $J_2$ is large compared with $J_1\sim Ck_1^{1+a}$ by a factor of $(k_1 r)^{-1-a}\gg 1$.

For $a<-1$, the integral in (C8) is divergent at the lower limit. Its leading-order asymptotic as $k_1 r \to 0$ is

(C10)\begin{equation} J_2 = \frac{2Ck_1^{1+a}}{3|1+a|}\{1+O[(k_1 r)^2]\},\quad a<{-}1,\end{equation}

which is the same size as $J_1\sim Ck_1^{1+a}$, and is independent of $r$.

In the particular case of $a=-1$, taking the leading-order asymptotic of (C8), we have

(C11)\begin{equation} J_2 = \frac{2C}{3}[\ln{(k_1 r)}+O(1)],\quad a={-}1.\end{equation}

Despite appearances, $J_2$ does not blow up as $k_1 r \to 0$, because $C\to 0$ as $k_1 r\to 0$ in order for the total energy to be finite, as (C3) shows. Nonetheless, $J_2$ dominates over $J_1\sim C$ in (C6), by a factor of $\ln (k_1 r)$.

C.1.2. Calculation of $J_3$

In $J_3$, we define $x = k/k_2$ and integrate by parts

(C12)\begin{align} J_3 & \equiv 2 \int^{\infty}_{k_2} \mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(k r)^3} \nonumber\\ & ={-} 2 k_2 \left[ \mathcal{E}(k_2 x) \frac{\cos(k_2 r x)+k_2 r x\sin(k_2 r x)}{(k_2 r x)^4} \right]^{\infty}_{1} \nonumber\\ & \quad + 2 k_2 \int^{\infty}_{1} \mathrm{d}\kern 0.06em x \left[\frac{\partial}{\partial x} \frac{\mathcal{E}(k_2 x)}{(k_2 r x)^3} \right] \frac{\cos(k_2 r x)+k_2 r x\sin(k_2 r x)}{k_2 r x}. \end{align}

The boundary term at $\infty$ is exponentially small, by assumption, while the remaining integral is small compared with the boundary term by a factor of $(k_2 r)^{-1} \ll 1$, by the Riemann–Lebesgue lemma. This leaves

(C13)\begin{align} J_3 & = 2 k_2 \mathcal{E}(k_2 ) \frac{\cos(k_2 r )+k_2 r\sin(k_2 r)}{(k_2 r)^4} + O[ (k_2 r)^{{-}1} ] \end{align}
(C14)\begin{align} & = 2 C r^{{-}a -1 } (k_2 r)^{a-3} [\cos(k_2 r )+k_2 r\sin(k_2 r)] + O[ (k_2 r)^{{-}1} ], \end{align}

which is manifestly small compared with $J_2$ for all $a<2$.

C.2. Case of $2 \leq a \leq 3$

In this range of $a$, it is more convenient to use integration by parts in (C7) before splitting the integration domain. We have

(C15)\begin{equation} 2 \int^{\infty}_{k_1}\mathrm{d} k \mathcal{E}(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3} = 2 C (k_1 r)^{a-1} r^{{-}1 - a} + 2\int^{\infty}_{k_1}\mathrm{d} k \left[\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{kr}\right]\frac{\sin (kr)}{kr^2}, \end{equation}

where the boundary term at $\infty$ vanishes by assumption.

C.2.1. Calculation of $J_2$

We take

(C16)\begin{align} J_2 & \equiv 2 C (k_1 r)^{a-1} r^{{-}1 - a} + 2\int^{k_2}_{k_1} \mathrm{d} k \left[ \frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{kr}\right]\frac{\sin (kr)}{kr^2} \nonumber\\ & = 2 C (k_1 r)^{a-1} r^{{-}1 - a} + 2 Cr^{{-}1-a}(a-1)\int^{k_2r}_{k_1r} {{\rm d} x} x^{a-3} \sin x. \end{align}

When $2\leq a<3$, the integral is convergent for $k_1 r \to 0$, $k_2 r \to \infty$, so

(C17)\begin{equation} J_2 ={-} 2 C r^{{-}a -1}(a-1)\varGamma(a-2)\sin \frac{a {\rm \pi}}{2},\end{equation}

to leading order, the boundary term being small by a factor of $(k_1 r)^{a-1}\ll 1$. As before, $J_2$ is large compared with $J_1\sim Ck_1^{1+a}$ by a factor of $(k_1 r)^{-1-a}\gg 1$.

If $a=3$, then (C16) does not converge as $k_2 r \to \infty$. Instead, we have

(C18)\begin{equation} J_2 = 4 Cr^{{-}4}[1-\cos(k_2 r)],\end{equation}

to leading order in $k_1 r\ll 1$.

C.2.2. Calculation of $J_3$

Taking $J_3$ to be the part of (C15) not included in $J_2$ as defined in (C16), integrating by parts and applying the Riemann–Lebesgue lemma as in (C12) gives

(C19)\begin{align} J_3 & \equiv 2\int^{\infty}_{k_2} \mathrm{d} k \left[\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{kr}\right]\frac{\sin (kr)}{kr^2} \nonumber\\ & = 2 (a-1) C r^{{-}1-a}(k_2 r)^{a-3} \cos (k_2r) + O[ (k_2 r)^{{-}1} ]. \end{align}

Comparison with (C17) shows that $J_3$ is small compared with $J_2$ for $a<3$. If $a=3$, then the leading-order part of $J_3$ cancels with the term proportional to $\cos (k_2 r)$ in (C18). The remaining term in (C18) is precisely (C17) with $a=3$, therefore we conclude that (C17) is valid for $2 \leq a \leq 3$, and provides the leading-order part of (C6) in this range.

C.3. Case of $a > 3$

The procedure for $a>3$ is similar to the one followed in § C.2: we continue to integrate (C15) by parts, split the resulting integral into pieces with $k< k_2$ and $k>k_2$, and use the Riemann–Lebesgue lemma to infer that the $k>k_2$ piece is subdominant. While a different number of integrations by parts will be required depending on the particular value of $a$, let us treat them all simultaneously, and integrate the integral appearing in (C15) by parts $n$ times. Using the complex representation of trigonometric terms for convenience, we have

(C20)\begin{align} 2\int^{\infty}_{k_1}\mathrm{d} k \left[\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{kr}\right]\frac{\sin (kr)}{kr^2} & = \frac{2}{r^3}\,\mathrm{Im} \left\{ \sum_{m=1}^n ({-}1)^{m-1}\left[\frac{{\rm e}^{{\rm i}kr}}{(ir)^m}\frac{\partial^{m-1}}{\partial k^{m-1}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k}\right]^{\infty}_{k_1}\right. \nonumber\\ & \quad \left.+({-}1)^n\int^{\infty}_{k_1} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{(ir^n)} \frac{\partial^{n}}{\partial k^{n}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k} \right\} \end{align}
(C21)\begin{align} & = 2\,\mathrm{Im} \left[ \sum_{m=1}^n i^{m}\,{\rm e}^{{\rm i} k_1 r} \frac{\varGamma (a-2)}{\varGamma (a-m-1)} (k_1 r)^{a-m-2} C r^{{-}a-1} \right. \nonumber\\ & \quad \left.+\frac{({-}1)^n}{r^3}\int^{\infty}_{k_1} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{(ir)^n} \frac{\partial^{n}}{\partial k^{n}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k} \right]. \end{align}

The last expression was obtained by simplifying the boundary term using the fact that $\mathcal {E}(k)= Ck^a$ around $k=k_1$.

C.3.1. Calculation of $J_2$

As before, we split (C15) into two components, $J_2$ and $J_3$. Formally, $J_2$ should be defined to include the $k=k_1$ boundary terms in (C15) and (C21). However, these are all small compared with $Cr^{-a-1}$ as long as we choose $n< a-2$. Doing so, we have

(C22)\begin{align} J_2 & = 2 \,\mathrm{Im}\left[ \frac{({-}1)^n}{r^3}\int^{k_2}_{k_1} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{(ir)^n} \frac{\partial^{n}}{\partial k^{n}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k}\right] \{ 1+O[(k_1 r)^{a-n-2}]\}\nonumber\\ & \simeq 2\frac{\varGamma(a-2)}{\varGamma(a-n-2)}(a-1)C r^{{-}1-a} \,\mathrm{Im} \left( i^n \int^{k_2r}_{k_1r}\mathrm{d}\kern 0.06em x \,{\rm e}^{{\rm i}x} x^{a-n-3} \right). \end{align}

Now, if $-1< a-n-3<0$ (which is consistent with our earlier choice of $n< a-2$), the integral in (C22) is convergent as we take $k_1 r \to 0$, $k_2 r \to \infty$. In that case, (C22) becomes the now-familiar

(C23)\begin{equation} J_2 ={-} 2 C r^{{-}a -1}(a-1)\varGamma(a-2)\sin \frac{a {\rm \pi}}{2}.\end{equation}

It follows that (C23) is valid for all non-integer $a>3$, since then it is always possible to choose $n$ such that $-1< a-n-3<0$.

If $a$ is an integer, we choose $n = a -3$ instead. In this case, (C22) becomes

(C24)\begin{equation} J_2 = 2(a-1)\varGamma(a-2)Cr^{{-}1-a} \,\mathrm{Im} \{ i^{a} [{\rm e}^{{\rm i}k_2r}-1+O(k_1 r)]\}.\end{equation}
C.3.2. Calculation of $J_3$

Taking $J_3$ to be the part of (C15) not included in $J_2$ and using (C21), we have

(C25)\begin{align} J_3 & \equiv 2 \,\mathrm{Im}\left[ \frac{({-}1)^n}{r^3}\int^{\infty}_{k_2} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{(ir)^n} \frac{\partial^{n}}{\partial k^{n}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k}\right] \nonumber\\ & = 2\,\mathrm{Im} (i^{n+1}\,{\rm e}^{{\rm i}k_2r})\frac{\varGamma(a-2)}{\varGamma(a-n-2)}(a-1)Cr^{{-}1-a}(k_2 r)^{a-n-3}\nonumber\\ &\quad + 2 \,\mathrm{Im}\left[ \frac{i^{n+1}}{r^3}\int^{\infty}_{k_2} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{r^{n+1}} \frac{\partial^{n+1}}{\partial k^{n+1}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k}\right]. \end{align}

In the case of non-integer $a$, our choice of $a-n-3<0$ ensures the boundary term in (C25) is small compared with $J_2$, while the integral is lower order still, by the Riemann–Lebesgue lemma. For integer $a$, our choice of $n=a-3$ means that the boundary term in (C25) cancels with the term that is proportional to $\mathrm {Im} (i^{a}\,{\rm e}^{{\rm i}k_2 r})$ in (C24). The remaining term in (C24) is precisely the right-hand side of (C23), which dominates over the remaining integral in (C25) (by the Riemann–Lebesgue lemma), provided that $a$ is odd.

However, special care must be taken when $a$ is an even integer greater than $3$. After the term proportional to $\mathrm {Im} (i^{a}\,{\rm e}^{{\rm i}k_2 r})$ in (C24) cancels with its partner in (C25), the other term in (C24) is proportional to $\mathrm {Im}(i^a)$, which vanishes if $a$ is even. This means that $J_2 = Cr^{-1-a}\times O(k_1 r)$. The leading-order non-vanishing term in the asymptotic expansion in $k_2 r\gg 1$ is contained within the integral in (C25). However, unlike for other values of $a$, this term cannot be extracted from (C25) by any number of integrations by parts, as the boundary terms so generated always vanish. This is because the desired leading-order contribution to $\chi (r)$ comes from the part of $\mathcal {E}(k)$ with $k>k_2$, where the form of $\mathcal {E}(k)$ is unknown. At best, we can constrain the dependence on $r$ by noting that, for $n=a-3$,

(C26)\begin{equation} \left| 2\, \mathrm{Im}\left[ \frac{i^{n+1}}{r^3}\int^{\infty}_{k_2} \mathrm{d} k \frac{{\rm e}^{{\rm i}kr}}{r^{n+1}} \frac{\partial^{n+1}}{\partial k^{n+1}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k}\right] \right| < 2r^{{-}1-a} \int^{\infty}_{k_2} \mathrm{d} k \left| \frac{\partial^{a-2}}{\partial k^{a-2}}\frac{1}{k}\frac{\partial}{\partial k}\frac{\mathcal{E}(k)}{k} \right|. \end{equation}

The integral on the right-hand side of (C26) is independent of $r$. It has the same dimensions as $C$, and typically will be $\sim C$, in which case the integral in (C25) can be large compared with the $O(k_1 r)$ part of $J_2$, and therefore it provides the dominant contribution to $\chi (r)$ in (C6). Thus, for even integer $a>3$, we are only able to conclude that $\chi (r)$ decays as $Cr^{-1-a}$ or faster as $r\to \infty$.

C.4. The leading-order correction in $k_1 r\ll 1$

The conclusion that $\chi (r\to \infty )\leq O(Cr^{-1-a})$ may fail if the contribution to $\chi (r)$ from $\mathcal {E}(k>k_2)$ is much smaller than the upper limit enforced by (C26), with the result that the strongest long-range correlations are instead determined by the $k< k_1$ part of the spectrum. An example of this is the superposition of a small-scale velocity field with exponentially decaying correlations and a second velocity field at much larger scales. On scales longer than those of the former field, correlations are dominated by the latter, so $\chi (r)\sim \mathrm {const.}$ if $r$ is small compared with the characteristic scale of the large-scale field.

A convenient way to obtain this result formally, without the need to keep track of higher-order terms in the expansion in $k_1 r$ above, is to define an auxiliary spectrum

(C27)\begin{equation} \mathscr{E}(k) = \begin{cases} C k^a & \text{if }k< k_1 ,\\ \mathcal{E}(k) & \text{otherwise}, \end{cases} \end{equation}

i.e. $\mathscr {E}(k)$ is the spectrum obtained by extending the power-law behaviour of $\mathcal {E}(k)$ at $k_1< k< k_2$ to $k< k_1$. Let us denote by $\chi _{\mathscr {E}}(r)$ the function obtained by replacing $\mathcal {E}(k)$ by $\mathscr {E}(k)$ in (C4) while retaining the definition of $u^2$. In this notation, $\chi (r)=\chi _{\mathcal {E}}(k)$ is the correlation function that we have been concerned with thus far. Then, $\chi _{\mathscr {E}}$ may be computed by simply setting $k_1 = 0$ in the calculation above:

(C28)\begin{equation} \chi_{\mathscr{E}}(r) = \chi_{\mathcal{E}}(r)|_{k_1=0}. \end{equation}

Now we define ${\rm \Delta} (k) = \mathcal {E}(k) - \mathscr {E}(k)$ to be the spectrum of the ‘excess energy’ contained at large scales (which can be negative). Owing to the linearity of (C4) in $\mathcal {E}$, we have

(C29)\begin{equation} \chi_{\mathcal{E}}(r) = \chi_{\mathscr{E} + \varDelta}(r) = \chi_{\mathscr{E}}(r) + \chi_{\varDelta}(r) = \chi_{\mathcal{E}}(r)|_{k_1=0} + \chi_{\varDelta}(r). \end{equation}

Thus, the finite-$k_1 r$ correction to the correlation function corresponding to $\mathcal {E}(k)$ is exactly the function obtained by replacing $\mathcal {E}(k)$ by ${\rm \Delta} (k)$ in (C4) (again, while retaining the definition of $u^2$). This is easily computed to leading order in $k_1 r$: since ${\rm \Delta} (k>k_1) = 0$,

(C30)\begin{equation} u^2 \chi_{\varDelta}(r) \equiv 2 \int^{k_1}_0\mathrm{d} k {\rm \Delta}(k) \frac{\sin(kr)-kr\cos(kr)}{(kr)^3} = \frac{2}{3} \int_0^{k_1} \mathrm{d} k {\rm \Delta}(k) + O[(k_1 r)^3].\end{equation}

As anticipated above, the leading-order correction is $\chi _{\varDelta }(r)\sim \mathrm {const.}$. The size of the correction depends on the amount of ‘excess energy’ contained at $k< k_1$, i.e. on the difference between the actual energy and the energy that would be present if the power law $\mathcal {E}(k)\propto Ck^a$ extended to $k< k_1$. We note that this excess energy can be negative, in which case $\chi _{\varDelta }(r)$ will also be negative, resulting in a tendency to introduce anti-correlations.

It should, however, be noted that this discussion is largely academic: while it is formally possible that the leading-order term in the $k_1 r$ expansion, (C30), will dominate $\chi (r)$, this is an artificial situation, because correlations do not fall off exponentially quickly in real turbulence (Batchelor & Proudman Reference Batchelor and Proudman1956). Instead, $\chi (r)$ typically decays as $r^{-6}$ in isotropic Batchelor turbulence (see Davidson Reference Davidson2013 for a discussion). In the general case, it should therefore be expected that the integral in (C25) provides the dominant part of $\chi (r)$.

C.5. Final expression

Assembling all the results derived in this appendix, and using (C3) to eliminate $C$, (C6) becomes (3.8) (in the case of $a = 2,4,6,\ldots$, we assume that the dominant correlations come from the part of $\mathcal {E}(k)$ with $k>k_2$; see the discussion in § C.4).

Appendix D. Evolution of mean square momentum in forced turbulence

In this appendix, we examine the evolution of $\langle \boldsymbol {P}_V^2 \rangle$ under the forced Navier–Stokes equation, as Davidson (Reference Davidson2015) did for decaying turbulence.

Ignoring viscous forces, the evolution of the total momentum contained within a volume $V$ is given by

(D1)\begin{equation} \frac{\mathrm{d} \boldsymbol{P}_V}{\mathrm{d} t} ={-}\int_{\partial V} \boldsymbol{u} (\boldsymbol{u} \boldsymbol{\cdot} \mathrm{d} \boldsymbol{S})-\int_{\partial V} p\mathrm{d} \boldsymbol{S} + \int_V \mathrm{d}^3\boldsymbol{x} \boldsymbol{f}.\end{equation}

The three terms appearing on the right-hand side are identified straightforwardly as the advection of momentum out of $V$, the net pressure force on $V$, and the net external forcing. Therefore,

(D2)\begin{equation} \frac{\mathrm{d} \boldsymbol{P}_V^2}{\mathrm{d} t} = 2\int_V \mathrm{d}^3 \boldsymbol{x}' \boldsymbol{u}' \boldsymbol{\cdot} \left[-\int_{\partial V} \boldsymbol{u} (\boldsymbol{u} \boldsymbol{\cdot} \mathrm{d} \boldsymbol{S})-\int_{\partial V} p\mathrm{d} \boldsymbol{S}+\int_V \mathrm{d}^3\boldsymbol{x} \boldsymbol{f}\right].\end{equation}

For simplicity, we shall assume that the correlation time of the forcing is short compared with the eddy-turnover time. In that case, (2.16) and (2.17) give

(D3)\begin{equation} \langle \boldsymbol{u} \,\boldsymbol{\cdot} \,\boldsymbol{f}' \rangle = \frac{1}{r^2}\frac{\partial}{\partial r}[r^3 H(t,r)],\end{equation}

while isotropy also demands (cf. (2.2))

(D4)\begin{equation} \langle u_i f_j'\rangle = \frac{1}{2r}[(r^2 H)' \delta_{ij} - H'(r) r_i r_j].\end{equation}

Taking an ensemble average of (D2), using (D3) and (D4), and restricting attention to spherical $V$ with radius $R$, one can show that

(D5)\begin{align} \frac{\mathrm{d} \langle \boldsymbol{P}_V^2 \rangle}{\mathrm{d} t} &= 4{\rm \pi}^2 R^2u^3 \int_0^{2R} \mathrm{d} r \left[1-\left(\frac{r}{2R}\right)^2\right] \frac{1}{r}\frac{\partial}{\partial r}(r^4 K) \nonumber\\ &\quad + 8{\rm \pi}^2 R^2 \int_0^{2R} \mathrm{d} r \left[1-\left(\frac{r}{2R}\right)^2\right]r^3 H. \end{align}

The derivation of this equation is closely analogous to the one presented in Davidson (Reference Davidson2015) for $\mathrm {d} \langle \boldsymbol {P}_V^2 \rangle / \mathrm {d} t$ in decaying turbulence, to which we refer the reader for details.

Equation (D5) shows that there are two relevant processes that can change the expectation value of the squared linear momentum contained in a volume of size $V$. The first, represented by the first term on the right-hand side, encodes the effect of the advection of momentum by the flow, the second the injection of momentum by the forcing. Both terms are at most $O(R^2)$ as $R\to \infty$, as long as $K(r\to \infty )=o(r^{-3})$ and $H(r\to \infty )=o(r^{-3})$. This makes sense, given that both effects are surface processes (as long as long-range correlations in the forcing are weak). This means that neither term can spontaneously generate $\langle \boldsymbol {P}_V^2 \rangle \propto R^3$ for arbitrarily large $R$, as is consistent with the conclusion of § 2.

In fact, the net effect of the two terms on the right-hand side of (D5) may be a scaling of $\mathrm {d} \langle \boldsymbol {P}_V^2 \rangle /\mathrm {d} t$ vs $R$ even weaker than $R^2$, because there can be partial cancellation between them, or else because $\partial (r^4 K)/\partial r$ may change sign at some $r$. Of course, some cancellation between the two terms is inevitable: the forcing cannot perpetually increase $\langle \boldsymbol {P}_V^2 \rangle$ unchecked, because, even if turbulence could maintain very short-range correlations and so not grow a $k^2$ spectrum, there would still be a cascade of energy to smaller scales (encoded in the first term) resulting in the destruction of the local structures. Nonetheless, it appears robustly the case that the net size of the right-hand side of (D5) must scale at most as $R^2$ for large $R$.

Let us then consider some $R\gg l$, and suppose that, indeed, the aggregate of the terms on the right-hand side of (D5) scales as $R^2$. In that case, on dimensional grounds,

(D6)\begin{equation} \frac{\mathrm{d} \langle \boldsymbol{P}_V^2 \rangle}{\mathrm{d} t} \sim u^3 R^2 l^3. \end{equation}

Suppose, roughly, that $\langle \boldsymbol {P}_V^2 \rangle$ increases according to (D6) until it saturates. Saturation occurs when $\langle \boldsymbol {P}_V^2 \rangle \sim u^2 R^3 l^3$ – this must be the case on dimensional grounds, because the $R$ dependence is fixed by (3.10) with $a=2$. We then find that the time taken for saturation at scale $R$ is $t_{sat}\sim R/u$. Equivalently, the scale $R_c$ at which the growth of $\langle \boldsymbol {P}_V^2 \rangle$ has just saturated at time $t$ is given by

(D7)\begin{equation} R_c \sim u t. \end{equation}

That this should be the strongest-allowed scaling of $R_c$ vs $t$ is intuitive: (D7) simply represents the limit on the growth of $R_c$ imposed by causality. At distances greater than $R_c \sim u t$, two points in the flow cannot have exchanged momentum, because there has not been enough time for local processes, namely, advection at speed $u$ (or else the cumulative effect of local forces acting on scales $\sim l$ with time scales $\sim l/u$), to act.

Appendix E. Derivation of the passive-momentum equations

In this appendix, we show how the mode-coupling equation (4.12) is obtained under the assumptions explained in § 4.2. Our goal is to compute the evolution of the spectrum $\mathcal {E}_{\boldsymbol {w}}(t, k)$ of the passive vector field $\boldsymbol {w}$ satisfying (4.6). To do this, we shall derive an evolution equation for the correlation function $\langle w_i (t, \boldsymbol {k}) w_j (t, \boldsymbol {k}')\rangle$.

E.1. Homogeneous and isotropic forms of relevant correlators

Let us first note the restrictions imposed by symmetries on the various correlators involved. Due to statistical homogeneity of $\boldsymbol {w}$, this function is restricted to satisfy

(E1)\begin{equation} \langle w_i (t, \boldsymbol{k}) w_j (t, \boldsymbol{k}')\rangle = (2{\rm \pi})^3 \delta ( \boldsymbol{k} + \boldsymbol{k}') \varPsi_{ij}(t, \boldsymbol{k}). \end{equation}

The form of the tensor $\varPsi _{ij}(t, \boldsymbol {k})$ is further restricted by isotropy and incompressibility

(E2)\begin{equation} \varPsi_{ij}(t, \boldsymbol{k}) = \tfrac{1}{2}\varPsi(t, k) \mathcal{P}_{ij} (\boldsymbol{k}), \end{equation}

where

(E3)\begin{equation} \mathcal{P}_{ij}(\boldsymbol{k}) = \delta_{ij} - \frac{k_i k_j}{k^2}, \end{equation}

is the usual projection operator onto the plane perpendicular to $\boldsymbol {k}$, and the isotropic function $\varPsi (t, k)$ is related to $\mathcal {E}_{\boldsymbol {w}}(t, k)$ via

(E4)\begin{equation} \varPsi(t, k) = 2{\rm \pi}^2 \frac{\mathcal{E}_{\boldsymbol{w}}(t, k)}{k^2}. \end{equation}

In order to compute (E1), we require similar correlation functions for $\boldsymbol {u}$ and $\boldsymbol {f}_{\boldsymbol {w}}$. In particular, we shall need $\langle u_i(t, \boldsymbol {k}) u_j (t', \boldsymbol {k}')\rangle$, whose general form is also restricted by homogeneity and statistical invariance of $\boldsymbol {u}$ in time:

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

As explained in § 4.2, we shall take $\boldsymbol {u}$ to have zero correlation time:

(E6)\begin{equation} \kappa_{ij}(\boldsymbol{k}, t-t') = \kappa_{ij}(\boldsymbol{k})\delta(t-t'). \end{equation}

Together, incompressibility and isotropy further imply

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

As explained in § 4.2, we shall assume $\boldsymbol {u}$ to have a single scale, so that the isotropic function $\kappa (k)$ is

(E8)\begin{equation} \kappa (k) = \kappa_0 \delta (k-k_f). \end{equation}

Finally, we shall require $\langle f_i(t, \boldsymbol {k}) f_j (t', \boldsymbol {k}')\rangle$, which satisfies

(E9)\begin{equation} \langle f_i(t, \boldsymbol{k}) f_j (t', \boldsymbol{k}')\rangle = \tfrac{1}{2}(2{\rm \pi})^3 \delta (\boldsymbol{k} + \boldsymbol{k}') \delta (t-t')\varPhi(k)\mathcal{P}_{ij}(\boldsymbol{k}),\end{equation}

for the same reasons as the other two correlation functions.

E.2. Derivation of the mode-coupling equation

Equation (4.6) reads, after dropping the subscript $\boldsymbol {w}$ from $p_{\boldsymbol {w}}$ and $\boldsymbol {f}_{\boldsymbol {w}}$ (there can be no confusion with the corresponding fields for $\boldsymbol {u}$, as $\boldsymbol {u}$ is now prescribed artificially),

(E10)\begin{equation} \partial_t w_i ={-} \partial_j (u_j w_i) - \partial_i p + \nu \nabla^2 w_i + f_i. \end{equation}

In Fourier space, this is

(E11)\begin{equation} \partial_t w_i(\boldsymbol{k}) + \nu k^2 w_i(\boldsymbol{k}) ={-} i k_l \mathcal{P}_{ip}(\boldsymbol{k})\int \frac{\mathrm{d}^3 \boldsymbol{k}'}{(2{\rm \pi})^3} u_l(\boldsymbol{k}') w_p(\boldsymbol{k} - \boldsymbol{k}') + f_i(\boldsymbol{k}), \end{equation}

where we have taken $\boldsymbol {f}$ to be solenoidal – there is no loss of generality here, as if one is interested in a non-solenoidal forcing, one can interpret $\boldsymbol {f}$ to be its solenoidal part. Then,

(E12)\begin{align} &\partial_t \langle w_i(\boldsymbol{k}) w_j(\boldsymbol{k}') \rangle + \nu \left( k^2+k'^2 \right) \langle w_i(\boldsymbol{k}) w_j(\boldsymbol{k}') \rangle \nonumber\\ &\quad ={-} i \int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} [ k_l \mathcal{P}_{ip}(\boldsymbol{k}) \langle u_l(\boldsymbol{k}'') w_p(\boldsymbol{k} - \boldsymbol{k}'') w_j(\boldsymbol{k}') \rangle\nonumber\\ &\quad + k'_l \mathcal{P}_{jp}(\boldsymbol{k}') \langle u_l(\boldsymbol{k}'') w_p(\boldsymbol{k}' - \boldsymbol{k}'') w_i(\boldsymbol{k}) \rangle]\nonumber\\ &\qquad +\langle f_i(\boldsymbol{k})w_j(\boldsymbol{k}')\rangle+\langle f_j(\boldsymbol{k}')w_i(\boldsymbol{k})\rangle. \end{align}

In order to simplify the correlators appearing on the right-hand side of (E12), we can make use of the zero-correlation-time assumption for $\boldsymbol {u}$ and $\boldsymbol {f}$. In the latter case, we note that integrating (E11) over time, multiplying by $f_i(\boldsymbol {k})$ and ensemble averaging yields

(E13)\begin{align} \langle w_i(\boldsymbol{k})f_j(\boldsymbol{k}')\rangle & = \int^t\mathrm{d} t' \left[ -\nu k^2 \langle w_i(\boldsymbol{k}, t')f_j(\boldsymbol{k}',t) \rangle \vphantom{\frac{1}{2}}\right.\nonumber\\ & \quad -\left. i k_l \mathcal{P}_{ip}(\boldsymbol{k})\int \frac{\mathrm{d}^3 \boldsymbol{k}'}{(2{\rm \pi})^3} \langle u_l(\boldsymbol{k}',t') w_p(\boldsymbol{k} - \boldsymbol{k}',t')f_j(\boldsymbol{k}',t)\rangle + \langle f_i(\boldsymbol{k},t')f_j(\boldsymbol{k}',t) \rangle\right] \nonumber\\ &= \frac{1}{2}(2{\rm \pi})^3 \delta (\boldsymbol{k} + \boldsymbol{k}') \varPhi(k)\mathcal{P}_{ij}(\boldsymbol{k}), \end{align}

where, in the second equality, we have taken the contributions of the first two correlators inside the $t'$ integral to vanish, as demanded by causality, and used (E9) to express the final correlator in terms of $\varPhi (k)$.

A similar strategy may be employed to treat the triple correlations appearing in (E12). Multiplying by $u_l(t,\boldsymbol {k}_3)$ the formal solution of the unaveraged version of (E12), which may be written as

(E14)\begin{align} w_i(\boldsymbol{k}_1) w_j(\boldsymbol{k}_2) &= \int^t \mathrm{d} t' \left\{- \nu ( k_1^2+k_2^2 ) w_i(\boldsymbol{k}_1, t') w_j(\boldsymbol{k}_2, t') \vphantom{\int^t}\right.\nonumber\\ &\quad - i \int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} [ k_{1n} \mathcal{P}_{iq}(\boldsymbol{k}_1) u_n(\boldsymbol{k}'', t') w_q(\boldsymbol{k}_1 - \boldsymbol{k}'', t') w_j(\boldsymbol{k}_2, t') \nonumber\\ &\quad + k_{2n} \mathcal{P}_{jq}(\boldsymbol{k}_2) u_n(\boldsymbol{k}'', t') w_q(\boldsymbol{k}_2 - \boldsymbol{k}'', t') w_i(\boldsymbol{k}_1, t') ]\nonumber\\ &\quad +\left.\vphantom{\int^t}f_i(\boldsymbol{k}_1,t')w_j(\boldsymbol{k}_2,t') +f_j(\boldsymbol{k}_2,t')w_i(\boldsymbol{k}_1,t')\right\}, \end{align}

ensemble averaging, and, finally, using the short-correlation-time assumption to split correlators by invoking causality, yields

(E15)\begin{align} &\langle w_i(t, \boldsymbol{k}_1) w_j(t, \boldsymbol{k}_2) u_l(t, \boldsymbol{k}_3) \rangle ={-} \frac{i}{2} (2 {\rm \pi})^3 \kappa_{ln}(\boldsymbol{k}_3) \delta (\boldsymbol{k}_1 + \boldsymbol{k}_2 + \boldsymbol{k}_3)\nonumber\\ &\quad \times[k_{1n}\mathcal{P}_{iq}(\boldsymbol{k}_1) \varPsi_{qj} ( \boldsymbol{k}_1 + \boldsymbol{k}_3 ) + k_{2n}\mathcal{P}_{jq}(\boldsymbol{k}_2) \varPsi_{iq} ( \boldsymbol{k}_2 + \boldsymbol{k}_3 )]. \end{align}

Using (E13) and (E15), (E12) becomes

(E16)\begin{align} &\partial_t \langle w_i(\boldsymbol{k}) w_j(\boldsymbol{k}') \rangle + \nu ( k^2+k'^2 ) \langle w_i(\boldsymbol{k}) w_j(\boldsymbol{k}') \rangle ={-}\frac{1}{2} \int \mathrm{d}^3 \boldsymbol{k}'' \kappa_{ln}(\boldsymbol{k}'') \delta (\boldsymbol{k} + \boldsymbol{k}') \nonumber\\ &\quad \times \{ k_l \mathcal{P}_{ip}(\boldsymbol{k}) [ (k_n-k''_n)\mathcal{P}_{pq}(\boldsymbol{k} - \boldsymbol{k}'') \varPsi_{qj} ( \boldsymbol{k} ) + k'_{n}\mathcal{P}_{jq}(\boldsymbol{k}') \varPsi_{pq} ( \boldsymbol{k}' + \boldsymbol{k}'' ) ] \nonumber\\ &\quad + k'_l \mathcal{P}_{jp}(\boldsymbol{k}') [ k_n \mathcal{P}_{iq}(\boldsymbol{k}) \varPsi_{qp} ( \boldsymbol{k} + \boldsymbol{k}'' ) + (k'_{n}-k''_{n})\mathcal{P}_{pq}(\boldsymbol{k}'-\boldsymbol{k}'') \varPsi_{iq} ( \boldsymbol{k}') ]\}\nonumber\\ & \quad + (2{\rm \pi})^3 \delta (\boldsymbol{k} + \boldsymbol{k}') \varPhi(k)\mathcal{P}_{ij}(\boldsymbol{k}). \end{align}

Integrating (E16) over $\boldsymbol {k}'$, and taking the trace (i.e. setting $i=j$ and summing over $i$), we obtain, after a small amount of algebra,

(E17)\begin{align} &\partial_t \varPsi(k) + 2 \nu k^2 \varPsi(k)\nonumber\\ & ={-}\frac{1}{2}\int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} \kappa_{ln}(\boldsymbol{k}'') k_l k_n \mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k} + \boldsymbol{k}'') [\varPsi(k) - \varPsi(|\boldsymbol{k} +\boldsymbol{k}''|)] + 2 \varPhi(k), \end{align}

where a useful step in the simplification of the integrand was noting that $\kappa _{ln }(\boldsymbol {k}'')k''_n \propto \mathcal {P}_{ln} (\boldsymbol {k}'')k''_n=0$.

Now, $\varPsi (k)$ may be brought outside the first integral in (E17), giving rise to a turbulent-viscosity term, $-\nu _T(k)k^2 \varPsi (k)$, where

(E18)\begin{align} 2 \nu_T(k) & \equiv \frac{1}{2}\int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} \kappa_{ln}(\boldsymbol{k}'') \frac{k_l k_n}{k^2} \mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k} + \boldsymbol{k}'') \nonumber\\ & = \frac{1}{(2{\rm \pi})^3} \int_0^{\infty} \mathrm{d} k'' k''^2 \kappa (k'') \int \mathrm{d}^2 \varOmega'' \left[1-\frac{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{k}'')^2}{k^2 k''^2}\right] \left[1 - \frac{(\boldsymbol{k} \times \boldsymbol{k}'')^2}{2 k^2 (\boldsymbol{k} + \boldsymbol{k}'')^2}\right]\nonumber\\ & \equiv \frac{1}{(2{\rm \pi})^2} \int_0^\infty \mathrm{d} k'' k''^2 \kappa (k'') G\left(\frac{k''}{k}\right). \end{align}

The function $G(x)$ may be obtained by choosing a spherical coordinate system about $\boldsymbol {k}$. Computing the angle integral then yields

(E19)\begin{equation} G(x) = \frac{1}{96 x^3} \left[ 2 x (3+ 53x^2 - 11x^4 +3x^6) - 3 (1-x^2)^4 \log\left|\frac{1+x}{1-x}\right| \right].\end{equation}

Finally, substituting $\kappa (k'')=\kappa _0 \delta (k''-k_f)$, we get

(E20)\begin{equation} 2 \nu_T(k) = \frac{1}{(2{\rm \pi})^2} \kappa_0 k_f^2 G\left(\frac{k_f}{k}\right). \end{equation}

The other term inside the integral in (E17), containing $\varPsi (|\boldsymbol {k}+\boldsymbol {k}''|)$, is handled in the following manner:

(E21)\begin{align} &\frac{1}{2}\int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} \kappa_{ln}(\boldsymbol{k}'') k_l k_n \mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k} + \boldsymbol{k}'') \varPsi ( | \boldsymbol{k} + \boldsymbol{k}'' |) \nonumber\\ &\quad = \frac{1}{2}\frac{1}{(2{\rm \pi})^3}\int_0^{\infty} \mathrm{d} k'' k''^2 \kappa (k'') \int \mathrm{d}^2 \varOmega'' k_l k_n\mathcal{P}_{ln}(\boldsymbol{k}'')\mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k} + \boldsymbol{k}'') \varPsi ( | \boldsymbol{k} + \boldsymbol{k}'' |)\nonumber\\ &\quad = \frac{1}{2}\frac{1}{(2{\rm \pi})^3} \int_0^{\infty} \mathrm{d} k'' k''\kappa (k'')\int^{2{\rm \pi}}_0 \mathrm{d} \varphi'' \int^{k+k''}_{|k-k''|} \mathrm{d} k' \frac{k'}{k} k_l k_n \mathcal{P}_{ln}(\boldsymbol{k}'') \mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k}')\varPsi ( k')\nonumber\\ &\quad = \frac{1}{(2{\rm \pi})^2} \int_0^{\infty} \mathrm{d} k'' k''\kappa (k'') \int^{k+k''}_{|k-k''|} \mathrm{d} k' \frac{k'}{k} K(k'', k', k) \varPsi ( k'), \end{align}

where, in the second equality, we set $\boldsymbol {k}' = \boldsymbol {k} + \boldsymbol {k}''$, chose the polar axis to be along $\boldsymbol {k}$, and changed variables from the polar angle $\theta ''$ to $k'=|\boldsymbol {k}'|$, noting that $k' \mathrm {d} k' = -k k'' \sin \theta ''\mathrm {d}\theta ''$. In the final equality, we have used $\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {k}''=(k'^2-k^2-k''^2)/2$ and $\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {k}' = (k'^2+k^2 - k''^2)/2$ to write the integration kernel $K(k'', k',k)$ in terms of the magnitudes of $k$, $k'$ and $k''$ only. Explicitly, it is

(E22)\begin{align} K(k'', k',k)&={-}\frac{(k^2-k''^2)^4}{32 k^2 k''^2}\frac{1}{k'^2}-\frac{(k^2-k''^2)^3}{8k^2 k''^2}+\frac{(5k^4+6k^2k''^2 - 3k''^4)}{16k^2 k''^2}k'^2\nonumber\\ &\quad - \frac{k^2 - k''^2}{8k^2k''^2}k'^4 - \frac{1}{32 k^2 k''^2}k'^6. \end{align}

Finally, taking $\kappa (k'')=\kappa _0 (k''-k_f)$ and integrating over $k''$, we find

(E23)\begin{align} &\frac{1}{2}\int \frac{\mathrm{d}^3 \boldsymbol{k}''}{(2{\rm \pi})^3} \kappa_{ln}(\boldsymbol{k}'') k_l k_n \mathcal{P}_{pq} (\boldsymbol{k}) \mathcal{P}_{pq} (\boldsymbol{k} + \boldsymbol{k}'') \varPsi ( | \boldsymbol{k} + \boldsymbol{k}'' |) \nonumber\\ &\quad =\frac{\kappa_0}{(2{\rm \pi})^2} k_f \int^{k+k_f}_{|k-k_f|} \mathrm{d} k' \frac{k'}{k} K(k_f, k', k) \varPsi ( k'). \end{align}

Using (E18) and (E23), (E17) becomes

(E24)\begin{equation} \partial_t \varPsi(k) + 2 [\nu+\nu_T(k)] k^2 \varPsi(k) = \frac{\kappa_0}{(2{\rm \pi})^2} k_f \int^{k+k_f}_{|k-k_f|} \mathrm{d} k' \frac{k'}{k} K(k', k) \varPsi ( k') + 2 \varPhi(k), \end{equation}

where we have suppressed the explicit dependence of $K(k_f, k', k)$ on $k_f$. Using ${\mathcal {E}(k)= k^2 \varPsi (k) / 2{\rm \pi} ^2}$, and identifying ${F=k^2\varPhi (k)/{\rm \pi} ^2}$ as the spectrum of the energy injection, the evolution equation for the spectrum of the passive vector becomes

(E25)\begin{equation} \partial_t \mathcal{E}(k) + 2 [\nu + \nu_T(k)] k^2 \mathcal{E}(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} ( k') + F(k),\end{equation}

which is (4.12), as promised. This is the analogue for the passive vector field $\boldsymbol {w}$ of the similar mode-coupling equations for the magnetic field (Kulsrud & Anderson Reference Kulsrud and Anderson1992) and the passive scalar (Schekochihin, Haynes & Cowley Reference Schekochihin, Haynes and Cowley2004).

E.3. Small-$k$ limit of the mode-coupling equation

In the present study, we require only the limit of (E25) with $k\ll k_f$. Noting that $\lim _{x\to \infty }G(x)=4/5$, and

(E26)\begin{equation} \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}

we find that (E25) becomes

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

where $\beta = \kappa _0 / 5{\rm \pi} ^2$ and we have assumed that $\nu _T\gg \nu$. This is (4.15).

Appendix F. Alternative derivation of the scaling of $\langle {P}_V^2\rangle$ vs $R$ in turbulence forced with long-range correlations

In this appendix, we show how the scalings (4.21), viz.,

(F1) \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}

which were derived in § 4.2 from the spectral evolution equation (4.16) for the passive vector field, can be obtained by instead considering momentum diffusion as the sum of many instances of a decaying passive vector field.

Between the times $t_{\textrm{inj}}$ and $t_{\textrm{inj}}+\mathrm {d} t_{\textrm{inj}}$, the forcing causes the spectral energy density to increase by ${\mathrm {d}\mathcal {E}_{\boldsymbol {w}} = Ck^b \mathrm {d} t_{\textrm{inj}}}$. As (4.16) is linear, we may consider the evolution of the spectral-energy-density increment $\mathrm {d}\mathcal {E}_{\boldsymbol {w}}$ in isolation from the rest of $\mathcal {E}_{\boldsymbol {w}}$. Assuming that $k\ll k_f$ and $b<4$, (4.16) implies that $\mathrm {d}\mathcal {E}_{\boldsymbol {w}}$ decays with time $t$ according to

(F2)\begin{equation} \mathrm{d}\mathcal{E}_{\boldsymbol{w}}(k, t) = C k^b \,\mathrm{d} t_{\textrm{inj}} \exp(-\beta k_f^2 k^2 {\rm \Delta} t),\end{equation}

where ${\rm \Delta} t = t-t_{\textrm{inj}}>0$. The energy-containing scale and total energy associated with $\mathrm {d} \mathcal {E}_{\boldsymbol {w}}$ are respectively ${\lambda \sim (\beta k_f^2 {\rm \Delta} t)^{1/2}}$ (note that $\lambda$ is not the same as $l$, the latter being the energy-containing scale associated with the full field $\boldsymbol {w}$, rather than solely the increment generated at $t_{\textrm{inj}}$) and

(F3)\begin{equation} \mathrm{d} E = C \,\mathrm{d} t_{\textrm{inj}} \int_0^{\infty}\mathrm{d} k k^b \exp(-\beta k_f^2 k^2 {\rm \Delta} t ) = \frac{1}{2} \varGamma \left( \frac{1+b}{2}\right)C \mathrm{d} t_{\textrm{inj}} (\beta k_f^2 {\rm \Delta} t)^{-(1+b)/2}. \end{equation}

Let us now determine the contribution of $\mathrm {d}\mathcal {E}_{\boldsymbol {w}}(k, t)$ to the mean square momentum $\mathrm {d} \langle \boldsymbol {P}^2_V \rangle$ contained within a spherical control volume $V$ with radius $R$. There are three different cases to consider:

  1. (i) If $R\ll \lambda$, then (3.10) gives $\mathrm {d} \langle \boldsymbol {P}^2_V \rangle \propto R^{6}\,\mathrm {d} E \propto R^{6} ({\rm \Delta} t)^{-(1+b)/2} \mathrm {d} t_{\textrm{inj}}$, as the decaying field has reached scales much larger than $R$, so there is little variation of $\boldsymbol {w}$ within $V$.

  2. (ii) If $R\gg \lambda$ and $b<3$, then (3.10) gives $\mathrm {d} \langle \boldsymbol {P}^2_V \rangle \propto R^{5-b} \mathrm {d} t_{\textrm{inj}}$, with no time dependence, owing to the fact that the momentum scaling is set by the power-law part of (F2) (see (3.9) and surrounding discussion), which is constant in time, a consequence of momentum conservation.

  3. (iii) If $R\gg \lambda$ and $3< b<4$, then (3.10) gives $\mathrm {d} \langle \boldsymbol {P}^2_V \rangle \propto R^2 c({\rm \Delta} t) \mathrm {d} t_{\textrm{inj}}$. An as-yet undetermined function $c({\rm \Delta} t)$ appears here because the integral (3.9) is dominated by the evolving contribution from the energy-containing scales, rather than from the invariant $k^b$ tail of (F2). In principle, we can determine $c({\rm \Delta} t)$ by substituting (F2) into our equation for $\chi (r)$ in terms of $\mathcal {E}(k)$, (3.7), and then evaluating the integral (3.4) exactly. However, some effort can be spared by using the self-similar form of (F2). Substituting $\mathcal {E}(k) = k^b g(k^2 {\rm \Delta} t)$ in (3.7), and changing the integration variable to $x=kr$, we have

    (F4)\begin{equation} u^2 \chi(r) = \frac{2}{r^{1+b}}\int_0^{\infty} \mathrm{d}\kern 0.06em x x^b g \left(\frac{x^2 t}{r^2}\right)\frac{\sin x - x \cos x}{x} \equiv \frac{2}{r^{1+b}} G\left(\frac{{\rm \Delta} t}{r^2}\right). \end{equation}
    Then, with a change of variables to $y = r'/\sqrt {{\rm \Delta} t}$, (3.4) becomes
    (F5)\begin{equation} \langle \boldsymbol{P}^2_V \rangle = 4{\rm \pi}^2 u^2 ({\rm \Delta} t)^{(3-b)/2}\int^{2R}_0 \mathrm{d} r r \int^{r/\sqrt{{\rm \Delta} t}}_0 \mathrm{d} y y^{2-b} G\left(\frac{1}{y}\right).\end{equation}
    For $b>3$, (3.10) demands that $\langle \boldsymbol {P}_{V}^2\rangle \propto R^2$ at any fixed time, so the integral over $y$ in (F5) must be independent of its upper limit. In this case, the result of the double integration is time independent, so we deduce $c({\rm \Delta} t) = ({\rm \Delta} t)^{(3-b)/2}$.

To summarise these results, we have established that

(F6) \begin{equation} \mathrm{d} \langle \boldsymbol{P}^2_V \rangle \propto \mathrm{d} t_{\textrm{inj}} \begin{cases} R^{6} ({\rm \Delta} t)^{-(1+b)/2} & \text{if }R\ll \lambda\sim (\beta k_f^2 {\rm \Delta} t)^{1/2}, \\ R^{5-b} & \text{if }R \gg \lambda \sim (\beta k_f^2 {\rm \Delta} t)^{1/2}\ \&\ b<3, \\ R^{2} ({\rm \Delta} t)^{(3-b)/2} & \text{if }R \gg \lambda\sim (\beta k_f^2 {\rm \Delta} t)^{1/2}\ \&\ 3< b<4. \\ \end{cases}\end{equation}

Under the passive-field assumption, and owing to the fact that the forcing at any given time is uncorrelated with the forcing at any other time, the total squared momentum contained within a volume of continually forced turbulence may be obtained as the sum of contributions from passive decays initialised continuously and uniformly in time. Let us consider a fixed sphere of radius $R$ and let $t_c = R^2/\beta k_f^2$ be the time at which the energy-containing scale of the decaying field initialised at $t_{\textrm{inj}}=0$ reaches the scale $R$. Then, for $t< t_c$, each part of the sum comes from the $R \gg \lambda$ part of (F6), so

(F7) \begin{equation} \langle \boldsymbol{P}^2_V \rangle = \int \mathrm{d} \langle \boldsymbol{P}^2_V \rangle = \int^{t}_0 \mathrm{d} t_{\textrm{inj}} \frac{\mathrm{d} \langle \boldsymbol{P}^2_V \rangle}{\mathrm{d} t_{\textrm{inj}}}\propto \begin{cases} R^{5-b} t & \text{if }b<3, \\ R^{2} t^{(5-b)/2} & \text{if }3< b<4. \\ \end{cases} \end{equation}

If, instead, $t>t_c$, then

(F8) \begin{equation} \langle \boldsymbol{P}^2_V \rangle = \int \mathrm{d} \langle \boldsymbol{P}^2_V \rangle = \int^{t}_{t-t_c} \mathrm{d} t_{\textrm{inj}} \frac{\mathrm{d} \langle \boldsymbol{P}^2_V \rangle}{\mathrm{d} t_{\textrm{inj}}} + \int^{t-t_c}_{0} \mathrm{d} t_{\textrm{inj}} \frac{\mathrm{d} \langle \boldsymbol{P}^2_V \rangle}{\mathrm{d} t_{\textrm{inj}}}.\end{equation}

The first integral encodes all the decays that have not yet reached the scale $R$, as ${t-t_{\textrm{inj}}< t_c}$ for them. Therefore, we may again substitute for $\mathrm {d} \langle \boldsymbol {P}^2_V \rangle$ using the $R \gg \lambda$ part of (F6), giving

(F9) \begin{equation} \int^{t}_{t-t_c} \mathrm{d} t_{\textrm{inj}} \frac{\mathrm{d} \langle \boldsymbol{P}^2_V \rangle}{\mathrm{d} t_{\textrm{inj}}} \propto \begin{cases} R^{5-b} t_c & \text{if }b<3, \\ R^{2} t_c^{(5-b)/2} & \text{if }3< b<4, \\ \end{cases} \propto R^{7-b}. \end{equation}

The second integral in (F8) encodes all the decays that were initialised at $t_{\textrm{inj}}< t-t_c$, so have reached a scale larger than $R$ at time $t$. It is

(F10) \begin{equation} \int^{t}_{t_c} \mathrm{d} t_{\textrm{inj}} \frac{\mathrm{d} \langle \boldsymbol{P}^2_V \rangle}{\mathrm{d} t_{\textrm{inj}}} \propto R^6 \int^{t-t_c}_{0} \mathrm{d} t_{\textrm{inj}} (t-t_{\textrm{inj}})^{-(1+b)/2}\propto R^6 t_c^{(1-b)/2} \propto R^{7-b}. \end{equation}

Thus, decays that have reached $\lambda >R$ and ones with $\lambda < R$ both contribute a term $\propto R^{7-b}$ to $\langle \boldsymbol {P}_V^2\rangle$, explaining the scaling found in (4.21). In particular, this calculation explains why forcing with different values of $3< b<4$ results in saturation with different scalings for $\langle \boldsymbol {P}_V^2\rangle$, despite both having $\langle \boldsymbol {P}_V^2\rangle \propto R^2$ at $R>R_c$: the total energy contained within an initially more diffuse blob of momentum decays more slowly than an initially less diffuse one (i.e. one with smaller $b$), even though the energy-containing scales of both grow at the same rate.

As an aside, we note that the saturated spectrum is also derivable directly from the ensemble of decaying states, (F2):

(F11) \begin{equation} \mathcal{E}_{\boldsymbol{w}}(k, t) = \int \mathrm{d} \mathcal{E}_{\boldsymbol{w}}(k,t) = \int^t_0 \mathrm{d} t_{\textrm{inj}} C k^b \exp(-\beta k_f^2 k^2 (t-t_{\textrm{inj}}) ) = [1- {\rm e}^{-\beta k^2 k_f^2 t }]\frac{C k^{b-2}}{\beta k_f^2}, \end{equation}

which is (4.17) (with the $\beta \mathcal {E}_{\boldsymbol {w}}(k_f)k^4$ term neglected).

Appendix G. Non-solenoidal forcing

In this appendix, we formalise the discussion of non-solenoidal forcing in § 4.3. The essential result is due to Saffman (Reference Saffman1967), which may be stated in the notation of present study as follows:

Theorem (Saffman)

Let $\boldsymbol {f}(\boldsymbol {x})$ be a random function of $\boldsymbol {x}$ that is statistically isotropic and homogeneous, and that has an analytic spectral tensor

(G1)\begin{equation} M_{ij}(\boldsymbol{k}, t, t') \equiv \int \mathrm{d}^3 \boldsymbol{r} \langle f_{\alpha}(t, \boldsymbol{x}) f_{\beta}(t', \boldsymbol{x} + \boldsymbol{r})\rangle \, {\rm e}^{-{\rm i}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{r}}. \end{equation}

Let $\boldsymbol {f}^{{(s)}}$ be the solenoidal part of $\boldsymbol {f}$, with spectral tensor $M^{{(s)}}_{ij}(\boldsymbol {k}, t, t')$. Then, the ‘Saffman integrals’ associated with $\boldsymbol {f}$ and $\boldsymbol {f}^{{(s)}}$, given by ${L_{\boldsymbol {f}}= \int ^t \mathrm {d} t' M_{ii}(\boldsymbol {0}, t, t')}$ and ${L_{\boldsymbol {f}^{{(s)}}}= \int ^t \mathrm {d} t' M^{{(s)}}_{ii}(\boldsymbol {0}, t, t')}$, respectively (cf. (4.23)), satisfy

(G2)\begin{equation} L_{\boldsymbol{f}^{{(s)}}}=\frac{2}{3}L_{\boldsymbol{f}}. \end{equation}

G.1. Proof of Saffman's theorem

A proof of Saffman's theorem, adapted from Saffman (Reference Saffman1967), is as follows (we are grateful to an anonymous referee for suggesting the following proof, which is more direct than the proof that we had originally). Because $\boldsymbol {f}$ is statistically isotropic and homogeneous, we may write

(G3)\begin{align} M_{ij}(\boldsymbol{k},t,t') & = a(k,t,t')P_{ij}(\boldsymbol{k}) + b(k,t,t')\frac{k_i k_j}{k^2} \nonumber\\ & =a(k,t,t')\delta_{ij} +[b(k,t,t') - a(k,t,t')]\frac{k_i k_j}{k^2}, \end{align}

from which it follows that

(G4)\begin{equation} M^{(s)}_{ij}(\boldsymbol{k},t,t') = a(k,t,t')P_{ij}(\boldsymbol{k}).\end{equation}

Since $M_{ij}(\boldsymbol {k},t,t')$ is analytic in $\boldsymbol {k}$ by assumption, it must be the case that ${a(0,t,t')=b(0,t,t')}$. It follows immediately that

(G5)\begin{equation} M^{(s)}_{ii}(\boldsymbol{0},t,t') = 2 a(0,t,t') = \tfrac{2}{3}M_{ii}(\boldsymbol{0},t,t').\end{equation}

Integrating (G5) over $t'$, we recover (G2), QED.

G.2. Consequences of Saffman's theorem

We make the following remarks regarding (G2):

  1. (i) The utility of Saffman's theorem for forced turbulence is as follows. Equations (2.15) and (2.16) together imply that

    (G6)\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) +L_{\boldsymbol{f}^{{(s)}}}, \end{equation}
    for a turbulence forced by a series of realisations of $\boldsymbol {f}$ that are delta-correlated in time. According to the argument presented in § 2.2, the term involving $K(r)$ vanishes (note that, if $L_{\boldsymbol {f}^{{(s)}}}\neq 0$, $K(r)=O(r^{-3})$ as $r\to \infty$, not $O(r^{-4})$ (as in (2.14)). This is because long-range correlations in $\boldsymbol {f}$ can propagate into $K$ via the correlators of $\boldsymbol {u}$ and $\boldsymbol {f}$ that are present in (2.13). Nonetheless, such a decay of $K(r\to \infty )$ is still sufficiently rapid for the term containing $K(r)$ in (G6) to be negligible). Thus, according to (G2),
    (G7)\begin{equation} \frac{\mathrm{d} L}{\mathrm{d} t} = \frac{2}{3}L_{\boldsymbol{f}} ,\end{equation}
    so the rate of growth of $L$ in forced turbulence is finite, and proportional to $L_{\boldsymbol {f}}$.
  2. (ii) Saffman's theorem requires that $M_{ij}$ be analytic. This condition is guaranteed by choosing $\boldsymbol {f}$ so that $\langle f_i(\boldsymbol {x}) f_j(\boldsymbol {x} +\boldsymbol {r})\rangle$ decays rapidly with distance – for example, we might take $\boldsymbol {f}$ to consist of an ensemble of local patches of uniformly directed force, whose magnitude decays exponentially away from their centre. However, it should be emphasised that Saffman's result does fail if $\langle f_i(\boldsymbol {x}) f_j(\boldsymbol {x} +\boldsymbol {r})\rangle$ decays slowly with $r$. An immediate example of this is the case where $\boldsymbol {f}$ is itself solenoidal, so $\boldsymbol {f} = \boldsymbol {f}^{{(s)}}$, and therefore naive application of (G2) would suggest $L_{\boldsymbol {f}}=0$. However, this does not mean that $L_{\boldsymbol {f}}=0$ in general for solenoidal forcing, as the argument applies only to functions $\boldsymbol {f}(\boldsymbol {x})$ with an analytic spectral tensor. Long-range correlations in $\boldsymbol {f}^{(s)}(\boldsymbol {x})$ are necessary for $L_{\boldsymbol {f}^{(s)}}\neq 0$, in which case its spectral tensor is not analytic (Saffman Reference Saffman1967).

  3. (iii) Saffman's theorem does not imply that any non-solenoidal forcing with an analytic spectral tensor will have $L_{\boldsymbol {f}^{{(s)}}}\neq 0$, and so induce a flow with $\mathcal {E}(k\to 0)\propto k^2$. This is because it remains possible that $L_{\boldsymbol {f}} = 0$ (i.e. $a(0,t,t')=0$ in (G5)). An example of this would be a ensemble of pairs of oppositely directed, local instantaneous impulses, separated by a small distance (cf. § 3.2). Another example, pertinent to numerical studies of forced turbulence, is forcing in a finite spectral band. As shown in Appendix B, long-range correlations decay arbitrarily rapidly for such a forcing, meaning that $M_{ij}$ is analytic – nonetheless, $L_{\boldsymbol {f}}$ is manifestly zero for such a $\boldsymbol {f}$, and, therefore, so is $L_{\boldsymbol {f}^{{(s)}}}$.

References

Adzhemyan, L.T., Antonov, N.V., Mazzino, A., Muratore-Ginanneschi, P. & Runov, A.V. 2001 a Pressure and intermittency in passive vector turbulence. Europhys. Lett. 55, 801.CrossRefGoogle Scholar
Adzhemyan, L.T., Antonov, N.V. & Runov, A.V. 2001 b Anomalous scaling, nonlocality, and anisotropy in a model of the passively advected vector field. Phys. Rev. E 64, 046310.CrossRefGoogle Scholar
Alexakis, A. & Biferale, L. 2018 Cascades and transitions in turbulent flows. Phys. Rep. 767–769, 1101.CrossRefGoogle Scholar
Alexakis, A. & Brachet, M.-E. 2019 On the thermal equilibrium state of large-scale flows. J. Fluid Mech. 872, 594625.CrossRefGoogle Scholar
Antonov, N.V., Hnatich, M., Honkonen, J. & Jurčišin, M. 2003 Turbulence with pressure: anomalous scaling of a passive vector field. Phys. Rev. E 68, 046306.CrossRefGoogle ScholarPubMed
Arponen, H. 2009 Anomalous scaling and anisotropy in models of passively advected vector fields. Phys. Rev. E 79, 056303.CrossRefGoogle ScholarPubMed
Batchelor, G.K. & Proudman, I. 1956 The large-scale structure of homogeneous turbulence. Phil. Trans. R. Soc. A 248, 369405.Google Scholar
Benzi, R., Biferale, L. & Toschi, F. 2001 Universality in passively advected hydrodynamic fields: the case of a passive vector with pressure. Eur. Phys. J. B 24, 125133.CrossRefGoogle Scholar
Birkhoff, G. 1954 Fourier synthesis of homogeneous turbulence. Commun. Pure Appl. Maths 7, 1944.CrossRefGoogle Scholar
Biskamp, D. & Müller, W.-C. 1999 Decay laws for three-dimensional magnetohydrodynamic turbulence. Phys. Rev. Lett. 83, 2195.CrossRefGoogle Scholar
Brandenburg, A. & Kahniashvili, T. 2017 Classes of hydrodynamic and magnetohydrodynamic turbulent decay. Phys. Rev. Lett. 118, 055102.CrossRefGoogle ScholarPubMed
Brandenburg, A., Kahniashvili, T. & Tevzadze, A.G. 2015 Nonhelical inverse transfer of a decaying turbulent magnetic field. Phys. Rev. Lett. 114, 075001.CrossRefGoogle ScholarPubMed
Brandenburg, A., Zhou, H. & Sharma, R. 2023 Batchelor, Saffman, and Kazantsev spectra in galactic small-scale dynamos. Mon. Not. R. Astron. Soc. 518, 3312.CrossRefGoogle Scholar
Cameron, A., Alexakis, A. & Brachet, M.-E. 2017 Effect of helicity on the correlation time of large scales in turbulent flows. Phys. Rev. Fluids 2, 114602.CrossRefGoogle Scholar
Dallas, V., Fauve, S. & Alexakis, A. 2015 Statistical equilibria of large scales in dissipative hydrodynamic turbulence. Phys. Rev. Lett. 115, 204501.CrossRefGoogle ScholarPubMed
Davidson, P.A. 2009 The role of angular momentum conservation in homogeneous turbulence. J. Fluid Mech. 632, 329358.CrossRefGoogle Scholar
Davidson, P.A. 2011 The minimum energy decay rate in quasi-isotropic grid turbulence. Phys. Fluids 23, 085108.CrossRefGoogle Scholar
Davidson, P.A. 2013 Turbulence in Rotating, Stratified and Electrically Conducting Fluids. Cambridge University Press.CrossRefGoogle Scholar
Davidson, P.A. 2015 Turbulence: An Introduction for Scientists and Engineers. Oxford University Press.CrossRefGoogle Scholar
Durrer, R. & Caprini, C. 2003 Primordial magnetic fields and causality. J. Cosmol. Astropart. Phys. 2003, 010.CrossRefGoogle Scholar
Eyink, G.L. & Thomson, D.J. 2000 Free decay of turbulence and breakdown of self-similarity. Phys. Fluids 12 (3), 477479.CrossRefGoogle Scholar
Falkovich, G., Gaweȩdzki, K. & Vergassola, M. 2001 Particles and fields in fluid turbulence. Rev. Mod. Phys. 73, 913.CrossRefGoogle Scholar
Galishnikova, A.K., Kunz, M.W. & Schekochihin, A.A. 2022 Tearing instability and current-sheet disruption in the turbulent dynamo. Phys. Rev. X 12, 041027.Google Scholar
Hatori, T. 1984 Kolmogorov-style argument for the decaying homogeneous MHD turbulence. J. Phys. Soc. Japan 53, 25392545.CrossRefGoogle Scholar
Hosking, D.N. & Schekochihin, A.A. 2021 Reconnection-controlled decay of magnetohydrodynamic turbulence and the role of invariants. Phys. Rev. X 11, 041005.Google Scholar
Hosking, D.N. & Schekochihin, A.A. 2022 Cosmic-void observations reconciled with primordial magnetogenesis. arXiv:2203.03573.CrossRefGoogle Scholar
Ishida, T., Davidson, P.A. & Kaneda, Y. 2006 On the decay of isotropic turbulence. J. Fluid Mech. 564, 455475.CrossRefGoogle Scholar
von Kármán, T. & Howarth, L. 1938 On the statistical theory of isotropic turbulence. Proc. R. Soc. Lond. A 164, 192215.CrossRefGoogle Scholar
Kazantsev, A.P. 1968 Enhancement of a magnetic field by a conducting fluid. Sov. Phys. JETP 26, 10311034.Google Scholar
Kolmogorov, A.N. 1941 a Dissipation of energy in locally isotropic turbulence. Dokl. Acad. Nauk SSSR 32, 1517.Google Scholar
Kolmogorov, A.N. 1941 b Local structure of turbulence in incompressible viscous fluid at very large Reynolds numbers. Dokl. Acad. Nauk SSSR 30, 301305.Google Scholar
Kraichnan, R.H. 1965 Inertial-range spectrum of hydromagnetic turbulence. Phys. Fluids 8, 13851387.CrossRefGoogle Scholar
Kraichnan, R.H. 1968 Small-scale structure of a scalar field convected by turbulence. Phys. Fluids 11, 945953.CrossRefGoogle Scholar
Kraichnan, R.H. 1973 Helical turbulence and absolute equilibrium. J. Fluid Mech. 59, 745752.CrossRefGoogle Scholar
Kraichnan, R.H. 1987 a An interpretation of the Yakhot-Orszag turbulence theory. Phys. Fluids 30, 24002405.CrossRefGoogle Scholar
Kraichnan, R.H. 1987 b Kolmogorov's constant and local interactions. Phys. Fluids 30, 1583.CrossRefGoogle Scholar
Kraichnan, R.H. 1994 Anomalous scaling of a randomly advected passive scalar. Phys. Rev. Lett. 72, 10161019.CrossRefGoogle ScholarPubMed
Kulsrud, R.M. & Anderson, S.W. 1992 The spectrum of random magnetic fields in the mean-field dynamo theory of the galactic magnetic field. Astrophys. J. 396, 606630.CrossRefGoogle Scholar
Landau, L.D. & Lifshitz, E.M. 1959 Fluid Mechanics. Pergamon.Google Scholar
Lee, T.D. 1952 On some statistical properties of hydrodynamical and magneto-hydrodynamical fields. Q. Appl. Maths 10, 6974.CrossRefGoogle Scholar
Lesieur, M. 2008 Turbulence in Fluids. Springer.CrossRefGoogle Scholar
Lesieur, M., Métais, O. & Comte, P. 2005 Large-Eddy Simulations of Turbulence. Cambridge University Press.CrossRefGoogle Scholar
Lesur, G. 2015 Snoopy: general purpose spectral solver. Astrophysics Source Code Library (ascl:1505.022).Google Scholar
Loitsyansky, L.G. 1939 Some basic laws for isotropic turbulent flow. Tsentralni Aerogidrodinamicheskii Institut USSR 440, p. 3.Google Scholar
Maron, J. & Blackman, E.G. 2002 Effect of fractional kinetic helicity on turbulent magnetic dynamo spectra. Astrophys. J. Lett. 566, L41.CrossRefGoogle Scholar
Michel, G., Pétrélis, F. & Fauve, S. 2017 Observation of thermal equilibrium in capillary wave turbulence. Phys. Rev. Lett. 118, 144502.CrossRefGoogle ScholarPubMed
Orszag, S.A. 1977 Lectures on the statistical theory of turbulence. In Fluid Dynamics. Les Houches Summer School, 1973 (ed. R. Balian & J.-L. Peube), p. 235. Gordon and Breach.Google Scholar
Panickacheril John, J., Donzis, D.A. & Sreenivasan, K.R. 2022 Laws of turbulence decay from direct numerical simulations. Phil. Trans. R. Soc. A 380, 20210089.CrossRefGoogle ScholarPubMed
Reppin, J. & Banerjee, R. 2017 Nonhelical turbulence and the inverse transfer of energy: a parameter study. Phys. Rev. E 96, 053105.CrossRefGoogle ScholarPubMed
Rincon, F. 2019 Dynamo theories. J. Plasma Phys. 85, 205850401.CrossRefGoogle Scholar
Saffman, P.G. 1967 The large-scale structure of homogeneous turbulence. J. Fluid Mech. 27, 581593.CrossRefGoogle Scholar
Schekochihin, A.A. 2022 MHD turbulence: a biased review. J. Plasma Phys. 88, 155880501.CrossRefGoogle Scholar
Schekochihin, A.A., Haynes, P.H. & Cowley, S.C. 2004 Diffusion of passive scalar in a finite-scale random flow. Phys. Rev. E 70, 046304.CrossRefGoogle Scholar
Schekochihin, A.A., Parker, J.T., Highcock, E.G., Dellar, P.J., Dorland, W. & Hammett, G.W. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma turbulence. J. Plasma Phys. 82, 905820212.CrossRefGoogle Scholar
Figure 0

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$.

Figure 1

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.

Figure 2

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$.

Figure 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$.

Figure 4

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.

Figure 5

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.

Figure 6

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.