2.1. The statistical closure problem

The fundamental PDE of many-body plasma physics is the Klimontovich equation (Klimontovich Reference Klimontovich and ter Harr1967), written here for zero magnetic field $\boldsymbol{B}$ :

where

is the phase-space microdensity of species $s$ , $\doteq$ denotes a definition, $\overline{n}$ is the mean density, $\boldsymbol{\partial }\doteq (q/m)\,\partial /\partial \boldsymbol{v}$ , $q$ and $m$ are charge and mass, and $\widetilde{\boldsymbol{E}}$ is the microscopic electric field. In electrostatics, to which this article is restricted for simplicity, it is sufficient to write $\widetilde{\boldsymbol{E}}=-\boldsymbol{{\rm\nabla}}\widetilde{{\it\phi}}$ and to determine $\widetilde{\boldsymbol{E}}$ or $\widetilde{{\it\phi}}$ from Poisson’s equation,

The functions $\widetilde{\boldsymbol{x}}(t)$ and $\widetilde{\boldsymbol{v}}(t)$ are the random particle trajectories that evolve from certain initial conditions. The quantities $\widetilde{\boldsymbol{x}}$ and $\widetilde{\boldsymbol{v}}$ , and hence $\widetilde{N}$ and $\widetilde{\boldsymbol{E}}$ , are random because it is assumed that at some initial time the phase-space locations of the particles are sampled from some ensemble, i.e., $\widetilde{N}$ describes the phase-space dynamics of one particular realization of that ensemble. An observable quantity is the ensemble or statistical average, denoted by $\langle \cdots \rangle$ , of a random variable. For example, the so-called one-particle distribution function isFootnote ¹² $f_{s}(\boldsymbol{x},\boldsymbol{v},t)\doteq \langle \widetilde{N}_{s}(\boldsymbol{x},\boldsymbol{v},t)\rangle$ . One is also interested in higher-order statistics of the fluctuation ${\it\delta}N\doteq \widetilde{N}-\langle \widetilde{N}\rangle$ , such as $\langle {\it\delta}N(t){\it\delta}N(t^{\prime })\rangle$ . In classical, near-equilibrium plasmas, that two-point correlation function is used to derive the Balescu–Lenard collision operator (usually approximated in practice by the Landau operatorFootnote ¹³ ) that describes the effects of particle discreteness in the limit of weak coupling ( ${\it\epsilon}_{\text{p}}\doteq 1/(n{\it\lambda}_{\text{D}}^{3})\ll 1$ , where ${\it\lambda}_{\text{D}}$ is the Debye lengthFootnote ¹⁴ ). But even in the continuum limit ${\it\epsilon}_{\text{p}}\rightarrow 0$ , turbulent fluctuations of the phase-space fluid can occur (Davidson Reference Davidson1972). The goal is to determine their effects on observable spectra and turbulent transport.Footnote ¹⁵

For neutral fluids, the paradigmatic PDE is the Navier–Stokes equation constrained by incompressibility:

Here ${\it\rho}_{\text{m}}$ is the mass density, ${\it\mu}$ is the kinematic viscosity (here assumed to be constant), and $\widetilde{\boldsymbol{f}}\vphantom{f}^{\text{ext}}$ is a solenoidal external random forcing.Footnote ¹⁶ The pressure $\widetilde{P}$ in (2.4a ) is determined in terms of $\widetilde{\boldsymbol{u}}$ by taking the divergence of (2.4a ), enforcing (2.4b ), and solving the resulting Poisson equation. The stochastic velocity $\widetilde{\boldsymbol{u}}$ is random either because the initial velocities are distributed randomly or because of the random forcing. Because a single realization varies erratically in space and time, one is interested in the ensemble average $\boldsymbol{u}(\boldsymbol{x},t)\doteq \langle \widetilde{\boldsymbol{u}}(\boldsymbol{x},t)\rangle$ and in higher-order statistics such as the two-point correlation function $\langle {\it\delta}\boldsymbol{u}(\boldsymbol{x},t){\it\delta}\boldsymbol{u}(\boldsymbol{x}^{\prime },t^{\prime })\rangle$ , the Fourier transform of which with respect to difference variables is the energy spectrum in wavenumber and frequency. The mean flow $\boldsymbol{u}(\boldsymbol{x},t)$ is analogous to the plasma $f_{s}(\boldsymbol{x},\boldsymbol{v},t)$ in the sense that both are the ensemble averages of random variables. (Unlike $f$ , however, $\boldsymbol{u}$ is not a PDF.)

Operationally, ensemble averaging is well defined. One is given some probability density function(al) of initial conditions and/or random forcing, and averages are taken with that PDF. A question is whether the result of the averaging has something useful to say about the observable macroscopic world. In particular, whereas Newton’s laws of motion are time-reversible, ensemble averaging (often together with certain asymptotic limiting procedures) typically introduces irreversibility (one manifestation of which one sees in the laboratory as turbulent transport). The use of ensembles and the origins of irreversibility are subtle topics.Footnote ¹⁷ A readable article on the subject as applied to many-body theory is by Lebowitz (Reference Lebowitz1993), who says (the italics are his):

Having thus warned the reader, I shall continue with the use of ensemble averages, which are technically convenient. In practice they can often be equated to some sort of time average, providing a link between the processing of simulation data and the analytical theory.

Both of equations (2.1) and (2.4) are quadratically nonlinear, so they are special cases of the generic quadratic PDE

where 1 stands for the complete set of independent variables,Footnote ¹⁸ the coupling coefficients $U_{i}$ are specified (e.g., see (3.10a ) and (3.10b )), $U_{3}$ can be assumed to be symmetric in its last two arguments, an integration–summation convention over repeated arguments is assumed, and $\widetilde{f}^{\text{ext}}$ is external random forcing. It is conventional to choose $\widetilde{f}^{\text{ext}}$ to be Gaussian for simplicity, so it is completely specified by its mean $\langle \,\widetilde{f}^{\text{ext}}\rangle$ and covariance $F^{\text{ext}}(1,1^{\prime })\doteq \langle {\it\delta}f^{\text{ext}}(1){\it\delta}f^{\text{ext}}(1^{\prime })\rangle$ . There are two basic approaches to determining the statistics of $\widetilde{{\it\psi}}$ : either formulate approximate coupled equations for a few low-order cumulants, or ambitiously try to find an equation for the entire PDFFootnote ¹⁹ of $\widetilde{{\it\psi}}$ . The cumulant approach is easier and traditional, and it is the one I shall describe here. An article that summarizes aspects of both approaches is by Kraichnan (Reference Kraichnan and Sirovich1991).

The simplest statistical quantity is the mean field, which exactly obeys

This already exhibits a closure problem: the mean field is determined by the covariance $C(2,3)\doteq \langle {\it\delta}{\it\psi}(2){\it\delta}{\it\psi}(3)\rangle$ , which is unknown. To obtain an equation for $C$ , subtract (2.6) from (2.5) to find that the fluctuations obey

Upon multiplying (2.7) by ${\it\delta}{\it\psi}(1^{\prime })$ and averaging, one finds that $C(1,1^{\prime })$ obeysFootnote ²⁰

where the three-point correlation function is

Quantities such as $\langle {\it\psi}\rangle$ , $C$ , or $T$ are cumulants, a central concept.Footnote ²¹ An important property is that multipoint cumulants of statistically independent quantities vanish. We see that $n$ -point cumulants are driven by $(n+1)$ -point cumulants; the cumulant hierarchy is unclosed. This is the classical statistical closure problem (Kraichnan Reference Kraichnan1962; Leslie Reference Leslie1973; Orszag Reference Orszag, Balian and Peube1977; McComb Reference McComb1990, Reference McComb2014; Krommes Reference Krommes2002, and references therein). To close the system at the covariance level, one must express the triplet correlation function $T$ and the cross-correlation of ${\it\delta}f^{\text{ext}}$ and ${\it\delta}{\it\psi}$ in terms of one-point (the mean field) and two-point quantities (such as $C$ ; a two-point response function $R$ will also be important).

Statistical closure would not be a problem if the PDF of $\widetilde{{\it\psi}}$ were Gaussian. General PDFs can usually be represented in terms of an infinite set of cumulants;Footnote ²² Gaussians are special because they possess just two cumulants – namely, the mean and the variance. A theorem states that the sum of two Gaussian variables is again Gaussian, so if the dynamics were linear (i.e., $U_{3}=0$ ), $\widetilde{{\it\psi}}$ would be Gaussian if it evolved from Gaussian initial conditions under Gaussian forcing. Unfortunately, nonlinearity confounds this result. The product of two Gaussian variables is not Gaussian, so even if the initial conditions were Gaussian, non-Gaussian statistics would develop after just one time step.Footnote ²³ Once non-Gaussianity arises, it is compounded during each subsequent time interval. After long times a statistically steady PDF can develop, but it will be non-Gaussian in a self-consistent way that is difficult to determine analytically.Footnote ²⁴ The inevitable development of non-Gaussian statistics in nonlinear PDEs is the essence of the closure problem.

Over roughly the last century, enormous efforts have been expended toward obtaining satisfactory statistical closures. A famous example is the DIA; it will be discussed in § 4. Some simpler approximations have enjoyed good successes for the description of the interactions of mean fields with turbulence; see § 6.3.

2.2. Interlude: Langevin equations

Although statistical closures for realistic nonlinear PDEs such as the Klimontovich equation are difficult and nontrivial, a good deal of the physics embodied in such equations can be understood with simple intuition. A good example is the use of Langevin equations (Lemons & Gythiel Reference Lemons and Gythiel1997; Gardiner Reference Gardiner2004) to describe the statistical properties of classical Brownian motion. I shall include a few words of introduction here so that later I can interpret certain exact and closure results in terms of generalized Langevin descriptions.

2.3. An exactly solvable model of passive advection

There are actually two levels of difficulty in cumulant-based statistical closure. Problems arise even for passive advection. (In plasmas, that is frequently called stochastic acceleration (Orszag Reference Orszag and Fox1969), implying no back-reaction of the particles on the electric field.) Self-consistency (in plasmas, determining the field from the particles rather than prescribing the field externally) adds more complications (polarization effects) and leads to the concept of a dielectric function ${\mathcal{D}}$ .

2.4. Statistical closure and the stochastic oscillator

It is instructive to try to recover the essence of these exact results from statistical closure approximations because exact solutions are not available for realistic, nonlinear PDEs. This topic has been treated extensively in the literature. An important early paper (the first part of which is quite pedagogical) is by Kraichnan (Reference Kraichnan1961), who considered aspects of the infinite- ${\it\tau}_{\text{ac}}$ problem; further discussion and references are given by Krommes (Reference Krommes, Galeev and Sudan1984) and, more thoroughly, by Krommes (Reference Krommes2002). Here I shall just give a brief introduction.

Consider the restricted problem of finding an approximation to the mean response function $g({\it\tau})$ for (2.24). The random response function $\widetilde{g}(t;t^{\prime })$ obeys (2.27). Instead of writing a formally exact solution (which would not be useful if the primitive amplitude equation were really a PDE), let us write the evolution equation for $g$ by averaging (2.27). One gets

To be consistent with (2.47), I shall write this as

or

Either of these forms is called a Dyson equationFootnote ⁴⁶ after the famous unification by Dyson (Reference Dyson1949) of various approaches to quantum electrodynamics.

The goal is to find an approximate formula for ${\it\Sigma}$ . One has exactly

To proceed, one must assess the importance of the (stochastically) nonlinear term ${\it\delta}{\it\omega}\,{\it\delta}g$ , and that depends on the values of the dimensionless parameters. For present purposes, if either ${\mathcal{K}}$ or ${\mathcal{R}}$ is small, it is appropriate to drop the nonlinear term, in which case the solution of (2.53) is

where $g_{0}({\it\tau})\doteq H({\it\tau})\exp (-\text{i}\overline{{\it\omega}}{\it\tau}-{\it\nu}{\it\tau})$ is the lowest-order Green’s function. Thus

where

This has the form of (2.52) with a particular approximation for ${\it\Sigma}$ . Since in the present case ${\it\Sigma}$ has been constructed with the unperturbed propagator $g_{0}$ but describes the nonlinearity, one may call this the quasi-linear approximation.Footnote ⁴⁷

(As a consistency check, for ${\mathcal{K}}\rightarrow 0$ one has ${\it\Sigma}({\it\tau})\rightarrow g_{0}({\it\tau})2k^{2}D{\it\delta}({\it\tau})=k^{2}D{\it\delta}({\it\tau})$ , in agreement with (2.32).)

For the quasi-linear approximation to be valid, the nonlinear term in the equation for the fluctuations must be smaller than the linear term. For simplicity, assume that $\overline{{\it\omega}}=0$ . Then the ratio of nonlinear to linear rates is ${\it\beta}^{2}{\it\tau}_{\text{ac}}^{({\it\Sigma}_{0})}/{\it\nu}$ , which becomes large as ${\it\nu}\rightarrow 0~({\mathcal{R}}\rightarrow \infty )$ . In this situation, quasi-linear theory is unjustifiable and advanced techniques must be used to effect a proper closure. As discussed by Kraichnan (Reference Kraichnan1961), a plausible approximation is the DIA. For the oscillator model, that amounts to replacing the $g_{0}$ in (2.56) by the as-yet-unknown $g$ : ${\it\Sigma}({\it\tau})\approx g({\it\tau})W({\it\tau})$ . The exact Dyson equation (2.52) is then used to solve for $g$ self-consistently.Footnote ⁴⁸ For ${\mathcal{K}}\rightarrow \infty$ one has ${\it\eta}={\it\beta}^{2}{\it\tau}_{\text{ac}}$ , and the balance $\partial _{{\it\tau}}\sim {\it\eta}$ self-consistently determines ${\it\tau}_{\text{ac}}$ to be $O({\it\beta}^{-1})$ . This is consistent with the exact solution; see (2.36a ). Note that in this strong-turbulence limit the Markovian approximation breaks down because ${\it\Sigma}({\it\tau})$ and $g({\it\tau})$ vary on the same time scale; the physics is intrinsically nonlocal in time. I shall give a derivation of the DIA in § 4.1.

2.5. Renormalization

The replacement of $g_{0}$ by $g$ is an example of renormalization. A thorough discussion of the technical meaning of renormalization and of all of its ramifications is unfortunately beyond the scope of this article. Heuristically, it implies the broadening of resonance functions due to random nonlinear effects; in perturbation theory, that broadening is represented as a sum of certain terms through all orders of the nonlinear coupling. In general, the broadening must be determined self-consistently. An elucidating collection of articles with historical perspectives was edited by Brown (Reference Brown1993), a beginner’s guide to methodology is by McComb (Reference McComb2004), and an advanced book is by Zinn-Justin (Reference Zinn-Justin1996). Some of the connections to quantum field theory are reviewed by Krommes (Reference Krommes2002). Renormalization has a deep connection to so-called anomalous scaling, frequently discussed in the theory of phase transitions (Goldenfeld Reference Goldenfeld1992). A brief discussion of anomalous scaling that uses the stochastic oscillator as an example is given in § 6.1.2 of Krommes (Reference Krommes2002). The basic point is that the naive quasi-linear result for ${\it\Sigma}$ scales with ${\it\beta}^{2}$ , but a self-consistent calculation of the nonlinear damping rate for large ${\mathcal{K}}$ and ${\mathcal{R}}$ leads to a scaling with ${\it\beta}^{1}$ ; here the power 1 is called an anomalous exponent. Introductory pedagogical discussion of renormalization and anomalous scaling can also be found in the article by Krommes (Reference Krommes2009), which provides some useful background for the present tutorial.

3.1. Prelude: dielectric polarization in electromagnetism

Ultimately I shall consider turbulent systems that are not necessarily electromagnetic in nature, such as the Navier–Stokes fluid. However, it is useful to begin on familiar footing by discussing electromagnetic systems, fully ionized plasmas being a prominent example.

3.2. Self-consistency and the dielectric function of a turbulent plasma

I shall show in § 4.2 how dielectric response can be defined for an arbitrary quadratically nonlinear PDE (2.5). However, for initial discussion it is best to be more concrete, so I shall first discuss the Klimontovich equation (2.1) as a particular, very important example. Here $\widetilde{\boldsymbol{E}}$ is linearly related to $\widetilde{N}$ through Poisson’s equation, so one has $\widetilde{\boldsymbol{E}}=\pmb{\mathbb{E}}\widetilde{N}$ for a particular linear (vector) operator $\pmb{\mathbb{E}}$ . In general, it is necessary to introduce a dielectric tensor that describes the response to external currents. However, it is simplest for pedagogical purposes to consider the electrostatic limit $\boldsymbol{E}=-\boldsymbol{{\rm\nabla}}{\it\phi}$ , in which case one can work with a dielectric scalar that describes the response to external charge or electrostatic potential. In electrostatics, the Green’s function solution of Poisson’s equation leads toFootnote ⁵⁵

Thus the Klimontovich equation can be written as

This is a self-consistent ( $\widetilde{\boldsymbol{E}}$ is determined by $\widetilde{N}$ , not specified externally), multiplicatively stochastic PDE. It has the form of (2.5) with coupling coefficients

(The symmetrization in (3.10b ) is a manifestation of the self-consistency.) To enquire about linear response to an external or free charge, introduce an external (statistically sharp) electric field ${\rm\Delta}\boldsymbol{E}^{\text{ext}}~(\widetilde{\boldsymbol{E}}\rightarrow \pmb{\mathbb{E}}\widetilde{N}+{\rm\Delta}\boldsymbol{E}^{\text{ext}})$ and consider the equation for first-order response ${\rm\Delta}\widetilde{N}$ :

The underlining here and subsequently identifies the self-consistent response term that would be absent for passive advection. Note that, although ${\rm\Delta}\boldsymbol{E}^{\text{ext}}$ is not random, ${\rm\Delta}\widetilde{N}$ is random because it evolves in the random background of $\widetilde{N}$ . That is, there are random coefficients $\widetilde{\boldsymbol{E}}$ and $\widetilde{N}$ on the left-hand side of (3.11), and ${\rm\Delta}\boldsymbol{E}^{\text{ext}}$ couples to the random $\widetilde{N}$ on the right-hand side. Ultimately we must calculate the mean response $\langle {\rm\Delta}\widetilde{N}\rangle$ , which will involve us with a closure problem.

The left-hand side of (3.11) is just the linearization of the Klimontovich equation (2.1). It can be solved by introducing a random infinitesimal response function $\widetilde{R}$ , which is just Green’s functionFootnote ⁵⁶ for the left-hand side of (3.11). It obeys

Thus, upon omitting all variables except the time for brevity, one finds

This result is analogous to (2.28); the right-hand side of (3.11) can be recognized as a particular kind of forcing. However, there are two crucial differences. First, the forcing is not purely external because it involves the background quantity $\widetilde{N}$ ; consequently, it is correlated to $\widetilde{R}$ and the neat factorizations involved in (2.29) and (2.35a ) no longer hold. Second, $\widetilde{R}$ is distinct from $\widetilde{g}$ because of the presence of the $\pmb{\mathbb{E}}$ term on the left-hand side of (3.12). That describes the self-consistent response to a small perturbation; it is absent in a passive problem.

Given $\widetilde{R}$ , one can calculate the mean internally produced field by averaging (3.13) and evaluating $\pmb{\mathbb{E}}\langle {\rm\Delta}\widetilde{N}\rangle$ . One can then add ${\rm\Delta}\boldsymbol{E}^{\text{ext}}$ to obtain the total mean field in the medium. Finally, ${\mathcal{D}}$ is defined from the relation

one findsFootnote ⁵⁷

(This symbolic notation is a bit tricky because ${\mathcal{D}}$ is a function of only spatial coordinates whereas $\pmb{\mathbb{E}}$ and $\widetilde{R}$ operate in phase space.Footnote ⁵⁸ ) This is as far as one can go without facing up to statistical closure, i.e., an approximate calculation of $\langle \widetilde{R}\,\widetilde{N}\rangle$ .

One way of proceeding, ultimately not fruitful, is to first express $\widetilde{R}$ in terms of $\widetilde{g}$ , then attempt to average. The Klimontovich $\widetilde{R}$ obeys (3.12) with the last term on the left-hand side omitted. Thus (3.12) can be written schematically as

Given $\pmb{\mathbb{E}}\widetilde{R}$ , the formal solution of (3.16) is

An equation for $\pmb{\mathbb{E}}\widetilde{R}$ follows by applying $\pmb{\mathbb{E}}$ to both sides of (3.17), then bringing the $\pmb{\mathbb{E}}\widetilde{R}$ term on the right-hand side to the left:

In the electrostatic approximation, this can be solved to find

where

Thus

The quantity $\widetilde{{\mathcal{D}}}$ is a kind of ‘random dielectric function’ (although that concept is not defined), and the form (3.21) emphasizes, through its second term, that self-consistency is crucial. However, how to calculate the correlation between $\widetilde{R}$ and $\widetilde{N}$ that is required for formula (3.15) is not apparent because of the wildly nonlinear way in which random variables are mixed throughout the last term of (3.21). Thus a different approach must be taken, as I now describe.

Instead of following the route of first solving for ${\rm\Delta}\widetilde{N}$ , then averaging, in traditional closures one averages (3.11) directly. Upon recalling that $f=\langle \widetilde{N}\rangle$ , one finds

The last term on the left-hand side must be evaluated by statistical closure. (A doubly underlined term indicates that at least part of the term contains self-consistent response that would be absent for passive advection.) The problem is essentially the same one that arises in closing the equation for $R$ itself. The most general way of proceeding is to use the method of sources introduced in the classical context by Martin et al. (Reference Martin, Siggia and Rose1973, henceforth MSR), motivated by prior developments in quantum field theory. I shall not give a full description of the MSR formalism here; see Krommes (Reference Krommes2002) for extensive discussion. Briefly, the procedure is to enquire about the response to an external, nonrandom source. By looking at infinitesimal variations of that source, one can obtain functional relations between the various statistical observables that can be usefully closed. For example, add an external sourceFootnote ⁵⁹ $\widehat{{\it\eta}}$ to the right-hand side of the Klimontovich equation (3.9). The random response function $\widetilde{R}$ is defined by the functional derivative (Beran Reference Beran1968)

and the mean response is

Thus one can obtain an equation for $R$ by functionally differentiating the mean equation,

to find

The key to usefully reducing the last term on the left-hand side is to recognize that, as $\widehat{{\it\eta}}$ changes, $f[\widehat{{\it\eta}}]$ changes in concert. One can therefore use a functional chain rule:

where (3.24) was used in the last step. This formally closes the equation for $R$ in the form of a Dyson equation (Dyson Reference Dyson1949; Krommes Reference Krommes2002, Reference Krommes2009):

whereFootnote ⁶⁰

is still unknown at this point. It can be evaluated in various statistical closures such as the DIA; see § 4.

In a seminal paper, DuBois & Espedal (Reference DuBois and Espedal1978) showed that one further manipulation puts the $R$ equation into a very suggestive form that is an obvious generalization of linear theory. Note that ${\it\Sigma}$ , whatever it is, will contain two kinds of terms: ones that operate directly on $R$ , such as $\ldots \boldsymbol{\partial }R$ ; and ones that apply the electric-field operator to $R~(\ldots \pmb{\mathbb{E}}R)$ . Call those terms ${\it\Sigma}^{g}R$ and ${\it\Sigma}^{f}R$ , respectively, and group the ${\it\Sigma}^{f}$ terms with the explicit $\boldsymbol{\partial }f\boldsymbol{\cdot }\pmb{\mathbb{E}}R$ term in (3.28), thereby changing $\boldsymbol{\partial }f$ to a fluctuation-dependent quantity $\boldsymbol{\partial }\overline{f}\doteq \boldsymbol{\partial }f+\boldsymbol{\partial }{\it\delta}\overline{f}$ with $\boldsymbol{\partial }{\it\delta}\overline{f}\boldsymbol{\cdot }\pmb{\mathbb{E}}\equiv {\it\Sigma}^{f}$ . Equation (3.28) then becomes

where the renormalized particle propagator $g$ obeys

Equation (3.30) now has the same form as (3.16), but importantly it contains only statistically averaged quantities. The method of solution is the same, however, so one obtains (cf. (3.20) and (3.21))

where ${\mathcal{D}}$ is the required dielectric function. The first term of (3.32a ) describes the propagation of a renormalized test particle; the second one represents the shielding cloud that arises as the test particle polarizes the medium. The formula (3.32b ) has the same form as the linearized Vlasov dielectric, but it contains a renormalized particle propagator $g$ and a fluctuation-dependent function $\boldsymbol{\partial }\overline{f}$ that formally replaces the derivative of the background distribution $f$ . In plasmas, where there is inevitably a nonvanishing mean field (namely $f$ ), formula (3.32b ) has a nontrivial limit in which all nonlinear terms are neglected, i.e., the linearized Vlasov ${\mathcal{D}}^{(0)}$ discussed in the textbooks and illustrated by formula (3.5). In other situations with no mean field, ${\mathcal{D}}$ is still nontrivial, but the corrections to the vacuum limit ${\mathcal{D}}=1$ are entirely nonlinear, being proportional to ${\it\delta}\overline{f}$ .Footnote ⁶¹

The somewhat mysterious correction ${\it\delta}\overline{f}$ is necessary for one to recover the proper form of the matrix elements in a reduction to weak-turbulence theory (Krommes Reference Krommes2002). It is also required for the satisfaction of constraints such as energy conservation. A special case was considered by Dupree & Tetreault (Reference Dupree and Tetreault1978) and also discussed by Diamond et al. (Reference Diamond, Itoh and Itoh2010, § 4.4.4).

It follows from (3.32a ) (and was an intermediate step of the solution methodology; cf. (3.19)) that

This crucial result is a statement of dielectric shielding. Its interpretation is that the field due to a particle streaming through the random medium (that streaming being appropriately renormalized by the effects contained in ${\it\Sigma}^{g}$ , such as turbulent diffusion) is shielded by the self-consistent response of that medium, which behaves on the average as a macroscopic dielectric.

One can generalize these results to define the dielectric function for an arbitrary turbulent system. I shall show how to do this in § 4.2 after I have given further background on closure by deriving the DIA in § 4.1.

3.3. Spectral balance and dielectric response

In § 2.3 I showed that, for a particular problem of passive advection, the two-time correlation function $C({\it\tau})$ or the fluctuation spectrum $\widehat{C}({\it\omega})$ can be written schematically as (cf. (2.43))

where $g$ is the mean Green’s function, T denotes transpose [ $g^{\text{T}}(1,2)\doteq g(2,1)$ ], and $F$ is the sum of the internal and external covariances of the incoherent noise. Although I shall not prove it here, it is plausible that a similar formula holds for self-consistent systems:

for some noise covariance $F=F^{\text{int}}+F^{\text{ext}}$ , where $F^{\text{int}}$ is determined by closure. (It is intimately related to ${\it\Sigma}$ by a generalization of (2.48).) The form $R^{-1}C=FR^{\text{T}}$ can also be viewed as a Dyson equation.Footnote ⁶² This equation follows in general from the MSR formalism.Footnote ⁶³ Closures such as the DIA discussed in § 4 provide specific formulas for ${\it\Sigma}$ and $F^{\text{int}}$ . Note that each of the functions in (3.35) depends on all of the independent variables of the problem. For the immediately following discussion, it is useful to imagine that one is dealing with a kinetic plasma description, for which the independent variables are space, velocity, species, and time. Electromagnetic fields are calculated by performing weighted integrations over velocity and summations over species. By doing so, we shall see how dielectric response is contained in the balance equation (3.35).

Because $R$ comprises two terms according to (3.32a ), there are various pieces to the formula (3.35), but those collapse when one calculates the electric-field fluctuation spectrum by applying $\pmb{\mathbb{E}}$ to each of the arguments 1 and $1^{\prime }$ of (3.35) and using (3.33):

Thus the appropriately propagated incoherent noise is shielded by the coherent dielectric polarization. The beauty and elegance of this representation, which is one form of a spectral balance equation, are hopefully apparent. This balance was already known to Kadomtsev & Petviashvili (Reference Kadomtsev and Petviashvili1963) in the context of weak-turbulence theory. Some physical interpretations of the balance equation for steady-state, isotropic, neutral-fluid turbulence were given by Kraichnan (Reference Kraichnan1964a ), although he did not explicitly introduce the concept of a dielectric function.

Equations (3.35) and (3.36) are generalized fluctuation–dissipation relations. They determine the fluctuation spectrum implicitly since it occurs in all of $g$ , ${\mathcal{D}}$ , and $F$ (or in $R$ and $F$ ).

The fact that the fluctuation balance can be represented as $C=RFR^{\text{T}}$ (3.35) implies that it is unnecessary to explicitly introduce either the renormalized propagator $g$ or the dielectric function ${\mathcal{D}}$ ; one only requires approximate formulas for ${\it\Sigma}$ and $F^{\text{int}}$ . Indeed, initial-value solutions of two-time statistical closures such as the DIA merely advance $R(t;t^{\prime })$ and $C(t,t^{\prime })$ . The purposes of stressing the decomposition (3.32a ) for $R$ and the formula (3.32b ) for ${\mathcal{D}}$ are to provide insight into the physical content of $R$ and to show that the effects of turbulence can be incorporated as a natural generalization of the familiar linear theory. Awareness of the relationships between $R$ , $g$ , and ${\mathcal{D}}$ is useful in understanding the content and limitations of various approximations, beginning with the approach taken in Dupree’s early work, discussed in the next section.

The principal results of this section are summarized in table 1.

3.4. The $D_{\bot }={\it\gamma}/k_{\bot }^{2}$ formula

In the last section I argued that one representation of the steady-state fluctuation spectrum, ${\mathcal{E}}(\boldsymbol{k},{\it\omega})\doteq \langle {\it\delta}\boldsymbol{E}\boldsymbol{\cdot }{\it\delta}\boldsymbol{E}\rangle (\boldsymbol{k},{\it\omega})$ , is given by

where ${\mathcal{N}}$ is the trace of the numerator of (3.36). In general, explicitly solving this equation for ${\mathcal{E}}$ is quite complicated because ${\mathcal{E}}$ is buried in both ${\mathcal{N}}$ and ${\mathcal{D}}$ . A straightforward though nontrivial method of attack would be to solve the kinetic DIA directly as an initial-value problem; if that were done, spectral functions would be calculated as subsidiary quantities. However, one can use (3.37) directly to interpret the approach taken by Dupree, which was very popular for a time and is still instructive. Let us enquire about the wavenumber spectrum

The ${\it\omega}$ dependence of both ${\mathcal{N}}$ and ${\mathcal{D}}$ contribute to this integral, so it cannot be done in general.Footnote ⁶⁴ But assume that the system supports a weakly damped normal mode, i.e., ${\mathcal{D}}(\boldsymbol{k},{\it\omega}_{\boldsymbol{k}})=0$ for complex frequency ${\it\omega}_{\boldsymbol{k}}={\it\Omega}_{\boldsymbol{k}}+\text{i}{\it\gamma}_{\boldsymbol{k}}$ with $|{\it\gamma}_{\boldsymbol{k}}/{\it\Omega}_{\boldsymbol{k}}|\ll 1$ . Here ${\it\gamma}_{\boldsymbol{k}}$ is the total growth rate, including all nonlinear corrections; crucially, as I remarked previously, it must be negative in order that a steady state can be achieved. (Therefore nonlinear contributions to ${\it\gamma}_{\boldsymbol{k}}$ are essential in order to overcompensate a positive ${\it\gamma}_{\boldsymbol{k}}^{\text{lin}}$ and balance the positive-definite nonlinear noise.) To enquire about the contribution of a given normal mode to the spectrum, one can expand

The right-hand side of (3.38) then involves

Let us assume that only this mode contributes to the fluctuation level. (This cannot really be correct since the modes are coupled nonlinearly.) Then the ${\it\omega}$ integration can be performed:

where I used the fact that ${\it\gamma}_{\boldsymbol{k}}<0$ . This is still an implicit equation for ${\mathcal{E}}(\boldsymbol{k})$ , which enters in both ${\mathcal{N}}$ and ${\mathcal{D}}^{\prime }$ . However, if the nonlinear noise is neglected altogether (see further discussion at the end of this section), one is led to the saturation condition

If one further writes

where ${\it\eta}_{\boldsymbol{k}}$ describesFootnote ⁶⁵ the nonlinear effects, then the saturation condition becomes

describing a balance between the nonlinear coherent damping and the linear drive. For the case of magnetized plasma, for which fluid elements move cross-field with $\boldsymbol{E}\times \boldsymbol{B}$ drifts, the approximations of Dupree (Reference Dupree1967) and Weinstock (Reference Weinstock1968, Reference Weinstock1969, Reference Weinstock1970) for ${\it\eta}_{\boldsymbol{k}}$ amounted to saying

for some $D_{\bot }$ that is a measure of cross-field transport. Not accidentally, this is the same formula for ${\it\eta}_{\boldsymbol{k}}$ that we found in the small- ${\it\tau}_{\text{ac}}$ passive oscillator discussed in § 2.3.

The virtue of (3.45) is that in conjunction with (3.44) it seems to directly determine the diffusion coefficient to be $D_{\bot }={\it\gamma}_{\boldsymbol{k}}^{\text{lin}}/k_{\bot }^{2}$ . An immediate objection is that this cannot literally be so because $D_{\bot }$ was assumed to be $\boldsymbol{k}$ -independent but ${\it\gamma}_{\boldsymbol{k}}^{\text{lin}}$ is not proportional to $k_{\bot }^{2}$ in general. (Furthermore, not all ${\it\gamma}_{\boldsymbol{k}}^{\text{lin}}$ need be positive.) The more general (3.44) does not necessarily have this problem because ${\it\eta}_{\boldsymbol{k}}$ is a functional of the wavenumber spectrum, which could in principle adjust so that (3.44) holds for all  $\boldsymbol{k}$ . However, the diffusive approximation is valid only at small $\boldsymbol{k}$ ,Footnote ⁶⁶ so it is pointless to worry about detailed $\boldsymbol{k}$ dependence. The saturation criterion is best evaluated at some typicalFootnote ⁶⁷ wavevector $\overline{\boldsymbol{k}}$ (one must have ${\it\gamma}_{\overline{\boldsymbol{k}}}^{\text{lin}}>0$ ); then a gross estimate is

This formula, essentially also given by Kadomtsev (p. 107),Footnote ⁶⁸ is crude but effective. It states that the fluctuation and transport levels scale with a typical linear growth rate, which can be tested. It is often reasonable.

This is not the place for a thorough review of the subsequent developments of RBT, which addressed various refinements and applications of the basic theme. Regarding its validity, there are really two issues: How does one adequately approximate the coherent damping  ${\it\eta}_{\boldsymbol{k}}$ ? And how important is the nonlinear noise  ${\mathcal{N}}_{\boldsymbol{k}}$ ? Note that there are important situations in which ${\mathcal{N}}$  can be neglected. In weak-turbulence theory, where the spectrum comprises weakly coupled waves, the lowest-order approximation to  ${\mathcal{N}}$ vanishes unless the resonance condition ${\it\Omega}_{\boldsymbol{k}}+{\it\Omega}_{\boldsymbol{p}}+{\it\Omega}_{\boldsymbol{q}}=0$ is satisfied for $\boldsymbol{k}+\boldsymbol{p}+\boldsymbol{q}=\mathbf{0}$ . Linear dispersion relations of ‘nondecay’ type (Kadomtsev Reference Kadomtsev1965; Sagdeev & Galeev Reference Sagdeev and Galeev1969) do not (by definition) obey that condition. In that case, nontrivial steady states require that (3.42) be satisfied (unless higher-order, e.g., four-wave, interactions are considered). The nonlinear corrections to  ${\mathcal{D}}$ in this case describe wave–wave–particle (and higher-order) interactions and can be systematically calculated by perturbation theory.Footnote ⁶⁹ Kadomtsev (Reference Kadomtsev1965) and later Sagdeev & Galeev (Reference Sagdeev and Galeev1969) and Sagdeev (Reference Sagdeev1979) discussed some instances of this situation; an example of a greatly refined analysis of such effects is the work of Horton & Choi (Reference Horton and Choi1979).

In later research, Dupree attempted to take the nonlinear noise into consideration. A difficulty is that at the kinetic level the incoherent noise is associated with several distinct kinds of physical processes: renormalized free-streaming particles, which enter into discussions by Dupree (Reference Dupree1972) and later authors (DuBois & Espedal Reference DuBois and Espedal1978; Boutros-Ghali & Dupree Reference Boutros-Ghali and Dupree1981) of phase-space granulationsFootnote ⁷⁰ in plasmas; and fluid mode-coupling effects that are not sensitive to the details of the velocity space. Those effects must not be neglected; failure to recognize this can easily lead to paradoxes or misconceptions (Krommes & Kim Reference Krommes and Kim1988). The mode coupling is exposed more transparently in statistical closures such as the DIA applied directly to fluid models (Weinstock & Williams Reference Weinstock and Williams1971), for which there has been much work (some of which is described below).

4.1. Heuristic derivation of the DIA

In the present section, I shall present a route to the DIA that seems to be modestly new. It is designed to emphasize the principal intuitive points that all beginners must grasp, which are that the statistics of turbulence are non-Gaussian and the goal of closure is to calculate those non-Gaussian corrections self-consistently. The method essentially provides an $\boldsymbol{x}$ -space version of Kraichnan’s original approach (Kraichnan Reference Kraichnan1959), which was anchored in $\boldsymbol{k}$ space.

If one were dealing with Gaussian statistics, the Furutsu–Novikov theorem (Novikov’s theorem for brevity) would be useful. That theorem (also known as Gaussian integration by parts) states that for a centred Gaussian random variable $\widetilde{{\it\psi}}^{\prime }$ (with covariance  $C$ ) and a functional  ${\mathcal{F}}$ of that variable, one has

where ${\it\delta}/{\it\delta}\widetilde{{\it\psi}}^{\prime }$ denotes functional differentiation and I use  $\widetilde{{\it\psi}}^{\prime }$ in this section as an alternative notation for the fluctuation $\widetilde{{\it\psi}}-\langle \widetilde{{\it\psi}}\rangle$ . See appendix B for the proof. Cook (Reference Cook1978) discussed the use of Novikov’s theorem to obtain closures for passive advection; the present work focuses on self-consistent problems. Thus the goal is to approximate the triplet correlation function  $T$ introduced in § 2; see (2.8) and (2.9). One can cast  $T$ into the form of the left-hand side of Novikov’s theorem:

where ${\mathcal{F}}[\widetilde{{\it\psi}}^{\prime }]=\widetilde{C}(2,3)$ with $\widetilde{C}(2,3)\doteq \widetilde{{\it\psi}}^{\prime }(2)\widetilde{{\it\psi}}^{\prime }(3)$ . However, one cannot trivially employ Novikov’s theorem to evaluate  $T$ because that theorem applies only to Gaussian statistics and, as I discussed in § 2.1, $\widetilde{{\it\psi}}^{\prime }$  becomes non-Gaussian already after the first time step. If one attempts to assert that $\widetilde{{\it\psi}}^{\prime }$ is Gaussian, then one finds the trivial result $T=0$ (third- and higher-order cumulants of a Gaussian vanish), consistent with the direct use of Novikov’s theorem:

because $\langle \widetilde{{\it\psi}}^{\prime }\rangle \equiv 0$ .

Non-Gaussian statistics are induced by the last $U_{3}$ term in (2.7). Let us write

where $\widetilde{{\it\psi}}_{\text{G}}^{\prime }$ is a Gaussian fluctuation with covariance $\langle \widetilde{{\it\psi}}_{\text{G}}^{\prime }(1)\widetilde{{\it\psi}}_{\text{G}}^{\prime }(1^{\prime })\rangle =C(1,1^{\prime })$ (this will later be determined self-consistently) and ${\rm\Delta}\widetilde{{\it\psi}}^{\prime }$  is a non-Gaussian correction. One therefore has

The first term vanishes because it is the third-order cumulant of a Gaussian process. The remaining three explicit terms have been written in a form suitable for application of Novikov’s theorem; the terms in square brackets define  ${\mathcal{F}}$ . (I shall actually evaluate the second explicit term as $\langle [{\rm\Delta}\widetilde{{\it\psi}}^{\prime }(3)\widetilde{{\it\psi}}_{\text{G}}^{\prime }(2)]\widetilde{{\it\psi}}_{\text{G}}^{\prime }(1^{\prime })\rangle$ .) The DIA is defined by the neglect of the terms of $O({\rm\Delta}\widetilde{{\it\psi}}^{\prime 2})$ . (Unfortunately, ${\rm\Delta}\widetilde{{\it\psi}}^{\prime }$  is not small in general, which leads to subtle deficiencies in the DIA; see § 5.1 for further discussion.) An equation for  ${\rm\Delta}\widetilde{{\it\psi}}^{\prime }$ follows from (2.7); correct to first order in  ${\rm\Delta}\widetilde{{\it\psi}}^{\prime }$ , it is

This equation can be represented more compactly by introducing the random infinitesimal response function

which according to (2.5) obeys

Here $\widetilde{R}$  is a functional of  $\widetilde{{\it\psi}}$ : $\widetilde{R}=\widetilde{R}[\widetilde{{\it\psi}}]$ . The quantity $\widetilde{R}[\langle {\it\psi}\rangle +\widetilde{{\it\psi}}_{\text{G}}^{\prime }]\equiv \widetilde{R}_{\text{G}}$ is Green’s function for the terms explicitly involving  ${\rm\Delta}\widetilde{{\it\psi}}^{\prime }$ in (4.6); thus (4.6) has the formal solution

This describes how non-Gaussian statistics are generated by the nonlinear coupling of Gaussian fluctuations.

One can now work out all of the $O({\rm\Delta}\widetilde{{\it\psi}}^{\prime })$ terms in (4.5) by approximating the consequences of Novikov’s theorem, assuming that all terms are functionals of $\widetilde{{\it\psi}}_{\text{G}}^{\prime }$ . One has

From (4.9), the functional derivative is approximately

Upon differentiating the identity

which leads to ${\rm\Delta}\widetilde{R}=-\widetilde{R}{\rm\Delta}(\widetilde{R}^{-1})\widetilde{R}$ , and noting from (4.8) that

one finds

For use in (4.10), one must average (4.11) using the result (4.14b ). Since all quantities have been expressed in terms of  $\widetilde{{\it\psi}}_{\text{G}}^{\prime }$ , which is Gaussian by assumption, the averages can be performed by using the functional generalization of results like

(Here ${\it\sigma}^{2}\doteq \langle \widetilde{{\it\psi}}^{\prime }\widetilde{{\it\psi}}^{\prime }\rangle$ ; see appendix B.) However, terms that bring in extra factors of  $U_{3}$ (any of the terms that arise from differentiating $\widetilde{R}^{-1}$ ) should be neglected because the role of  $U_{3}$ is to induce allegedly small non-Gaussian statistics and we have already (for better or worse) neglected terms of higher than first order in that correction. This amounts to saying that one should retain only the first term of (4.15) and should approximate $\widetilde{R}_{\text{G}}\approx \langle \widetilde{R}\rangle \equiv R$ . Only the first line of (4.11) then survives the average, and one finds a contribution to $U_{3}T/2$ of the form $-{\it\Sigma}(1;\overline{1})C(\overline{1},1^{\prime })/2$ , where

The second and third terms of $O({\rm\Delta}\widetilde{{\it\psi}}^{\prime })$ in (4.5) can be worked out in a similar fashion. By symmetry, the second term is equal to the first. The third term gives a distinct contribution to the term $U_{3}T/2$ in (2.8) of the form $F^{\text{int}}(1,\overline{1})R(1^{\prime };\overline{1})$ , where

This function can be interpreted as the covariance of an internal emission or noise term arising from the nonlinearity.

Next, one must evaluate the cross-correlation $\langle \,\widetilde{f}^{\text{ext}\prime }(1)\widetilde{{\it\psi}}^{\prime }(1^{\prime })\rangle$ required in (2.8). This can be done exactly if one assumes that  $\widetilde{{\it\psi}}$ evolves from Gaussian initial conditions $\widetilde{{\it\psi}}_{0}$ that are uncorrelated with the Gaussian $\widetilde{f}^{\text{ext}}$ . From Novikov’s theorem, one has

(Here the notation  $\underline{1}$ denotes the set of all independent variables except time.) By assumption one has $\langle \widetilde{{\it\psi}}_{0}\,\widetilde{f}^{\text{ext}\prime }\rangle =0$ , and by definition one has $\langle {\it\delta}\widetilde{{\it\psi}}(1^{\prime })/{\it\delta}\widetilde{f}^{\text{ext}}(\overline{1})\rangle =R(1^{\prime };\overline{1})$ . ThusFootnote ⁷¹

This leads to the covariance equation

Finally, one can obtain a closed equation for the mean response function  $R$ by noting from (4.8) that $R$ obeys

If one approximates $\widetilde{{\it\psi}}^{\prime }\approx \widetilde{{\it\psi}}_{\text{G}}^{\prime }$ , the last term can be evaluated by Novikov’s theorem:

Thus, upon recalling (4.16), one finds that the last term on the left-hand side of (4.21) reduces to ${\it\Sigma}(1;\overline{1})R(\overline{1};1^{\prime })$ ; therefore

This result shows that the operator on the left-hand side of (4.20) is just  $R^{-1}$ . Upon formally solving (4.20), one is then led to the symmetrical form

or concisely $C=RFR^{\text{T}}$ .

This closure, consisting of (2.6) for the mean field, equations (4.20) or (4.24) for the covariance, and (4.23) for the response function, along with the definitions (4.16) for  ${\it\Sigma}$ and (4.17) for $F^{\text{int}}$ , is precisely the DIA. A thorough discussion and alternative derivation of the DIA for Vlasov turbulence was given by DuBois & Pesme (Reference DuBois and Pesme1985).

For homogeneous turbulence, the DIA equations are often written in $\boldsymbol{k}$ space. If the quadratic PDE for a fluid description (no velocity space) is written asFootnote ⁷²

then in the DIA

and

where $M\equiv U_{3}$ . If ${\it\psi}_{\boldsymbol{k}}$ has been normalized in such a way that $M_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}+M_{\boldsymbol{p}\boldsymbol{q}\boldsymbol{k}}+M_{\boldsymbol{q}\boldsymbol{k}\boldsymbol{p}}=0$ , then one can verify that

This is a statement of ‘energy’ conservation by the nonlinear terms: $\sum _{\boldsymbol{k}}\!\left.C_{\boldsymbol{k}}(t,t)\right|_{\text{n.l.}}=0$ . (The proof can as easily be conducted in $\boldsymbol{x}$ space.) This result is inherited from the primitive amplitude equation; it is an essential property that any credible closure should satisfy.

The pros and cons of the DIA have been discussed at great length (Krommes Reference Krommes2002, and references therein). Significantly, it is known that the DIA is realizable; it provides the exact description of the second-order statistics for certain stochastic models. Some general discussion and original references are given by Kraichnan (Reference Kraichnan and Sirovich1991). Although the first demonstration of realizability involved a so-called random-coupling model (Kraichnan Reference Kraichnan1961), a more intuitive generalized Langevin representation was later displayed by Kraichnan (Reference Kraichnan1970) and Leith (Reference Leith1971). For homogeneous statistics, that model is

where $\widehat{{\it\Sigma}}_{\boldsymbol{k}}$  has the DIA form (4.26),

and $\widetilde{{\it\xi}}_{\boldsymbol{p}}$ is an auxiliary random variable having covariance  $C_{\boldsymbol{p}}$ . It is left as an instructive exercise for the reader to verify that the equations for the response function and covariance that follow from (4.29) reproduce the DIA.

In addition to realizability (which guarantees a positive-semidefinite covariance), the DIA properly reduces to weak-turbulence theory at the level of wave–wave–particle or three-wave coupling (but omits some terms associated with higher-order processesFootnote ⁷³ ). Its self-consistent, nonlinearly energy-conserving determination of  $R$ and  $C$ is a definite plus, and it has enjoyed considerable successes (Kraichnan Reference Kraichnan1964a ; Krommes Reference Krommes2002). However, the approximation described in the present section, and others like it, merely define algorithms; for arbitrarily strong turbulence, they cannot be systematically defended because one does not have a controlled estimate of the size of the non-Gaussian effects. A well-known qualitative issue related to the lack of random Galilean invariance is discussed in § 5.1. The DIA cannot deal properly with mixtures of turbulent motions and coherent structures. And even for regimes in which the DIA is expected to be good, it is quite complicated. Alternative, simpler approximations are mentioned later in the article.

4.2. Perturbative algorithm for the dielectric function of a general turbulent medium in the direct-interaction approximation

As I discussed after (3.29), the essence of the $R{-}g{-}{\mathcal{D}}$ factorization obtained by DuBois and Espedal was the grouping of the nonlinear terms in the equation for the response function into two parts. In their words, they ‘temporarily consider [  $f$ ] and [ $\boldsymbol{E}$ ] as independent functional variables’ and replace ${\it\delta}/{\it\delta}f$ in (3.29) by

However, if one tries the same approach for the Navier–Stokes equation, the nonlinear term of which is essentially (modulo the constraint of incompressibility) $\langle {\it\delta}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\delta}\boldsymbol{u}\rangle$ , confusion arises because the two  $\boldsymbol{u}$ ’s are obviously not functionally independent. In fact, the assumption of functional independence is problematic even for the plasma case because $\boldsymbol{E}$ and  $f$ are related by $\boldsymbol{E}=\pmb{\mathbb{E}}f$ . Operationally, what is really being done is keeping track of which term is responsible for the advection, then organizing the terms in such a way that the self-consistent response is split off. This leads one to a formula for the dielectric of a general (quadratically nonlinear) system. I shall illustrate the procedure in the DIA.

The approach will be to calculate nonlinear corrections to $R$ perturbatively, then to heuristically renormalize by replacing certain zeroth- or second-order quantities by their exact counterparts. This is one of the standard elementary procedures for obtaining the DIA itself (Kraichnan Reference Kraichnan1977; DuBois & Pesme Reference DuBois and Pesme1985); where the methodology differs is that the symmetrized form of the mode-coupling coefficient will not be used here. Working with the unsymmetrized form allows one to distinguish terms associated with advection and with self-consistent response; this leads to a procedure that generalizes the approach of DuBois and Espedal.

For conciseness, I shall use a subscript notation that looks like a Fourier representation. In fact, however, the notation works in either $\boldsymbol{x}$ space or $\boldsymbol{k}$ space. A subscript  $k$ could be a combination of discrete index (say for a Cartesian velocity component) and continuous variable (either  $\boldsymbol{k}$ or  $\boldsymbol{x}$ , and  $t$ ). Again, sums over repeated indices are assumed. Let the dynamical equation be written

The hatted index identifies the advection term. Here $\widehat{M}_{}$ is not symmetric in its last two indices. Let me again underline terms that do not exist for passive advection. Then the random response function obeys

At first order, one has

Upon averaging (4.33), one finds that through second order $R$  obeys

Substitution of (4.34) gives

In the last expression, the underlined quantities identify terms in which a hatted index is associated with  $R$ . Applying a hatted index to  $R$ is the generalization of applying the $\pmb{\mathbb{E}}$ operator in the original manipulations of DuBois and Espedal; i.e., it is related to the last term in (4.31). I shall group terms involving a hatted index on  $R$ separately, and I shall also interchange indices so that the names of the indices on each  $C$ and  $R$ match. Furthermore, one can heuristically renormalize by dropping the superscripts on all quantities. Then one obtains

where

and

(Here $M_{k\,p\,q}\doteq \widehat{M}_{k\,p\,q}+\widehat{M}_{k\,q\,p}$ and redundant hats were dropped in each of the final expressions.) The quantity  $g$ is a generalized ‘particle’ or passive fluid-element propagator, while ${\it\Sigma}^{f}$ generalizes the plasma $\boldsymbol{\partial }{\it\delta}\overline{f}\boldsymbol{\cdot }\pmb{\mathbb{E}}$ .

These general formulas work for both kinetic problems, where the indices are six-dimensional (plus time), and fluid problems, where there is no velocity variable. When the latter case is applied to a scalar field, the operations become purely multiplicative; no summations or integrations occur. Then one can define a dielectric function simply by the ratio of particle to total response:

where ${\it\Sigma}\doteq {\it\Sigma}^{g}+{\it\Sigma}^{f}$ and formula (4.41) holds for the case of zero mean field. An important example of this situation is the guiding-centre plasma model (a system of charged rods moving across a very large magnetic field with $\boldsymbol{E}\times \boldsymbol{B}$ motion). In a most elegant paper, Taylor (Reference Taylor1974) found the proper form of  ${\mathcal{D}}$ for the guiding-centre model ( $g_{0}^{-1}=-\text{i}{\it\omega}$ ) in thermal equilibrium. He used a linear response formalism that applied only to Gibbsian thermal equilibrium; subsequently Krommes & Similon (Reference Krommes and Similon1980) showed how Taylor’s result follows from the DIA. The present formalism recovers those results.

It should be noted that self-consistency is not only responsible for dielectric shielding, it also affects the form of the renormalized particle propagator  $g$ . Whereas passive advection leads to a generalized diffusion operator for  ${\it\Sigma}^{g}$ , polarization modifies that operator, producing among other things a smaller effective diffusion coefficient. This was discussed by Krommes and Similon for the guiding-centre model; the general polarization effect is described by the second term in (4.39a ). (The first term in that equation would be the sole one arising in passive advection because it involves a  $C$ with two hatted indices, i.e., the covariance of the advecting velocity.)

The method of construction guarantees that the formulas (3.32b ) and (3.33) for the dielectric function and shielded response function hold provided that one can identify in  $\widehat{M}_{}$ the analogue of the  $\pmb{\mathbb{E}}$ and $\boldsymbol{\partial }$ operators; cf. (3.10b ). This appears to be most useful for cases with a two-level description – for example, involving a velocity space that, with the $\widehat{p}$ index in (4.32), must be integrated over, as in Klimontovich theory. For the plasma case, the results are recorded in table 2. For fluid cases without a phase-space variable (such as the guiding-centre plasma or the Navier–Stokes equation), ${\mathcal{D}}$  can still be defined although the implications of this are not fully understood. It is noteworthy that dielectric response is mentioned nowhere in the most recent book on homogeneous turbulence (McComb Reference McComb2014). Indeed, it is unnecessary to explicitly introduce  ${\mathcal{D}}$ since the response function  $R$ carries the relevant information. Nevertheless, it is clear that nonlinear dielectric response operates in Navier–Stokes turbulence. It is not excluded that new perspectives on certain nonlinear processes in fluid turbulence may be obtained by expressing them in terms of the nonlinear  ${\mathcal{D}}$ .

The key DIA formulas obtained in this section are summarized in table 2.

5.1. Kadomtsev and the DIA

The most general mathematical formulation of the last statement was presumably not known to Kadomtsev, although it is probable that he appreciated it at least intuitively. What is clear is that Kadomtsev understood an important deficiency of the DIA,Footnote ⁷⁴ which I shall now describe.

It was already known to Kraichnan in 1959 that the DIA did not predict the famous $k^{-5/3}$ inertial-range spectrum of Navier–Stokes turbulence that follows from the $\boldsymbol{x}$ -space self-similarity considerations of Kolmogorov (Reference Kolmogorov1941). (For some careful discussion, see Frisch (Reference Frisch1995).) Instead, the DIA predicts $k^{-3/2}$ , and for a time Kraichnan entertained the possibility that Kolmogorov was incorrect. However, the tidal-basin measurements of Grant, Stewart & Moilliet (Reference Grant, Stewart and Moilliet1962) confirmed Kolmogorov’s prediction, so Kraichnan was led to reconsider the structure of the DIA. In 1964 he explained (Kraichnan Reference Kraichnan1964c ) that the difficulty is that the DIA is not invariant to random Galilean transformations. Ordinary Galilean invariance implies that, under the addition of an infinite-wavelength flow, short-scale eddies should be carried along unchanged. In fact, the DIA is Galilean invariant in that respect (McComb Reference McComb2014). But if one considers an ensemble of such flows with random direction and strength, phase-mixing effects arise as a consequence of ensemble averaging. Those are captured correctly by the two-time  $R$ , as in the stochastic-oscillator model. However, because of the intimate way in which the DIA couples $R(t;t^{\prime })$ and  $C(t,t^{\prime })$ to $C(t,t)$ (discussed in mathematical detail by Kraichnan (Reference Kraichnan1964c ); see also footnote 33 on page 16), that phase mixing introduces spurious distortion into the small-scale eddies (see figure 4), which shows up as an incorrect inertial-range spectrum.

Figure 4. Illustration of difficulties with random Galilean invariance. In the top panel, an infinite-wavelength flow (no shear) translates an eddy unchanged. In the bottom panel, the average over an ensemble of such flows produces in the DIA a spurious distortion of a typical eddy, violating random Galilean invariance.

Also in 1964, Kadomtsev’s review appeared in Russian. It contains, after some discussion about the differences between resonant and adiabatic wave interactions and the role of wavepackets, the following observation (p. 55):

Kadomtsev’s discussion was not as crisply technical as Kraichnan’s, but certainly he understood the intuitive essence of the issue. Because of the publication dates, writing styles, and lack of cross-references, it seems clear that the analyses were independent.

The issue of random Galilean invariance has received much attention and fostered many attempts at a cure. In general, it has been understood that the problem cannot be cured within a systematically renormalized Eulerian formulation (which includes the MSR formalism). A heuristic Eulerian approximation that does respect random Galilean invariance is the eddy-damped quasi-normal Markovian closure (Orszag Reference Orszag1970), and various more or less successful Lagrangian schemes have also been discussed (Kraichnan Reference Kraichnan1965; McComb Reference McComb2014). But it should be emphasized that in many problems of practical interest the issue is only of secondary importance. In particular, turbulent transport is usually dominated by the long-wavelength, energy-containing fluctuations, in which case it is insensitive to the precise form of the inertial range. Although there are situations (for example, the intermittent mixing of contaminants) in which one cares about detailed inertial-range structure, they are well beyond the scope of this article; see Falkovich, Gawȩdzki & Vergassola (Reference Falkovich, Gawȩdzki and Vergassola2001). Fusion plasma physicists should not be overly concerned with inertial-range spectra, although the shapes of those spectra sometimes figure in discussions of closures for large-eddy simulations (see footnote 9 on page 6).

5.2. Drift-wave saturation and the mixing-length formula

The DIA for Vlasov turbulence was first discussed by Orszag & Kraichnan (Reference Orszag and Kraichnan1967); further contributions were made by DuBois & Espedal (Reference DuBois and Espedal1978), Krommes (Reference Krommes1978), and DuBois & Pesme (Reference DuBois and Pesme1985). However, kinetic DIA calculations proved mostly too hard to perform at the time. Plasma physicists slowly latched on to the use of the DIA for fluid models (Krommes Reference Krommes1982; Waltz Reference Waltz1983; Sudan & Pfirsch Reference Sudan and Pfirsch1985), though even there it is rather complicated for anisotropic turbulence. That difficulty fostered the development of DIA-based Markovian closures for plasma turbulence (Waltz Reference Waltz1983; Bowman Reference Bowman1992; Bowman et al. Reference Bowman, Krommes and Ottaviani1993; Bowman & Krommes Reference Bowman and Krommes1997), and those can be used to quantify some important heuristic drift-wave scalings discussed by Kadomtsev. Here I shall give a brief introduction. Definitions of the symbols used in the following discussion are given in appendix A.

5.3. The impact of zonal flows on the Hasegawa–Mima equation

Zonal potentials are ones with $k_{y}=k_{\text{z}}=0$ (in slab geometry) or with no variation within a flux surface (in toroidal geometry). The resulting $\boldsymbol{E}\times \boldsymbol{B}$ drifts (mostly poloidal) are called zonal flows; they are important in multiple contexts.Footnote ⁸⁸ In fusion plasmas, for example, they play an important role in the regulation of the levels of drift-wave turbulence and transport, either by shearing of turbulent eddies (Diamond et al. Reference Diamond, Itoh, Itoh and Hahm2005) or by catalysing the coupling of unstable fluctuations to stable normal modes (Hatch et al. Reference Hatch, Terry, Jenko, Merz, Pueschel, Nevins and Wang2011b ).

Zonal flows are strongly non-adiabatic: because they have $k_{\Vert }=0$ , the adiabatic ordering ${\it\omega}/k_{\Vert }v_{\text{t}e}\ll 1$ is violated. Dorland & Hammett (Reference Dorland and Hammett1993) emphasized that Poisson’s equation must be modified as a consequence, giving rise to an enhanced response for the zonal modes. In the simplest model, the electron response is taken to vanish altogether for zonal modes. This can be implemented by writing ${\it\delta}n_{e,\boldsymbol{k}}/n_{e}={\it\alpha}_{\boldsymbol{k}}{\it\delta}{\it\varphi}_{\boldsymbol{k}}$ , where ${\it\alpha}_{\boldsymbol{k}}$ vanishes for zonal modes and is equal to 1 otherwise. A consequence is that the derivation of the Hasegawa–Mima equation must be reconsidered; one finds instead the modified Hasegawa–Mima equation

This equation has been used in various studies related to the physics of zonal flows (Krommes & Kim Reference Krommes and Kim2000; Parker & Krommes Reference Parker and Krommes2013); it is further discussed in § 6.2.

5.4. Further developments in closure theory

For other facets of closure theory, I point the reader to the discussion and references in Krommes (Reference Krommes2002) and McComb (Reference McComb2014). The latter focuses on homogeneous turbulence. In the next section, I shall describe some aspects of the theory of inhomogeneous turbulence.

6.1. Reynolds stress

Kadomtsev (p. 59) wrote,

In the last paragraph, Kadomtsev was referring to the modelling of the so-called Reynolds stress by means of ‘the mixing-length concept introduced by Prandtl (Reference Prandtl1925)’. In its original interpretation, Reynolds stress  ${\bf\tau}$ described the contribution of fluctuations in the fluid velocity to the mean Navier–Stokes momentum equation:

whereFootnote ⁸⁹ ${\bf\tau}\doteq {\it\rho}_{\text{m}}\langle {\it\delta}\boldsymbol{u}\,{\it\delta}\boldsymbol{u}\rangle$ . Prandtl’s estimate for the size of a fluctuating velocity component in a free jet was ${\it\delta}u\approx \ell |\partial \langle u\rangle /\partial y|$ , where $u$  is the longitudinal velocity, $y$  is the transverse coordinate, and $\ell$  is the transverse mixing length. Kadomtsev (p. 59) reviewed a solution by Tollmien of the steady-state system of momentum and continuity equations that results when one takes $\ell =cx$ , where $c$  is a fitting constant; agreement with experiment is excellent when $c$  is chosen appropriately.

Note that it is only the divergence of the stress tensor that contributes to the mean momentum equation. For homogeneous turbulence, that divergence vanishes, so inhomogeneity is essential.Footnote ⁹⁰ As an example, let us discuss the off-diagonal Reynolds stress due to $\boldsymbol{E}\times \boldsymbol{B}$ flows in a two-dimensional model with constant magnetic field. One has ${\it\tau}_{xy}\propto \langle {\it\delta}V_{E,x}{\it\delta}V_{E,y}\rangle$ with $\boldsymbol{V}_{E}\propto \,\widehat{\boldsymbol{b}}\times \boldsymbol{{\rm\nabla}}\widetilde{{\it\phi}}$ and $\widetilde{{\it\phi}}$  being the random electrostatic potential. Contours of constant  $\widetilde{{\it\phi}}$ are the streamlines of the velocity. Consider an elliptical ‘eddy’, of major axis  $a$ and minor axis  $b$ , inclined at random angle  $\widetilde{{\it\theta}}$ from the $x$  direction:

where $\widetilde{c}\equiv \cos \widetilde{{\it\theta}}$ and $\widetilde{s}\equiv \sin \widetilde{{\it\theta}}$ . It is straightforward to calculate $\widetilde{V}_{E,x}\propto -\partial _{\text{y}}\widetilde{{\it\phi}}$ , $\widetilde{V}_{E,y}\propto \partial _{\text{x}}\widetilde{{\it\phi}}$ , and the product $\widetilde{{\it\tau}}_{xy}\doteq {\it\delta}\widetilde{V}_{E,x}{\it\delta}\widetilde{V}_{E,y}$ . I assume that $y$  is a direction of symmetry.Footnote ⁹¹ The purely $x$ -dependent part of the result (which after averaging gives the contribution to $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\tau}$ ) is

For the case of circular eddies ( $a^{2}=b^{2}=r^{2}$ ), this vanishes.Footnote ⁹² Even for elliptical eddies, ${\it\tau}_{xy}\doteq \langle \widetilde{{\it\tau}}_{xy}\rangle$ vanishes by symmetry when it is averaged over a symmetrical distribution in  ${\it\theta}$ . But if $\widetilde{{\it\theta}}$  is distributed according to the PDF  $f({\it\theta})=\,f_{\text{s}}({\it\theta})+{\it\epsilon}f_{\text{a}}({\it\theta})$ , where s and a denote symmetric and antisymmetric parts, then the antisymmetric part contributes and one finds ${\it\tau}_{xy}\propto {\it\epsilon}{\it\Delta}x^{2}$ , where ${\it\Delta}\doteq b^{-2}-a^{-2}$ . This example shows the importance of asymmetry in producing nonvanishing Reynolds stress, and it suggests that toroidal devices with up–down asymmetry may, in the presence of microturbulence, have larger Reynolds stresses. Those may be important in limiting the size of the turbulence and in the generation of intrinsic rotation.Footnote ⁹³

It is well known that the deformation tensor $\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ of a general flow field  $\boldsymbol{u}(\boldsymbol{x})$ can be decomposed into a symmetric part (the rate-of-strain tensor  $\unicode[STIX]{x1D64E}$ , associated with both stretching at constant volume and possibly compression) and an antisymmetric part (associated with vorticity and leading to rigid-body rotation). Therefore tilted elliptical eddies, hence Reynolds stress, can arise by subjecting a circular eddy to an inhomogeneous velocity field. One source of inhomogeneity is the generation of zonal flows due to the (divergence of the) Reynolds stress. As an example, consider a slab zonal flow with constant shear:Footnote ⁹⁴ $\boldsymbol{u}={\it\Omega}_{0}x\widehat{\boldsymbol{y}}$ , with $\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}={\it\Omega}_{0}\widehat{\boldsymbol{z}}$ and $\unicode[STIX]{x1D64E}={\it\Omega}_{0}(\widehat{\boldsymbol{x}}\,\widehat{\boldsymbol{y}}+\widehat{\boldsymbol{y}}\,\widehat{\boldsymbol{x}})/2$ (the principal axes of which are $(1,\hspace{2.22198pt}\pm 1)^{\text{T}}/\sqrt{2}$  ). Tilt arises from the vorticity, and ellipticity (at constant area) arises from the straining. Thus one has a classic problem of the chicken and the egg. Inhomogeneous Reynolds stress drives zonal flows, and such stress can arise from the modulation of homogeneous turbulence by zonal flows; a feedback loop is indicated. But given homogeneous turbulence, there is no divergence of the Reynolds stress, hence no zonal-flow generation (see footnote 90), hence no inhomogeneity, hence no effect due to Reynolds stress. However, such homogeneous turbulence can be unstable to the generation of inhomogeneous flows, as I shall describe in the next section. Thus inhomogeneous states can result in which zonal flows and microturbulence are coupled self-consistently.

6.2. Infinitesimal response, symmetry breaking, and zonostrophic instability

A basic question is this: Under what circumstances is a state of steady homogeneous turbulence unstable to the symmetry-breaking generation of inhomogeneity? This can be addressed by linearizing around the homogeneous state (allowing for inhomogeneous perturbations) and enquiring about the onset of a zero eigenvalue that solves a certain dispersion relation. I shall show, expanding on the original beautiful paper of Carnevale & Martin (Reference Carnevale and Martin1982), that that dispersion relation is intimately related to a certain approximation to the inverse of the infinitesimal response function  $R$ for the homogeneous turbulent state. In the theory of zonal-flow generation, the instability of homogeneous turbulence is called the zonostrophic instability (Srinivasan & Young Reference Srinivasan and Young2012).

Mention of a zero eigenvalue in the context of homogeneous turbulence might be confusing. Recall that in earlier sections we saw how the nonlinear dielectric function participates in the spectral balance of the turbulence by shielding nonlinear noise, which I argued is important and cannot generally be neglected. Also recall that in steady state the zeros of ${\mathcal{D}}(\boldsymbol{k},{\it\omega})$ or $R^{-1}(\boldsymbol{k},{\it\omega})$ are in the lower half of the ${\it\omega}$  plane, not on the real line. Now the dielectric function (or inverse response function) referred to here relates internal incoherent noise to the covariance of the fluctuations (for homogeneous turbulence, there is no mean flow). But if one admits the possibility of inhomogeneity, there will arise an effective dielectric tensor that operates on an extended vector containing not only fluctuations but also mean fields. At the symmetry-breaking bifurcation of homogeneous turbulence to the generation of mean fields, that tensor will develop a zero eigenvalue.

If one writes $Z\equiv \langle {\it\psi}\rangle$ , (2.6) for zero forcing becomes

The linearization of this equation is

In particular, I shall linearize around a homogeneous state for which it is assumed that $Z=0$ ; the terms that vanish with this assumption are underlined. The (generally inhomogeneous) fluctuation spectrum is assumed to obey the Dyson equation

Linear perturbations to the spectrum obey

Doubly underlined terms are $O(C\,{\rm\Delta}C)$ ; they describe fluctuation–fluctuation interactions. Let us neglect those and focus on the coupling between the mean fields and the turbulence. Because ${\rm\Delta}R=-R{\rm\Delta}(R^{-1})R$ and

one has

The doubly underlined term contains pieces of order $C\,{\rm\Delta}R$ and $R\,{\rm\Delta}C$ ; those can be ignored because when multiplied by  $F$ , as required for (6.7), they contribute effects of order $C^{2}$ or  $C\,{\rm\Delta}C$ . The solution to (6.7) is then, after some permutation of indices,

(The last term, which ensures symmetry, comes from the first term on the right-hand side in (6.9) after using steady-state balance to replace the  $F$ in $F{\rm\Delta}R$ with $R^{-1}C(R^{-1})^{\text{T}}$ .) When this solution is inserted into (6.5), one finds

where

A consistent (fluctuation–dissipation) approximation (Carnevale & Martin Reference Carnevale and Martin1982) is to replace the two-time dependence of  $C(3,\overline{3})$ with that of $R_{0}(3;\overline{3})$ . Because the coefficients of  ${\rm\Delta}Z$ are all calculated in the homogeneous state, one can Fourier-transform as usual. If one follows the common convention of using  $\boldsymbol{q}$ to refer to the wavevector of the mean field, one then arrives at the dispersion relation

where $R_{1}$ is the first renormalized approximation to $R$ . (In (6.8), replace  ${\it\Sigma}$ by  ${\it\Sigma}_{0}$ and  $R$ by  $R_{1}$ .) The solution determines the eigenvalue  ${\it\omega}$ as a functional of the level  $C$ of the homogeneous turbulence (as well as  $\boldsymbol{q}$ and the various parameters): ${\it\omega}={\it\omega}_{\boldsymbol{q}}[C]$ . Setting ${\it\omega}_{\boldsymbol{q}}[C]=0$ determines the neutral curve $C=C(\boldsymbol{q})$ ; larger values of  $C$ are unstable to the generation of inhomogeneous mean fields (see figure 5). This process is analogous to the general phenomenon of pattern formation, as discussed by Parker & Krommes (Reference Parker and Krommes2013).Footnote ⁹⁵ In the context of the barotropic vorticity equation, the zonostrophic dispersion relation was analysed in detail in a beautiful paper by Srinivasan & Young (Reference Srinivasan and Young2012). Further results have been reported by Parker & Krommes (Reference Parker and Krommes2014, Reference Parker, Krommes, Galperin and Read2015).

Figure 5. Cartoon of a neutral curve $C=C(q)$ for zonostrophic instability. Below the neutral curve, the homogeneous turbulent state is stable. Above it, the homogeneous turbulence is unstable to the generation of inhomogeneous mean fields. Instability first sets in at a critical wavenumber  $q_{c}$ . For given  $C$ in the unstable region, zonal-flow equilibria with a continuum of  $q$ values are allowed. Further analysis must be done to examine the stability of those solutions (Parker & Krommes Reference Parker and Krommes2014).

This zonostrophic instability is a generalization of the often-discussed modulational instability (Dewar & Abdullatif Reference Dewar, Abdullatif, Denier and Frederiksen2006; Connaughton et al. Reference Connaughton, Nadiga, Nazarenko and Quinn2010), in which one considers the instability of a fixed ‘pump’ wave at wavevector  $\boldsymbol{K}$ to the generation of fluctuations at wavevector  $\boldsymbol{Q}$ with sidebands at $\boldsymbol{P}_{\pm }=-\boldsymbol{K}\pm \boldsymbol{Q}$ (see figure 6). The relationship between the zonostrophic and modulational instabilities was first demonstrated by Carnevale & Martin (Reference Carnevale and Martin1982), though their result was not widely appreciated. It has recently been discussed in more detail by Parker & Krommes (Reference Parker, Krommes, Galperin and Read2015), whose work was further generalized by Bakas, Constantinou & Ioannou (Reference Bakas, Constantinou and Ioannou2015, appendix C).

Figure 6. Wavevectors for the modulational instability of a fixed pump wave  $\boldsymbol{K}$ to fluctuations at  $\boldsymbol{Q}$ ; the sidebands are $\boldsymbol{P}_{\pm }\doteq -(\boldsymbol{K}\pm \boldsymbol{Q})$ . The mode coupling obeys $\boldsymbol{K}+\boldsymbol{P}_{\pm }\pm \boldsymbol{Q}=\mathbf{0}$ . This process can be extracted from the general dispersion relation for zonostrophic instability by inserting the background spectrum $C_{\boldsymbol{k}}\propto {\it\delta}(\boldsymbol{k}-\boldsymbol{K})+{\it\delta}(\boldsymbol{k}+\boldsymbol{K})$ .

6.3. The CE2 and S3T closures

The dispersion relation (6.13) also follows from a superficially distinct closure theory, the so-called CE2 (second-order cumulant) closure (see Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011, Srinivasan & Young Reference Srinivasan and Young2012, and references therein). In this approximation, the small-scale turbulence is represented by merely an external white-noise forcing; otherwise, ‘eddy–eddy interactions’ are neglected and the retained interactions are between only the mean fields and the turbulence. It is left to the reader to learn more about this closure (see the discussion below as well as various articles in Galperin & Read (Reference Galperin and Read2015)) and to ponder why the dispersion relations should be equivalent.

In the general notation of (2.6)–(2.8), the CE2 closure omits the terms of $O({\it\delta}{\it\psi}^{2})$ in (2.7) or the third-order cumulant  $T$ in (2.8) (those terms describe the eddy–eddy interactions) while keeping (2.6) for the mean field intact. With white-noise forcing, (2.8) closes in terms of $\langle {\it\psi}\rangle (t)$ , $C(t,t)$ , and the known strength of the forcing; the truncated form of (2.7) ensures that this approximation is realizable. Tobias et al. (Reference Tobias, Dagon and Marston2011), Tobias & Marston (Reference Tobias and Marston2013), and Marston, Qi & Tobias (Reference Marston, Qi, Tobias, Galperin and Read2015) have referred to the numerical solutions of such truncated cumulant systems as ‘direct statistical simulations’.

Realizability is an important property of CE2. While one can define higher-order approximations CE $n$ for $n>2$ , none of those are realizable.Footnote ⁹⁶ Marston (Reference Marston2012) and Marston et al. (Reference Marston, Qi, Tobias, Galperin and Read2015) have explored several phenomenologically modified versions of CE3 that appear to pragmatically restore realizability. The presence of eddy–eddy interactions in those closures leads to results that are superior to those from CE2 in some situations. However, a more systematic approach is desirable.

Studies of CE $n$ closures (by that name) were preceded by the research of Farrell and coworkers on the so-called stochastic structural stability theory (SSST or S3T), the fundamental paper being that of Farrell & Ioannou (Reference Farrell and Ioannou2003); see also Farrell & Ioannou (Reference Farrell, Ioannou, Galperin and Read2015) for an overview. One version of the theory is mathematically identical to CE2, a philosophical difference being the interpretation of the random forcing. In strict CE2, the forcing is viewed as being entirely extrinsic – for example, due to modes that are otherwise unaccounted for. In S3T, the forcing is frequently said to represent the effects of the eddy–eddy interactions as well as any extrinsic noise. That would be in accord with a Langevin representation of the turbulence provided that an energy-conserving coherent damping due to the eddy–eddy terms was also included. Sometimes that has been done, as in the work of DelSole (Reference DelSole2001) and Farrell & Ioannou (Reference Farrell and Ioannou2009b ); other times it has been omitted, as in the application by Farrell & Ioannou (Reference Farrell and Ioannou2009a ) of S3T to the Hasegawa–Wakatani system. Omission of that damping fosters uncertainty about some of the details reported in the latter work (without, however, vitiating the principal message).

Farrell & Ioannou (Reference Farrell, Ioannou, Galperin and Read2015) have reviewed some calculations using S3T and have also embedded that approximation into a more general framework that they call statistical state dynamics. The basic philosophy is that working with an ensemble (even one that is characterized by only the first few cumulants) provides information about cooperative and self-consistent behaviour that is superior to that gained from studies of individual realizations. A technical justification for a simple stochastic parametrization of the eddy–eddy interactions involves the non-normalityFootnote ⁹⁷ of the linear operator in the presence of mean fields; consult Farrell & Ioannou (Reference Farrell, Ioannou, Galperin and Read2015, and references therein) for further discussion of the implications of non-normality.

In the original applications of these types of closures, coupled equations for mean fields and eddy covariance were obtained by applying a zonal average to the original PDE. However, Bakas et al. (Reference Bakas, Constantinou and Ioannou2015) have generalized the procedure to show that one can consider the onset of both zonal and nonzonal perturbations; many detailed results about barotropic beta-plane turbulence are contained in that beautiful paper, which has implications for drift-wave physics as well.

There are many opportunities for further research. For example, some current work focuses on applying CE2-type ideas to the study of dynamo action in accretion disks (Squire & Bhattacharjee Reference Squire and Bhattacharjee2015). Also, the kinds of inhomogeneous states described in this section, which involve steady or quasi-steady zonal jets, do not appear to be relevant for microturbulence in tokamaks; see footnote 90 on page 52 as well as challenge 6 in § 8.

7.1. Some ideas discussed by Kadomtsev

7.2. Some ideas discussed by Dupree and Weinstock

7.3. Some miscellaneous topics related to analytical methods

The following subsections mention representative topics that have emerged subsequent to the early works by Kadomtsev, Dupree, and Weinstock.