text:
# Non-white noise and a multiple-rate Markovian closure theory for turbulence

## I Introduction

Our derivation builds on and closely follows the work by Bowman, Krommes, and Ottaviani (we will frequently refer to this paper as BKO), on realizable Markovian closures derived from Kraichnan’s direct-interaction-approximation (the DIA). The DIA is based on a renormalized perturbation theory and gives an integro-differential set of equations to determine the two-time correlation function. The DIA involves time integrals over the past history of the system, which can be computationally expensive. Markovian approximations give a simpler set of differential equations that involve only information from the present time. They approximate two-time information in the correlation function and in the renormalized Green’s function by a decorrelation rate parameter. The structure of the equations we derive here is similar to the realizable Markovian closure (RMC) of BKO, but with extensions such as replacing a single decorrelation rate parameter with several different nonlinear rate parameters, to allow for a more accurate model of effects such as non-white noise. (As will be discussed more below, the RMC does include more non-white-noise effects than one might think at first.)

The basic issue studied in the present paper can be illustrated by a simple Langevin equation (which will be discussed in more detail in the next section):

$$\left(\frac{}{t}+\eta \right)\psi (t)=f(t),$$
(1)
where $`\eta `$ is the decay rate and $`f`$ is a random forcing or stirring term (also known as noise). As is well known, if $`f`$ is white noise, then the decorrelation rate for $`\psi `$ is given by $`\eta `$, so that in a statistical steady state the two-time correlation function $`\psi (t)\psi ^{}(t^{})=C_0\mathrm{exp}(\eta |tt^{}|)`$ (assuming constant real $`\eta `$ here). However, if $`f(t)`$ varies slowly compared to the 1/$`\eta `$ time scale, then the solution to the Langevin equation is just $`\psi (t)f(t)/\eta `$, and the decorrelation rate for $`\psi `$ is instead given by the decorrelation rate for $`f`$. Note that the Green’s function (the response to a perturbation at time $`t^{}`$) is still $`\mathrm{exp}(\eta (tt^{}))`$. Previous Markovian closures employed some variant of an ansatz, based on the fluctuation–dissipation theorem, that the two-time correlation function and the Green’s function were proportional to each other. This is a rigorous result for a system in thermal equilibrium, but may not necessarily apply to a turbulent system. The purpose of the present paper is to explore an extended Markovian closure, which we will call the Multiple-Rate Markovian Closure (MRMC), that allows the decorrelation rate of $`\psi `$ to differ from the decay rate $`\eta `$.

In practice, the corrections due to non-white-noise effects may be quantitatively modest, as the decorrelation rate for the turbulent noise $`f`$ that is driving $`\psi `$ at a particular wave number $`𝒌`$ is often comparable to or greater than the nonlinear damping rate $`\eta `$ at that $`𝒌`$. This is because the turbulent noise driving mode $`𝒌`$ arises from the nonlinear beating of other modes $`𝒑`$ and $`𝒒`$ such that $`𝒌=𝒑+𝒒`$. Thus $`|𝒑|`$ or $`|𝒒|`$ has to be comparable to or larger than $`|𝒌|`$, and will thus have comparable or larger decorrelation rates, since the decay rate $`\eta `$ is usually an increasing function of $`|𝒌|`$. Furthermore, there are some offsetting effects due to the time-history integrals in the DIA’s generalized Langevin equation that might further reduce the difference between the decay rate and the decorrelation rate. Indeed, past comparisons of the RMC with the full DIA or with the full nonlinear dynamics have generally found fairly good agreement in many cases, including two-field Hasegawa–Wakatani drift-wave turbulence and galactic dynamo MHD turbulence. Some of the results in this paper help to give a deeper insight into why this agreement is often fairly good, despite the arguments of the previous paragraph, i.e., why the fluctuation–dissipation ansatz is often a reasonable approximation even out of thermal equilibrium. But there may be some regimes where the differences are important and the improvements suggested here would be welcome. These might include include plasma cases where the wave dynamics can make $`\eta `$ vary strongly with the direction of $`𝒌`$ in some cases (with strong Landau damping in some directions and strong instabilities in other directions, for example), or non-steady-state cases involving zonal flows exhibiting predator-prey dynamics.

Markovian closures such as the test-field model (TFM) or Orszag’s eddy-damped quasinormal Markovian (EDQNM) closure have been extensively used to study turbulence in incompressible fluids and plasmas. The introduction of BKO provides useful discussions of the background of the DIA and Markovian closures, and we will add just a few remarks here (there are also many reviews on these topics, such as Refs. ). The RMC developed in BKO is similar to the EDQNM, but has features that ensure realizability even in the presence of the linear wave phenomena exhibited by plasmas (e.g. drift waves) and rotating planetary flows (e.g. Rossby waves). “Realizability” is a property of a statistical closure approximation that ensures that, even though it is only an approximate solution of the original equations, it is an exact solution to some other underlying stochastic equation, such as a Langevin equation. The absence of realizability can cause serious physical and numerical problems, such as the prediction of negative or even divergent energies. The RMC reduces to the DIA-based version of the EDQNM in a statistical steady state, so in some cases the issue of realizability is only important in the transient phase as a steady state is approached or in freely decaying turbulence. Realizability may also be important in certain cases of recent interest among fusion researchers where oscillations may occur between various parts of the spectrum (such as predator–prey type oscillations between drift waves and zonal flows) where a simple statistical steady state might not exist, or where one is interested in the transient dynamics. Unlike some Markovian models that differ only in the transient dynamics, the Multiple-Rate Markovian Closure presented here could also alter the steady-state spectrum.

Our results apply to a Markovian approximation of the DIA for a generic one-field system with a quadratic nonlinearity. They are immediately applicable to some simple drift-wave plasma turbulence problems, Rossby-wave problems, or two-dimensional hydrodynamics. Future work could extend this approach to multiple fields, similar to the covariant multifield RMC of BKO or their later realizable test-field model. Multiple field equations can get computationally difficult (with the compute time scaling as $`n^6`$, where $`n`$ is the number of fields), though two-field studies have been done and disparate scale approximations or other approximations might make them more tractable. In addition to their direct use in studying turbulence in a wide range of systems, the Markovian closures discussed here might also be useful in deriving subgrid turbulence models for computer simulations.

While our formulation is general and potentially applicable to a wide range of nonlinear problems involving Markovian approximations, we were motivated by some recent problems of interest in plasma physics and fusion energy research, such as zonal flows. Initial analytic work elucidating the essentials of nonlinear zonal flow generation used weak-turbulence approximations or secondary-instability analysis. Recent interesting work by Krommes and Kim uses a Markovian statistical theory to extend the study of zonal flows to the strong turbulence regime. An important question is why the strong generation of zonal flows seen near marginal stability is not as important in stronger instability regimes (i.e., why is the Dimits nonlinear shift finite?). A strong turbulence theory is needed to study this. An alternative approach, which has been fruitful in providing the main answers to the finite Dimits shift question, is to analyze the secondary and tertiary instabilities involved in the generation and breakup of zonal flows. That work suggests that a complete strong-turbulence Markovian model of this problem would also need multi-field and geometrical effects (involving at least the potential and temperature fields, along with certain neoclassical effects in toroidal magnetic field geometry).

Based on the reasoning immediately following Eq. (1) above, one might think that the assumption that the two-time correlation function and the Green’s function are proportional to each other is rigorous only in the limit of white noise (which has an infinite decorrelation rate). The Realizable Markovian Closure has been shown to correspond exactly to a simple Langevin equation (where the effects of the turbulence appear in nonlinear damping and nonlinear noise terms), for which this might appear to be the implication. However, the mapping from statistically averaged equations (such as Markovian closures) back to a stochastic equation for which it is the solution, is not necessarily unique. In particular, the full DIA corresponds to a generalized Langevin equation (Eq. (49) below), in which the damping term $`\eta \psi (t)`$ in the simple Langevin equation is replaced by a time-history integral operator. As we will find, it is then possible for the two-time correlation function and the Green’s function to be proportional to each other even when the noise has a finite correlation time. This allows the fluctuation–dissipation theorem (which is rigorous in thermal equilibrium) to be satisfied without requiring the noise to be white (since the noise is not necessarily white in thermal equilibrium). Thus the fluctuation–dissipation ansatz of BKO is a less restrictive assumption than one might have at first thought. \[It should be noted that previous Markovian models account implicitly for at least some non-white-noise effects. For example, in the calculation of the triad interaction time $`\theta _{𝒌𝒑𝒒}=1/(\eta _𝒌+\eta _𝒑+\eta _𝒒)`$ for three-wave interactions, finite values of the assumed noise decorrelation rate $`\eta _𝒑+\eta _𝒒`$ are used.\]

Nevertheless, there is still no reason in a general situation that the two-time correlation function and the Green’s function be constrained to be proportional to each other. As described elsewhere, there may be regimes where the resulting differences between the decorrelation rate and the decay rate are significant.

The outline of this paper is as follows. Sec. (II) presents some of the essential ideas of this paper for a very simple Langevin equation, and includes a section motivating the choice of the limit operator introduced by BKO to ensure realizability. Sec. (III) presents a more detailed calculation of the non-white Markovian model for a simple Langevin equation, including the effects of complex damping rates (to represent the wave frequency) and the issue of Galilean invariance. The model is compared with exact results in the steady-state limit, and then an extension to time-dependent Langevin statistics is presented (along with, in Appendix A, an alternative proof of realizability for this case). Sec. (IV) presents the notation of the full many-mode nonlinear equations we will solve and summarizes the direct interaction approximation (DIA), which is our starting point. Sec. (V) summarizes how the non-white Markovian approach is derived for the steady-state limit (with further details given in Appendix B), while Sec. (VI) presents the full non-white Markovian approximation for the time-dependent DIA. Sec. (VII) discusses some important properties of these equations, including the limits of thermal equilibrium and inertial range scaling, and some difficulties due to the lack of random Galilean invariance in the DIA and described in Appendix C. The conclusions include some suggestions for future research.

## II Simple examples based on the Langevin equation

Here we expand upon the analogy given in Sec. (I) using a simple Langevin equation, which provides a useful paradigm for understanding the essential ideas considered in this paper. Since realizable Markovian closure approximations to the DIA can be shown to correspond exactly to an underlying set of coupled Langevin equations, the analogy is quite relevant. In this section we will consider heuristic arguments based on some simple scalings; later sections will be more rigorous.

Consider the simple Langevin equation

$$\left(\frac{}{t}+\eta (t)\right)\psi =f^{}(t),$$
(2)
where $`\eta `$ is a damping rate and $`f^{}`$ is a random forcing or stirring term (also known as “noise”). \[Here we now force Eq. (1) with the complex conjugate of $`f`$, for consistency with the form of the equations used later for a generic quadratically nonlinear equation.\] The statistics of the noise are given by a specified two-time correlation function $`C_f(t,t^{})=f(t)f^{}(t^{})`$. \[In the white-noise limit, $`C_f(t,t^{})=2D\delta (tt^{})`$, and the power spectrum of the Fourier-transform of $`f(t)`$ is independent of frequency, and is thus called a “white” spectrum.\] The Langevin equation is used to model many kinds of systems exhibiting random-walk or Brownian motion features. Here we can think of $`\psi `$ as the complex amplitude of one component of the turbulence with a specified Fourier wave number $`𝒌`$. Note that $`\eta `$ may be complex (representing both damping and wave-like motions) and represents both linear and nonlinear (renormalized) damping or frequency shifts due to interactions with other modes. The random forcing $`f^{}`$ represents nonlinear driving by other modes beating together to drive this mode.

The response function (or Green’s function or propagator) for this equation satisfies

$$\left(\frac{}{t}+\eta \right)R(t,t^{})=\delta (tt^{}),$$
(3)
which easily yields $`R(t,t^{})=\mathrm{exp}(\eta (tt^{}))H(tt^{})`$ if $`\eta `$ is independent of time, where $`H(t)`$ is the Heaviside step function. The solution to the Langevin equation is just $`\psi (t)=_0^t𝑑\overline{t}R(t,\overline{t})f^{}(\overline{t})`$ (for the initial condition $`\psi (0)=0`$). It is then straightforward to demonstrate the standard result that, if $`f`$ is white noise and the long-time statistical steady-state limit is considered, then the correlation function for $`\psi `$ is

$`C(t,t^{})\psi (t)\psi ^{}(t^{})=C_0\mathrm{exp}(\eta (tt^{}))`$

for $`t>t^{}`$, where $`C_0=2D/(\eta +\eta ^{})`$ (we emphasize definitions with the notation $``$). \[For $`t<t^{}`$, one can use the symmetry condition $`C(t,t^{})=C^{}(t^{},t)`$.\] That is, the decorrelation rate for $`\psi `$ is just $`\eta `$. This is equivalent to the assumption in a broad class of Markovian models that the decorrelation rate for $`\psi `$ is the same as the decay rate of the response function.

However, consider the opposite of the white-noise limit, where $`f`$ varies slowly in time compared to the $`1/\eta `$ time scale. Then the solution of Eq. (2) is approximately $`\psi (t)=f^{}(t)/\eta `$, and the decorrelation rate for $`\psi `$ will be the same as the decorrelation rate for $`f^{}`$. In this limit, the assumption in many Markovian models that the decorrelation rate is $`\eta `$ is not valid.

Denote the decorrelation rate for $`f^{}`$ as $`\eta _f^{}`$, and the decorrelation rate for $`\psi `$ as $`\eta _C`$. Then one might guess that a simple Padé-type formula that roughly interpolates between the white-noise limit $`\eta _f\eta `$ and the opposite “red-noise” limit $`\eta _f\eta `$ would be something like $`1/\eta _C1/\eta +1/\eta _f^{}`$, or

$$\eta _C=\frac{\eta \eta _f^{}}{\eta +\eta _f^{}}.$$
(4)
In fact, we will discover in the next section that more detailed calculations give similar results in the limit of real $`\eta `$ and $`\eta _f`$, though the formulas are more complicated in the presence of wave behavior with complex $`\eta `$ and $`\eta _f`$.

We note that in many cases of interest, the noise decorrelation rate $`\eta _f`$ turns out to be of comparable magnitude to $`\eta `$ (for example, if the dominant interactions involve modes of comparable scale). In this case, while the white-noise approximation is not rigorously valid, the corrections to the decorrelation rate considered in this paper might turn out to be quantitatively modest, $`50\%`$. Furthermore, in the case of the full DIA and its corresponding generalized Langevin equation, we will find additional corrections that can, in some cases, offset the effects in Eq. (4) and cause $`\eta _C`$ to be closer to $`\eta `$.

Before going on to the more detailed results in the next section, we consider the meaning of an operator introduced in the BKO derivations in order to preserve realizability in the time-dependent case, where $`\eta (t)`$ varies in time and may be negative (transiently), representing an instability. \[In order for a meaningful long-time steady-state limit to exist, the net $`\eta `$ (which is the sum of linear and nonlinear terms) must eventually go positive to provide a sink for the noise term. But it is important to preserve realizability during the transient times when $`\eta `$ may be negative.\] Based on arguments about symmetry and the steady-state fluctuation–dissipation theorem, they initially proposed a time-dependent ansatz of the form

$$C(t,t^{})=C^{1/2}(t)C^{1/2}(t^{})\mathrm{exp}\left(_t^{}^t𝑑\overline{t}\eta (\overline{t})\right)$$
(5)
(for $`t>t^{}`$), where $`C(t)C(t,t)`$ is the equal-time covariance. Later in their derivation, they state that in order to ensure realizability, $`\eta (\overline{t})`$ in this expression had to be replaced with $`𝒫(\eta (\overline{t}))`$, where the operator $`𝒫(\eta )=Re\eta H(Re\eta )+iIm\eta `$ prevents the real part of the effective $`\eta `$ in Eq. (5) from going negative.

Physically this makes sense for the following reasons. Consider Eq. (2) with white noise $`f`$ (thus ignoring the non-white-noise effects). Then in a normal statistical steady state where $`\eta (t)`$ is constant and $`Re\eta >0`$, Eq. (5) properly reproduces the usual result $`C(t,t^{})=C_0\mathrm{exp}(\eta |tt^{}|)`$. However, if $`Re\eta <0`$ (which it might do at least transiently in the full turbulent system considered later), then Eq. (2) can’t reach a steady state, and the solution is eventually $`\psi (t)=\psi _0\mathrm{exp}(\eta t)=|\psi (t)|\mathrm{exp}(iIm\eta t)`$, after an initial transient phase. Thus $`C(t,t^{})=C^{1/2}(t)C^{1/2}(t^{})\mathrm{exp}(iIm\eta (tt^{}))`$, in agreement with and providing an additional intuitive argument for BKO’s modified form of Eq. (5), including the $`𝒫(\eta )`$ operator. \[There may be an initial phase where the noise term $`f`$ in Eq. (2) dominates and causes $`C(t)`$ to grow linearly in time, $`C(t)=\psi (t)\psi ^{}(t)=2Dt`$. But eventually the unstable $`\eta \psi `$ term will become large enough to dominate and lead to exponential growth of $`\psi `$.\]

The model we will introduce below replaces $`\eta `$ in Eq. (5) with a separate parameter $`\eta _C`$, and develops a formula to relate $`\eta _C`$ to other parameters in the problem such as $`\eta `$ and $`\eta _f`$. In the white-noise limit, the formula for $`\eta _C`$ automatically reproduces the effects of the $`𝒫`$ limiting operator, as will be described in the next section and in Appendix (A). But numerical investigation of non-white noise with wave dynamics ($`Im\eta 0`$ or $`Im\eta _f0`$) uncovered cases where the $`𝒫`$ limiting operator is still needed to ensure realizability. This will be explained at the end of Sec. (III C).

We considered naming the method described in this paper the Non-White Markovian Closure since, for the simple Langevin equation considered here and in the next section, the decorrelation rate and the decay rate are equal only in the white-noise limit, and this approach allows these rates to differ. \[Alternatively, to emphasize the flexibility of this method one might have called it the Colored-Noise Markovian Closure since instead of being restricted to white-noise (a uniform spectrum) we can allow a noise spectrum of width $`\delta \omega Re\eta _f`$ peaked near an arbitrary frequency $`\omega Im\eta _f`$. In other words, this closure can model spectra with a range of possible colors.\] However, as we will discuss further, while a simple Langevin equation is sometimes used to demonstrate realizability of Markovian approximations, the DIA is actually based on a generalized Langevin equation involving a non-local time-history integral (compare Eq. (2) with Eq. (49)). Because non-white fluctuations enter not only by making the noise term non-white but also by affecting this time-history integral, it is possible for the decay-rate and the decorrelation rate to be equal even in some cases where the noise is not white (as indeed is the case in thermal equilibrium where the fluctuation–dissipation theorem must hold but the noise is not necessarily white). We thus favor the name Multiple-Rate Markovian Closure (MRMC), to emphasize that the method developed here is a generalization of the previous Realizable Markovian Closure (RMC) to allow for multiple rates (i.e., separate decay and decorrelation rates).

## III Detailed demonstration of the Multiple-Rate Markovian method with the Langevin equation

In this section, we demonstrate the Multiple-Rate Markovian approach starting with a simple Langevin equation. The steps in the derivation are quite similar to the steps that will be taken in the following sections for the case of the more complete DIA for more complicated nonlinear problems, and thus help build insight and familiarity. In this section, we will be introducing various approximations that may seem unnecessary for the simple Langevin problem, which can be solved exactly in many cases (for simple forms of the noise correlation function). But these are the same approximations that will be used later in deriving Markovian approximations to the DIA, and so it is useful to be able to test their accuracy in the Langevin case.

Our starting point is the Langevin Eq. (2), but we allow $`\eta (t)`$ to be a function of time, so that the solution to Eq. (3) for the response function is

$$R(t,t^{})=\mathrm{exp}\left(_t^{}^t𝑑\overline{t}\eta (\overline{t})\right)H(tt^{})$$
(6)
(instead of the solution given immediately after Eq. (3), which assumes that $`\eta `$ is independent of time). The solution to the Langevin equation is

$$\psi (t)=R(t,0)\psi (0)+_0^t𝑑\overline{t}R(t,\overline{t})f^{}(\overline{t}).$$
(7)
In principle it is possible to calculate directly two-time statistics like $`C(t,t^{})=\psi (t)\psi ^{}(t^{})`$ from this, but in practice it is often convenient to consider instead the differential equation for $`C(t,t^{})/t`$, which from Eq. (2) and Eq. (7) is

$`\left({\displaystyle \frac{}{t}}+\eta \right)C(t,t^{})`$ $`=`$ $`f^{}(t)\psi ^{}(t^{})`$ (8)
$`=`$ $`{\displaystyle _0^t^{}}𝑑\overline{t}R^{}(t^{},\overline{t})C_f^{}(t,\overline{t}),`$ (9)
where the noise correlation function is defined as $`C_f(t,t^{})=(f(t)f^{}(t^{})`$, and we have assumed that the initial condition $`\psi (0)`$ has a random phase. This equation is the analog of the DIA equations for the two-time correlation function (compare with Eq. (46) and Eqs. (IV B)), but with an integral only over the noise and no nonlinear modification of the damping term.

We define the equal-time correlation function $`C(t)`$ in terms of the two-time correlation function $`C(t,t^{})`$ as $`C(t)=C(t,t)=\psi (t)\psi ^{}(t^{})`$ (note that these two functions are distinguished only by the number of arguments). Then

$$\frac{C(t)}{t}+2Re\eta C(t)=2Re_0^t𝑑\overline{t}R^{}(t,\overline{t})C_f^{}(t,\overline{t}).$$
(10)
This is the analog of the DIA equal-time covariance equation, Eq. (IV B).

### A Langevin statistics in the steady-state limit

Consider the steady-state limit where $`t,t^{}\mathrm{}`$ (but with finite time separation $`tt^{}`$), and assume the noise correlation function has the simple form $`C_f(t,t^{})=C_{f0}\mathrm{exp}[\eta _f(tt^{})]`$ for $`t>t^{}`$. In this section we assume $`\eta `$ and $`\eta _f`$ are time-independent constants. The response function reduces back to its steady-state form $`R(t,t^{})=\mathrm{exp}[\eta (tt^{})]H(tt^{})`$. Then Eq. (10) in steady state gives

$$C_0\underset{t\mathrm{}}{lim}C(t)=C_{f0}\frac{Re(\eta +\eta _f)}{Re(\eta )(\eta +\eta _f)(\eta ^{}+\eta _f^{})}.$$
(11)
Writing $`\eta =\nu +i\omega `$ and $`\eta _f=\nu _f+i\omega _f`$ in terms of their real and imaginary components, and denoting the frequency mismatch $`\mathrm{\Delta }\omega =\omega +\omega _f`$ (remember, because the complex conjugate $`f^{}`$ is used as the forcing term, resonance occurs when $`Im(\eta )=Im(\eta _f^{})`$) this can be written as

$$C_0=\frac{C_{f0}}{\nu }\frac{(\nu +\nu _f)}{(\nu +\nu _f)^2+(\mathrm{\Delta }\omega )^2}.$$
(12)
This has a familiar Lorentzian form characteristic of resonances.

To find the two-time correlation function, the time integral in Eq. (9) can be evaluated for $`t>t^{}`$ to give

$$\left(\frac{}{t}+\eta \right)C(t,t^{})=\frac{C_{f0}}{\eta ^{}+\eta _f^{}}\mathrm{exp}[\eta _f^{}(tt^{})].$$
(13)
With the steady-state boundary condition $`C(t=t^{},t^{})=C_0`$, this can be solved to give

$`C(t,t^{})`$ $`=`$ $`C_0\left[1{\displaystyle \frac{Re(\eta )(\eta +\eta _f)}{Re(\eta +\eta _f)(\eta \eta _f^{})}}\right]\mathrm{exp}[\eta (tt^{})]`$ (14)
$`+`$ $`C_0{\displaystyle \frac{Re(\eta )(\eta +\eta _f)}{Re(\eta +\eta _f)(\eta \eta _f^{})}}\mathrm{exp}[\eta _f^{}(tt^{})].`$ (15)
In the white-noise limit, $`|\eta _f||\eta |`$, this reduces to the standard simple result $`C(t,t^{})=C_0\mathrm{exp}[\eta (tt^{})]`$. But in the more general case of non-white noise, the two-time correlation function is more complicated. \[Despite the apparent singularity in the denominator, it is cancelled by the exponentials so that $`C(t,t^{})`$ is well-behaved in the limit $`\eta \eta _f^{}`$.\] In the context of the turbulent interaction of many modes, $`C_f(t,t^{})`$ and thus $`C(t,t^{})`$ may be very complicated functions. Even if the noise correlation function has a simple exponential dependence $`C_f(t,t^{})\mathrm{exp}[\eta _f(tt^{})]`$, we see that the resulting correlation function for $`\psi `$ is more complicated.

Consider the task of fitting this complicated $`C(t,t^{})`$ with a simpler model of the form

$$C_{\mathrm{mod}}(t,t^{})=C_0\mathrm{exp}[\eta _C(tt^{})]$$
(16)
(for $`t>t^{}`$). One way to define the effective decorrelation rate $`\eta _C`$ might be based on the area under the time integral,

$$_{\mathrm{}^{}}^t𝑑t^{}C_{\mathrm{mod}}(t,t^{})=\frac{C_0}{\eta _C}=_{\mathrm{}^{}}^t𝑑t^{}C(t,t^{}).$$
(17)
This can be evaluated either by directly substituting Eq. (15), or by taking a time average of Eq. (13); the same answer results either way. It turns out that in the later versions of this calculation it is easier to determine $`\eta _C`$ by integrating the governing differential equation over time. Operating on Eq. (13) with $`_{\mathrm{}}^t𝑑t^{}`$ and using

$$_{\mathrm{}}^t𝑑t^{}\frac{C(t,t^{})}{t}=\frac{}{t}_{\mathrm{}}^t𝑑t^{}C(t,t^{})C(t,t),$$
(18)
we find

$$\frac{1}{\eta _C}=\frac{1}{\eta }+\frac{Re(\eta )(\eta +\eta _f)}{Re(\eta +\eta _f)\eta \eta _f^{}}.$$
(19)
This recovers the white-noise limit $`\eta _f\eta `$ and the red-noise limit $`\eta _f\eta `$ discussed in Sec. II. In the limit of real $`\eta `$ and real $`\eta _f`$ it simplifies to the Padé approximation $`\eta _C=\eta \eta _f/(\eta +\eta _f)`$ also suggested in the introduction. However, there is a problem with Eq. (19) related to Galilean invariance. Suppose we make the substitutions $`\psi =\widehat{\psi }\mathrm{exp}[i\omega _2t]`$ and $`f^{}=\widehat{f}^{}\mathrm{exp}[i\omega _2t]`$ into the Langevin Eq. (2). Then it can be written as

$$\left(\frac{}{t}+\widehat{\eta }\right)\widehat{\psi }=\widehat{f}^{}(t),$$
(20)
where $`\widehat{\eta }=\eta +i\omega _2`$, and the results should be the same if written in terms of the transformed variables. In particular, the correlation function should transform as $`\widehat{\psi }(t)\widehat{\psi }^{}(t^{})=\mathrm{exp}[i\omega _2(tt^{})]\psi (t)\psi ^{}(t^{})=\mathrm{exp}[i\omega _2(tt^{})]C(t,t^{})`$. Thus the decorrelation rate $`\widehat{\eta }_C`$ for $`\widehat{\psi }`$ should be related to the decorrelation rate $`\eta _C`$ for $`\psi `$ by $`\widehat{\eta }_C=\eta _C+i\omega _2`$. The decorrelation rate for the transformed noise term $`\widehat{f}^{}`$ also transforms as $`\widehat{\eta }_f^{}=\eta _f^{}+i\omega _2`$. In the case of fluid or plasma turbulence where $`\psi `$ represents the amplitude of a Fourier mode $`\mathrm{exp}[i𝒌𝒙]`$ and $`f^{}`$ represents the amplitude of two modes with wave numbers $`𝒑`$ and $`𝒒`$ beating together to drive the $`𝒌`$ mode (so $`𝒑+𝒒=𝒌`$), these transformations correspond to a Galilean transformation to a moving frame $`𝒙=𝒙_0+𝒗t`$, with $`\omega _2=𝒌𝒗`$.

So all results should be independent of $`\omega _2`$ under the transformation $`\eta =\widehat{\eta }i\omega _2`$, $`\eta _f^{}=\widehat{\eta }_f^{}i\omega _2`$, (thus $`\eta _f=\widehat{\eta }_f+i\omega _2`$), $`\eta _C=\widehat{\eta }_Ci\omega _2`$. Eq. (12) satisfies this, but Eq. (19) fails this test. This problem and its solution is described in the review paper by Krommes, who shows it is related to other differences in various previous Markovian closures. The problem can be traced to the definition of Eq. (17), which doesn’t satisfy the invariance for general forms of $`C(t,t^{})`$. For example, we could have multiplied the integrand in Eq. (17) by an arbitrary weight function (such as $`\mathrm{exp}[i\omega _2(tt^{})]`$) before taking the time average, and the results would have changed. The way to fix this problem is to do the time average in a natural frame of reference for $`\psi `$ that accounts for its frequency dependence. This leads us to the definition:

$$\frac{C_0^2}{\eta _C+\eta _C^{}}_{\mathrm{}}^t𝑑t^{}C_{\mathrm{mod}}^{}(t,t^{})C(t,t^{}).$$
(21)
This corresponds to fitting $`C_{\mathrm{mod}}(t,t^{})`$ to $`C(t,t^{})`$ by requiring that both effectively have the same projection onto the function $`C_{\mathrm{mod}}(t,t^{})`$. \[As Krommes points out, using the invariant definition Eq. (21) instead of Eq. (17) is a non-trivial point needed to ensure realizability and avoid spurious nonphysical solutions in some cases.\]

Operating on Eq. (13) with $`_{\mathrm{}}^t𝑑t^{}C_{\mathrm{mod}}^{}(t,t^{})`$, using a generalization of Eq. (18), and doing a little rearranging yields

$$\eta _C=\eta \frac{C_{f0}(\eta _C+\eta _C^{})}{C_0(\eta ^{}+\eta _f^{})(\eta _C^{}+\eta _f^{})}.$$
(22)
This is properly invariant to the transformation described in the previous paragraph. Solving for $`\eta _C`$ while leaving $`\eta _C^{}`$ on the other side of the equation, eventually leads to

$$\eta _C=\frac{\eta \eta _f^{}Re(\eta +\eta _f)+i\eta _C^{}Im(\eta \eta _f^{})}{(\eta +\eta _C^{})Re(\eta +\eta _f)+(\eta _f^{}+\eta ^{})Re(\eta _f)}.$$
(23)
If we consider the limit where $`\eta `$, $`\eta _f`$, and thus $`\eta _C`$ are all real, this simplifies to the form

$$\eta _C=\frac{\eta \eta _f}{\eta +\eta _f+\eta _C}.$$
(24)
This is similar to (but more accurate than) the rough interpolation formula Eq. (4) suggested in the introduction. This kind of recursive definition, with $`\eta _C`$ appearing on both sides, is a common feature of the steady-state limit of theories based on the renormalized DIA equations, and can be solved in practice by iteration, or by considering the time-dependent versions of the theories. In Eq. (24) with real coefficients, one can easily solve this equation for $`\eta _C`$, but the solution is much more difficult in the case of complex coefficients in Eq. (23). The resulting calculation is laborious, so we used the symbolic algebra package Maple to solve for $`\eta _C`$ with complex coefficients. Looking at the real and imaginary parts of Eq. (23) separately eventually leads to a quadratic equation and a linear equation to determine the real and imaginary parts of $`\eta _C`$. Unfortunately it takes 16 lines of code to write down the resulting closed-form solution (though perhaps there are common subexpressions that would simplify it). (Maple worksheets that show this calculation and check other main results in this paper are available online.) This is tedious for humans but easy to evaluate in Fortran, C, or other computer language. On the other hand, this is only helpful for the simple Langevin problem anyway since direct solution is not really practical for the full nonlinear problem considered by the DIA, where the noise term of the Langevin equation is replaced by a sum over many modes. In many cases of interest, the noise decorrelation rate $`\eta _f`$ turns out to be comparable in magnitude to $`\eta `$, so iteration of Eq. (23) usually converges quickly. (However, there are limits where convergence is slow, such as some strongly non-resonant cases where $`Re\eta `$ is very close to $`Re\eta _f`$ and both are very small compared to $`Im(\eta _f+\eta )`$.) The other option is to consider the time-dependent problem, the topic of the subsection after next, which effectively performs an iteration in time as a steady state is approached.

### B Comparison of the Multiple-Rate model with exact Langevin result

Figs. (1-6) provide a comparison of the exact and model results for various parameters. The exact Langevin solution for $`C(t,t^{})/C_0`$ is given by Eq. (15). The curves labeled “model” are for the Multiple-Rate Markovian model $`C_{\mathrm{mod}}(t,t^{})/C_0=\mathrm{exp}(\eta _C|tt^{}|)`$, where $`\eta _C`$ is obtained by solving Eq. (23). The curves labeled “wn” are the results for a simple white-noise assumption $`C(t,t^{})/C_0=\mathrm{exp}(\eta |tt^{}|)`$. The plots show both the real and imaginary parts of $`C(t,t^{})`$, except when $`C(t,^{}t)`$ is purely real.

The results are shown in Figs. (1-6) for a variety of parameters. We choose $`\eta =1`$ as a standard normalization in all cases. Only the frequency mismatch ($`\mathrm{\Delta }\omega =Im(\eta \eta _f^{})`$) between the oscillator and the random driving term and the relative decorrelation rate ($`Re(\eta _f)/Re(\eta )`$) can matter. Thus we choose a frame of reference where $`Im\eta =0`$ and any frequency mismatch is reflected in the value of the noise frequency $`Im\eta _f`$.

These comparisons show that the non-white-noise Multiple-Rate model for $`\eta _C`$ does fairly well in most cases. All formulas of course agree well in the white-noise limit of $`Re\eta _fRe\eta `$. The errors of the white-noise model are particularly large in the “red-noise limit” $`Re\eta _fRe\eta `$, though they are noticeable even if $`Re\eta _fRe\eta `$. The white-noise model has a purely real correlation function in all cases, thus missing the frequency shifts that arise when $`Im\eta _f0`$, while the multiple-rate model does a fairly good job of capturing the real and imaginary parts of $`C(t,t^{})`$ in most cases. The most challenging case for even the multiple-rate model is depicted in Fig. (5), where there is a large frequency mismatch but comparable decorrelation rates, $`Re\eta _fRe\eta `$. However, Eq. (12) shows that the amplitude, $`C_02C_{f0}/(\mathrm{\Delta }\omega )^22C_{f0}/(Im(\eta \eta _f^{}))^2`$, will be small in this strongly non-resonant case, and perhaps does not matter much compared to resonant interactions in realistic many-mode turbulence cases. Strongly non-resonant cases are easier to model with disparate values of $`Re\eta _f`$ and $`Re\eta `$, as shown in Fig. (6) and Fig. (2), because interference effects are less important. To do better for the non-resonant case with $`Re\eta _fRe\eta `$ would probably require a more elaborate two-exponential model than Eq. (16), to allow for the constructive and destructive interference effects represented in Fig. (5). Of course, for the simple Langevin case of this section, such a model could exactly reproduce Eq. (15), although for more complicated cases it would again become a model to be fit to the true $`C(t,t^{})`$ dynamics. (Another approach, which might improve the long-time fit a bit, might be to use $`C_{\mathrm{mod}}^{}(t,t^{})(tt^{})`$ as the weight function in Eq. (21) instead of just $`C_{\mathrm{mod}}^{}(t,t^{})`$.)

### C Time-dependent Langevin statistics

We now return our attention to the more general Langevin problem with time-dependent $`\eta (t)`$ and time-varying statistics for the noise term $`f^{}(t)`$. That is, for generality, we also allow the noise amplitude (given by the equal-time covariance $`C_f(t)C_f(t,t)`$) and the noise decorrelation rate to vary in time. Our choice of a self-consistent model for $`C_f(t,t^{})`$ to accomplish this is motivated by BKO’s demonstration that the following form is a realizable correlation function:

$$C_f(t,\overline{t})=C_f^{1/2}(t)\mathrm{exp}\left[_{\overline{t}}^t𝑑t^{\prime \prime }\eta _f(t^{\prime \prime })\right]C_f^{1/2}(\overline{t})$$
(25)
(for $`t\overline{t}`$). \[BKO show this is realizable as long as $`Re(\eta _f(t))0`$ almost everywhere.\] Using this expression, Eq. (10) can be written as

$$\frac{C(t)}{t}+2Re\eta (t)C(t)=2ReC_f^{1/2}(t)\mathrm{\Theta }^{}(t),$$
(26)
where

$$\mathrm{\Theta }(t)_0^t𝑑\overline{t}R(t,\overline{t})\mathrm{exp}\left[_{\overline{t}}^t𝑑t^{\prime \prime }\eta _f(t^{\prime \prime })\right]C_f^{1/2}(\overline{t}).$$
(27)
Taking the time derivative of this expression, and using Eq. (6), leads to

$$\frac{\mathrm{\Theta }(t)}{t}=[\eta (t)+\eta _f(t)]\mathrm{\Theta }(t)+C_f^{1/2}(t),$$
(28)
which is more convenient to use in a time-dependent calculation than Eq. (27). The initial condition is $`\mathrm{\Theta }(0)=0`$. Eq. (26) and Eq. (28) can be used to determine the equal-time covariance $`C(t)`$, but how can we determine the decorrelation rate $`\eta _C`$ from the two-time correlation function $`C(t,t^{})`$? \[In the full nonlinear equations used for the DIA, $`\psi `$ for one mode appears in noise terms for other modes, and so we would like to know the decorrelation rate as well as the amplitude $`C(t)`$.\] Even in the steady-state limit of the previous section, we found that the full two-time correlation function $`C(t,t^{})`$ had a more complicated form than a simple exponential, and so we fit a simpler model $`C_{\mathrm{mod}}(t,t^{})`$ to it in order to determine an effective decorrelation rate $`\eta _C`$.

We follow a similar procedure here. We again use BKO’s form for a realizable time-dependent two-time correlation function to provide a model of $`C(t,t^{})`$,

$$C_{\mathrm{mod}}(t,t^{})=C^{1/2}(t)\mathrm{exp}\left[_t^{}^t𝑑t^{\prime \prime }\eta _C(t^{\prime \prime })\right]C^{1/2}(t^{})$$
(29)
(for $`tt^{}`$). Consider the integral

$$A(t)=_0^t𝑑t^{}C_{\mathrm{mod}}^{}(t,t^{})C(t,t^{}).$$
(30)
This is the time-dependent analog of Eq. (21). Rather than try to use this to determine $`\eta _C`$ directly, it is more convenient to again take time derivatives. If $`C(t,t^{})`$ in Eq. (30) is replaced with $`C_{\mathrm{mod}}(t,t^{})`$ of Eq. (29), then

$$\frac{A(t)}{t}=C^2(t)[\eta _C(t)+\eta _C^{}(t)]A+\frac{1}{C(t)}\frac{C(t)}{t}A.$$
(31)
If we instead calculate $`A/t`$ with the full $`C(t,t^{})`$ in Eq. (30), and use Eq. (9) to evaluate $`C(t,t^{})/t`$, then

$`{\displaystyle \frac{A(t)}{t}}`$ $`=`$ $`C^2(t)[\eta (t)+\eta _C^{}(t)]A+{\displaystyle \frac{1}{2C(t)}}{\displaystyle \frac{C(t)}{t}}A`$ (33)
$`+\mathrm{\Theta }_3^{}(t)C^{1/2}(t)C_f^{1/2}(t),`$
where

$$\mathrm{\Theta }_3^{}(t)=_0^t𝑑t^{}\frac{C_{\mathrm{mod}}^{}(t,t^{})}{C^{1/2}(t)}_0^t^{}𝑑\overline{t}R^{}(t^{},\overline{t})\frac{C_f^{}(t,\overline{t})}{C_f^{1/2}(t)}.$$
(34)
Taking the time derivative of this, and using Eq. (25) for $`C_f(t,\overline{t})`$, gives

$$\frac{\mathrm{\Theta }_3(t)}{t}=C^{1/2}(t)\mathrm{\Theta }(t)[\eta _C(t)+\eta _f(t)]\mathrm{\Theta }_3(t).$$
(35)
Equating Eq. (31) and Eq. (33), one can then solve for the effective decorrelation rate $`\eta _C`$. Using Eq. (26) to eliminate the $`C(t)/t`$ term, the result is

$`\eta _C(t)`$ $`=`$ $`𝒫(\eta (t)Re\eta (t)+{\displaystyle \frac{C_f^{1/2}(t)Re\mathrm{\Theta }(t)}{C(t)}}`$ (37)
$`{\displaystyle \frac{\mathrm{\Theta }_3^{}(t)C^{1/2}(t)C_f^{1/2}(t)}{A(t)}}),`$
where we have added the $`𝒫`$ operator to enforce realizability for the reasons discussed below. Here $`𝒫(z)=z`$ if $`Rez0`$ and $`𝒫(z)=iImz`$ if $`Rez<0`$. Substituting Eq. (26) into Eq. (31) gives

$`{\displaystyle \frac{A(t)}{t}}`$ $`=`$ $`C^2(t)2Re(\eta (t)+\eta _C(t))A(t)`$ (39)
$`+\mathrm{\hspace{0.17em}2}Re\mathrm{\Theta }(t){\displaystyle \frac{C_f^{1/2}(t)}{C(t)}}A(t).`$
Eqs. (26), (28), and (35-39) provide a complete set of equations that can be integrated forward in time. They comprise a Markovian closure theory (including non-white noise effects) for the time-dependent Langevin equation. The relevant initial conditions are discussed below. This set of equations can be used to determine the amplitude $`C(t)`$ and the effective decorrelation rate $`\eta _C(t)`$ used to model the two-time behavior $`C(t,t^{})`$.

In a normal long-time statistical steady state, where $`\eta `$, $`\eta _f`$ and $`C_f`$ are constants (and $`Re(\eta )>0`$ and $`Re(\eta _f)>0`$), then one can show that the second and third terms on the right-hand side of Eq. (37) cancel and that it reproduces the steady-state result for $`\eta _C`$ in Eq. (22).

Consider the behavior of these equations in an unstable case, with $`Re\eta =\gamma <0`$. For simplicity, assume the coefficients $`\eta `$, $`\eta _f`$ and $`C_f`$ are all constant in time, with $`Re(\eta _f)>0`$. Then one can show that $`C(t)`$ eventually grows as $`\mathrm{exp}(2\gamma t)`$, while $`\mathrm{\Theta }(t)\mathrm{exp}((\gamma \eta _f)t)`$ grows more slowly, so that the third term on the right-hand side of Eq. (37) vanishes. The fourth term on the right-hand side of Eq. (37) also vanishes because $`\mathrm{\Theta }_3\mathrm{exp}(2\gamma t)`$ while $`A\mathrm{exp}(4\gamma t)`$. In this limit, $`\eta _C=\eta Re(\eta )`$.

Thus with constant coefficients, the two cases of positive or negative $`Re\eta `$ will, at least in the long-time limit, naturally reproduce the limiting operator $`𝒫(\eta )=Re\eta H(Re\eta )+iIm\eta `$, which was introduced by BKO to preserve realizability for the assumed form of $`C(t,t^{})`$ in Eq. (29). In the white-noise limit $`\eta _f\eta `$, it is straightforward to show that realizability is ensured for all time, not just in the long-time limit (see also Appendix (A)). These results might suggest that the $`𝒫`$ operator in Eq. (37) is not needed, if its argument always has a positive real part anyway. However, by numerically integrating Eqs. (26), (28), and (35-39), we have found cases where this is not true and the $`𝒫`$ operator is needed in Eq. (37) to enforce the realizability condition $`Re\eta _C0`$. \[Without the $`𝒫`$ operator, $`Re\eta _C`$ will transiently go negative in some strongly non-resonant cases such as $`\eta =1`$ and $`\eta _f=0.25+16i`$.\] Eqs. (26, 28) are an exact system of equations for the equal time covariance $`C(t)`$ for Langevin dynamics, which ensures that $`C(t)`$ is always positive. But according to Theorem 2 of BKO (and Appendix A of the present paper), $`Re\eta _C0`$ is necessary for $`C_{\mathrm{mod}}(t,t^{})`$ as given by Eq. (29) to be a realizable two-time correlation function. This may be important if $`C_{\mathrm{mod}}(t,t^{})`$ is in turn used in a noise term driving some other Fourier mode.

Formally, the initial conditions for this system of equations require some care to handle an apparent singularity, but in practice this should not be a problem. With a finite initial $`\psi (0)`$ in Eq. (7), the initial conditions for the Markovian closure equations are $`\mathrm{\Theta }(0)=0`$, $`A(0)=0`$, $`\mathrm{\Theta }_3(0)=0`$, and $`C(0)=C_1`$. For short times, we then have $`C(t)C_1`$, $`\mathrm{\Theta }(t)C_f^{1/2}t`$. If $`\eta _C`$ is finite, then for short times we also have $`A(t)=C_1^2t`$ and $`\mathrm{\Theta }_3=(C_1C_f)^{1/2}t^2/2`$. It follows from Eq. (37) that $`\eta _C=\eta Re\eta +tC_f/(2C_1)`$ for short times, which is a consistent solution that is finite and continuous, resolving the $`0/0`$ ambiguity in the last term of Eq. (37). In a numerical code, it is convenient to use the initial conditions $`\mathrm{\Theta }(0)=0`$, $`\mathrm{\Theta }_3(0)=0`$ (thus assuming the initial noise $`C_f=0`$), $`C(0)=C_1`$, and $`A(0)=C_1^2\mathrm{\Delta }t`$, where $`\mathrm{\Delta }t`$ is a time step smaller than any other relevant time scales in the problem.

## IV Formulation of the full nonlinear problem and statistical closures

In this section we provide background on the general form of the nonlinear problem we are considering and on the general theory of statistical closures. In particular we will write down Kraichnan’s direct-interaction approximation, which is the starting point of our calculation. This section borrows heavily from the BKO paper (including some of their wording), but is provided for completeness to define our starting point.

### A The fundamental nonlinear stochastic process

Consider a quadratically nonlinear equation, written in Fourier space, for some variable $`\psi _𝒌`$:

$$(\frac{}{t}+\nu _𝒌)\psi _𝒌(t)=\frac{1}{2}\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}M_{𝒌𝒑𝒒}\psi _𝒑^{}(t)\psi _𝒒^{}(t).$$
(40)
Here the time-independent coefficients of linear “damping” $`\nu _𝒌`$ and mode-coupling $`M_{𝒌𝒑𝒒}`$ may be complex. Given random initial conditions, we seek ensemble-averaged (or, if the system is ergodic, time-averaged) moments of $`\psi _𝒌(t)`$, taking for simplicity the mean value of $`\psi _𝒌`$ to be zero.

Many important nonlinear problems can be represented in this form with a simple quadratic nonlinearity. For example, the two-dimensional Navier–Stokes equation for neutral fluid turbulence can be written in this form, where $`\psi `$ represents the stream function such that the velocity $`𝒗=\widehat{𝒛}\times \psi `$, and $`M_{𝒌𝒑𝒒}=\widehat{𝐳}𝒑\times 𝒒(q^2p^2)/k^2`$. Other examples include Charney’s barotropic vorticity equation for planetary fluid flow, and a class of two-dimensional plasma drift wave turbulence problems (such as the Hasegawa–Mima equation or the Terry–Horton equation). Some three-dimensional one-field plasma turbulence problems can also be written in this form since the dominant $`\stackrel{}{E}\times \stackrel{}{B}`$ nonlinearity acts only in two dimensions perpendicular to the magnetic field. The three-dimensional Navier–Stokes equations and general multi-field plasma turbulence equations can also be written in the form of Eq. (40) if $`\psi _𝒌`$ is considered as a vector and $`\nu _𝒌`$ and $`M_{𝒌𝒑𝒒}`$ become matrices or tensors. In fact, BKO consider covariant multiple-field formulations of the DIA and Markovian closures. Here we will focus on the one-field case, where $`\psi _𝒌`$ is a scalar amplitude for mode $`𝒌`$.

For each $`𝒌`$ in Eq. (40), the summation on the right-hand-side involves a sum over all possible $`𝒑`$ and $`𝒒`$ that satisfy the three-wave interaction $`𝒌+𝒑+𝒒=0`$ (this is sometimes expressed as $`𝒌=𝒑+𝒒`$, but the reality conditions $`\psi _𝒌=\psi _𝒌^{}`$ has been used to rearrange it). Without any loss of generality one may assume the symmetry

$$M_{𝒌𝒑𝒒}=M_{𝒌𝒒𝒑}.$$
(41)
Another important symmetry possessed by many such systems is

$$\sigma _𝒌M_{𝒌𝒑𝒒}+\sigma _𝒑M_{𝒑𝒒𝒌}+\sigma _𝒒M_{𝒒𝒌𝒑}=0$$
(42)
for some time-independent nonrandom real quantity $`\sigma _𝒌`$. \[See Refs.  and  for the relation between this symmetry and the Manley-Rowe relations for wave actions.\] Equation (42) is easily shown to imply that the nonlinear terms of Eq. (40) conserve the ensemble-averaged total generalized energy $`E\frac{1}{2}_𝒌\sigma _𝒌\left|\psi _𝒌(t)\right|^2`$. \[The nonlinear terms also conserve the generalized energy in each individual realization, although we will be focusing on ensemble-averaged quantities, where $`\mathrm{}`$ denotes ensemble-averaging.\] For some problems, Eq. (42) may be satisfied by more than one choice of $`\sigma _𝒌`$; this implies the existence of more than one nonlinear invariant. For example, in the case of two-dimensional hydrodynamics, Eq. (42) is satisfied for both $`\sigma _𝒌=k^2`$ and $`\sigma _𝒌=k^4`$, which correspond to the conservation of energy and enstrophy, respectively.

We define the two-time correlation function $`C_𝒌(t,t^{})\psi _𝒌(t)\psi _𝒌^{}(t^{})`$ and the equal-time correlation function $`C_𝒌(t)C_𝒌(t,t)`$ (note that the two functions are distinguished only by the number of arguments), so that $`E=\frac{1}{2}_𝒌\sigma _𝒌C_𝒌(t)`$. In stationary turbulence, the two-time correlation function depends on only the difference of its time arguments: $`C_𝒌(t,t^{})𝒞_𝒌(tt^{})`$. The renormalized infinitesimal response function (nonlinear Green’s function) $`R_𝒌(t,t^{})`$ is the ensemble-averaged infinitesimal response to a source function $`S_𝒌(t)`$ added to the right-hand side of Eq. (40) for mode $`𝒌`$ alone. As a functional derivative,

$$R_𝒌(t,t^{})\frac{\delta \psi _𝒌(t)}{\delta S_𝒌(t^{})}|_{S_𝒌=0}.$$
(43)
We adopt the convention that the equal-time response function $`R_𝒌(t,t)`$ evaluates to $`1/2`$ \[although $`lim_{ϵ0+}`$ $`R_𝒌(t+ϵ,t)=1`$\].

### B Statistical closures; the direct-interaction approximation

The starting point of our derivation will be the equations of Kraichnan’s direct-interaction approximation (DIA), as given in Eqs. (6-7) of BKO, and reproduced below as Eqs. (IV B-IV B).

The general form of a statistical closure in the absence of mean fields is

$`({\displaystyle \frac{}{t}}+\nu _𝒌)C_𝒌(t,t^{})`$ $`+`$ $`{\displaystyle _0^t}𝑑\overline{t}\mathrm{\Sigma }_𝒌(t,\overline{t})C_𝒌(\overline{t},t^{})`$ (45)
$`=`$ $`{\displaystyle _0^t^{}}𝑑\overline{t}_𝒌(t,\overline{t})R_𝒌^{}(t^{},\overline{t}),`$ (46)
$`({\displaystyle \frac{}{t}}+\nu _𝒌)R_𝒌(t,t^{})`$ $`+`$ $`{\displaystyle _t^{}^t}𝑑\overline{t}\mathrm{\Sigma }_𝒌(t,\overline{t})R_𝒌(\overline{t},t^{})`$ (47)
$`=`$ $`\delta \left(tt^{}\right).`$ (48)
While these equations (with the expressions for $`\mathrm{\Sigma }_𝒌`$ and $`_𝒌`$ given below) are an approximate statistical solution to Eq. (40), they are the exact statistical solution to a generalized Langevin equation

$$\left(\frac{}{t}+\nu _𝒌\right)\psi _𝒌(t)+_0^t𝑑\overline{t}\mathrm{\Sigma }_𝒌(t,\overline{t})\psi _𝒌(\overline{t})=f_𝒌(t),$$
(49)
where $`\mathrm{\Sigma }_𝒌`$ is the kernel of a non-local damping/propagation operator, and $`_𝒌(t,\overline{t})=f_k(t)f_k^{}(\overline{t})`$. These equations specify an initial-value problem for which $`t=0`$ is the initial time.

The original nonlinearity in Eq. (40) gives rise to two types of terms in Eqs. (IV B): those describing nonlinear damping ($`\mathrm{\Sigma }_𝒌`$) and one modeling nonlinear noise ($`_𝒌`$). The nonlinear damping and noise in Eqs. (IV B) are determined on the basis of fully nonlinear statistics.

The direct-interaction approximation provides specific approximate forms for $`\mathrm{\Sigma }_𝒌`$ and $`_𝒌`$:

$$\mathrm{\Sigma }_𝒌(t,\overline{t})=\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}R_𝒑^{}(t,\overline{t})C_𝒒^{}(t,\overline{t}),$$
(51)
$$_𝒌(t,\overline{t})=\frac{1}{2}\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}\left|M_{𝒌𝒑𝒒}\right|^2C_𝒑^{}(t,\overline{t})C_𝒒^{}(t,\overline{t}).$$
(52)
These renormalized forms can be obtained from the formal perturbation series by retaining only selected terms. While there are infinitely many ways of obtaining a renormalized expression, Kraichnan has shown that most of the resulting closed systems of equations lead to physically unacceptable solutions. For example, they might predict the physically impossible situation of a negative value for $`C_𝒌(t,t)`$ (i.e., a negative energy)! Such behavior cannot occur in the DIA or other realizable closures.

The DIA also conserves all of the same generalized energies ($`\frac{1}{2}_𝒌\sigma _𝒌\left|\psi _𝒌(t)\right|^2`$) that are conserved by the primitive dynamics. To show this important property, it is useful to write the equal-time covariance equation in the form

$$\frac{}{t}C_𝒌(t)+2ReN_𝒌(t)=2ReF_𝒌(t),$$
(54)
where

$$N_𝒌(t)\nu _𝒌C_𝒌(t)\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\overline{\mathrm{\Theta }}_{𝒑𝒒𝒌}^{}(t),$$
(55)
$$F_𝒌(t)\frac{1}{2}\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}\left|M_{𝒌𝒑𝒒}\right|^2\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}^{}(t),$$
(56)
$$\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}(t)_{t_0}^t𝑑\overline{t}R_𝒌(t,\overline{t})C_𝒑(t,\overline{t})C_𝒒(t,\overline{t}),$$
(57)
given initial conditions at the time $`t=t_0`$ (unless otherwise stated, we will take $`t_0=0`$). As shown in BKO, the symmetries (41) and (42) ensure that Eq. (54) conserves all quadratic nonlinear invariants of the form $`E\frac{1}{2}_𝒌\sigma _𝒌C_𝒌(t)`$ in the dissipationless case where $`Re\nu _𝒌=0`$. The Markovian closures that BKO developed, and that we extend here, preserve the structure of Eqs. (IV B) and so have all of the same quadratic nonlinear conservation properties as the original equations. \[One can show that $`F_𝒌`$ is always real, so the $`Re`$ operation on $`F_𝒌`$ in Eq. (54) is redundant.\]

The DIA equations (IV B) and (IV B) provide a closed set of equations, but are fairly complicated because they involve convolutions over two-time functions. Their general numerical solution requires $`𝒪(N_t^3)`$ operations, or $`𝒪(N_t^2)`$ operations in steady state. As described in BKO and Krommes, a Markovian approximation seeks to simplify this complexity by parameterizing the two-time functions in terms of a single decorrelation rate. Our approach here is essentially to generalize this to allow several rate parameters to be used, to allow the decorrelation rate for $`C_𝒌(t,t^{})`$ to differ from the decay rate for $`R_𝒌(t,t^{})`$.

## V Response functions in a statistical steady state

Markovian models provide approximations that can simplify the integrals in Eqs. (IV B). For insight, we will first investigate the long-time limit where a statistical steady-state should be reached, so that the two-time correlation function $`C(t,t^{})`$ and response function $`R(t,t^{})`$ can depend only on the time difference $`tt^{}`$. In a statistical steady state, all of the Markovian models in BKO use a simple exponential behavior for $`C_𝒌(t,t^{})`$ and $`R_𝒌(t,t^{})`$. Here we will assume the model forms

$$R_{\mathrm{mod},𝒌}(t,t^{})=\mathrm{exp}\text{(}\eta _𝒌(tt^{})\text{)}H(tt^{})$$
(58)
and

$$C_{\mathrm{mod},𝒌}(t,t^{})\{\begin{array}{cc}C_{0𝒌}\mathrm{exp}\text{(}\eta _{C𝒌}(tt^{})\text{)}\hfill & \text{for }tt^{}\text{,}\hfill \\ C_{0𝒌}\mathrm{exp}\text{(}\eta _{C𝒌}^{}(tt^{})\text{)}\hfill & \text{for }t<t^{}\text{.}\hfill \end{array}$$
(59)
Note that $`\eta _𝒌`$ is the decay rate for the infinitesimal response function $`R_𝒌`$, while $`\eta _{C𝒌}`$ is the decorrelation rate for $`C_𝒌(t,t^{})`$.

Inserting Eq. (51) into Eq. (48) and using the exponential forms of Eq. (58) and Eq. (59) in the integrals yields

$`({\displaystyle \frac{}{t}}+\nu _𝒌)R_𝒌(t,t^{})=\delta \left(tt^{}\right)`$ (60)
$`+`$ $`{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_{0𝒒}H(tt^{})`$ (62)
$`\times {\displaystyle _t^{}^t}d\overline{t}\mathrm{exp}\text{(}(\eta _𝒑^{}+\eta _{C𝒒}^{})(t\overline{t})\eta _𝒌(\overline{t}t^{})\text{)}.`$
Evaluating the integral gives

$`({\displaystyle \frac{}{t}}+\nu _𝒌)R_𝒌(t,t^{})=\delta \left(tt^{}\right)`$ (63)
$`+`$ $`{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_{0𝒒}}{\eta _𝒑^{}+\eta _{C𝒒}^{}\eta _𝒌}}H(t,t^{})`$ (64)
$`\times `$ $`\left[\mathrm{exp}\text{(}\eta _𝒌(tt^{})\text{)}\mathrm{exp}\text{(}(\eta _𝒑^{}+\eta _{C𝒒}^{})(tt^{})\text{)}\right].`$ (65)
The solution to this equation for $`t>t^{}`$ is

$`R_𝒌(t,t^{})=\mathrm{exp}\text{(}\nu _𝒌(tt^{})\text{)}`$ (66)
$`+{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_{0𝒒}}{\eta _𝒑^{}+\eta _{C𝒒}^{}\eta _𝒌}}`$ (67)
$`\times [{\displaystyle \frac{\mathrm{exp}\text{(}\nu _𝒌(tt^{})\text{)}\mathrm{exp}\text{(}\eta _𝒌(tt^{})\text{)}}{\eta _𝒌\nu _𝒌}}`$ (68)
$`{\displaystyle \frac{\mathrm{exp}\text{(}\nu _𝒌(tt^{})\text{)}\mathrm{exp}\text{(}(\eta _𝒑^{}+\eta _{C𝒒}^{})(tt^{})\text{)}}{\eta _𝒑^{}+\eta _{C𝒒}^{}\nu _𝒌}}].`$ (69)
Clearly this is not strictly consistent with the simple exponential form for $`R_𝒌`$ assumed in Eq. (58) and used to evaluate the integrals in Eq. (48). We will instead fit the model Eq. (58) to Eq. (69), in the same way that we did in the Langevin case for Eq. (21). Requiring that both Eq. (58) and the full Eq. (69) give the same weighted average over time (where $`R_{\mathrm{mod},𝒌}^{}`$ is used as the weight to ensure invariance to frequency shifts) gives

$`{\displaystyle \frac{1}{\eta _𝒌+\eta _𝒌^{}}}{\displaystyle _t^{}^{\mathrm{}}}𝑑tR_{\mathrm{mod},𝒌}^{}(t,t^{})R_𝒌(t,t^{}).`$ (70)
Inserting Eq. (69) on the right-hand side, and carrying out a few lines of algebra, the result is

$`{\displaystyle \frac{1}{\eta _𝒌+\eta _𝒌^{}}}={\displaystyle \frac{1}{\nu _𝒌+\eta _𝒌^{}}}`$ (71)
$`+{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_{0𝒒}}{(\nu _𝒌+\eta _𝒌^{})(\eta _𝒌+\eta _𝒌^{})(\eta _𝒑^{}+\eta _{C𝒒}^{}+\eta _𝒌^{})}}.`$ (72)
A little rearranging leads to

$$\eta _𝒌\nu _𝒌\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}\frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_{0𝒒}}{\eta _𝒌^{}+\eta _𝒑^{}+\eta _{C𝒒}^{}}.$$
(73)
Note that this has a similar form to the steady-state decay rate in the DIA-based EDQNM, such as in Eq. (39b) of BKO (but with their $`\eta _𝒒^{}`$ replaced by $`\eta _{C𝒒}^{}`$).

One can go through a similar calculation of $`C_𝒌(t,t^{})`$, and calculate its weighted time average to determine the decorrelation rate $`\eta _{C𝒌}`$. We will not do so now, as one can instead just take the steady-state limit of the results in the next section. Eq. (73) can also be obtained from the steady-state limit of the results in the next section, and so provides a useful cross-check.

We note that there is some flexibility in the choice of weighting in Eq. (70). We could use $`C_{\mathrm{mod},𝒌}^{}(t,t^{})`$ as the weight instead of $`R_{\mathrm{mod},𝒌}^{}(t,t^{})`$. Either choice preserves Galilean invariance. Using this alternate weight, Eq. (70) becomes

$`{\displaystyle \frac{C_{𝒌0}}{\eta _𝒌+\eta _{C𝒌}^{}}}{\displaystyle _t^{}^{\mathrm{}}}𝑑tC_{\mathrm{mod},𝒌}^{}(t,t^{})R_𝒌(t,t^{})`$ (74)
and the resulting expression for $`\eta _𝒌`$ is like Eq. (73) but with $`\eta _𝒌^{}`$ on the right-hand side of Eq. (73) replaced by $`\eta _{C𝒌}^{}`$, which would automatically agree with the steady-state $`\overline{\eta }_k`$ to be defined in Eq. (110). But it turns out that the main steady-state results of Sec. (VII) hold with either choice of weights, and it seems more symmetric and makes more sense as a standard fitting procedure to use $`R_{\mathrm{mod},𝒌}^{}`$ as the weight for integrating $`R_𝒌`$ in Eq. (70). This raises the question of whether to use $`C_{\mathrm{mod},𝒌}^{}`$ or $`R_{\mathrm{mod},𝒌}^{}`$ as the weight function for time averages of $`C_𝒌(t,t^{})`$, as we will do in the next section. We can resolve this ambiguity by going back to the steady-state Langevin problem of Sec. (III A). If one tries to use $`R^{}(t,t^{})`$ as the weight in Eq. (21), so that it becomes

$$\frac{C_0}{\eta _C+\eta ^{}}_{\mathrm{}}^t𝑑t^{}\mathrm{exp}(\eta ^{}(tt^{}))C(t,t^{}),$$
(75)
then one can go through the same steps used to derive Eq. (24) and find that in the limit of real coefficients it gives $`\eta _C=\eta \eta _f/(2\eta +\eta _f)`$. In the red-noise limit $`\eta _f\eta `$, this gives $`\eta _C=\eta _f/2`$, which is a factor of 2 off from the correct result ($`\eta _C=\eta _f`$) for the red noise limit. Thus, we will use $`C_{\mathrm{mod},𝒌}^{}(t,t^{})`$ as the weight for taking time-averages of $`C_𝒌(t,t^{})`$ and use $`R_{\mathrm{mod},𝒌}^{}(t,t^{})`$ for time-averaging $`R_𝒌(t,t^{})`$. The weighting choices might be reconsidered in a multi-field generalization of a Markovian closure, where the requirement of covariance may impose constraints on the choice of the weight functions, but it seems that the symmetric choices made here are most likely to generalize well.

## VI Time-Dependent Multiple-Rate Markovian Closure

Applying these techniques in a straightforward way to the time-dependent DIA equations leads to the Multiple-Rate Markovian Closure (MRMC) equations. The two-time correlation function is modeled with the realizable form

$$C_{\mathrm{mod},𝒌}(t,t^{})=C_𝒌^{1/2}(t)C_𝒌^{1/2}(t^{})\mathrm{exp}\left(_t^{}^t𝑑\overline{t}\eta _{C𝒌}(\overline{t})\right)$$
(76)
(for $`t>t^{}`$, with $`C_{\mathrm{mod},𝒌}(t,t^{})=C_{\mathrm{mod},𝒌}^{}(t^{},t)`$ for $`t<t^{}`$), and the response function is modeled as

$$R_{\mathrm{mod},𝒌}(t,t^{})=\mathrm{exp}\left(_t^{}^t𝑑\overline{t}\eta _𝒌(\overline{t})\right)H(tt^{}).$$
(77)
Denoting $`\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}(t)=\mathrm{\Theta }_{𝒌𝒑𝒒}(t)C_p^{1/2}(t)C_q^{1/2}(t)`$, and inserting Eqs. (76-77) into Eq. (57), we can write the equal-time DIA covariance equations of Eq. (IV B) as

$$\frac{}{t}C_𝒌(t)+2Re\overline{\eta }_𝒌(t)C_𝒌(t)=2F_𝒌(t),$$
(79)
$$\overline{\eta }_𝒌\nu _𝒌\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\mathrm{\Theta }_{𝒑𝒒𝒌}^{}(t)C_𝒒^{1/2}(t)C_𝒌^{1/2}(t),$$
(80)
$$F_𝒌\frac{1}{2}\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}\left|M_{𝒌𝒑𝒒}\right|^2\mathrm{\Theta }_{𝒌𝒑𝒒}^{}(t)C_𝒑^{1/2}(t)C_𝒒^{1/2}(t),$$
(81)
$$\frac{}{t}\mathrm{\Theta }_{𝒌𝒑𝒒}+(\eta _𝒌+\eta _{C𝒑}+\eta _{C𝒒})\mathrm{\Theta }_{𝒌𝒑𝒒}=C_p^{1/2}(t)C_q^{1/2}(t),$$
(82)
$$\mathrm{\Theta }_{𝒌𝒑𝒒}(0)=0.$$
(83)
This is very similar to the Bowman–Krommes–Ottaviani Realizable Markovian Closure (RMC) (as given by Eqs. (66a–e) of BKO), but with the replacement of the single decay/decorrelation rate of the RMC with three different rates in these equations. \[Other Markovian models, such as the EDQNM closure, also use a single decorrelation rate parameter.\] If in Eq. (82) we replace $`\eta _𝒌=\overline{\eta }_𝒌`$, $`\eta _{C𝒑}=𝒫(\overline{\eta }_𝒑)`$, and $`\eta _{C𝒒}=𝒫(\overline{\eta }_𝒒)`$, then these equations become identical to the RMC.

To summarize the three different rates used here:

* $`\overline{\eta }_𝒌`$ is the nonlinear energy damping rate for the wave energy equation for the equal-time covariance $`C_k(t)`$ in Eq. (79), and is defined in Eq. (80);
* $`\eta _𝒌`$ is the decay rate for the infinitesimal response function $`R_k(t,t^{})`$ in Eq. (77), and is defined in Eq. (91);
* and $`\eta _{C𝒌}`$ is the decorrelation rate for $`C_k(t,t^{})`$ in Eq. (76), and is defined in Eq. (107).

To determine $`\eta _𝒌(t)`$ and $`\eta _{C𝒌}(t)`$, we follow a similar procedure as we did for the time-dependent Langevin equation in Sec. (III C). Define $`A_𝒌(t)`$ as the following weighted time-average of $`R_𝒌`$

$$A_𝒌(t)=_0^t𝑑t^{}R_{\mathrm{mod},𝒌}^{}(t,t^{})R_𝒌(t,t^{}).$$
(84)
If $`R_𝒌(t,t^{})=R_{\mathrm{mod},𝒌}(t,t^{})`$ as given by Eq. (77), then

$$\frac{A_𝒌}{t}=1(\eta _𝒌^{}+\eta _𝒌)A_𝒌,$$
(85)
while if $`R_𝒌(t,t^{})`$ satisfies Eqs.(48,51), then

$`{\displaystyle \frac{A_𝒌}{t}}`$ $`=`$ $`1(\eta _𝒌^{}+\nu _𝒌)A_𝒌`$ (87)
$`+{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒒^{1/2}\mathrm{\Theta }_{1,𝒑𝒒𝒌}^{},`$
where

$`\mathrm{\Theta }_{1,𝒑𝒒𝒌}^{}(t)=`$ (88)
$`{\displaystyle _0^t}𝑑t^{}R_{\mathrm{mod},𝒌}^{}(t,t^{}){\displaystyle _t^{}^t}𝑑\overline{t}{\displaystyle \frac{C_𝒒^{}(t,\overline{t})}{C_𝒒^{1/2}(t)}}R_𝒑^{}(t,\overline{t})R_𝒌(\overline{t},t^{}).`$ (89)
It is often more convenient to work with the differential version of this, which, after using Eqs. (76-77) to replace $`C_𝒒(t,t^{})`$ and $`R_𝒑(t,t^{})`$ with their model forms, is

$$\frac{\mathrm{\Theta }_{1,𝒑𝒒𝒌}^{}}{t}=(\eta _𝒌^{}+\eta _{C𝒒}^{}+\eta _𝒑^{})\mathrm{\Theta }_{1,𝒑𝒒𝒌}^{}+C_𝒒^{1/2}A_k(t)$$
(90)
(with the initial condition $`\mathrm{\Theta }_{1,𝒌𝒑𝒒}(0)=0)`$). Requiring that Eq. (85) and Eq. (87) be equivalent determines $`\eta _𝒌`$ to be

$$\eta _𝒌=\nu _𝒌\frac{1}{A_𝒌}\underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒒^{1/2}\mathrm{\Theta }_{1,𝒑𝒒𝒌}^{}.$$
(91)
The calculation of $`\eta _{C𝒌}`$ proceeds in a similar way. $`A_{C𝒌}(t)`$ is defined as a weighted time integral of $`C_𝒌(t,t^{})`$:

$$A_{C𝒌}(t)=_0^t𝑑t^{}C_{\mathrm{mod},𝒌}^{}(t,t^{})C_𝒌(t,t^{}).$$
(92)
If $`C_𝒌(t,t^{})`$ in this integral is replaced by $`C_{\mathrm{mod},𝒌}(t,t^{})`$ as given by Eq. (76), then

$$\frac{A_{C𝒌}}{t}=C_𝒌^2(t)+\frac{1}{C_𝒌(t)}\frac{C_𝒌(t)}{t}A_{C𝒌}(\eta _{C𝒌}^{}+\eta _{C𝒌})A_{C𝒌},$$
(93)
(where we make the time dependence of $`C_𝒌(t)`$ explicit to distinguish it from the two-time $`C_𝒌(t,t^{})`$). If the exact dynamics for $`C_𝒌(t,t^{})`$ given by Eqs. (46,IV B) are used, then

$`{\displaystyle \frac{A_{C𝒌}}{t}}`$ $`=`$ $`C_𝒌^2(t)+{\displaystyle \frac{1}{2C_𝒌(t)}}{\displaystyle \frac{C_𝒌(t)}{t}}A_{C𝒌}(\eta _{C𝒌}^{}+\nu _𝒌)A_{C𝒌}`$ (94)
$`+`$ $`{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒌^{1/2}C_𝒒^{1/2}\mathrm{\Theta }_{2,𝒑𝒒𝒌}^{}`$ (95)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}|M_{𝒌𝒑𝒒}|^2C_𝒌^{1/2}C_𝒑^{1/2}C_𝒒^{1/2}\mathrm{\Theta }_{3,𝒌𝒑𝒒}^{},`$ (96)
where

$`\mathrm{\Theta }_{2,𝒑𝒒𝒌}^{}(t)=`$ (97)
$`{\displaystyle _0^t}𝑑t^{}{\displaystyle \frac{C_{\mathrm{mod},𝒌}^{}(t,t^{})}{C_𝒌^{1/2}(t)}}{\displaystyle _0^t}𝑑\overline{t}{\displaystyle \frac{C_𝒒^{}(t,\overline{t})}{C_𝒒^{1/2}(t)}}R_𝒑^{}(t,\overline{t})C_𝒌(\overline{t},t^{})`$ (98)
and

$`\mathrm{\Theta }_{3,𝒌𝒑𝒒}^{}(t)=`$ (99)
$`{\displaystyle _0^t}𝑑t^{}{\displaystyle \frac{C_{\mathrm{mod},𝒌}^{}(t,t^{})}{C_𝒌^{1/2}(t)}}{\displaystyle _0^t^{}}𝑑\overline{t}{\displaystyle \frac{C_𝒒^{}(t,\overline{t})}{C_𝒒^{1/2}(t)}}{\displaystyle \frac{C_𝒑^{}(t,\overline{t})}{C_𝒑^{1/2}(t)}}R_𝒌^{}(t^{},\overline{t}).`$ (100)
Using Eqs. (76-77), the differential versions of these are

$`{\displaystyle \frac{\mathrm{\Theta }_{2,𝒑𝒒𝒌}^{}}{t}}`$ $`=`$ $`(\eta _{C𝒌}^{}+\eta _{C𝒒}^{}+\eta _𝒑^{})\mathrm{\Theta }_{2,𝒑𝒒𝒌}^{}`$ (102)
$`+C_𝒌(t)\mathrm{\Theta }_{𝒑𝒒𝒌}^{}(t)+{\displaystyle \frac{C_𝒒^{1/2}(t)}{C_𝒌^{1/2}(t)}}A_{C𝒌}(t)`$
and

$`{\displaystyle \frac{\mathrm{\Theta }_{3,𝒌𝒑𝒒}^{}}{t}}`$ $`=`$ $`(\eta _{C𝒌}^{}+\eta _{C𝒒}^{}+\eta _{C𝒑}^{})\mathrm{\Theta }_{3,𝒌𝒑𝒒}^{}`$ (104)
$`+C_𝒌^{1/2}(t)\mathrm{\Theta }_{𝒌𝒑𝒒}^{}(t).`$
The quantity $`\eta _{C𝒌}`$ is then determined by equating Eq. (93) and Eq. (96), yielding

$`\eta _{C𝒌}=\nu _𝒌+{\displaystyle \frac{1}{2C_𝒌(t)}}{\displaystyle \frac{C_𝒌(t)}{t}}`$ (105)
$`{\displaystyle \frac{1}{A_{C𝒌}}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒌^{1/2}C_𝒒^{1/2}\mathrm{\Theta }_{2,𝒑𝒒𝒌}^{}`$ (106)
$`{\displaystyle \frac{1}{2A_{C𝒌}}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}|M_{𝒌𝒑𝒒}|^2C_𝒌^{1/2}C_𝒑^{1/2}C_𝒒^{1/2}\mathrm{\Theta }_{3,𝒌𝒑𝒒}^{},`$ (107)
where Eq. (79) could be used to eliminate $`C_𝒌(t)/t`$. As in Eq. (37) for the case of the time-dependent Langevin equation, while there are effects in this equation that will tend to give $`Re\eta _{C𝒌}0`$, it may be necessary to modify this equation to enforce realizability in all cases. This is done by replacing this equation, of the form $`\eta _{C𝒌}=\mathrm{RHS}`$, with $`\eta _{C𝒌}=𝒫(\mathrm{RHS})`$. Note that it is only $`Re\eta _{C𝒌}0`$ that is needed for realizability, while $`Re\eta _𝒌`$ can transiently go negative (as it does in two-dimensional hydrodynamics because of the inverse cascade, or in some plasma problems where the zonal flows may become nonlinearly unstable). This is similar to BKO’s treatment.

The complete set of equations that constitutes the Multiple-Rate Markovian Closure (MRMC) are Eqs. (77) for the equal-time covariance $`C_𝒌(t)`$ and related quantities, Eqs. (85,90,91) for quantities related to the response function, and Eqs. (93,102-107) for quantities related to the two-time correlation function. The MRMC extends the RMC to make less restrictive assumptions and include additional effects, but at the expense of a few new parameters. In addition to replacing the single decay/decorrelation rate of the RMC with 3 different rates, $`\eta _𝒌`$, $`\overline{\eta }_𝒌`$, and $`\eta _{C𝒌}`$, it also replaces the single triad interaction time of the RMC with 4 different triad interaction times, $`\mathrm{\Theta }_{𝒌𝒑𝒒}`$, $`\mathrm{\Theta }_{1,𝒌𝒑𝒒}`$, $`\mathrm{\Theta }_{2,𝒌𝒑𝒒}`$, and $`\mathrm{\Theta }_{3,𝒌𝒑𝒒}`$. Each of these triad interaction times has a different weighting of response functions and two-time correlation functions. While this increases the complexity some, the overall computational scaling of this system is still $`𝒪(N_t)`$, a significant improvement over the $`𝒪(N_t^2)`$ or $`𝒪(N_t^3)`$ scaling of the DIA.

## VII Properties of the Multiple-Rate Markovian Closure

In a steady-state limit, Eq. (91) simplifies to

$`\eta _𝒌=\nu _𝒌{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒒}{\eta _𝒌^{}+\eta _𝒑^{}+\eta _{C𝒒}^{}}}.`$ (108)
The steady-state balance $`Re\overline{\eta }_𝒌C_k=F_k`$ from Eq. (79) simplifies to

$`C_𝒌Re\left[\nu _𝒌{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒒}{\eta _{C𝒌}^{}+\eta _𝒑^{}+\eta _{C𝒒}^{}}}\right]`$ (109)
$`={\displaystyle \frac{1}{2}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{|M_{𝒌𝒑𝒒}|^2C_𝒑C_𝒒}{\eta _𝒌^{}+\eta _{C𝒑}^{}+\eta _{C𝒒}^{}}}.`$ (110)
\[Note the subtle notational differences: the expression for $`\eta _𝒌`$ becomes the expression for $`\overline{\eta }_𝒌`$ if $`\eta _𝒌^{}`$ on the RHS of Eq. (108) is replaced by $`\eta _{C𝒌}^{}`$.\] Finally, Eq. (107) reduces to

$`\eta _{C𝒌}\eta _𝒌(\eta _{C𝒌}+\eta _{C𝒌}^{})[{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}C_𝒒}{(\eta _{C𝒌}^{}+\eta _𝒑^{}+\eta _{C𝒒}^{})^2}}`$ (111)
$`+{\displaystyle \frac{1}{2C_𝒌}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{|M_{𝒌𝒑𝒒}|^2C_𝒑C_𝒒}{(\eta _{C𝒌}^{}+\eta _{C𝒑}^{}+\eta _{C𝒒}^{})(\eta _𝒌^{}+\eta _{C𝒑}^{}+\eta _{C𝒒}^{})}}].`$ (112)
(113)
Thus the decorrelation rate $`\eta _{C𝒌}`$ equals the response function decay rate $`\eta _𝒌`$ plus the two correction terms in brackets. For the simple steady-state non-wave case with real and positive $`\eta `$’s, the second correction term will cause $`\eta _{C𝒌}`$ to decrease (as expected for non-white noise), while the first term will usually have an offsetting opposite sign and cause $`\eta _{C𝒌}`$ to increase. The origin of these two terms can be traced back to the DIA Eq. (46). The second correction term corresponds to the usual effects of non-white noise (related to the integral involving $`_𝒌(t,\overline{t})`$ in Eq. (46)), but the first correction term in Eq. (113) is related to the time-history integral involving the renormalized propagator $`\mathrm{\Sigma }_𝒌(t,\overline{t})`$ in Eq. (46). Thus non-white fluctuations in other modes $`C_𝒒(t,\overline{t})`$ not only change the noise term for the $`𝒌`$ mode, but also change the effective damping from the time-history integral, broadening the width of $`\mathrm{\Sigma }_𝒌(t,\overline{t})`$ in time (if the fluctuations $`C_𝒒`$ were treated as white noise, then Eq. (51) would give $`\mathrm{\Sigma }_𝒌(t,\overline{t})\delta (t\overline{t})`$).

An important property to demonstrate is that in thermal equilibrium it is possible for these two terms to cancel exactly. Then the decorrelation rate and response function decay rate are equivalent, $`\eta _{C𝒌}=\eta _𝒌`$, and the fluctuation–dissipation theorem is satisfied. To demonstrate that this is true, we assume the result ($`\eta _{C𝒌}=\eta _𝒌`$) to simplify some of the equations and then show that this is a self-consistent assumption. (Note also that if $`\eta _{C𝒌}=\eta _𝒌`$, then $`\overline{\eta }_𝒌=\eta _𝒌`$ also.) Splitting the first summation in brackets in Eq. (113) into two equal parts and interchanging the $`𝒑`$ and $`𝒒`$ labels for one of these parts (i.e., using an identity of the form $`\mathrm{\Sigma }G_{𝒌𝒑𝒒}=\mathrm{\Sigma }G_{𝒌𝒑𝒒}/2+\mathrm{\Sigma }G_{𝒌𝒒𝒑}/2`$), the terms in brackets in Eq. (113) can be written as

$`{\displaystyle \frac{1}{2C_𝒌}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}{\displaystyle \frac{M_{𝒌𝒑𝒒}}{(\eta _𝒌^{}+\eta _𝒑^{}+\eta _𝒒^{})^2}}`$ (114)
$`\times (M_{𝒑𝒒𝒌}^{}C_𝒒C_𝒌+M_{𝒒𝒑𝒌}^{}C_𝒑C_𝒌+M_{𝒌𝒑𝒒}^{}C_pC_q).`$ (115)
In thermal equilibrium, the spectrum $`C_𝒌`$ is given by equipartition among modes of a generalized energy-like conserved quantity. Consider an equipartition spectrum of the form $`C_𝒌=1/\lambda _𝒌`$, where $`\lambda _𝒌=_i\alpha ^{(i)}\sigma _𝒌^{(i)}`$, the $`\sigma _𝒌^{(i)}`$ are the coefficients in Eq. (42) (related to the quadratic invariants), and $`\alpha ^{(i)}`$ are determined by the initial conditions. Substituting $`C_𝒌=1/\lambda _𝒌`$, $`C_𝒑=1/\lambda _𝒑`$, and $`C_𝒒=1/\lambda _𝒒`$ into Eq. (115), and using Eqs. (41-42), one can show that Eq. (115) indeed vanishes, so that Eq. (113) simplifies to $`\eta _{C𝒌}=\eta _𝒌`$. The proof that $`C_𝒌=1/\lambda _𝒌`$ is a solution of the steady-state Eq. (110) proceeds in a similar way, interchanging the $`𝒑`$ and $`𝒒`$ labels for half of the summation on the left-hand side of Eq. (110), and noting that $`Re\nu _𝒌=0`$ in an isolated thermal system, etc. Rigorously, this only shows that the equipartition spectrum $`C_𝒌=1/\lambda _𝒌`$ is an equilibrium solution. This paper doesn’t demonstrate that it is a stable equilibrium or that mixing dynamics will necessarily relax to this state. For a discussion of the Gibbs-type $`H`$ theorem that leads to this result, see Appendix H of Ref.  and Refs. and . It is significant to note that the thermal equilibrium result holds even if the number of modes is small, and it does not assume that the noise spectrum is white. This is unlike a simple Langevin equation of the form of Eq. (2) (which has a local damping term in contrast to the time-history integral of Eq. (49)), where the two-time correlation function and the infinitesimal response function are proportional only if the noise is white.

We next estimate the importance of the correction terms in the decorrelation rate for an inertial range of a turbulent steady state, such as in two-dimensional hydrodynamics. Typically most of the energy is at long wavelengths ($`C_𝒒`$ is peaked at sufficiently low $`𝒒`$), so that the dominant contributions to the sums in Eq. (113) come from long wavelengths: $`|𝒒||𝒌|`$ in the first sum, and $`|𝒒||𝒌|`$ or $`|𝒑|=|𝒌+𝒒||𝒌|`$ in the second sum. This means that one can approximate the denominators in the sums of Eq. (113) using, for example, $`(\eta _{C𝒌}+\eta _𝒑+\eta _{C𝒒})(\eta _{C𝒌}+\eta _𝒌)`$ (since $`\eta _𝒌`$ and $`\eta _{C𝒌}`$ are typically increasing functions of $`𝒌`$). Similar approximations give $`(\overline{\eta }_𝒌\nu _k)(\eta _𝒌\nu _k)2\eta _𝒌^{}/(\eta _{C𝒌}^{}+\eta _𝒌^{})`$. Using the steady-state relation $`F_𝒌=Re\overline{\eta }_𝒌C_k`$ from Eq. (110) and the disparate scale approximations to rewrite the second sum in Eq. (113) in terms of $`\overline{\eta }_𝒌`$, and allowing finite dissipation but ignoring wave dynamics (so that $`\nu _k`$ and the various $`\eta `$ coefficients are real), one can show that Eq. (113) simplifies in this disparate scale limit to

$`\eta _{C𝒌}\eta _𝒌\nu _𝒌{\displaystyle \frac{2\eta _𝒌(\eta _𝒌\nu _𝒌)(\eta _𝒌\eta _{C𝒌})}{(\eta _{C𝒌}+\eta _𝒌)^2}}.`$ (116)
This gives a cubic equation for $`\eta _{C𝒌}`$. For $`\nu _k=0`$ the roots are $`\eta _{C𝒌}=\eta _𝒌`$ and $`\eta _{C𝒌}=(1\pm \sqrt{2})\eta _𝒌`$. Our speculation is that $`\eta _{C𝒌}=\eta _𝒌`$ will be the usual case in a steady-state inertial range. (This appears reasonable, but it might require numerical simulations to test it more definitively.) The other roots are probably unstable equilibria, so that any perturbation away from it would eventually approach the stable root, or may only be relevant in transient inverse-cascade cases where $`Re\eta _𝒌<0`$ ($`Re\eta _{C𝒌}0`$ being required to satisfy realizability).

Thus the two correction terms in Eq. (113) again exactly cancel each other (assuming the root choice made above), leading to $`\eta _{C𝒌}=\eta _𝒌`$ and the result that non-white-noise corrections are asymptotically unimportant in a wide inertial range ($`𝒌`$ large compared to the long-wavelength energy-containing wave number scale $`𝒌_\mathrm{𝟎}`$). However, this may be an artifact of the problem that the underlying DIA, on which the MRMC is based, does not satisfy random Galilean invariance. As is well known, the reason the DIA predicts a slightly different spectrum ($`E(k)k^{3/2}`$ in the energy cascade inertial range) than the Kolmogorov result ($`E(k)k^{5/3}`$) is because of this lost random Galilean invariance. \[The standard definitions for two-dimensional hydrodynamics use $`E(k)k^3C_{|𝒌|}`$ when $`\psi _𝒌`$ represents the stream function, so that the total energy is $`𝑑kE(k)`$, a one-dimensional integral over the magnitude of $`𝒌`$.\] The magnitude of this discrepancy between the DIA and dimensionally self-similar predictions is calculated for a general equation of the form Eq. (40) in Appendix C.

The underlying reason for this failure of the DIA is that the nonlinear damping and noise terms (the left- and right-hand sides of Eq. (110)) are dominated by contributions from the energy at long wavelengths. A random-Galilean invariant theory should depend only on the shear of longer wavelength modes (as $`\eta _𝒌`$ does in Orszag’s phenomenological EQDNM) and the most energetically significant interactions should occur among comparable scales ($`|𝒒||𝒑||𝒌|`$). Then the disparate scale approximations that led to Eq. (116) would no longer be valid. In such a case, it would seem unlikely that the two terms in Eq. (113) would still exactly cancel, and there would probably be some difference between the decorrelation rate $`\eta _{C𝒌}`$ and the decay rate $`\eta _𝒌`$. It would therefore be interesting to try to apply the techniques developed here (for allowing multiple rates) to other starting equations that respect random Galilean invariance, such as the Lagrangian-history DIA, test-field model, or renormalization group methods.

A regime where the correction terms might not cancel each other, and the differences between $`\eta _{C𝒌}`$ and $`\eta _𝒌`$ might be significant, even with the DIA’s overemphasis of long-wavelength contributions to the eddy turnover rate, is in ITG/drift-wave plasma turbulence, where the spectrum can often be anisotropic and have strong wave effects. That is, $`\nu _𝒌`$ can be complex, with unstable modes in some directions and damped modes in others, so that $`\nu _𝒌`$ and $`C_𝒌`$ vary strongly with the direction of $`𝒌`$. Some plasma cases have a reduced range of relevant nonlinearly interacting scales, and the simplifications of disparate scales in an inertial range used to derive Eq. (116) are not appropriate. The corrections might also be important in non-steady-state transient cases (such as zonal flows with predator–prey dynamics) or in other regimes where interactions between comparable scales dominate. Evaluating the difference between the decorrelation rate $`\eta _{C𝒌}`$ and the decay rate $`\eta _𝒌`$ in more general cases such as these probably requires a numerical treatment.

Finally, it is useful to demonstrate that the Multiple-Rate Markovian Closure approximation preserves realizability, which turns out to require one additional constraint. The MRMC equations (77) have the underlying Langevin equation

$$\frac{\psi _𝒌}{t}+\overline{\eta }_𝒌(t)\psi _𝒌(t)=f_𝒌^{}(t),$$
(117)
where $`\overline{\eta }_𝒌`$ is given by Eq. (80). The statistics that $`f_𝒌`$ must satisfy can be found by comparing the solution for such a Langevin equation, given by Eq. (10), with Eqs. (77), finding the constraint $`Re_0^t𝑑\overline{t}\overline{R}^{}(t,\overline{t})C_f^{}(t,\overline{t})=F_k`$, where $`F_k`$ is given by Eq. (81) and $`\overline{R}(t,\overline{t})=\mathrm{exp}(_{\overline{t}}^t𝑑t^{\prime \prime }\overline{\eta }_𝒌(t^{\prime \prime }))`$ is the propagator for Eq. (117). Using an integral form for $`\mathrm{\Theta }_{𝒌𝒑𝒒}`$ (similar to Eq. (57)),

$`\mathrm{\Theta }_{𝒌𝒑𝒒}(t)={\displaystyle _0^t}𝑑\overline{t}R_{\mathrm{mod},𝒌}(t,\overline{t}){\displaystyle \frac{C_{\mathrm{mod},𝒑}(t,\overline{t})C_{\mathrm{mod},𝒒}(t,\overline{t})}{C_𝒑^{1/2}(t)C_𝒒^{1/2}(t)}},`$

we find that if the two-time statistics of $`f_𝒌`$ satisfy

$`C_f(t,\overline{t})=\mathrm{exp}\left[{\displaystyle _{\overline{t}}^t}𝑑t^{\prime \prime }\left(\eta _k(t^{\prime \prime })\overline{\eta }_𝒌(t^{\prime \prime })\right)\right]`$
$`\times {\displaystyle \frac{1}{2}}{\displaystyle \underset{𝒌+𝒑+𝒒=\mathrm{𝟎}}{}}|M_{𝒌𝒑𝒒}|^2C_{\mathrm{mod},𝒑}(t,\overline{t})C_{\mathrm{mod},𝒒}(t,\overline{t})`$
(for $`t>\overline{t}`$), then the MRMC is the statistical solution of Eq. (117). As shown in Theorem 1 of Appendix B of BKO (and as can be inferred from considering the statistics of $`f(t)=g(t)h(t)`$, where $`g`$ and $`h`$ are statistically independent), a product of realizable correlation functions is also a realizable correlation function. $`C_{\mathrm{mod},𝒑}(t,t^{})`$ and $`C_{\mathrm{mod},𝒒}(t,t^{})`$ are individually realizable because $`Re\eta _{C𝒌}>0`$ for all $`𝒌`$. So in order to guarantee realizability of $`C_f(t,\overline{t})`$, we need to impose the additional condition that $`Re\eta _𝒌Re\overline{\eta }_𝒌`$. This constraint seems physically reasonable. The parameter $`\eta _𝒌`$ measures the decay rate for the ensemble averaged response $`\delta \psi _𝒌(t)`$, which can decay either as energy is nonlinear transferred out of mode $`𝒌`$ or as the energy that is in $`\delta \psi _𝒌`$ becomes randomly phased. The quantity $`\overline{\eta }_𝒌`$ used in Eq. (77) measures only the rate at which net energy (regardless of phase) is transferred out of mode $`𝒌`$ into other modes, so it would seem reasonable that $`\overline{\eta }_𝒌\eta _𝒌`$ will naturally result.

## VIII Conclusions

In summary, we have demonstrated a method for extending Markovian approximations of the DIA, to allow the decorrelation rate for fluctuations to differ from the decay rate for the infinitesimal response function (the renormalized Green’s function or nonlinear propagator). This can give a more accurate treatment of various effects such as non-white-noise forcing terms. In practice, the corrections to the decorrelation rate are modest, at least in isotropic non-wave cases, since the decorrelation rate of the noise is usually comparable to, if not much larger than, the decay rate for the response function. For example, if $`\eta _f=\eta `$ in the simple Langevin example of Eq. (24), then the decorrelation rate is $`60`$% lower than its white-noise value. Furthermore, the Multiple-Rate Markovian Closure Eq. (113) for the full DIA contains an offsetting term that can increase $`\eta _{C𝒌}`$, so the net result is less clear. This is because the DIA is related to a generalized Langevin equation Eq. (49), where non-white fluctuations modify not only the noise term (which tends to reduce the decorrelation rate) but also modify the renormalized propagator in the time-history integral (which tends to increase the decorrelation rate).

We have demonstrated that these two terms in fact exactly cancel each other as they should in thermal equilibrium where the fluctuation–dissipation theorem applies. We have also found another case, that of a wide inertial range with no waves, where it is possible for these two corrections to offset each other exactly, so that the decorrelation rate and the decay rates become equal. However, this may be an artifact of the loss of Galilean invariance in the Eulerian DIA, where modes in the inertial range nonlinearly interact predominantly with long wavelength modes. Thus it would be interesting to try to apply the techniques developed in this paper to other renormalized statistical theories, in which the dominant nonlinear interactions in an inertial range are between comparable scales instead of disparate scales and which properly reproduce Kolmogorov’s $`E(k)k^{5/3}`$ inertial-range energy spectrum instead of the Eulerian DIA’s $`E(k)k^{3/2}`$. Single-rate Markovian approximations have been applied in the past to other renormalized statistical theories and white-noise assumptions have also been employed in renormalization group calculations of turbulence. An interesting question is whether there is some way to generalize such calculations to allow for multiple rates and non-white noise as considered here. Another question is whether multiple-rate extensions might modify subgrid turbulence models. \[Such corrections would probably be important only at short scales near the transition from resolved to unresolved scales.\]

Even in the context of an Eulerian DIA-based theory, there may be some regimes where the multiple-rate corrections in this paper may be important and warrant further investigation. These might include cases where non-steady-state dynamics are important (i.e., predator–prey oscillations between different parts of the spectrum, such as between drift waves and zonal flows), or where interactions between comparable $`|𝒌|`$ scales are more important, such as might occur in anisotropic plasma turbulence with wave dynamics and with instability growth rates or Landau damping rates that vary strongly with the magnitude and direction of the wavenumber. One could test whether these corrections are important or negligible in various regimes by looking at 3-mode coupling cases, or by numerically comparing with the DIA or direct numerical simulations.

The complete set of equations that constitutes the Multiple-Rate Markovian Closure (MRMC) are summarized in the final paragraph of Sec. (VI). The MRMC extends the Realizable Markovian Closure (RMC) of BKO to allow various nonlinear rates and interaction times to differ. The single decay/decorrelation rate of the RMC is replaced with 3 different rates, $`\eta _𝒌`$ (the response function decay rate), $`\eta _{C𝒌}`$ (the decorrelation rate for the two-time correlation function), and $`\overline{\eta }_𝒌`$ (the energy damping rate). The triad interaction time of the RMC is replaced with 4 different triad interaction times with various weightings of decorrelation and decay rates. While this increases the complexity of the equations somewhat, the main computational advantages of a local-in-time Markovian closure relative to the non-local-in-time DIA are retained.

###### Acknowledgements.

We thank Prof. John A. Krommes for many helpful discussions sharing his insights into the DIA and Markovian approximations. In particular, we thank him for pointing out how to get the symmetric form of Eq. (21), which helps preserve important invariance properties. This work was supported by U.S. Department of Energy Contract No. DE–AC02–76CHO3073 and by the Natural Sciences and Engineering Research Council of Canada.

## A Realizability of a particular two-point correlation function

In theorem 2 of their Appendix B, Bowman, Krommes, and Ottaviani show one way to prove that a two-point correlation function of the form of Eq. (29) is “realizable” (if $`Re\eta _C(t)>0`$ is satisfied almost everywhere). Realizability means that this two-point correlation function is the exact solution to some underlying stochastic problem, such as a Langevin equation. In the absence of realizability, non-physical difficulties can sometimes develop, such as the predicted energy $`C(t)=C(t,t)`$ going negative or diverging. Here we present an alternate proof that Eq. (29) is realizable.

Consider the standard Langevin equation with time-dependent coefficients, but in the white-noise limit $`f(t)f^{}(t^{})=2D(t)\delta (tt^{})`$. Then Eq. (10) simplifies to

$$\frac{C(t)}{t}+2Re\eta (t)C(t)=2D(t),$$
(A1)
while the equation for the two-time correlation function, Eq. (9), becomes just $`C(t,t^{})/t+\eta (t)C(t,t^{})=0`$ for $`t>t^{}`$, with the boundary condition $`C(t^{},t^{})=C(t^{})`$. Taking the time derivative of Eq. (29) gives

$$\left(\frac{}{t}+\eta _C(t)\right)C_{\mathrm{mod}}(t,t^{})=\frac{1}{2C(t)}\frac{C(t)}{t}C_{\mathrm{mod}}(t,t^{}).$$
(A2)
If $`C(t,t^{})=C_{\mathrm{mod}}(t,t^{})`$, then these last two equations give $`\eta =\eta _C(C(t)/t)/(2C(t))`$. Using Eq. (A1), this becomes $`\eta _C(t)=\eta (t)Re\eta (t)+D(t)/C(t)`$. It is interesting to note that this ensures $`Re\eta _C0`$ even if $`Re\eta <0`$. These equations can be rearranged to give $`Im\eta (t)=Im\eta _C(t)`$, $`D(t)=C(t)Re\eta _C(t)`$, and $`Re\eta (t)=Re\eta _C(t)(\mathrm{log}C(t)/t)/2`$. Thus, given any 3 arbitrary functions $`C(t)0`$, $`Re\eta _C(t)`$, and $`Im\eta _C(t)`$ that determine the model Eq. (29), it is possible to find a white-noise Langevin equation for which it is the exact solution (as long as $`Re\eta _C0`$ so that $`D0`$). Conversely, for any arbitrary complex $`\eta (t)`$ and real $`D(t)0`$ that specify a white-noise Langevin problem, one can find a corresponding solution of the form Eq. (29).

It is interesting to note that $`C(t,t^{})=`$ $`C(t^{})\mathrm{exp}[_t^{}^t𝑑t^{\prime \prime }\eta (t^{\prime \prime })]`$ (for $`t>t^{}`$) is also an exact solution for this same white-noise Langevin problem. This form is valid for arbitrary $`\eta (t)`$ (even $`Re\eta <0`$). However, BKO and references therein indicate that this fails to preserve realizability when used in the context of Markovian approximations to the DIA, so they instead use Eq. (29).

On a related topic, BKO showed that their realizable Markovian closure (RMC), as given by their Eqs. (66a-e), has an underlying Langevin representation given by their Eq. (67) with a two-time noise correlation function $`f_k(t)f_k^{}(t^{})`$ of the form of their Eq. (64), which is not necessarily white noise. However, other two-time noise correlation functions can also give the same equal-time statistics equivalent to their Eq. (66a). This requires $`F_k(t)=Re_0^t𝑑\overline{t}f_k(t)f_k^{}(\overline{t})R_k^{}(t,\overline{t})`$, where $`F_k(t)`$ is the noise term in their Eq. (66a). For a case where $`F_k(t)`$ is always positive, then the RMC is also equivalent to a Langevin representation with white-noise, $`f_k(t)f_k^{}(t^{})=2D_f(t)\delta (tt^{})`$, where $`D_f(t)`$ $`=\frac{1}{2}Re_{𝒌+𝒑+𝒒=0}|M_{𝒌𝒑𝒒}|^2\mathrm{\Theta }_{𝒌𝒑𝒒}C_p^{1/2}C_q^{1/2}`$. While both white and non-white noise can give the same equal-time equations for $`C(t)`$, they will give different results for the two-time correlation function $`C(t,t^{})`$. However, there can be cases where $`F_𝒌(t)<0`$, for which a realizable Langevin representation must use non-white noise, as in their Eq. (64). \[Note that while $`C(t)0`$ is a fundamental requirement preserved by a realizable theory, the “triad interaction time” $`Re\mathrm{\Theta }_{𝒌𝒑𝒒}`$ may go negative. An example, similar to Eq. (47) of BKO, can be constructed for the realizable $`\mathrm{\Theta }_{𝒌𝒑𝒒}`$ of Eq. (66d) of BKO in the limit of constant $`C_𝒑`$ and $`C_𝒒`$ with $`\eta _𝒌=\eta _𝒑+\eta _𝒒=\rho +ia`$.\]

## B Fitting models to the two-time correlation function

Conceptually the process of fitting an exponential model of decorrelation to the actual two-time correlation function seems straightforward. But as described in Sec. (III A) and Sec. (V), there are various choices one could make in the weights used to fit the models. Galilean invariance imposes some constraints, but does not completely constrain the problem. In this appendix we further describe some options and our choices.

Consider the following measure of the error between the actual two-time correlation function and a model correlation function:

$$S(t)=_0^t𝑑t^{}|C(t,t^{})C_{\mathrm{mod}}(t,t^{})|^2.$$
(B1)
We will assume $`C_{\mathrm{mod}}(t,^{}t)`$ is of the form of Eq. (29). The equal time correlation function $`C(t)=C(t,t)`$ is already specified, so our task is to choose $`\eta _C(t)`$ in Eq. (29) in such a way as to minimize the squared error $`S`$. We want to stay in a Markovian framework, where $`\eta _C(t)`$ depends on parameters only from the present time. We assume that $`\eta _C(t^{})`$ for times $`t^{}<t`$ has already been chosen optimally. But we can choose $`\eta _C(t)`$ at the present time so that the extrapolation of $`S(t)`$ into the future is minimized. That is, we want to minimize $`S/t`$, which, after using Eq. (9) for $`C(t,t^{})/t`$ and Eq. (29) to evaluate $`C_{\mathrm{mod}}(t,t^{})/t`$, is

$`{\displaystyle \frac{S}{t}}`$ $`=`$ $`2{\displaystyle _0^t}dt^{}[\eta (t)C(t,t^{})+{\displaystyle _0^t^{}}d\overline{t}R^{}(t^{},\overline{t})C_f^{}(t,\overline{t})`$ (B4)
$`{\displaystyle \frac{1}{2C(t)}}{\displaystyle \frac{C(t)}{t}}C_{\mathrm{mod}}(t,t^{})+\eta _C(t)C_{\mathrm{mod}}(t,t^{})]`$
$`\times (C^{}(t,t^{})C_{\mathrm{mod}}^{}(t,t^{}))+\mathrm{c}.\mathrm{c}.,`$
where $`\mathrm{c}.\mathrm{c}.`$ indicates the complex conjugate of the previous expression. Separately minimizing $`S/t`$ with respect to the real part $`\eta _{Cr}`$ and imaginary part $`\eta _{Ci}`$ of $`\eta _C(t)`$ (i.e., set $`(S/t)/\eta _{Cr}=0`$, and then $`(S/t)/\eta _{Ci}=0`$) leads to the requirement that

$$_0^t𝑑t^{}C_{\mathrm{mod}}^{}(t,t^{})C_{\mathrm{mod}}(t,t^{})=_0^t𝑑t^{}C_{\mathrm{mod}}^{}(t,t^{})C(t,t^{}).$$
(B5)
\[Note that when evaluating derivatives of Eq. (B4) with respect to $`\eta _{Cr}`$ and $`\eta _{Ci}`$, it is only the explicit appearance of $`\eta _C(t)`$ in Eq. (B4) that is important. The parameter $`\eta _C(t)`$ also appears implicitly via the definition of $`C_{\mathrm{mod}}(t,t^{})`$, but there it has an impact on the integral defining $`S/t`$ only through a set of measure zero, and so can be neglected as long as $`\eta _C(t)`$ is bounded.\]

In the steady-state limit, Eq. (B5) is equivalent to Eq. (21). For a time-dependent case, consider Eq. (B5) as providing a constraint of the form $`A_{\mathrm{mod}}(t)=A(t)`$. Assuming that this has already been satisfied for earlier times, we want it to remain satisfied for future times, i.e., we need to require that $`A_{\mathrm{mod}}/t=A/t`$. This is precisely what we are doing when we set Eq. (31) and Eq. (33) to be equal, and it leads to a formula for $`\eta _C(t)`$ at the present time that minimizes the errors as time advances.

The same procedures as described here are used in fitting a model response function $`R_{\mathrm{mod}}(t,t^{})`$ of the form of Eq. (77) to the actual response function, leading to the constraint

$$_0^t𝑑t^{}R_{\mathrm{mod}}^{}(t,t^{})R_{\mathrm{mod}}(t,t^{})=_0^t𝑑t^{}R_{\mathrm{mod}}^{}(t,t^{})R(t,t^{})$$
(B6)
As mentioned at the end of Sec. (V), $`R_{\mathrm{mod}}^{}`$ in this expression could be replaced with $`C_{\mathrm{mod}}^{}`$ and one would still get an expression defining $`\eta `$ that was Galilean invariant. However, Eq. (B6) seems to make more sense as a least-squares best fit of $`R_{\mathrm{mod}}`$ to $`R`$, and that is the choice we have made.

But consider Eq. (B1) in the steady-state limit where $`\eta _C`$ is a constant and $`C(t,t^{})`$ depends only on $`tt^{}`$,

$$S_0=_{\mathrm{}}^t𝑑t^{}|C(t,t^{})C_0e^{\eta _C(tt^{})}|^2.$$
(B7)
It is straightforward to show that choosing $`\eta _C`$ to minimize the total squared error $`S_0`$ leads to the condition

$`{\displaystyle _{\mathrm{}}^t}𝑑t^{}e^{\eta _C^{}(tt^{})}C_0e^{\eta _C(tt^{})}(tt^{})`$ (B8)
$`=`$ $`{\displaystyle _0^t}𝑑t^{}e^{\eta _C^{}(tt^{})}C(t,t^{})(tt^{}).`$ (B9)
Note that this differs from Eq. (B5) by an additional factor of $`(tt^{})`$, which weights errors at larger time separation more strongly. Including an extra weighting factor of $`(tt^{})`$ in Eq. (21) might help to refine the model, particularly for cases such as in Fig. (5), where the short time behavior is reasonable but the long-time fit needs improvement.

It is perhaps not surprising that optimizing a constant $`\eta _C`$ to minimize the global error $`S_0`$ gives a somewhat different result than optimizing $`\eta _C(t)`$ to minimize the local error $`S/t`$. In order for the time-dependent fitting procedures to reproduce this steady-state result, one could modify Eq. (B1) by multiplying the integrand by a factor of $`(tt^{})`$. Working through the derivation, one finds that the integrands in Eq. (B5) would be modified to also have an additional weighting factor of $`(tt^{})`$. Thus one might be able to improve the results in this paper some by including an extra weighting of $`(tt^{})`$ in the appropriate places, Eq. (30), Eq. (70), Eq. (84), and Eq. (92), and working through the derivations to see the modified results. While such modifications could lead to an improved model, and would be interesting for future work, one should realize that the dynamics are complicated and no choice of weights is perfect. For example, what one really wants is a best fit model for the triad interaction times which are weighted by interactions between three modes as given in Eq. (57), not necessarily best fits for the decorrelation rates of just individual modes. Probably a higher priority for future work is to use a starting set of equations that satisfy random Galilean invariance, so that interactions with large scales are not overemphasized as they are in the Eulerian DIA.

## C Inertial-range scaling of DIA-based closures

Here we determine steady-state self-similar inertial-range solutions in $`d`$ dimensions to closures of the form (IV B) in an unbounded domain (so that $`_{𝒌+𝒑+𝒒=\mathrm{𝟎}}_{\mathrm{\Delta }_𝒌}𝑑𝒑𝑑𝒒𝑑𝒑𝑑𝒒\delta (𝒌+𝒑+𝒒)`$), taking the initial time $`t_0=\mathrm{}`$. This extends previous derivations in the literature to self-similar spectra consistent with generic DIA-based closures (IV B) of the quadratically nonlinear equation (40), arising from the cascade of a generalized energy $`\frac{1}{2}_𝒌\sigma _𝒌\left|\psi _𝒌(t)\right|^2`$. Assuming self-similar scalings of the mode-coupling and statistical variables, our derivation requires only the additional condition (C10), which is somewhat weaker than statistical stationarity.

The turbulence could be forced with a linear instability, incorporated with dissipation into the linear coefficient $`\nu _𝒌`$, or else a random force could be added to the right-hand side of Eq. (54). By definition, both the external forcing and dissipation $`\nu _𝒌`$ vanish in the inertial range. The symmetry (42) then implies that the nonlinear terms in Eq. (54), weighted by $`\sigma _𝒌`$, must balance. It is convenient to define

$`S_𝒌\sigma _𝒌Re(F_𝒌N_𝒌)`$
$`=`$ $`{\displaystyle \frac{1}{2}}Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒\sigma _𝒌M_{𝒌𝒑𝒒}M_{𝒌𝒑𝒒}^{}\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}`$
$`+Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒\sigma _𝒌M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\overline{\mathrm{\Theta }}_{𝒑𝒒𝒌}^{}`$
$`=`$ $`{\displaystyle \frac{1}{2}}Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒M_{𝒌𝒑𝒒}(\sigma _𝒑M_{𝒑𝒒𝒌}^{}+\sigma _𝒒M_{𝒒𝒌𝒑}^{})\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}`$
$`+Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒\sigma _𝒌M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\overline{\mathrm{\Theta }}_{𝒑𝒒𝒌}^{}`$
$`=`$ $`Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒\sigma _𝒑M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}`$
$`+Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒\sigma _𝒌M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}\overline{\mathrm{\Theta }}_{𝒑𝒒𝒌}^{}`$
$`=`$ $`Re{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}(\sigma _𝒌\overline{\mathrm{\Theta }}_{𝒑𝒒𝒌}^{}\sigma _𝒑\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}).`$
We seek self-similar solutions of the DIA that obey the scalings (for $`\lambda >0`$)

$$M_{\lambda 𝒌,\lambda 𝒑,\lambda 𝒒}=\lambda ^mM_{𝒌𝒑𝒒},$$
(C2)
$$\sigma _{\lambda 𝒌}=\lambda ^s\sigma _𝒌,$$
(C3)
$$R_{\lambda 𝒌}(t,t^{})=R_𝒌(t,t\lambda ^{\mathrm{}}(tt^{})),$$
(C4)
$$C_{\lambda 𝒌}(t,t^{})=\lambda ^nC_𝒌(t,t\lambda ^{\mathrm{}}(tt^{})),$$
(C5)
so that, upon making the change of variables $`\overline{s}t\lambda ^{\mathrm{}}(t\overline{t})`$ in Eq. (57),

$$\overline{\mathrm{\Theta }}_{\lambda 𝒌,\lambda 𝒑,\lambda 𝒒}=\lambda ^{\mathrm{}+2n}\overline{\mathrm{\Theta }}_{𝒌𝒑𝒒}.$$
(C6)
Once we have determined suitable values of the scaling exponent $`n`$, we may compute the wavenumber exponent $`\beta `$ for the energy spectrum $`E(k)ϵ^\alpha k^\beta `$. If the total energy $`E`$ is related to the correlation function $`C_𝒌`$ of the fundamental variable $`\psi `$ by $`E=𝑑𝒌k^\gamma C_𝒌=𝑑kE(k)`$, then $`\beta =d1+\gamma +n`$.

Following Ref , we will use the change of variables $`z=k^2/p`$, $`w=kq/p`$ to determine values of the exponents $`\mathrm{}`$ and $`n`$ for which the angular average $`S(k)`$ of $`S_𝒌`$ vanishes. In terms of the scaling factor $`\lambda =k/z`$ we note that $`k=\lambda z`$, $`p=\lambda k`$, and $`q=\lambda w`$. Letting $`𝒛=z\widehat{𝒑}`$ and $`𝒘=w\widehat{𝒒}`$, we may then express $`d𝒑d𝒒=\lambda ^{3d}d𝒛d𝒘`$ and $`\delta (𝒌+𝒑+𝒒)=\lambda ^d\delta (z\widehat{𝒌}+k\widehat{𝒑}+𝒘)`$. Hence, upon interchanging $`\widehat{𝒑}`$ and $`\widehat{𝒌}`$ in the integration, we deduce

$`S(k)`$ $``$ $`{\displaystyle 𝑑\widehat{𝒌}S_𝒌}=Re{\displaystyle 𝑑\widehat{𝒌}_{\mathrm{\Delta }_𝒌}𝑑𝒛𝑑𝒘\lambda ^{3dd+2m+s+\mathrm{}+2n}}`$
$`\times M_{𝒛,𝒌,𝒘}M_{𝒌,𝒘,𝒛}^{}(\sigma _𝒛\overline{\mathrm{\Theta }}_{𝒌,𝒘,𝒛}^{}\sigma _𝒌\overline{\mathrm{\Theta }}_{𝒛,𝒌,𝒘})`$
$`=`$ $`Re{\displaystyle 𝑑\widehat{𝒌}_{\mathrm{\Delta }_𝒌}𝑑𝒛𝑑𝒘\lambda ^{2d+2m+s+\mathrm{}+2n}}`$
$`\times M_{𝒌,𝒛,𝒘}^{}M_{𝒛,𝒘,𝒌}(\sigma _𝒌\overline{\mathrm{\Theta }}_{𝒛,𝒘,𝒌}\sigma _𝒛\overline{\mathrm{\Theta }}_{𝒌,𝒛,𝒘}^{})`$
$`=`$ $`S(k),`$
provided that

$$2d+2m+s+\mathrm{}+2n=0.$$
(C7)
The condition (C7) guarantees that the angle-averaged nonlinear terms in Eq. (54) will balance in a steady state and lead to an inertial range.

The exponent $`\mathrm{}`$ can be determined by integrating the DIA response function equation

$`{\displaystyle \frac{}{t}}R_𝒌(t,t^{}){\displaystyle _{\mathrm{}}^t}𝑑\overline{t}{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}`$ (C8)
$`\times R_𝒑^{}(t,\overline{t})C_𝒒^{}(t,\overline{t})R_𝒌(\overline{t},t^{})=\delta (tt^{}),`$ (C9)
over all $`t^{}`$, using the steady-state condition

$$\underset{t\mathrm{}}{lim}\frac{}{t}_{\mathrm{}}^{\mathrm{}}𝑑t^{}R(t,t^{})=0.$$
(C10)
One obtains

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\overline{t}{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}`$ (C11)
$`\times R_𝒑^{}(\mathrm{},\overline{t})C_𝒒^{}(\mathrm{},\overline{t}){\displaystyle _{\mathrm{}}^{\mathrm{}}}dt^{}R_𝒌(\overline{t},t^{})=1.`$ (C12)
Upon replacing $`𝒌`$ by $`\lambda 𝒌`$ (for any constant $`\lambda `$) and exploiting the self-similar scalings given in Eqs. (C), we make the change of variable $`s^{}=\overline{t}\lambda ^{\mathrm{}}(\overline{t}t^{})`$ to obtain

$`\lambda ^{d+2m+\mathrm{}+n}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\overline{t}{\displaystyle _{\mathrm{\Delta }_𝒌}}𝑑𝒑𝑑𝒒M_{𝒌𝒑𝒒}M_{𝒑𝒒𝒌}^{}`$
$`\times R_𝒑^{}(\mathrm{},\overline{s})C_𝒒^{}(\mathrm{},\overline{s}){\displaystyle _{\mathrm{}}^{\mathrm{}}}ds^{}R_𝒌(\overline{t},s^{})=1,`$
where $`\overline{s}t\lambda ^{\mathrm{}}(t\overline{t})`$. The integral over $`\overline{t}`$ is dominated by contributions from large $`\overline{t}`$, for which the integral over $`s^{}`$ asymptotically approaches a constant (with respect to $`\overline{t}`$), according to Eq. (C10). Hence, after making a final change of variables from $`\overline{t}`$ to $`\overline{s}`$, we see that the balance expressed in Eq. (C12) is recovered if

$$\lambda ^{d+2m+2\mathrm{}+n}=1,$$
(C13)
from which we conclude that $`\mathrm{}=(d+n)/2m`$. If one inserts this result into Eq. (C7), one obtains the Kolmogorov scalings

$$\mathrm{}=\frac{1}{3}s\frac{2}{3}m,$$
(C15)
$$n=d\frac{2}{3}(m+s),$$
(C16)
$$\beta =\gamma 1\frac{2}{3}(m+s).$$
(C17)
Alternatively, one could adopt instead of Eq. (C10) the stronger condition of statistical stationarity, $`R_𝒌(t,t^{})=_𝒌(tt^{})`$ and $`C_𝒌(t,t^{})=𝒞_𝒌(tt^{})`$. Equation (C13) is then readily seen to follow directly from Eq. (C9). In either case we have only shown that Eq. (C13) is a necessary condition for self-similar solutions of the form (C) to exist. In order that these solutions actually satisfy Eq. (C9), it is also necessary at the very least that the wavenumber integral in Eq. (C9) converges.

Unfortunately, the scaling expressed in Eq. (C13) often leads to a divergence of the $`q`$ integral in Eq. (C9), preventing self-similar solutions from existing. Typically, the mode-coupling coefficients $`M_{𝒌,𝒌𝒒,𝒒}`$ asymptotically approach a constant as $`q`$ goes to zero while $`𝒌`$ is held fixed. Upon performing the $`𝒑`$ integration in Eq. (C9), we then see that the $`q`$ integrand will scale like $`q^{d1}C_𝒒^{}(t,\overline{t})`$ for small $`q`$. If $`C_𝒒`$ asymptotically scales as $`q^n`$, then the integrand will scale like $`q^{d1+n}`$. But Eq. (C16) implies that $`d1+n=12(m+s)/3`$. Normally $`m+s>0`$ (see Table I); in these cases there would be a divergence of the $`q`$ integral in Eq. (C9) if self-similar solutions really were to exist.

This divergence indicates that the dominant contributions to the eddy-turnover time come from the energy spectrum at large scales, where self-similarity no longer holds. (For this reason, the DIA is not invariant to random Galilean transformations.) The actual value of the scaling $`\mathrm{}`$ that appears in the DIA response must be calculated by taking into account that $`C_𝒒`$ does not actually behave as $`q^n`$ for small $`q`$. The DIA equations apply to the case of zero mean flow, where the energy spectrum goes to zero at small wavenumbers. This means that the integration in Eq. (C9) must be effectively cut off at some fixed large scale wavenumber $`k_0`$. The introduction of this cutoff wavenumber removes the divergence in the integral, but it also changes the above scaling argument. Since the dominant contribution to Eq. (C9) still comes from small $`q`$, we need to identify the scaling of the mode-coupling coefficients with $`k`$ for $`qk`$, $`M_{\lambda 𝒌,\lambda 𝒌,\lambda 𝒒}=\lambda ^m^{}M_{𝒌,𝒌,𝒒}.`$ Since the lower wavenumber limit is now fixed, no self-similar scaling in $`𝒒`$ can be made; the scaling with $`k`$ for small $`q`$ then leads to $`\lambda ^{2\mathrm{}+2m^{}}=1`$. Hence for the DIA equations the actual scalings of the response function, correlation function, and energy spectrum are given by

$$\mathrm{}_{\mathrm{DIA}}=m^{},$$
(C19)
$$n_{\mathrm{DIA}}=dm+\frac{m^{}s}{2},$$
(C20)
$$\beta _{DIA}=\gamma 1m+\frac{m^{}s}{2}.$$
(C21)
In Table I we compare the scalings in Eqs. (C13) with the anomalous DIA scalings given by Eq. (C). The scalings given by Eq. (C13) are consistent with Kolmogorov’s dimensional analysis. We emphasize that these scalings would have also been obtained for the DIA equations (they too are dimensionally consistent) had the wavenumber integral in Eq. (C12) converged.