1. Introduction
A common problem in fluid dynamics, and indeed nonlinear physics in general, is to obtain statistical descriptions for quantities whose master dynamical equations are known. One approach is to accumulate statistics from numerical solutions of the master equations; for example, one may solve the Navier–Stokes equations for a turbulent flow and gather statistics via, for example, time averaging, ensemble averaging or averaging over a spatial coordinate. An alternative approach, termed direct statistical simulation (DSS) seeks solutions of closed-form equations for the statistics themselves that are derived from the dynamical equations; these may be equations that describe the probability density function of the system or the low-order statistics thereof (Allawala & Marston Reference Allawala and Marston2016). In many contexts, particularly in fluid dynamics, these statistical approaches are more favourable computationally; they can alleviate the stringent requirements of spatial resolution and sample sizes. The underlying expectation in this approach, however, is that solutions of the statistical equations will match those obtained from solutions of the dynamical equation.
Here, we consider the specific case of quasi-linear theory, a type of mean-field approximation that has been widely applied to a range of physical systems in fluids (and indeed plasmas) (Malkus Reference Malkus1954; Vedenov, Velikhov & Sagdeev Reference Vedenov, Velikhov and Sagdeev1961; Herring Reference Herring1963). In recent years, quasi-linear model (QL) equations have been used to simulate zonal jets – strong east–west mean flows – found in many astrophysical and geophysical situations (Farrell & Ioannou Reference Farrell and Ioannou2007; Srinivasan & Young Reference Srinivasan and Young2012; Galperin & Read Reference Galperin and Read2019) as well as to simulate flows in rotating dynamos (Calkins et al. Reference Calkins, Julien, Tobias and Aurnou2015), plasmas (Parker & Krommes Reference Parker and Krommes2014), wall-bounded shear flows (Farrell et al. Reference Farrell, Ioannou, Jimenez, Constantinou, Lozano-Duran and Nikolaidis2016) and convection (O'Connor, Lecoanet & Anders Reference O'Connor, Lecoanet and Anders2021). The simplest non-trivial statistical closure – a cumulant expansion closed at second order, denoted CE2 – naturally follows the assumption of the QL by retaining only up to second-order cumulants, and has enabled efficient DSS (Marston, Conover & Schneider Reference Marston, Conover and Schneider2008; Marston Reference Marston2010; Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011; Marston, Qi & Tobias Reference Marston, Qi and Tobias2019). Also see Farrell & Ioannou (Reference Farrell and Ioannou2007, Reference Farrell and Ioannou2013) for the closely related closure stochastic structural stability theory, where stochastic forcing replaces the neglected cumulants.
It is often supposed (usually implicitly) that the solution of the statistical equations should mirror those of the dynamical system from which they are derived if no approximation or truncation is made. In this paper, we demonstrate that this is not always the case for QL dynamics in a range of fluid problems, solved via numerical computation. We associate the discrepancy between the statistical theory (CE2) and a given realisation of the QL dynamics with a ‘rank instability’ of the CE2 system that is not available to QL dynamics.
This paper is organised as follows. In § 2, the fluid models used as a test bed are described, and the numerical framework is set up. In § 3, results are presented for deterministic and stochastically forced systems, showing the potential discrepancy between the QL and CE2. Theory for the origin of the discrepancy, including an exactly solvable illustration, is given in § 4, and conclusions are drawn in § 5.
2. Set-up of the models and numerical framework
2.1. The fully nonlinear model
We conduct two-dimensional simulations of a rotating, incompressible fluid with velocity 
${\boldsymbol u}$ on a doubly periodic 
$\beta$-plane 
$[0,2{\rm \pi} ]^2$. The time evolution of the relative vorticity 
$\zeta \equiv \hat {z}\boldsymbol {\cdot } (\boldsymbol {\nabla } \times {\boldsymbol u})$ is given by
\begin{equation} \partial_t \zeta = \beta\,\partial_x \psi - \kappa\zeta + \nu\,\nabla^2 \zeta + J[\psi,\zeta] + F, \end{equation}
where the Jacobian is 
$J[\psi,\zeta ] \equiv \partial _x\psi \,\partial _y \zeta - \partial _x \zeta \,\partial _y \psi$, and the streamfunction 
$\psi$ is given by 
$\psi \equiv \nabla ^{-2}\zeta$. Gradients of rotation are included via the 
$\beta$ term, ‘bottom’ friction via the 
$\kappa$ term, and viscosity via the 
$\nu$ term. We adopt three different forcing models to enable the comparison of the QL and CE2 with different dynamics. The first two forcings considered are deterministic. In the first, the flow is made to relax to an unstable jet profile (Marston et al. Reference Marston, Conover and Schneider2008); see § 3.1.1. In the second, a steady, two-scale Kolmogorov-type forcing is used (Tobias & Marston Reference Tobias and Marston2017); see § 3.1.2. The third forcing that we consider is stochastic in nature (Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2016); see § 3.2.
2.2. The quasi-linear model
The QL of the fully nonlinear system is obtained by decomposing the variables into a mean part and a fluctuating part, i.e. we write
\begin{equation} \zeta = \bar{\zeta}+ \zeta^\prime, \quad \psi = \bar{\psi} + \psi^\prime, \end{equation}where an overbar indicates a zonal average given by
\begin{equation} \bar{f}(y)= \frac{1}{2 {\rm \pi}}\int_0^{2 {\rm \pi}} f(x,y) \, {{\rm d}\kern0.06em x}, \quad f^\prime(x,y) = f(x,y) - \bar{f}(y). \end{equation}It is then possible to derive the QL by considering the interactions between the means and the fluctuations so that
\begin{equation} \partial_t \bar{\zeta} ={-} \kappa \bar{\zeta} + \nu\,\partial_{yy} \bar{\zeta} + \overline{J[\psi^\prime,\zeta^\prime]} + \bar{F} \end{equation}and
\begin{align} \partial_t \zeta^\prime = \beta\,\partial_x \psi^\prime - \kappa\zeta^\prime + \nu\,\nabla^2 \zeta^\prime + J[\bar{\psi},\zeta^\prime] + J[{\psi^\prime},\bar{\zeta}] + (J[\psi^\prime,\zeta^\prime] - \overline{J[\psi^\prime,\zeta^\prime]})+ F^\prime, \end{align}and then removing the fluctuation–fluctuation (or ‘eddy–eddy’) interactions (in round brackets) to yield the QL
$$\begin{gather} \partial_t \bar{\zeta} ={-} \kappa \bar{\zeta} + \nu\,\partial_{yy} \bar{\zeta} + \overline{J[\psi^\prime,\zeta^\prime]} + \bar{F}, \end{gather}$$
$$\begin{gather}\partial_t \zeta^\prime = \beta\,\partial_x \psi^\prime - \kappa\zeta^\prime + \nu\,\nabla^2 \zeta^\prime + J[\bar{\psi},\zeta^\prime] + J[{\psi^\prime},\bar{\zeta}] + F^\prime. \end{gather}$$It is this system that will be compared with the cumulant expansion (CE2) described below.
2.3. The cumulant expansion at second order (CE2)
The CE2 system is derived by first defining the equal-time cumulants. The first cumulant for the vorticity is given by 
$c_{\zeta }(y) = \bar {\zeta }$, whilst the second cumulants are defined by
\begin{equation} c_{\zeta\zeta}(\xi,y_1,y_2) = \overline{\zeta^\prime(x_1,y_1)\,\zeta^\prime(x_1-\xi,y_2)}, \end{equation}
where 
$\xi = x_1 - x_2$ as the system is translationally invariant in 
$x$. Here, 
$c_{\zeta \zeta }$ can be directly related to 
$c_{\zeta \psi }$ and 
$c_{\psi \zeta }$ by differentiation; for example, 
$c_{\zeta \psi } = \nabla ^2_1 c_{\psi \psi }$, with the differential operators 
$\nabla ^2_1 = \partial _\xi ^2 + \partial _{y_1}^2$ and 
$\nabla ^2_2 = \partial _\xi ^2 + \partial _{y_2}^2$ (for more detail, see Tobias et al. Reference Tobias, Dagon and Marston2011). We can then proceed by writing the equations in terms of 
$\zeta$ cumulants as follows.
With these definitions, the dynamical equations for the first cumulant 
$c_{\zeta }$ may be found immediately from (2.6):
\begin{equation} \dfrac{\partial c_{\zeta}}{\partial t} ={-}\kappa c_{\zeta} +\nu\,\dfrac{\partial^2 c_{\zeta}}{\partial y_1^2} - \left(\dfrac{\partial}{\partial y_1} + \dfrac{\partial}{\partial y_2}\right) \left.\dfrac{\partial c_{\psi \zeta}}{\partial \xi}\right|_{\xi \to 0}^{y_1 \to y_2} + \bar{F}. \end{equation}
Likewise, the equation of motion for the second cumulant 
$c_{\zeta \zeta }$ is obtained by multiplying (2.7) by 
$\zeta$ followed by a zonal average:
\begin{align} \dfrac{\partial c_{\zeta \zeta}}{\partial t} &= \dfrac{\partial c_{\psi}}{\partial y_1}\,\dfrac{\partial c_{\zeta \zeta}}{\partial \xi} - \left(\dfrac{\partial c_{\zeta}}{\partial y_1}-\beta \right) \dfrac{\partial c_{\psi \zeta}}{\partial \xi} - \dfrac{\partial c_{\psi}}{\partial y_2}\,\dfrac{\partial c_{\zeta \zeta}}{\partial \xi} + \left(\dfrac{\partial c_{\zeta}}{\partial y_2}-\beta \right) \dfrac{\partial c_{\zeta \psi}}{\partial \xi}\nonumber\\ &\quad +\nu (\nabla_1^2 + \nabla_2^2) c_{\zeta \zeta} - \kappa c_{\zeta \zeta} + \varGamma. \end{align}
Here, 
$\varGamma$ is given by the covariance of the forcing, i.e. 
$\varGamma (\xi, y_1, y_2) = $ 
$\overline{F^\prime ({x_1},y_1)\,F^\prime (x_1-\xi,y_2)}$, with the zonal average over 
$x_1$ combined with a short time average over the fast stochastic fluctuations. The QL defined by (2.6) and (2.7) can be transformed into wavevector space in a basis of Fourier harmonics. Similarly, the CE2 system defined by (2.9) and (2.10) can be solved in Fourier space by time-evolving 
$\hat {c}_\zeta (n)$ and 
$\hat {c}_{\zeta \zeta } (m, n_1, n_2)$ for the first and second cumulants, respectively. Here, 
$m$ is the wavenumber in the zonal (
$x$) direction, and 
$n$, 
$n_1$ and 
$n_2$ are wavenumbers in the non-zonal (
$\kern0.09em y$) direction. For convenience in the following, we define the matrix 
$C^{(m)}$ to be the zonal decomposition of the second cumulant: 
$C^{(m)} = \hat {c}_{\zeta \zeta } (m, n_1, n_2) = $
$({1}/{(2{\rm \pi} )^3})\int _{0}^{2{\rm \pi} }\int _{0}^{2{\rm \pi} }\int _{0}^{2{\rm \pi} }c_{\zeta \zeta }(\xi,y_1,y_2) \,{\rm e}^{-{\rm i} m\xi - {\rm i} n_1y_1 - {\rm i}n_2y_2}\,{\rm d}\xi \,{\rm d}y_1\,{\rm d}y_2$ for 
$m \in [1,M]$, where 
$M$ is the spectral cutoff in the zonal wavenumber.
2.4. Forcing and numerical scheme
We first examine cases that demonstrate divergences between QL and CE2 in deterministic systems, before adopting a stochastic forcing model (Srinivasan & Young Reference Srinivasan and Young2012; Farrell & Ioannou Reference Farrell and Ioannou2013; Tobias & Marston Reference Tobias and Marston2013; Constantinou et al. Reference Constantinou, Farrell and Ioannou2016; Marston et al. Reference Marston, Qi and Tobias2019). This enables us to highlight that disagreements between QL and CE2 may also appear in this commonly utilised scenario.
The spectral solver ZonalFlow.jl written in the Julia programming language (Bezanson et al. Reference Bezanson, Edelman, Karpinski and Shah2017) and available online (Nivarti, Marston & Tobias Reference Nivarti, Marston and Tobias2021) is used to obtain QL and CE2 solutions of (2.1). Time-stepping algorithms are imported from the well-tested ecosystem of the DifferentialEquations.jl package (Rackauckas & Nie Reference Rackauckas and Nie2017); here, we use the explicit 
$5/4$ Runge–Kutta method of Dormand and Prince, with a fixed time step for the deterministic cases, and an SRIW1 method of order 
$1.5$ for the stochastic cases. For the purpose of comparing QL and CE2 solutions, the key is to use the same spatial resolution for each model. We use small numbers of spectral modes – here, 
$M \times N = 8 \times 8$, 
$16 \times 16$ and 
$12\times 20$ – for the Kolmogorov flow, pointjet and stochastic cases, respectively. This allows us to demonstrate clearly the differences between QL and CE2. We have verified that the solution behaviour remains qualitatively similar for higher resolutions.
3. Results
3.1. Deterministic forcing
3.1.1. Relaxation to a pointjet
We begin with the case of relaxation to a pointjet. Here, the forcing is chosen to be 
$F(y) = - ({\varXi }/{\tau }) \tanh [({{\rm \pi} - y})/{\Delta y}]$. The 
$\beta$-plane is equatorial, i.e. 
$\beta = 2\varOmega \cos (0^\circ ) = 4{\rm \pi}$, where the period of rotation is 
$\varOmega = 2 {\rm \pi}$. Following Marston et al. (Reference Marston, Conover and Schneider2008) and Allawala, Tobias & Marston (Reference Allawala, Tobias and Marston2020), the viscosity is 
$\nu = 0$, and the relaxation time scale 
$\tau$ is set equal to the friction time scale 
$\kappa ^{-1}$. The jet strength is 
$\varXi = \varOmega = 2 {\rm \pi}$, and the jet width is 
$\Delta y = 0.1$. For the results shown, 
$\tau = 20$ days, though we stress that the comparisons hold for a wider range of 
$\tau$.
Figure 1 shows both QL and CE2 solutions obtained for this case. Figures 1(a–c) show the time evolution of the energy in zonal modes 
$m$, given by 
$E_m = \sum _n \hat {\zeta }_{m,n}^*\hat {\zeta }_{m,n}/(m^2 + n^2)$, where 
$\hat {\zeta }_{m,n}$ is the relevant Fourier coefficient of vorticity 
$\zeta (x,y)$, which can be calculated in the QL using 
$\zeta = \bar {\zeta }+ \zeta ^\prime$. In CE2, 
$E_m$ can be calculated as 
$E_m = \sum _n |\hat {c}_\zeta (n)|^2/n^2$ for 
$m = 0$, and 
$E_m = \sum _n C^{(m)}(n,n) /(m^2 + n^2)$ for 
$m \ge 1$. The QL solution, initialised with a random noise of strength 
$10^{-6}$, predicts a fixed point (FP) with significant energy in the zonal mean (
$m = 0$) mode and energy in two further non-zero zonal wavenumbers (
$m = 1$ and 
$m = 4$). We construct two different initial conditions (ICs) for CE2; in the first, the cumulants are constructed using the QL IC (figure 1b). For the second IC, we initialise the CE2 evolution using the stable QL FP (figure 1c) (Marston et al. Reference Marston, Qi and Tobias2019). Evidently, for this case, CE2 solutions for both ICs agree with the QL FP as expected.
Figure 1. Energy in zonal wavenumbers 
$E_m$ for unit-rank initialisation in the pointjet case up to a spin-up of 
$1000$ days, for (a) QL, (b) CE2 initialised with the QL IC, and (c) CE2 initialised with the QL fixed point solution (QL FP). Rank of cumulant submatrices (d) 
$C^{(1)}$ and (e) 
$C^{(4)}$ as found in the two CE2 solutions.
Since both ICs are constructed from a specific QL solution, the cumulant submatrices 
$C^{(m)}$ have rank unity at initialisation for all 
$m \leq M$. Numerically, we calculate this rank as the number of eigenvalues of the second cumulant larger than a small cutoff 
$10^{-12}$ of the order of the initial power. Notably, the ranks of 
$C^{(m)}$ for 
$m = 1,4$ remain unity for both CE2 solutions for all 
$t \geq 0$, as shown in figures 1(d,e), respectively, indicating that the statistical solutions are consistent with those of a single realisation of the QL system for the pointjet solution. Thus this is an example where the QL and CE2 are consistent. This agreement persists even if the second cumulant is initialised to be full rank in CE2 by assigning a value 
$10^{-6}$ to diagonal entries in 
$C^{(m)}$ for all 
$m \leq M$. The first cumulant is initialised at zero. Figure 2(a) shows 
$E_m$, with the zonal mean (
$m = 0$) and first harmonic (
$m = 1$) rapidly attaining values close to the QL time mean solution. The second harmonic (
$m = 4$), however, goes through an initially prolonged period of decay, unlike the unity rank initialisation in figure 1(b). After approximately 
$t = 150$ days, the 
$m = 4$ mode finally drops to unit rank (the figure 2(b) inset shows evolution of the largest two eigenvalues for 
$C^{(4)}$). This prompts the QL time mean solution to become stable once again. Thus, despite an initial full rank perturbation, CE2 is able to ‘find’ agreement with the QL by virtue of the stability of the unit rank solution; the solution drops rank until its rank is unity, and thereafter the agreement between CE2 and the QL is ensured. We stress that all the DSS solutions remain realisable, with positive eigenvalues for the second cumulant.
Figure 2. (a) Energy in zonal wavenumbers 
$E_m$. (b) Evolution of the rank of cumulant submatrix 
$C^{(4)}$ with full-rank initialisation for the pointjet case (inset shows evolution of its two largest eigenvalues, 
$\lambda _0$ and 
$\lambda _1$). The stability of the unit rank solution is consistent with the agreement of the CE2 and QL (figure 1a).
3.1.2. Kolmogorov forcing
Now we give what we believe is the first example of disagreement. Figure 3 shows solutions for a case with Kolmogorov flow forcing 
$F(y) = {-}\cos y -8\cos 4y$, where we choose 
$\nu = 0.02$ and 
$\beta = \kappa = 0$ in (2.1), which yields non-trivial dynamics in the resulting system (Tobias & Marston Reference Tobias and Marston2017). The QL (figure 3a) predicts a long-time solution with non-zero mean and two harmonics, i.e. 
$m = 1$ (blue) and 
$m = 2$ (orange). As for the pointjet case, CE2 simulations are performed with unit rank initialisations of the second cumulant using the QL IC (figure 3b) and the QL endpoint (EP) solution (figure 3c). Here, CE2 finds a single harmonic 
$m = 1$ in both cases, which disagrees with the QL (figure 3(e) shows the rank of 
$C^{(2)}$ collapsing to zero). The disagreement is caused by the increase in rank of 
$C^{(1)}$ from unity to nearly fully rank, as shown in figure 3(d), indicative of a ‘rank instability’.
Figure 3. Energy in zonal wavenumbers for the Kolmogorov flow case with unity rank initialisation, for (a) QL and (b) CE2 initialised with the QL initial solution (QL IC), and (c) CE2 initialised with the QL endpoint solution (QL EP). Rank of second cumulant submatrices (d) 
$C^{(1)}$ and (e) 
$C^{(2)}$ as found in the CE2 solutions.
To probe this disagreement, we re-run the CE2 simulation initialised using the QL EP solution with three separate modifications.
First, in case I (figures 4a,b), energy is removed from the zonal wavenumber 
$m = 2$ in the IC so that the submatrix 
$C^{(2)}$ has zero rank initially with the rest of the ICs unchanged. Further, a unit rank random noise (power 
$10^{-12}$) is added to wavenumbers with 
$m > 2$. The resulting CE2 solution (figure 4a) disagrees with the QL EP (figure 3a), the rank instability appearing again in 
$m =1$; see figure 4(b). Eigenvalue spectra of 
$C^{(1)}$ (inset of figure 4b) confirm that the additional eigenvalues contributing to its higher rank are non-zero. Hence the disagreement observed between the CE2 solution and the QL solution in case I is attributable to a rank instability in 
$m=1$.
Figure 4. (a,c,e) Energy 
$E_m$ in zonal wavenumbers for CE2 cases I, II and III. (b,f) Rank of cumulant submatrices 
$C^{(m)}$ with 
$m = 1,2$. (d) The QL solution for case II. The rank instability is isolated in case I, which disagrees with the QL EP of figure 3(a). Case II suppresses the rank instability by locking the base state, resulting in consistent QL and CE2 solutions. Case III recovers the QL EP using eigenvalue truncation (insets in (b) and (f) show spectra at the indicated time).
Second, in case II (figures 4c,d) starting with the QL EP solution, the 
$m=0$ and 
$m=1$ components are frozen, thereby ensuring their unit ranks, while the 
$m = 2$ harmonic is assigned a small full rank perturbation. As previously, for CE2, the 
$m=2$ component decays away (figure 4c). Repeating the same experiment for the QL (figure 4d) also now sees the 
$m=2$ component decay in agreement with CE2. This indicates that if the same base state is considered for both systems, then the zonal stability characteristics are the same.
Third, in case III (figures 4e,f), the QL EP is again used as the IC, but during time integration, only the largest eigenvalue is retained for each second cumulant submatrix 
$C^{(m)}$ for all 
$m \leq M$. This ensures that unit rank is maintained at all times for the entire second cumulant (eigenvalue spectra in the inset of figure 4(f) confirm that the initially observed departures are trivial). The resulting CE2 solution for case III (figure 4e) is now in agreement with the QL solution (figure 3a). Thus the CE2 solutions in cases II and III, though exhibiting different zonal stability characteristics, agree with their corresponding QL solutions if the rank of the second cumulant is artificially kept at unity; it is the rank instability in 
$m=1$ that causes the observed divergence between CE2 and QL. As the rank of the 
$m=1$ harmonic grows beyond unity, there are more channels available for wave energy to flow, depriving the 
$m=2$ harmonic of any energy. As a check on our simulations, we have reproduced the rank and zonal instabilities on the sphere with independently developed code that uses an adaptive RK4 time integration algorithm (Plummer, Marston & Tobias Reference Plummer, Marston and Tobias2019; Marston & Tobias Reference Marston and Tobias2023). The Kolmogorov forcing employed is 
$F(\theta ) = -P_2(\cos {\theta }) - 8 P_8(\cos {\theta })$, where 
$P_\ell$ are Legendre polynomials, and 
$\theta$ is the co-latitude.
3.2. Stochastic forcing
The above examples indicate that there may be disagreement between the numerical solutions of the QL and CE2 systems. Relying on the determinism of the cases considered thus far, we have linked the disagreement to (1) a rank instability permitted only in CE2, and (2) a possible subsequent zonal mode instability arising in the QL but not in CE2. Our example shows both causes of divergence, but the rank instability is the key that may or may not trigger a subsequent difference in the zonal instability characteristics of the two systems.
We now show that divergences arising from rank instabilities are not limited to deterministic systems. Adopting the stochastic forcing strategy detailed in Constantinou (Reference Constantinou2015) and Constantinou et al. (Reference Constantinou, Farrell and Ioannou2016), the forcing term in (2.1) becomes 
$F = \sqrt {\varepsilon }\,\eta (x,y)$, where 
$\varepsilon$ is the energy injection rate, and 
${\eta }$ is the stochastic noise with variance 
$Q(x,y)$. Hence 
$\bar {F} = \bar {\eta } = 0$, and the forcing term 
$F^\prime = F = \sqrt {\varepsilon }\,\eta (x,y)$ in (2.4) drives the QL. The CE2 system is driven by the forcing term in (2.10), which becomes 
$\varGamma = \varepsilon \,Q ( x, y)$, where 
$Q(x, y) = Q(x_1 - x_2, y_1 - y_2) = \overline {\hat {\eta }(x_1,y_1)\,\hat {\eta }(x_2,y_2)}$ due to the statistical homogeneity of the Gaussian noise 
$\hat {\eta }$. Following Constantinou (Reference Constantinou2015) (see H.4 there), the covariance 
$Q(x, y)$ is specified on the Fourier domain by the forcing spectrum 
$\hat {Q} (m,n) = c(m)^2 d^2 \,{\rm e}^{-n^2d^2}$, which is non-zero for zonal wavenumbers 
$m \in [m_f, m_f + \delta m)$, where 
$d = 0.1$, and 
$c(m)$ is such that the net energy injection rate is 
$\varepsilon$.
Figure 5 shows the time-dependent zonal energy obtained when the stochastic forcing described above is applied to two zonal modes, 
$m = 8,9$ (i.e. 
$k_f = 8$ and 
$\delta k = 2$). The QL initialised with random noise (figure 5a) predicts three unstable zonal modes 
$m = 3,4,5$ and a stronger mean (
$m = 0$, black line). Alongside the forced zonal modes (
$m = 8,9$), these are the only significantly energetic zonal modes at 
$t = 1000$ days. The CE2 system initialised with the same IC (figure 5b) diverges from the QL solution, predicting a significantly weaker 
$m = 5$ mode. In contrast to the QL, the 
$m = 6$ mode is among the most energetic modes in CE2. This feature is also reproduced when CE2 is initialised with the QL EP solution (figure 5c). Once again, the complete picture emerges when ranks achieved by second cumulant submatrices are plotted for the diverging solutions. Figure 5(d) compares the ranks of cumulant submatrices for the two CE2 solutions against the unit ranks seen in the QL. The forced modes 
$m = 8,9$ in CE2 take on full ranks as prescribed by the method of forcing, while the remaining modes depart significantly from unity. Once again, we note a departure in the CE2 second cumulant rank from unity in tandem with a divergence between QL and CE2 solutions.
Figure 5. Energy in zonal wavenumbers for the stochastically forced case with unit rank and full rank initialisation, for (a) QL and (b) CE2 initialised using the QL initial solution (QL IC), and (c) CE2 initialised using the QL endpoint solution (QL EP). (d) The corresponding rank of second cumulant submatrices.
We now modify the CE2 solution obtained above (as diverged from the QL) to include a full rank initialisation of the second cumulant (power 
$10^{-4}$) coupled with a weakly energetic first cumulant (power 
$10^{-8}$). The zonal energy evolution in figure 6(a) indicates that such a CE2 solution returns to the CE2 solutions obtained with unity rank initialisation (figures 5b,c), i.e. with a relatively weak 
$m = 5$ mode (purple) and significant energy in 
$m = 6$ mode (brown). Hence CE2 with full rank initialisation also differs from the corresponding QL solution (figure 5a). Figure 6(b) shows that the 
$m = 5$ mode, though considerably weaker than in the QL, initially evolves to a unity rank up until 
$t = 500$ days – but the 
$m = 6$ mode retains a high rank, leading it to be considerably more energetic than in the QL. These results confirm that the correspondence between the QL state and CE2 statistical dynamics is not restored in a system with stochastic (full rank) noise and/or full rank initialisation.
Figure 6. (a) Energy in zonal wavenumbers 
$E_m$. (b) Evolution of the rank of cumulant submatrices 
$C^{(5)}$ and 
$C^{(6)}$ with full rank initialisation for the stochastically forced case. Despite the initialisation and forcing being full rank, the CE2 solution does not correspond to the QL solution (figure 5a).
4. Theory
A CE2 solution can stay coincident with a QL solution only if each second cumulant across the same excited zonal wavenumbers stays unit rank. The possibility of a rank instability, which converts a unit rank second cumulant into one that is multi-rank, is most clearly illustrated with a simple idealisation of the stability problem. (Nonlinear feedback on the mean, which would control these instabilities, is ignored here; for a more detailed discussion see Appendix A of Markeviciute & Kerswell (Reference Markeviciute and Kerswell2023).) Let the linear stability problem for the QL problem at a given zonal wavenumber be
\begin{equation} \partial_t \boldsymbol{u} = \boldsymbol{A} \boldsymbol{u}, \end{equation}
where 
$\boldsymbol {A}$ is an 
$N \times N$ real matrix, and 
$N$ measures the discretisation in 
$y$. If 
$\lambda _i$ and 
$\boldsymbol {e}(i)$ are the 
$1\leq i \leq N$ eigenvalues and corresponding right eigenvectors of 
$\boldsymbol {A}$, then there will be as many unstable directions as eigenvalues for which 
${\rm Re}(\lambda _i)>0$. The corresponding cumulant perturbation equation for 
$C_{ij}:=u_i u_j$ is
\begin{equation} \partial_t \boldsymbol{C} = \boldsymbol{A} \boldsymbol{C}+\boldsymbol{C} \boldsymbol{A}^{\rm T}, \end{equation}
which has a corresponding matrix eigenvalue problem with 
$N^2$ eigenvalues. It is straightforward to show that these are made up of all possible (ordered) pairings 
$\lambda _i+\lambda _j$ with corresponding eigenmatrix 
$\boldsymbol {C}(i,j)=\boldsymbol {e}(i)\,\boldsymbol {e}(j)$. (Building in the symmetry 
$C_{ij}=C_{ji}$, as is commonly done, reduces the eigenvalue count to 
$N(N+1)/2$ by not counting 
$\lambda _i+\lambda _j$ and 
$\lambda _j+\lambda _i$ separately.) The key observation is that if 
$\lambda _i+\lambda _j$ has a positive real part with 
$i \neq j$, then there is a new instability direction excitable in CE2 that is not available to the QL computation because the corresponding eigenvector has no counterpart there, i.e. it is uninterpretable in the QL. This is always the case when 
$\lambda _i$ and 
$\lambda _j$ have positive real parts, but the associated degrees of freedom are already excited separately due to the individual unstable eigenvalues, and may not be dynamically important. The situation is different, however, if one unstable eigenvalue partnered with a stable eigenvalue makes the pair unstable in CE2. This is significant because it could lead to stable degrees of freedom in the QL being stimulated in CE2. In practice, numerical errors will inevitably excite these new modes – or rank instabilities – in CE2, increasing the rank of the cumulant above unity.
A simple example makes this phenomenon clear. Consider the QL 
$u_1'=2u_1$ and 
$u_2'=-u_2$, which has one stable and one unstable direction. The second cumulant is given by
\begin{equation} C =\begin{pmatrix} u_1\\ u_2 \end{pmatrix}\otimes (u_1 \quad u_2)= \begin{pmatrix} u_1^2 & u_1 u_2\\ u_2 u_1 & u_2^2 \end{pmatrix}, \end{equation}
and it has one eigenvector with non-zero eigenvalue 
$u_1^2 + u_2^2$,
\begin{equation} \begin{pmatrix} u_1\\ u_2 \end{pmatrix}, \end{equation}and a second eigenvector with zero eigenvalue,
\begin{equation} \begin{pmatrix} -u_2\\ u_1 \end{pmatrix}. \end{equation}
The second cumulant therefore has a rank 
$1$, reflecting the fact that it is constructed from the outer product of a single vector with itself (see (4.3)).
The equivalent CE2 system where 
$C_{11} := u_1^2$, 
$C_{12} = C_{21} := u_1 u_2$ and 
$C_{22}: = u_2^2$ (and the symmetry 
$C_{12} = C_{21}$ is built in), i.e.
\begin{equation} \partial_t \begin{pmatrix} C_{11} \\ C_{12} \\ C_{22} \end{pmatrix}= \begin{pmatrix} 4 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -2 \end{pmatrix} \begin{pmatrix} C_{11} \\ C_{12} \\ C_{22} \end{pmatrix}, \end{equation}
has two unstable directions and one stable direction. It is possible to have a perturbation of form 
$(0\ 1\ 0)^{\rm T}$ in the extra unstable direction, but this perturbation has no equivalent in the QL problem as it is inconsistent with any one choice of 
$u_1$ and 
$u_2$. In other words, 
$C_{12}^2$ need not equal 
$C_{11}\,C_{22}$ as it must for the QL (4.3).
Importantly, a rank instability in CE2 will generally cause the mean flow to diverge away from its QL equivalent, as seen in the numerical experiments. If this effect is sufficiently strong, it is possible that one system may have different zonal stability properties to the other. In particular, one may be unstable to a new zonal wavenumber, whereas the other may not, at the same parameter values (e.g. see figure 3). This zonal instability would further increase the divergence between the two systems.
Finally, it is worth remarking that the divergence of QL and CE2 dynamics can occur even if ensemble rather than spatial averaging is used. The only prerequisite is that the QL dynamics are not full ranked (so there are damped degrees of freedom), giving ‘headroom’ for a rank instability in CE2. This is, for example, always the case for a properly resolved QL representation of a dissipative system. To be specific, a QL consisting of 
$N$ ordinary differential equations that tracks an 
$n$-dimensional attractor will typically have 
$n \ll N$ for a well-resolved computation (otherwise 
$N$ needs to be increased). As a result, an ensemble of all possible ICs in 
${\mathbb {R}}^N$ will evolve over time in the QL to give a second cumulant of rank 
$O(n)< N$. This opens up the possibility of a rank instability occurring in CE2 so that the propagated cumulant has a higher rank than that in the QL computation, thus producing divergence between the systems. Adding stochastic forcing will not change the situation unless the QL dynamics is forced to have a full ranked second cumulant and all singular values are significant (i.e. the forcing is strong). The latter condition excludes a situation where QL is weakly forced so that the singular values are partitioned into significant values relevant for the inherent dynamical attractor, and much weaker values that reflect the passive response to the forcing in otherwise damped degrees of freedom.
5. Conclusion
Our work indicates that care needs to be taken in interpreting the results of the statistical representation of QL. The instabilities described above, which can be triggered either by numerical errors, or more physically, by any of the ubiquitous sources of noise such as thermal fluctuations (Bandak et al. Reference Bandak, Eyink, Mailybaev and Goldenfeld2021) present in real systems can lead to solutions emerging for the statistical representation different to those that describe any single realisation of a QL flow. Perhaps implicitly recognising this issue, researchers have often initialised CE2 with a full rank second cumulant (Plummer et al. Reference Plummer, Marston and Tobias2019) – often termed a ‘maximum ignorance’ initial condition – and accepted that the resulting statistical description may not agree with that of any realisation of the underlying QL theory. The instabilities described above indicate that care must be taken in comparing direct statistical simulation (DSS) with the results from single realisations of a QL. We note that both the QL and CE2 systems are approximations to the full underlying dynamics, and it is not a priori obvious which choice better approximates the statistics of the fully nonlinear system. Higher-order truncations of cumulant expansions (such as CE2.5 and CE3) lead naturally to a higher-rank second cumulant, owing to the feedback of the higher cumulants on the evolution of the second cumulant (see e.g. Marston et al. Reference Marston, Qi and Tobias2019). We conclude by stressing that these results are in no way a criticism of either DSS at CE2 or investigations utilising a single realisation of a QL. We simply emphasise that these types of investigation can only elucidate different processes. Realisations of a QL model are useful for determining which nonlinear processes may be important in the dynamics of a turbulent flow, whilst statistical models such as CE2 yield information about the expected (in a statistical sense) response of a system.
Acknowledgements
We acknowledge support by funding from the European Union Horizon 2020 research and innovation programme (grant agreement no. D5S-DLV-786780). J.B.M. and S.M.T. are supported in part by a grant from the Simons Foundation (grant no. 662962, GF). J.B.M. acknowledges funding support from the US Department of Energy (grant no. DE-SC0024572).
Declaration of interests
The authors report no conflict of interest.
Open access
