## 1. Introduction

This paper revisits the seminal work of Malkus (Reference Malkus1956) that attempted to build a theory of shear turbulence. This theory was based upon maximising the momentum transport (or, equivalently, dissipation rate) achieved by the flow amongst all those with a marginally stable mean profile. Malkus clearly had a statistical form of marginal stability in mind but, to make progress, had to resort to specifying marginality with respect to the then 50-year-old Orr–Sommerfeld (OS) equation (Orr Reference Orr1907; Sommerfeld Reference Sommerfeld1908). So posed, his marginal stability idea was quickly repudiated (Reynolds & Tiederman Reference Reynolds and Tiederman1967 and, more recently, Iyer *et al.* Reference Iyer, Witherden, Chernyshenko and Vincent2019) although further studies showed it could be made to work if an anisotropic eddy viscosity model was used (Reynolds & Hussain Reference Reynolds and Hussain1972; Malkus Reference Malkus1979; Sen & Veeravalli Reference Sen and Veeravalli2000; Sen *et al.* Reference Sen, Veeravalli, Carpenter and Joshi2007). The concept of statistical stability was, however, central to Malkus's thinking (remaining so throughout his career, e.g. Malkus Reference Malkus1996 and Reference Malkus2003) and is clearly different from stability as viewed within the context of the governing Navier–Stokes equations (epitomised by the celebrated OS equation). For example, it is fairly uncontentious to assert that the turbulent attractor in, say, pressure-driven channel flow at a high enough Reynolds number and large enough domain has stationary statistics (defined by averaging over one or more homogeneous directions or in an ensemble sense) and so within the partial differential equations (PDEs) that govern how these statistics evolve, the realised turbulence is a stable fixed point – i.e. turbulence is statistically stable to infinitesimal perturbations of the statistics. This is in contrast to the time-dependent turbulent attractor as viewed in the Navier–Stokes equations where adding a small disturbance may well show that disturbance grows and never decays yet the original statistics still recover (an example is shown herein). The difference, of course, is that a flow disturbance can have a component along the turbulent trajectory and, hence, acts as a time shift: this part of the disturbance never decays to zero but does not affect the statistics (see, e.g. Nikitin Reference Nikitin2018). For his theory, Malkus wanted a statistical stability criterion based only on the lowest-order statistic – the mean flow. Ideally more statistical information needs to be incorporated and it is our objective here to attempt this. At the very least, doing so should improve the now-standard approach of carrying out linear stability analysis of the mean profile of time-dependent flows in an attempt to understand observed coherent structures (see, e.g. Crighton & Gaster Reference Crighton and Gaster1976; Gaster, Kit & Wygnanski Reference Gaster, Kit and Wygnanski1985; Roshko Reference Roshko1993; Barkley Reference Barkley2006; Lesshafft *et al.* Reference Lesshafft, Huerre, Sagaut and Terracol2006; Sipp & Lebedev Reference Sipp and Lebedev2007; Akervik *et al.* Reference Akervik, Ehrenstein, Gallaire and Henningson2008; Sipp *et al.* Reference Sipp, Marquet, Meliga and Barbagallo2010; Mantic-Lugo, Arratia & Gallaire Reference Mantic-Lugo, Arratia and Gallaire2014; Beneddine *et al.* Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016; Lefauve *et al.* Reference Lefauve, Partridge, Zhou, Dalziel, Caulfield and Linden2018).

The motivation for this work comes from two different directions. The first is the general approach of applying linear analysis around the mean profile of a time varying, possibly turbulent flow to deduce information about the likely fluctuations seen. Initially, this took the form of linear stability analysis stimulated by Malkus's work that tends to work well in free-shear turbulent flows like jets where inviscid (inflectional) instabilities dominate (Crighton & Gaster Reference Crighton and Gaster1976; Gaster *et al.* Reference Gaster, Kit and Wygnanski1985; Roshko Reference Roshko1993), but less well in wall-bounded situations such as channel flow where viscosity can be important (although there have been successes, e.g. Barkley Reference Barkley2006; Sipp & Lebedev Reference Sipp and Lebedev2007; Beneddine *et al.* Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016). Driven by the fact that shear flow mean profiles tend to be linearly stable, this approach subsequently diversified into non-modal analysis (Butler & Farrell Reference Butler and Farrell1993; del Alamo & Jimenez Reference del Alamo and Jimenez2006; Cossu, Pujals & Depardon Reference Cossu, Pujals and Depardon2009) and input–output or resolvent analyses (Chernyshenko & Baig Reference Chernyshenko and Baig2005; Jovanovic & Bamieh Reference Jovanovic and Bamieh2005; Hwang & Cossu Reference Hwang and Cossu2010*a*,Reference Hwang and Cossu*b*; McKeon & Sharma Reference McKeon and Sharma2010; Moarref & Jovanovic Reference Moarref and Jovanovic2012; Blesbois *et al.* Reference Blesbois, Chernyshenko, Touber and Leschziner2013; Sharma & McKeon Reference Sharma and McKeon2013; Schmidt *et al.* Reference Schmidt, Towne, Rigas, Colonius and Bres2018). The resolvent approach has been particularly illuminating in showing exactly how far linear analysis can go given the mean flow profile (McKeon Reference McKeon2017; Jovanovic Reference Jovanovic2021). Predictions can be made of the dominant fluctuation response at a given frequency and wavevector that typically resonates with observations at least up to amplitude and phase. In some sense this solves ‘half’ the problem of turbulence – given a mean flow, linear analysis around this can extract the dominant fluctuation structures – and refocuses attention on the ‘other’ half – predicting the mean profile. Identifying the required amplitudes and phases of all the fluctuation fields to support the observed mean profile, however, may approach the difficulty of ‘just’ solving the Navier–Stokes equations. One alternative is to appeal to some simpler ‘organising’ physical principle as epitomised in Malkus (Reference Malkus1956) and then the idea of statistical stability seems a key concept.

The second motivation is the resurgence of interest in dealing with statistical quantities directly through a cumulant expansion of time-varying flows (Hopf Reference Hopf1952; Orszag Reference Orszag1977; Frisch Reference Frisch1995). Generating the evolution equations for cumulants immediately highlights the closure problem of turbulence in which the time derivative of a $n$th-order cumulant requires knowledge of the $(n+1)$th cumulant and so forth. Present-day computational power has, however, renewed interest in the pursuit of simple cumulant-discard closures to see how well they do in modelling flows in a variety of contexts, for example, in atmospheric dynamics (Farrell & Ioannou Reference Farrell and Ioannou2003; Marston, Conover & Schneider Reference Marston, Conover and Schneider2008; Srinivasan & Young Reference Srinivasan and Young2012; Tobias & Marston Reference Tobias and Marston2013; Parker & Krommes Reference Parker and Krommes2014), astrophysics (Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011), plasmas (Farrell & Ioannou Reference Farrell and Ioannou2009; Parker & Krommes Reference Parker and Krommes2013) and wall-bounded shear flows (Farrell & Ioannou Reference Farrell and Ioannou2012; Constantinou *et al.* Reference Constantinou, Lozano-Duran, Nikolaidis, Farrell, Ioannou and Jimenez2014; Farrell *et al.* Reference Farrell, Ioannou, Jimenez, Constantinou, Lozano-Duran and Nikolaidis2016). The most popular closure ignores third- and higher-order cumulants – commonly referred to as CE2 – and has the nice property of being exact for a quasilinear (or historically a ‘mean field theory’) version of the Navier–Stokes equations (e.g. see the recent review by Marston & Tobias Reference Marston and Tobias2022). Significantly for our purposes here, CE2 and higher closures (CE$n$ where $n \geq 3$ is the highest-order cumulant retained) present the most natural framework in which to extend Malkus's idea of statistical stability. In what follows, the focus will be on using CE2 and its relationship to the quasilinear (QL) version of the Navier–Stokes equations to develop an improved version of the usual OS analysis of the mean profile; see figure 1. The overall objective is to develop a way to judge whether a state is statistically stable or not using a subset of its statistics.

It is worth emphasizing that the approach taken here is perfectly general and not confined to shear flows or even fluid mechanics. The key idea given a physical system described by a PDE is to generate evolution equations – ‘statistical’ equations – for the first few statistics of the solution. If the solution of the original PDE is spatiotemporally complicated but has steady statistics, the premise pursued here is that a better way to assess the stability of this solution is to look at the linear stability within the statistical equations rather than the original PDEs. The former is an approximation given that a finite number of the system's statistics are considered, but yields a spectral problem built upon a steady state in statistical space. The latter has to rely on a costly ensemble of simulations where a distribution of small perturbations are each added to the solution and their evolution monitored (Iyer *et al.* Reference Iyer, Witherden, Chernyshenko and Vincent2019; Kashyap, Duguet & Dauchot Reference Kashyap, Duguet and Dauchot2022). A good complementary example of where this approach would be useful is convection (Malkus Reference Malkus1954; Wen *et al.* Reference Wen, Ding, Chini and Kerswell2022).

The plan of this paper is to illustrate the analysis within the context of channel flow described in § 2. The Reynolds number is assumed high enough for the computational box used that the flow can be taken as approximately statistically stationary. An averaging procedure will be assumed below such that the mean flow can be assumed to only depend on the cross-stream variable $y$, i.e. $\boldsymbol {U}=U(y)\boldsymbol {\hat {x}}$. This could be ensemble averaging or averaging over the streamwise ($x$) and spanwise ($z$) directions. Using the latter spatially averaging approach, Malkus pointedly only chose a spanwise average so that his mean profile $\boldsymbol {U}=U(x,y) \boldsymbol {\hat {x}}$ could depend on the cross-stream and streamwise variables. As a result, his OS analysis targeted the stability of a streamwise-independent mean profile $U(y)\boldsymbol {\hat {x}}$ to streamwise-dependent mean flow disturbances $\delta \boldsymbol {U}(y) \exp ({\rm i}k(x-ct))$ rather than fluctuation fields defined as having non-vanishing spanwise dependence. Contrarily, there are growing arguments to only streamwise average to retain the spanwise structure of the mean flow, i.e. $\boldsymbol {U}=U(y,z)\boldsymbol {\hat {x}}$ (e.g. see table 1 of Lozano-Duran *et al.* (Reference Lozano-Duran, Constantinou, Nikolaidis and Karp2021) for a sample list of relevant works). Despite this, the focus here is on the simplest mean flow definition for a statistically steady state, $\boldsymbol {U}=U(y)\boldsymbol {\hat {x}}$, given the central role this plays in resolvent flow analysis but, there is no doubt, that extending the mean flow definition is clearly an important direction to extend the approach discussed here. Standard OS analysis is recalled in § 2.1.

Section 2.2 then introduces the cumulant expansion approach (Hopf Reference Hopf1952; Orszag Reference Orszag1977; Frisch Reference Frisch1995) and the hierarchy of evolution equations for these statistical quantities. Solving these equations to (hopefully) reach a steady state is an appealingly direct way to estimate the properties of statistically steady turbulent flows since it avoids having to average across large direct numerical simulation (DNS)-generated datasets. Here, however, the emphasis is on the concept of stability in this statistical framework and its relationship to (dynamic) stability within the Navier–Stokes equations (see § 2.3) not on the accuracy of suitable closures of the cumulant equations in capturing the reference flow state. That is, our strategy is to apply a statistical stability criterion from a statistical closure to a solution of the Navier–Stokes equations to approximate its statistical stability there. Our particular focus will be on the simplest non-trivial cumulant-discard scheme CE2 given the rapidly increasing dimensionality of the approach: CE$n$ works with cumulants up to order $n$ that, before exploiting any symmetries, is typically a $3n$–$2$ rank tensor (three spatial coordinates per field reduced by two averaging directions). Note that CE2 is exact for the quasilinearized Navier–Stokes equations – or QL equations – and translating what a statistically steady state in CE2 means for the QL equations is a crucial step discussed in § 2.4. Appendix A presents an equally important discussion on how the stability predictions within the CE2 and QL systems are related. Quasilinearization (Vedenov, Velikhov & Sagdeev Reference Vedenov, Velikhov and Sagdeev1961; Herring Reference Herring1963, Reference Herring1964) has enjoyed a resurgence of interest recently (see, e.g. Hernandez & Hwang Reference Hernandez and Hwang2020; Skitka, Marston & Fox-Kemper Reference Skitka, Marston and Fox-Kemper2020; O'Connor, Lecoanet & Anders Reference O'Connor, Lecoanet and Anders2021; Marston & Tobias Reference Marston and Tobias2022) given its accessibility, and together with its ‘generalized’ elaboration (Marston, Chini & Tobias Reference Marston, Chini and Tobias2016) in which the definition of what constitutes a mean is extended, has the ability to focus on different parts of the nonlinearity in the Navier–Stokes equations (see, e.g. Hernandez, Yang & Hwang Reference Hernandez, Yang and Hwang2021).

Section 3.1 then introduces an extended Orr–Sommerfeld analysis – or ‘EOS’ analysis – based on translating the statistical stability problem in CE2 back to the dynamical equations. Intriguingly, this is not the same as just working within the QL approximation (as Appendix A makes clear). Applying EOS analysis carries a substantial overhead so we consider a reduced (practical) version referred to as ‘minimally extended Orr–Sommerfeld analysis’ – or ‘mEOS’ analysis – in § 3.2 that is almost as cheap to apply as OS analysis.

Interestingly, applying the same strategy of mapping a cumulant-based system back to the underlying dynamical equations can only go one level higher in sophistication and requires a jump directly to CE$\infty$. This transforms EOS equations into the familiar linearised Navier–Stokes (LNS) equations albeit based around the steady base state derived by assuming stationary second rank cumulants. This ‘infinitely extended OS (iEOS)’ analysis is described in § 3.2. Section 4 then explores the performance of OS, EOS and mEOS analyses on four different turbulent states realised in two-dimensional (2-D) channel flow that is used as an approximation of a statistically steady flow. The limitations of the analysis are discussed in § 5 followed by a summary and final thoughts in § 6.

## 2. Formulation: channel flow

For context in this work, we consider a channel flow $\boldsymbol {u}^*(\boldsymbol {x}^*,t^*)$ of a fluid with density $\rho ^*$ and kinematic viscosity $\nu ^*$ between two parallel plates at $y^*=\pm h^*$ across which a time-dependent pressure gradient $9\rho ^* U^{*2}G(t^*)/4h^*\,\boldsymbol {\hat {x}}$ is imposed such that the bulk flow

is kept constant (unstarred/starred quantities are dimensionless/dimensional and periodicity is imposed across the spanwise domain $z^* \in h^*[-L_z,L_z]$). Non- dimensionalizing the Navier–Stokes equations using $h^*$, $3U^*/2$ (so Reynolds numbers based on the bulk speed and the laminar centreline speed $U^{*c}=3U^*/2$ correspond for the laminar parabolic flow) and $\rho ^*$ leads to

with $({1}/{4L_z})\int ^{L_z}_{-L_z}\int ^1_{-1} u\, \textrm {d}y\,\textrm {d}z=2/3$, where $\boldsymbol {u}=\boldsymbol {u}^*/U^*=u \boldsymbol {\hat {x}}+v \boldsymbol {\hat {y}}+w \boldsymbol {\hat {z}}$, $t=3U^*/(2h^*) t^*$ and $Re:=3U^* h^*/2\nu ^*$. We also impose streamwise periodicity of the flow over $x^* \in h^*[-L_x,L_x]$ so the non-dimensionalised flow domain is $(x,y,z) \in [-L_x,L_x] \times [-1,1] \times [-L_z,L_z]$ with non-slip boundary conditions on the plates at $y=\pm 1$ (fundamental wavenumbers in $x$ and $z$ are labelled $\alpha :={\rm \pi} /L_x$ and $\beta :={\rm \pi} /L_z$, respectively). In terms of an averaging procedure, a number of choices present themselves: ensemble averaging, spatial averaging and time averaging (or even a combination thereof) that all should be equivalent for a statistically stationary system in a large enough domain. However, in what follows we want to treat a numerical experiment in a finite domain with finite $Re$ and then ensemble averaging is the most natural choice as will come clear below. This averaging process is denoted by an overbar, $\overline {(\boldsymbol {\cdot })}$, and then the flow can be decomposed into a mean $U(y,t) \boldsymbol {\hat {x}}:=\overline {\boldsymbol {u}(\boldsymbol {x},t)}$ and fluctuation part $\boldsymbol {\tilde {u}}:=\boldsymbol {u}-\overline {\boldsymbol {u}}=\tilde {u} \boldsymbol {\hat {x}}+\tilde {v} \boldsymbol {\hat {y}} + \tilde {w} \boldsymbol {\hat {z}}$ (due to symmetry, a vanishing mean spanwise component is assumed so $(\overline {\tilde {v} \tilde {w}} )_y=0$).

The Navier–Stokes equations can be similarly decomposed into a mean part,

and a fluctuation part,

which is incompressible

where subscripts denote derivatives (e.g. $U_y:=\textrm {d}U/{\textrm {d}y}$).

### 2.1. Orr–Sommerfeld stability analysis

Given a possibly turbulent flow $(U,\boldsymbol {\tilde {u}})$, the ‘standard’ linear stability analysis is to consider small (also known as infinitesimal) perturbations $(0, \delta \boldsymbol {\tilde {u}})$ to a base state $(U,\boldsymbol {0})$ where the fluctuation field is ignored and the mean flow $U$ is assumed steady. As a result, only the (2.4) and (2.5) need be perturbed (and, hence, linearised) and since $U=U(y)$, the ensuing eigenvalue calculation is parameterised by a streamwise and spanwise wavenumber. Squire's theorem Squire (Reference Squire1933) is usually invoked to focus the search for instability to spanwise-independent perturbations and leads to the celebrated OS equation (Orr Reference Orr1907; Sommerfeld Reference Sommerfeld1908). In primitive variables, $\tilde {w}$ then decouples from $\tilde {u}$ and $\tilde {v}$ and can be ignored, leaving the reduced eigenproblem

where $(\tilde {u},\tilde {v},\tilde {p}) \propto \exp ({\textrm {i}m \alpha (x-ct)})$ and $c:=c_r+ic_i$ is the (complex) eigenvalue. The OS equation is then reached by defining a streamfunction $\psi$ (so $\tilde {u}=\psi _y$ and $\tilde {v}=-\textrm {i}m \alpha \psi$) and eliminating the pressure so that

In what follows, we actually work with the primitive variable formulation (2.6)–(2.8) as it is clearer to interpret the origin of new terms added below, the eigenvalue problem is better conditioned (two equations with second-order operators as opposed to one with a fourth-order operator) and it is easily extended to three dimensions if needed (for three-dimensional (3-D) disturbances, the OS equation must be augmented with Squire's equation for the cross-stream vorticity). Nevertheless, we refer to this general approach of doing linear stability around a turbulent mean as ‘OS analysis’ in recognition of its conception in Malkus's work.

The (matrix) size of the eigenproblem (2.6)–(2.8) is only $3N_y \times 3N_y$ for each streamwise wavenumber $m$ and so needs to be repeated $N_x$ times ($N_x,N_y$ represent the streamwise and wall-normal truncations). When the mean profile is known from simulations or experiments, for typical resolutions, this is an easily accessible procedure that, through the structure of the most unstable eigenvectors, can shed some light on the dominant structures of the turbulent flow.

A popular extension to the standard OS analysis involves using an eddy viscosity $E(y)$ instead of the molecular viscosity. An eddy viscosity can be determined self-consistently from the Reynolds stress needed to sustain the turbulent mean profile (e.g. Reynolds & Tiederman Reference Reynolds and Tiederman1967) using the steady version of (2.3),

Then, the eigenproblem can be modified to

We will use this to assess the performance of our EOS analysis in § 4.2.

### 2.2. Statistics: cumulants

In this section we consider a statistical framework for the flow by working with the equal-time cumulants of the flow (Hopf Reference Hopf1952; Orszag Reference Orszag1977; Frisch Reference Frisch1995). Even within this framework, we specialise further to exclusively consider equal-$x$-and-$z$ cumulants that are the subset of cumulants which influence the mean flow. (In fact, only equal-$y$ cumulants are needed in the mean flow equation – see (2.23) below. However, the evolution equations for these are not available without also solving for ‘non-equal’ $y$ cumulants.) The first cumulant is the mean $U(y,t)$. The second cumulant is the symmetric matrix

where we introduce the notation $C_{ij}(1,2):= \overline {[\boldsymbol {\tilde {u}}(x,y_1,z,t)]_i [\boldsymbol {\tilde {u}}(x,y_2,z,t)]_j}$ (here $[\boldsymbol {\tilde {u}}]_1=\tilde {u}$, $[\boldsymbol {\tilde {u}}]_2=\tilde {v}$ and $[\boldsymbol {\tilde {u}}]_3=\tilde {w}$) to emphasize the $y$ arguments and de-emphasize the implicit time dependence, and the third cumulant is the third-order tensor

These correspond to the second and third central moments of the flow, respectively. Cumulants and central moments, however, diverge at fourth order and beyond, e.g.

We will not go this high in the cumulant expansion used here but just note that the $n$th-order cumulant is $n$-dimensional in space so that storage when doing computations becomes prohibitive very quickly. Hence, the onus is on applying some sort of closure as soon as possible.

To derive evolution equations for the cumulants, we introduce a double Fourier series representation of the flow

where $m,n \in {\mathbb {Z}}$. Clearly $\boldsymbol {\tilde {u}}^{-m-n}=\boldsymbol {\tilde {u}}^{*mn}$ (the complex conjugate of $\boldsymbol {\tilde {u}}^{mn}$) for a real flow but it will be clearer not to build this into the notation in anticipation of deriving perturbation equations later. Hence, we write

(note, e.g. $C_{ij}^{mn}(1,2)=C_{ji}^{-m-n}(2,1)$). Equations to evolve the cumulants are obtained by temporally differentiating their definitions in (2.19) and (2.20) and using (2.4). For example, for the second-order cumulant,

where $\partial _i:= \partial _{y_i}$, $U(i)=U(y_i)$ for $i=1,2$ and $C_{i4}^{mn}(1,2):= [\boldsymbol {\tilde {u}}^{mn} (y_1,t)]_i \tilde {p}^{-m-n}(y_2,t)=:C_{4i}^{-m-n}(2,1)$ are three extra ‘velocity–pressure’ cumulants that get generated. Incompressibility conditions give the required three extra matrix constraints

along with the equation for the mean equation (2.3)

to close the system. The infamous closure problem of the Navier–Stokes equation is immediately evident here in that the evolution equation for the second-order cumulant depends on the third-order cumulant, a pattern that continues for higher-order cumulants so the system never closes. A popular (lowest) closure – commonly called CE2 – is to simply ignore the third-order cumulant which is equivalent to ignoring the fluctuation-fluctuation term (last bracketed term on the right-hand side of (2.4)). This is the QL approximation or sometimes referred to as the mean field theory (e.g. Vedenov *et al.* Reference Vedenov, Velikhov and Sagdeev1961; Herring Reference Herring1963, Reference Herring1964),

The defining feature of this approximation is that (2.25) is linear in $\boldsymbol {\tilde {u}}$ so that fluctuations with different wavenumbers are only coupled in the mean flow equation (2.24). This linearity also means that any fluctuation field (parametrised by streamwise and spanwise wavenumbers) can not be consistently in the stable manifold of $U$ as it varies with time. In particular, if $U$ is steady, only marginally stable fluctuation fields (typically with a temporal frequency) can be non-vanishing. Malkus argued for this model (and its marginal stability implications) on the basis that the fluctuation-fluctuation nonlinear term was only stabilising. This would be reasonable if bifurcations from unidirectional shear flows were always supercritical but, some decades later, subcriticality was realised the more generic situation (Kerswell Reference Kerswell2005; Eckhardt, Schneider & Westerweel Reference Eckhardt, Schneider and Westerweel2007; Graham & Floryan Reference Graham and Floryan2021; Kawahara, Uhlmann & van Veen Reference Kawahara, Uhlmann and van Veen2012).

Applying a closure at next order so $C^{(4)}$ is some assumed function of the lower cumulants or simply ignored (termed CE3) is less straightforward as the ensuing positive definiteness of the second cumulant is not automatic (Marston, Qi & Tobias Reference Marston, Qi and Tobias2019). This difficulty explains the popularity of the lower-order CE2 approximation where, for example, the existence and stability of steady solutions has recently been investigated for ordinary differential equation (ODE) systems (Li, Marston & Tobias Reference Li, Marston and Tobias2021; Li *et al.* Reference Li, Marston, Saxena and Tobias2022).

### 2.3. Approximations to statistical stability

The approach here is to consider the stability within the cumulant framework as this presents a natural way to assess stability of the flow statistics. Ideally, this should be attempted for a cumulant expansion which is high enough order to show a robustness against including even higher-order cumulants. However, the rate at which the dimensionality of this procedure explodes means that only second- and perhaps third-order closures are currently practical. As a result, we focus on CE2 here and identify a clear way to progress to higher order (see § 3.2).

Here CE2 is the statistical equations (2.21)–(2.23) with the third-order cumulants in (2.21) set to zero. The corresponding equations for perturbations $(\delta U, \delta C_{ij}^{mn} )$ upon a base statistical state $(U,C_{ij}^{mn})$ – hereafter referred to as the $\delta$CE2 problem – are

An equivalent equation arises in non-modal stability theory (e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Jovanovic & Bamieh Reference Jovanovic and Bamieh2005). If the pressure gradient is kept fixed $\delta G=0$, otherwise the constant volume flux condition $\int ^1_{-1} \delta U {\textrm {d}\kern 0.05em y}=0$ is the extra constraint required. Crucially the ansatz $(\delta U, \delta C_{ij}^{mn}) \propto \textrm {e}^{ \lambda t}$ is possible if the base state is independent of time – in other words, the base flow is statistically steady.

Even in the cheapest 2-D situation, the CE2 stability problem requires handling $5 N_x$ correlation matrices of size $N_y^2$ all linked through the mean equation. This leads to a matrix of size $(N_y+5 N_x N_y^2)^2$ or $\approx 25N_x^2 N_y^4$ elements which is out of reach even for modest resolutions. Given this, we explore a route to potentially still capture the essence of the CE2 statistical stability approximation without the considerable cost. To do this, we discuss a connection back to the equations of motion that plausibly retains the steadiness of the stability problem.

### 2.4. Steady statistics and simplications

A statistically steady base flow has steady mean $U$ and cumulants $\boldsymbol {C}$ with their Fourier components $\boldsymbol {C}^{mn}(1,2)$ also steady under ensemble averaging. This ensures that the associated stability problem (2.27)–(2.29) has temporally constant coefficients and is therefore a (conceptually at least) simple eigenvalue problem. In what follows below, we will not attempt to solve this directly as it is too unwieldy. Instead, a smaller, more practical QL stability problem is sought as a good proxy for it. The key in doing this is identifying a suitable base velocity field around which to develop a QL-type stability problem. A straightforward approach is to find the ‘best’ rank-1 approximation of each $\boldsymbol {C}^{mn}$ and use the associated velocity field, $\boldsymbol {\tilde {u}}^{mn}_0(y)$, as representative of the base flow. This can be accomplished by minimising the Frobenius matrix norm of the difference,

where, e.g. $\tilde {u}^{mn}_{(1:N)}:=[ \tilde {u}^{mn}(y_1)\ \tilde {u}^{mn}(y_2)\ \ldots \tilde {u}^{mn}(y_N)]^\textrm {T}$. Since $\boldsymbol {C}^{mn}=[\boldsymbol {C}^{mn}]^H$ (the Hermitian conjugate) and positive definite, the required $\boldsymbol {\tilde {u}}^{mn}_0$ is the leading right eigenvector of $\boldsymbol {C}^{mn}$ associated with the largest eigenvalue $\sigma _1$ scaled so that $|\boldsymbol {\tilde {u}}^{mn}_0|=\sqrt {\sigma _1}$.

We now discuss how the QL stability problem based upon the $\boldsymbol {\tilde {u}}^{mn}_0$ relates to the CE2 stability problem based on $\boldsymbol {C}^{mn}=\boldsymbol {\tilde {u}}^{mn}_0 \otimes \boldsymbol {\tilde {u}}^{*mn}_0$. Possible disturbances partition into two types: type A where the disturbance has energy in Fourier pairings not excited in the base flow, and type B where the disturbance has energy in Fourier pairings that are a subset of those present in the base flow, i.e. $\boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn}$ is only non-zero if $\boldsymbol {\tilde {u}}_0^{mn}$ is (see Appendix A for more detail). The former (type A) case is straightforward (since $\boldsymbol {\delta } \boldsymbol {C}=\boldsymbol {0}$) so we focus here on the type B situation. The QL stability problem – hereafter the $\delta$QL problem – is

where

(as ${\mathbb {L}}^{mn}$ is affine in $U$). This has temporally constant coefficients and, therefore, admits an eigenfunction of the form

The key point is that this has a direct equivalent in the CE2 problem (2.27)–(2.29) where the same eigenvalue $\lambda$ has the corresponding eigenfunction $\boldsymbol {C}^{mn}:= \boldsymbol {\tilde {u}}_0^{mn} (y_1,t)\otimes \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{-m-n}(y_2,t)+ \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn} (y_1,t)\otimes \boldsymbol {\tilde {u}}_0^{-m-n}(y_2,t) =\widehat {\boldsymbol {C}}^{mn}(y_1,y_2)\textrm {e}^{\lambda t}$. Given this, instability within the much more tractable $\delta$QL problem is therefore sufficient to conclude statistical instability within $\delta$CE2. Going further to claim that the stability within the $\delta$QL problem also implies stability within the $\delta$CE2 problem seems very likely but is not assured (see Appendix A). With this one caveat, we nevertheless assume that the stability of the (much) smaller QL system is a proxy for the statistical stability of the CE2 system: in particular, stability in the $\delta$QL system is taken to imply statistical stability within $\delta$CE2. This allows us to generate an extended version of OS, or EOS analysis but, before pursuing this in the next section, we make a remark and an observation.

The remark is that certain notational liberties have been taken here to keep the discussion as clear as possible. For example, the operator ${\mathbb {L}}^{mn}$ strictly maps an incompressible flow field to another that involves a supplementary scalar field (the pressure) entering into the definition. As is well known, this is determined by imposing incompressibility. An implicitly incompressible representation for the velocity field could be used (e.g. reducing the problem down to just using wall-normal velocity and vorticity) to avoid this wrinkle, but staying with primitive variables makes the various manipulations as clear as possible.

The observation is that the special form of the nonlinearity in the QL and CE2 formulations means that different wavenumber pairings only interact through the mean flow equation. As a result, eigenvalues can be sought within any subset of the wavenumbers possible and these are still valid within the full system of wavenumbers, i.e. this is not a truncation merely a subclass of disturbances. In particular, for CE2, the simplest perturbation of the base state $(U,\boldsymbol {C})$ just consists of perturbations in one wavenumber pairing and the mean flow,

Even this leads to an eigenvalue matrix calculation of size $(N_y+9N_y^2) \times (N_y+9N_y^2)$ in three dimensions. The equivalent QL perturbation is

and requires an eigenvalue matrix calculation of size $(N_y+8N_y) \times (N_y+8N_y)$ that is much smaller (or $(N_y+6N_y)\times (N_y+6N_y)$ in two dimensions – see Appendix B).

## 3. Extended OS analysis

We refer to the $\delta$QL problem that includes some information on the second-order flow statistics of the base flow as EOS stability analysis. The approach is as follows.

(i) Estimate the base flow statistics $(U, \boldsymbol {C})$.

(ii) Approximate each base cumulant component tensor $\boldsymbol {C}^{mn}$ as rank 1 to obtain a representative physical base field $\boldsymbol {\tilde {u}}^{mn}_0(y)$.

(iii) Solve the EOS eigenvalue problem for small perturbations $(\delta U, \,\boldsymbol {\delta } \boldsymbol {\tilde {u}}, \,\delta \tilde {p})$ that is

(3.1)\begin{equation} \partial_t\delta U=\frac{1}{Re}\partial^2_y \delta U+\delta G - \sum_{m,n} \partial_y \left(\delta \tilde{u}^{mn} \tilde{v}^{{-}m-n}_0 +\tilde{u}_0^{mn}\delta \tilde{v}^{{-}m-n} \right) \end{equation}coupled with the fluctuation equation for every Fourier mode $(m,n)$,(3.2)\begin{gather} \partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{mn} = \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}-\boldsymbol{\nabla} \delta \tilde{p}^{mn} - {\rm i}m U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} - \delta \tilde{v}^{mn}U_y\boldsymbol{\hat{x}} -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{mn} -\tilde{v}_0^{mn} \delta U_y \boldsymbol{\hat{x}}, \end{gather}(3.3)\begin{gather}{\rm i}m \alpha\, \delta\tilde{u}^{mn}+\partial_y\, \delta \tilde{v}^{mn}+{\rm i}n\beta\, \delta \tilde{w}^{mn} = 0 \end{gather}(In (3.2), the last two terms depend on $\delta U$ and will be discussed later).

This eigenvalue problem has size $(N_y+4 N_x N_y N_z)^2 \approx 16N_x^2N_y^2N_z^2$ that is probably impractical for all but the smallest systems since all the wavenumber pairings are coupled through the mean equation. A natural way to simplify the calculation is to only include a targeted subset of the wavenumber pairings and, going further, to only consider one wavenumber pairing. This latter approximation is obviously the most extreme but is also closest in spirit to the original OS analysis – we call this minimally extended OS, or mEOS, analysis.

### 3.1. Minimally EOS equations

Here, the sum over wavenumber pairings in (3.1) is removed leaving the mean flow equation (3.1) coupled with just one Fourier mode equation (3.2). This eigenvalue problem focuses on the coupling between an individual Fourier mode and the mean flow and its size is now $(5 \times N_y)^2 = 25 N_y^2$, so comparable to OS analysis. Just as for OS analysis, it needs to be repeated for each Fourier wavenumber of interest.

### 3.2. Infinitely EOS equations

Before going on to test EOS and mEOS analyses on some statistically steady flows, it is worth briefly discussing how far this approach of using the cumulant framework to further motivate more sophisticated versions of OS analysis can go. In fact, it turns out that only one further enhancement is possible and then only for CE$\infty$. This is because any intermediate closure will lead to an equation that when ‘unwrapped’ at the highest level contradicts those obtained at lower cumulant orders. To see this, consider CE3, the next closure after CE2 in which $\boldsymbol {C}^{(4)}$ is ignored. Unwrapping the statistical equation for $\boldsymbol {C}$ will recover the Navier–Stokes equations whereas unravelling the $\boldsymbol {C}^{(3)}$ equations will not as the nonlinear term leading to $\boldsymbol {C}^{(4)}$ has been dropped (any other closure will also suffer this inconsistency). The only way to avoid this is: (1) to only unwrap one cumulant equation – so CE2 which leads to EOS, or (2) ensure that all unwrapped equations are the same, which means no closure and CE$\infty$. In this latter case, the Navier–Stokes equations re-emerge that, under perturbation around a steady base state, lead to the linearized Navier–Stokes equations as the ultimate extension of OS analysis, i.e.

and, for each wavenumber pair $(m,n)$,

This could be called an iEOS analysis as there is no further extension possible with the cumulant expansion framework considered here. It is worth remarking that the extra term added in the iEOS analysis would get dropped in any minimal version where the individual wavenumber pairings are considered separately, i.e. minimalizing iEOS is just mEOS. There is, however, considerable potential in avoiding this drastic reduction in favour of retaining small subsets of interacting wavenumber pairings such as triads that interact through the newly present term. This or iEOS will not be pursued further here to avoid overcomplicating the discussion.

## 4. Application

We now test the extended stability approaches, EOS and mEOS, upon turbulent states obtained by DNS of 2-D channel flow. At high enough $Re$, the expectation is that the flow will be statistically stationary (Malkus Reference Malkus1956). The simulations are performed using the open-source PDE solver Dedalus (Burns *et al.* Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). The 2-D flow is simulated using a vorticity–streamfunction formulation so that $(u,v):=(\psi _y,-\psi _x)$ and

where the flow is driven by a streamwise body force $f(y,t)\boldsymbol {\hat {x}}$ so that the volume flux is fixed. Here $Re$ is defined as in § 2.1 so that the conditions

are imposed with the latter two reflecting the presence of non-slip walls. Following earlier work (Falkovich & Vladimirova Reference Falkovich and Vladimirova2018), the length of the channel $2L_x$ is set to four times its height ($L_x=4$), $Re=36\,300$ and, computationally, $1024$ Fourier modes are used to discretize in the $x$ direction and $256$ Chebyshev modes in $y$ direction (see Markeviciute & Kerswell (Reference Markeviciute and Kerswell2021) for details).

The streamwise body force is defined in terms of a profile function $g(y)$ as

where setting $g=0$ recovers the usual $y$-independent applied pressure gradient $-G(t)$. For this situation, it is already known that there is bistability with two statistically steady states possible in a 2-D channel flow at $Re=36\,300$ (Markeviciute & Kerswell Reference Markeviciute and Kerswell2021): a state that is statistically symmetric about the channel midplane – the symmetric state S – and another that is statistically asymmetric – the asymmetric state A; see figure 2(*a*). Choosing non-zero $g(y)$ is a device to diversify the test states available with two extra examples generated by using profiles $g_1:=(1-y^2)^6$ (state F1) and $g_2:=\cos \tfrac {3}{2} {\rm \pi}y^2$ (state F2).

### 4.1. Statistically stable test states

A spatial average over the streamwise direction and temporal averaging was used as a proxy for the ensemble-averaged statistics discussed in deriving EOS equations (these are equivalent for a statistically steady flow over a large enough domain and long enough time). The mean profiles for states S, A, F1 and F2 are shown in figure 2 along with their streamwise root-mean-squared velocity profiles, all obtained by time averaging over $10^4$ time units. Their respective power spectra are shown in figure 3(*a*). While the dominant wavenumber is the same for all four states ($k_d=2\alpha ={\rm \pi} /2$), significant differences are seen at neighbouring streamwise wavenumbers $2k_d$ and $3k_d$ for example. Typical temporal variations of the fluctuation field energy $E_{fluct}$ compared with the total energy of the flow $E$ – figure 3(*b*) – show the desired wide variety of mean fluctuation energy and fluctuation amplitudes. In particular, the amplitude of the fluctuations in state F2 are $\approx 60\,\%$ of the mean and only $\approx 30\,\%$ for state A.

For each state, the second-order cumulant matrices were computed every unit of time and then averaged over a period of $10^3$. During the process of calculating EOS and mEOS eigenvalues, correlation matrices time averaged over different time windows, $t = 500, 1000, 1500$, were used as well as correlation matrices time averaged with different time steps (0.5 and 1) for one time window $t=500$. The qualitative results were all the same, with only minor quantitative differences that did not affect the relative positions of eigenvalues obtained by different methods. The time-averaged cumulant matrix was then diagonalised and the flow field corresponding to the leading (real) eigenvalue used to typify the base fluctuation field. The implications of the time averaging to the rank of the matrix are discussed in § 5.1.

The EOS problem – (3.1) and (3.2) – and the mEOS problem were posed as generalized matrix eigenvalue problems of matrix size $\approx ( 4 N_x N_y)^2$ and $(5N_y)^2$, respectively. For EOS analysis, $N_x=20$ streamwise wavenumbers were used as a balance between accurate eigenvalue convergence and computational accessibility; see, e.g. figure 4. Wall-normal resolution of the eigenvalue problem was also varied to ensure eigenvalue convergence in the mEOS case. Doubling the wall-normal resolution from the DNS resolution of $256$ Chebyshev modes to $512$ produced less than a $10^{-8}$ relative error in the eigenvalues, so all stability calculations were done using the DNS wall-normal resolution.

### 4.2. Comparisons

Here we compare the eigenvalues obtained by standard (OS), extended (EOS) and mEOS stability analyses on the four different turbulent states. Since all our test cases are presumed statistically steady, for a model to be good at predicting statistical stability of the turbulent state, we expect all the eigenvalues to be stable, i.e. to have a negative real part. For the symmetric (figure 5*a*) and asymmetric (figure 5*b*) states, we observe no significant difference in the leading eigenvalues between the three stability models. While all three models predict the asymmetric state to be statistically stable, they also all predict statistical instability for the symmetric state. For these test states, the leading eigenvalue is associated with the first streamwise wavenumber (marked blue in the plot). We observe some changes for the other eigenvalues which we do not examine further as they are all stable and do not affect characterization of the statistical stability of the turbulent states.

Since all streamwise wavenumbers are coupled in the extended model, the leading wavenumber for the eigenvalue is assigned by examining the power spectrum of the corresponding eigenvector. Example power spectra for leading eigenvalues are shown in figure 6 where it is apparent that the mean velocity component of an eigenvector is at least an order of magnitude smaller than the leading wavenumber contribution.

The leading eigenvalue for the body-forced states states F1 (figure 5*c*) and F2 (figure 5*d*) is associated with the second streamwise wavenumber (marked red in the plot). For both F1 and F2 states, this eigenvalue is unstable in the standard stability analysis but becomes stable in EOS analysis. The stabilisation effect is also captured by mEOS analysis but is not so strong. While the statistical stability of the F1 state is predicted correctly by the extended model, the improvement is only partial for the F2 state. In addition to the unstable leading eigenvalue that is stabilised by the extended model, there are two subsequent unstable eigenvalues corresponding to the first and third streamwise wavenumbers (marked blue and green on the plot). Similar to the symmetric state case, these two unstable eigenvalues see no improvement towards stability when the extended stability models are applied.

We also repeat the standard OS analysis with eddy viscosity as described in § 2.1 (marked using squares in figure 5). While eddy viscosity has a general stabilising effect on the OS eigenvalues, this effect is rarely enough to stabilise the unstable OS eigenvalues. Significant improvement is only seen for the symmetric S state (top left plot) where the unstable eigenvalue corresponding to $m=1$ is stabilised, while the unstable eigenvalue corresponding to $m=3$ for the F2 state (bottom right plot) moves to just ${O} (10^{-4})$ under the instability line. However, the most unstable eigenvalues corresponding to $m=2$ for the states F1 and F2 where significant improvement is seen using the EOS analyses, are not notably affected by the eddy viscosity. Note, to avoid the singularity in the expression for eddy viscosity (2.10) for the asymmetric state, the mean velocity profile $U$ was symmetrised by $\tfrac {1}{2}(U(y)+U(-y))$ and the Reynolds stress $\overline { \tilde {u} \tilde {v}}$ was anti-symmetrised by $\tfrac {1}{2}( \overline { \tilde {u} \tilde {v}}(y) \unicode{x2013} \overline { \tilde {u} \tilde {v}}(-y) )$.

So, in summary, for each unstable eigenvalue observed in our test cases, the standard and extended models either agree on their statistical stability prediction or the extended model shifts the prediction towards statistical stability. The stabilising effect is qualitatively captured by both EOS and mEOS analyses with the mean velocity components making a non-zero power contribution to the eigenvectors of the extended model. We now examine why the extended stability models sometimes perform better than OS analysis and sometimes make very little improvement.

### 4.3. Eigenvalue perturbation analysis

In an attempt to rationalise the improvement (stabilisation) or not in mEOS of the leading eigenvalues, we carry out an eigenvalue perturbation analysis. While this neglects the interactions between different streamwise wavenumbers, it includes the minimal improvement that our models offer: adding the mean velocity perturbation into the stability consideration. This leads to an additional equation that governs the evolution of the mean velocity perturbation, and two extra terms – advection and shear – which affect the fluctuation perturbations. While formally limited to only small effects, the perturbation analysis provides a framework to study the structures of the mean velocity perturbations, base fluctuations and fluctuation perturbations and how their interactions affect the eigenvalues.

The unperturbed mEOS eigenvalue problem is the OS problem ((3.2) for one Fourier pairing without the $\delta U$ dependent terms) and (3.1) (we are imagining a small number $\varepsilon$ inserted in front of the $\delta U$ dependent terms and developing a perturbation expansion in this but actually $\varepsilon =1$). The latter is passive: it sets $\delta U$ but there is no feedback to the OS equation. Discretized, it reads

with matrices $A,B$, eigenvalue $\sigma$ and right eigenvector $\boldsymbol {y_R}:=[\delta U\ \boldsymbol {\delta } \boldsymbol {\tilde {u}}^{m0}]^\textrm {T}$. The the last two $\delta U$ terms in (3.2) are then treated as perturbations of $A$ that causes $(\sigma,\boldsymbol {y_R}) \rightarrow (\sigma +\delta \sigma, \,\boldsymbol {y_R}+\boldsymbol {\delta y_R})$ so that

where ${\delta A_A\boldsymbol {y_R}}$ corresponds to the advection term ${{-\textrm {i}m\alpha \boldsymbol {\tilde {u}}_0^{m0}\delta U_R }}$ and $\delta A_S \boldsymbol{y_R}$ to the shear term $-\tilde {v}_0^{m0} \partial _y \delta U_R \boldsymbol {\hat {x}}$. By taking the inner product using the corresponding left eigenvector of the unperturbed system $\boldsymbol {y_L}$ (${\boldsymbol y_L}^{{\dagger} } A=\sigma {\boldsymbol y_L}^{{\dagger} } B$), and considering first-order terms only, an expression then follows for the first-order perturbation to the eigenvalue $\delta \sigma$,

We present the two most interesting cases of the eigenvalue perturbation analysis. In the first case, we perform eigenvalue perturbation analysis on the leading eigenvalue of the standard OS analysis for the F2 state, which is seen to stabilise in mEOS and EOS analysis. The second case corresponds to the leading eigenvalue of the symmetric state, which shows no significant difference between the standard and extended approaches. The eigenvalues as well as perturbation analysis predictions are summarised in table 1. The left and right eigenvectors and perturbed terms are shown in figures 7 and 8.

For the body-force-driven F2 case, the advection term has a small and positive contribution to the real part of $\delta \sigma$ while the shear term has a large negative contribution to the real part of $\delta \sigma$, as summarised in table 1. This latter prediction is only qualitatively correct as the shear term is not a small perturbation giving a much larger predicted decay rate of $-0.169$ than the actual value of $-0.023$ obtained from mEOS analysis. The perturbation analysis, however, indicates that the shear in the perturbed mean field is causing the stabilisation. Looking in more detail at the fields in figure 7, it is clear that the active regions of the left and right eigenvectors coincide. In particular, the shear in $\delta U_R$ coincides with where the left eigenvector is significant and is of opposite phase giving the large stabilising effect. In contrast, the resulting advection term is not only an order of magnitude smaller than the shear term, but also has the same sign contribution as the left eigenvector, yielding a small destabilising contribution to the eigenvalue.

For the symmetric S case, the perturbation analysis shows both advection and shear terms make small destabilising contributions to the eigenvalue commensurate with mEOS analysis although there is still an error due to the finiteness of the perturbation. Examining the structure of the left and right eigenvectors shown in figure 8 indicates that they are ill-matched. The left eigenvector is concentrated much closer to the channel walls than the right eigenvector. Moreover, even though the right eigenvector component $\delta U$ has some structure near the wall, the interaction of this component with the perturbed operators moves this structure towards the centre of the wall for the shear term, avoiding any significant interaction with the left eigenvector. The advection term has some overlap with the left eigenvector, but it is also an order of magnitude smaller than the shear term, overall resulting in small destabilising contributions. It is also worth remarking that even at a resolution of $2048$ wall-normal Chebyshev modes, the eigenvectors possess the same small-scale structure seen in the figures.

The ‘take home’ message from this section is that the mean velocity perturbation, even if an order of magnitude smaller than the other components of the eigenvector (see figure 6), can interact with the base flow fluctuation fields in a way that produces a strong stabilising influence.

## 5. Model limitations

In this section we discuss the various approximations made in deriving EOS and mEOS which include time averaging the statistics, approximating correlation matrices as rank 1 and neglecting higher-order statistics. As already mentioned, the channel flow states considered above were assumed statistically steady. Making this assumption is required to apply Malkus's OS analysis in the first place and so subsequently EOS and mEOS analysis, but is only an approximation. Time averaging the statistics in a small computational box is a commonly used approach to get a better estimate of the steady statistics hypothetically generated in an infinite computational box.

### 5.1. Rank 1 approximation of correlation matrices

Time averaging the correlation matrix increases its rank from 1 up to potentially its full dimension; see figure 9 for singular values $\sigma _i$ of correlation matrices time averaged over 1000 time units. This is not necessarily a reflection that the statistics are non-stationary as a fluctuation field with two frequencies will give a rank-2 correlation matrix under averaging. For $m=2$, which is the dominant streamwise wavenumber in all of the test states, a gap between the leading and the second singular values is larger than an order of magnitude indicating that a rank-1 approximation of the time-averaged correlation matrix may be reasonable. On the other hand, time-averaged correlation matrices for $m \in \{1,3\}$ do not show a significant gap between the leading and second singular values. In this case, it is much more likely that some important fluctuation field features might not be accurately represented by the rank-1 representation of the correlation matrix and, as a result, the approach is less justified.

In an attempt to explore a more accurate representation of $\boldsymbol {C}^{mn}=\sum _{i=1} \sigma _i \,\boldsymbol {\hat {e}_i} \otimes \, \boldsymbol {\hat {e}_i}$ (where $\boldsymbol {\hat {e}_i}$ is the $i$th normalised eigenvector and $\sigma _i$ the largest eigenvalue or singular value of $\boldsymbol {C}^{mn}$ as it is Hermitian) beyond just the rank-1 approximation $\boldsymbol {C}^{mn} \approx \sigma _1 \,\boldsymbol {\hat {e}_1} \otimes \, \boldsymbol {\hat {e}_1}$ ($\boldsymbol {\tilde {u}}^{mn}_0=\sqrt {\sigma _1} \boldsymbol {\hat {e}_1}$), a weighted average of the leading $N$ fields,

was also considered for the F2 state where both $m=1$ and $m=2$ wavenumbers have a non-rank-1 time-averaged correlation matrix $C^{m0}$. However, including $N=1,2,5$ or even $10$ eigenvectors into the expansion showed no significant difference in the leading eigenvalues (less than $1\,\%$ difference in the absolute value of the eigenvalue). Understanding the effect of a rank-1 approximation really requires time stepping the full statistical equations which is beyond the scope of this work.

### 5.2. Time-averaged statistics

To understand the effect of time averaging the statistics in the EOS stability models, numerical experiments were performed in which the perturbations to the turbulent base flow were exposed to a time-varying base flow. This experiment was first performed for the symmetric (S) state where the unstable OS eigenvalue remains unstable in EOS and mEOS analyses. From the DNS data, 10 turbulent flow histories were taken of length $\Delta t=500$ separated by at least 1000 time units to avoid statistical dependence. Each was taken as a turbulent, time-dependent base flow. Two random perturbations with energy $E_{pert} = 10^{-5} E_{base}$ of the base flow were then added to each of the 10 base flow histories (making a total of 20 cases) and then their evolution computed over a $\Delta t=500$ time window using the QL equations. For $(\delta U, \boldsymbol {\delta }\boldsymbol {\tilde {u}}^{m0})$, these read

so that the perturbation-base flow fluctuation interactions are not included in the perturbation fluctuation equation (compare with (5.4) below where they are retained). Averaging the perturbation energy growth over the 20 different simulations revealed an unstable mode with exponential growth after initial transients decayed. The growth rate of this unstable mode was found to be $\sigma =0.00263$ which is within $5\,\%$ of the growth rate predicted by the EOS analysis; see figure 10. The most unstable eigenvectors of the EOS and mEOS analyses agree and are observed as the growing structure in this QL experiment; see figure 11. We conclude that time dependence of the base flow is not enough to recover the statistical stability of the symmetric state S.

### 5.3. Importance of higher-order statistics

To evaluate the need of flow statistics beyond the second order to recover statistical stability of the symmetric (S) state, we repeated the numerical experiment described above but this time the perturbation fields were time stepped using the linearised Navier Stokes (LNS) equations where perturbation-base flow fluctuation interactions are now included and the base flow is time varying,

We repeat the experiment with initial $E_{pert} = 10^{-5} E_{base}$ and $E_{pert} = 10^{-7} E_{base}$ obtaining the same qualitative results. When the energy growth of the perturbation field is time averaged over 20 different runs, the perturbation growth is seen to saturate after an initial transient growth; see figure 10(*a*). This suggests that the higher-order statistics ignored in QL equations are important here for recovering statistical stability of the symmetric (S) turbulent state.

### 5.4. Nonlinear effects

Extended OS stability analysis applied to the body-force F2 state shows a significant stabilisation of the eigenvalue corresponding to the second streamwise wavenumber, but shows no significant difference for the other two unstable eigenvalues of the standard OS approach. Repeating the numerical experiments described above applied to the F2 state shows that both the QL equations and the LNS equations yield rapid exponential growth of the perturbations applied to the time-dependent base flow; see figure 10(*b*). Periodically renormalising the growing perturbation in the latter LNS case and continuing the simulation shows sustained growth over a period of $O(1000)$ time units (not shown). This behaviour is in contrast to what was found for state S where the perturbation saturated (see figure 10*a*). On the face of it, this experiment seems to show that state F2 is linearly unstable but this does not necessarily follow. Time stepping a perturbation using the LNS equations (or in the tangent space of a turbulent attractor) can continually sample the stretching or contracting dynamics along the attractor which can dominate the anticipated decaying behaviour perpendicular to the attractor. It seems stretching along the attractor dominates for state F2 but not for state S. A comparison of the leading growing structures in figure 12 is consistent with this: the QL and LNS flows bear no similarity with the predictions of the mEOS equations. These LNS calculations for the state F2 highlight perfectly the difficulty in assessing the (linear) statistical stability using the Navier–Stokes equations touched on in the introduction.

Repeating the experiment (with $E_{pert} = 10^{-5} E_{base}$) but now time stepping the perturbation field with the full Navier–Stokes equations, we see the perturbation energy saturate after an initial transient; see figure 10(*b*). Importantly, and as anticipated, the statistics of the turbulent state F2 recover albeit slowly. Figure 13 shows the relative difference between the correlation matrices averaged over a rolling time window between the perturbed and unperturbed F2 state. The statistical differences ebb away after an initial transient of growth. Since F2 is statistically linearly stable, this initial growth indicates that the statistical evolutionary operator linearised around the state F2 fixed point has to be non-normal. Finally, figure 14 shows the evolution of the streamwise perturbation velocity field. We see that the highly localised structure visible at $t=10$ spreads in the channel at $t=20$ and reaches its saturated state at $t=30$ that does not change qualitatively even at much longer times $t=200$. It takes much longer for the statistics to recover suggesting the adjustment along the attractor happens much quicker than normal to it; see figure 13.

## 6. Discussion

### 6.1. Summary

We first summarise what has been done in the paper. The motivating objective has been to develop a theoretical approach which can assess the statistical stability of a turbulent state and goes beyond just examining the mean flow as Malkus did in the 1950s (Malkus Reference Malkus1956). That such an approach should exist in some form is predicated on the fact that the statistics of a disturbed turbulent flow always return to their undisturbed values if the disturbance is not too large regardless of whether it decays or not. That is, a turbulent flow is assumed linearly statistically stable (note this does not preclude the possibility of the statistical disturbance initially increasing before ultimately decaying away as the linear statistical operator may be non-normal). The way forward is to include more statistical information than just the mean of the flow in the analysis and the simplest such is the second rank cumulant or two-point equal-time correlations of the flow fluctuations. To do this, we have proposed to examine the statistical stability of the turbulent flow from within a cumulant equation framework truncated to just consider these statistical quantities (called CE2). Even in CE2, however, the ensuing statistical spectral problem is so unwieldy that a substitute has had to be sought. Here the well-known connection between CE2 and the QL equations can be exploited and we have argued that a suitably designed QL spectral problem defined in § 3 can act as a good proxy to detect statistical instability. The key in designing this QL problem is to translate a statistically steady base flow in CE2 into a representative steady base flow in the QL problem with no pretence that this base flow satisfies the QL equations itself (the QL equations offer an approximation of uncertain accuracy to the Navier–Stokes equations after all).

The issue of how the statistical spectral problem in CE2 is related to the QL spectral problem is an intriguing one that seems not to have been discussed before. It is important because it underpins the tacit assumption that a CE2 simulation starting from given initial conditions should exactly shadow (statistically) the equivalent QL simulation given one is derivable from the other (e.g. Marston *et al.* Reference Marston, Conover and Schneider2008; Tobias *et al.* Reference Tobias, Dagon and Marston2011; Tobias & Marston Reference Tobias and Marston2013). However, this ignores the difference in dimensionality between the two formulations and the possibility of triggering unstable ‘non-physical’ perturbations in the larger CE2 problem (meaning perturbations that cannot exist in the QL setting). In Appendix A we have tried to clarify the different types of instability that exist in these formulations – type A and type B – and how they are related. In particular, for type A instabilities, the spectral QL problem captures the most unstable disturbance possible in the CE2 system and so is an accurate proxy for CE2: CE2 can only be type A unstable if the QL problem is. For type B, the situation is less certain. While it remains the case that any type B instability in the QL setting is mirrored in CE2 and CE2 will likely have further type B instabilities, it is not clear whether CE2 can have a type B instability and the QL problem not. Given the as-yet unreported consequence of this – the CE2 equations could suffer a bifurcation whereas the QL equations not – it is tempting to assume that this is not possible and then that the QL analysis is also a good proxy for type B perturbations. With this one caveat, probing the QL spectral problem – named here EOS analysis – then appears the best way to examine the statistical stability of a flow where the two-point correlations as well as the mean profile of the flow are incorporated.

Despite improving upon the CE2 spectral problem, EOS analysis is still unfortunately quite costly to implement so a reduced version was proposed in § 3.1 – mEOS analysis – which has a similar cost to the usual OS analysis of the flow mean and is the most straightforward extension of doing linear stability analysis only around the mean flow. Just as in OS analysis, one chooses a particular Fourier mode $(m,n)$ to analyse but now extra information about the corresponding base Fourier mode is included. Minimally EOS analysis was found to capture the qualitative effect of EOS analysis with either both producing an improved estimate of statistical stability (e.g. the $m=2$ mode for state F1; see figure 5*c*) or offering no correction of the OS eigenvalue (e.g. the $m=1$ mode for state F1). Reassuringly, both approaches never made a poorer statistical stability prediction than OS analysis at least over the turbulent channel flow states tested. Minimally EOS analysis works comparatively well because the structure of the key eigenmodes in the EOS problem tend to consist of one dominant (primary) fluctuation field interacting with a mean flow perturbation (see figure 6). This is presumably because the comparatively small size of the cumulants or base fluctuation fields relative to the mean profile (e.g. indicated by the singular values in figure 9) means that only the few fluctuation fields close to marginality on the mean flow play a role. This may well change with increasing $Re$ but the relative size of the cumulants should also increase possibly acting to deprioritise this near-marginal stability property. Clearly higher $Re$ flows warrant further study.

### 6.2. Implications

By incorporating two-point fluctuation correlations, EOS analysis and its slimmed down version mEOS analysis, represent an attempt to generate a better way to assess the statistical stability of a turbulent flow than just considering the mean profile. Malkus (Reference Malkus1956) wanted to use this to constrain the set of realisable turbulent mean profiles and it still remains an interesting requirement given the paucity of predictive alternatives. The concept of marginal stability still periodically rears its head in the literature to post-rationalise turbulent mean profiles (see, e.g. Brauckmann & Eckhardt (Reference Brauckmann and Eckhardt2017) and references therein) and EOS analysis represents a finer tool for this. As noted in the introduction, the mean flow is increasingly taken as the starting point for resolvent flow analysis. To move this towards a full wall-bounded turbulence theory, a way has to be found to close the cycle between the fluctuation field and the mean. Just maybe, requiring statistical stability with respect to the dominant fluctuation fields could help constrain this.

A byproduct of this work is that it has also suggested an improved approach for predicting coherent structures in turbulent flows. Coherence here means weakly damped or growing fluctuation and mean disturbances that emerge from the linear operator built around the base flow rather than, say, proper orthogonal decomposition modes derived from the two-point correlations that have no direct relation to the underlying equations. Initial work (Crighton & Gaster Reference Crighton and Gaster1976; Gaster *et al.* Reference Gaster, Kit and Wygnanski1985; Roshko Reference Roshko1993) linearising around the turbulent mean proved successful in free-shear flows but less so when viscosity plays a leading role. Admittedly, EOS or mEOS analysis needs extra (second cumulant) information about the base flow that adds cost to the analysis but the potential for better predictions is clear. A popular (cheap) alternative approach has been to apply OS analysis around the mean profile that incorporates an eddy viscosity to acknowledge the presence of the base fluctuations (Reynolds & Hussain Reference Reynolds and Hussain1972; Sen & Veeravalli Reference Sen and Veeravalli2000; del Alamo & Jimenez Reference del Alamo and Jimenez2006; Sen *et al.* Reference Sen, Veeravalli, Carpenter and Joshi2007; Cossu *et al.* Reference Cossu, Pujals and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010*a*,Reference Hwang and Cossu*b*). This eddy viscosity is taken to be that which is needed to sustain the turbulent mean profile and so makes no attempt to incorporate fluctuation information for each excited wavenumber of the base flow as in EOS analysis. In the well-known triple expansion of Hussain & Reynolds (Hussain & Reynolds Reference Hussain and Reynolds1970), the coherent wave part ($\tilde {f}$ in their equation (1.1)) is the equivalent of $\boldsymbol {\delta } \boldsymbol {\tilde {u}}^{mn}$ that there is excited by a wavemaker and then studied as it travels downstream in a turbulent flow. Orr–Sommerfeld analysis was found unable to reconcile what was seen experimentally but EOS analysis should do better.

The cumulant framework examined here to motivate EOS analysis was used only to the first non-trivial level where only the first and second cumulants are retained. Whether this represents a good approximation or not to modelling the turbulent flow itself is a question not considered here. Rather the focus is on the associated spectral problem characterising the statistical stability within this framework taking the turbulent base flow as a given from DNS or experiment. Clearly considering higher cumulants beyond the second would be better but the computational cost of handling the ensuing spectral problem becomes prohibitive. Interestingly, ‘projecting down’ to an underlying Navier–Stokes-type problem fails immediately unless no closure at all is performed, which means CE$\infty$ is considered. This gives what we have called the iEOS analysis that is as costly as EOS analysis and unfortunately does not have a ‘minimal’ version. Again, only targeting certain subsets of fluctuation fields – e.g. key triads suggested by resolvent flow analysis – may yet make this tractable.

### 6.3. Future perspectives

In terms of future work, an obvious focus is to improve the main approximation in EOS analysis that is deriving a representative base flow. Here, we have used the leading singular vector of the second cumulant for each excited (vector) wavenumber but other strategies involving information from higher-order cumulants are certainly possible (leaving aside the cost of computing these$\ldots$). This would seem sensible in the case of iEOS analysis that is a proxy for CE$\infty$ where all the cumulants are retained and the only other cumulant system beyond CE2 that allows a Navier–Stokes-type proxy. Pragmatically, a steady base flow has been sought that leads to a (conceptually) simple eigenvalue problem. Trying to go beyond this profoundly complicates matters.

One interesting future direction hinted at in the introduction is to relax the averaging procedure to just being over the streamwise direction or, more generally, redefine what constitutes the ‘mean’. A key feature we have been able to exploit is the connection between CE2 and the QL equations that still holds if the mean is more flexibly defined to include some wavenumber subset of the fluctuation field. The then ‘generalized’ QL equations (Marston *et al.* Reference Marston, Chini and Tobias2016) have the advantage of including more of the nonlinearity in the Navier–Stokes equations and can interpolate between the QL equations and the full Navier–Stokes equations depending on how the mean is defined. A minimal mean extension would be to include one spanwise wavenumber field whereas including all fluctuation fields with no streamwise variation is equivalent to a mean defined by only streamwise averaging. There are many theoretical reasons for considering the latter, where the mean is $U(y,z)\boldsymbol {\hat {x}}$ rather than just $U(y)\boldsymbol {\hat {x}}$, in wall-bounded flows (e.g. the SSP/VWI mechanism Waleffe Reference Waleffe1997; Hall & Sherwin Reference Hall and Sherwin2010), but capturing it experimentally is a very time-consuming task and using it reduces the predictive power of resolvent flow analysis. Even in computations, the spanwise homogeneity needs to be broken otherwise there is no reason to suppose the mean flow should not be spanwise invariant (and then the spanwise structure of the mean profile needs to be divorced from the source of spanwise inhomogeneity). Nevertheless, this choice is starting to receive attention at least at the level of doing OS analysis on $U(y,z)\boldsymbol {\hat {x}}$ (see Lozano-Duran *et al.* (Reference Lozano-Duran, Constantinou, Nikolaidis and Karp2021) and references herein) and this analysis can be extended as envisaged here.

Many challenges remain. Statistical stability is an important fundamental feature of turbulent flows yet is difficult to assess and, hence, to exploit. Hopefully, the work reported here represents a useful but admittedly small step forward.

## Acknowledgement

We thank one referee for a very helpful report that improved the presentation of the paper.

## Funding

V.K.M. acknowledges financial support from EPSRC through a studentship.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Relationship between dynamic stability in QL and statistical stability in CE2

In this appendix we argue that the dynamical stability problem in the QL equations and statistical stability in CE2 based on $\boldsymbol {\tilde {u}}_0^{mn}(y)$ can be used to predict statistical instability within CE2 based upon $\boldsymbol {C}^{mn}(1,2)= \boldsymbol {\tilde {u}}_0^{mn} (y_1) \otimes \boldsymbol {\tilde {u}}_0^{-m-n} (y_2)$. This is because any instability present in the QL equations has an equivalent in CE2. We have been unable to establish the converse however plausible: a statistical instability in CE2 does not necessarily imply a dynamical instability in the QL equations. Despite this lack of complete equivalence, the main text takes the stability problem in the QL equations to be a good proxy of that in CE2.

There are two types of instability in the QL and CE2 equations. To identify them, we define two sets: $\varLambda$ that contains all the Fourier modes that contain energy in the base state $\boldsymbol {\tilde {u}}_0$, and $\delta \varLambda$ that does the same for the perturbation field $\boldsymbol {\delta } \boldsymbol {\tilde {u}}$, i.e.

*a*,

*b*)\begin{equation} \varLambda:=\left\{(m,n) \left| \int^1_{{-}1} |\,\boldsymbol{\tilde{u}}_0^{mn}|^2 \, {{\rm d} y} > 0 \right. \right\} \quad \text{and} \quad \delta \varLambda:=\left\{(m,n) \left| \int^1_{{-}1} |\,\boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}\,|^2 \, {{\rm d} y} > 0 \right. \right\}. \end{equation}

There are then two types of eigenfunctions for the linear stability problem in the QL problem (the $\delta$QL problem). These are as follows.

(i) Type A: $\delta \varLambda \subseteq \varLambda '$, the complement of $\varLambda$ that is the set of all possible unexcited Fourier modes in $\boldsymbol {\tilde {u}}_0$.

(ii) Type B: $\delta \varLambda \subseteq \varLambda$ when there is an exact equivalent eigenfunction for $\delta$CE2 (the statistical stability problem in CE2).

Type A is the most straightforward to understand as there is no feedback of the fluctuation field onto the mean in the $\delta$QL problem (the feedback is quadratic in the fluctuation disturbance). In this case, the stability problem simplifies to an OS-type problem on the base mean profile. As a result, a type A eigenfunction $\delta \boldsymbol {\tilde {u}}_A$ in the $\delta$QL problem consists of just one wavenumber pair $(\hat {m},\hat {n})$. The situation is similar in CE2 with the corresponding eigenmatrix $\delta C_A=|\delta \boldsymbol {\tilde {u}}_A|^2$ but with two notable differences. The first is the corresponding eigenvalue in CE2 is $2\mathrm {Re}(\lambda _A)$, where $\lambda _A$ is the QL eigenvalue ($\partial _t |\delta \boldsymbol {\tilde {u}}_A|^2=(\lambda _A \boldsymbol {\tilde {u}}_A)\boldsymbol {\tilde {u}}_A^*+\boldsymbol {\tilde {u}}_A(\lambda ^*_A \boldsymbol {\tilde {u}}_A^*)=(\lambda _A+\lambda _A^*)|\delta \boldsymbol {\tilde {u}}_A|^2$). The second is that a mean flow perturbation is forced,

(in the simpler constant-pressure gradient situation), which in turn forces perturbation components across all wavenumber pairings in $\varLambda$ since

where ${\mathcal {L}}(\lambda _A)$ is the spatial operator defined by the system (2.27) and (2.28) with $\partial _t$ replaced by $\lambda _A$. Importantly, if a type A QL eigenvalue is unstable, i.e. $\mathrm {Re}({\lambda _A})>0$, then so is the corresponding CE2 eigenfunction.

The eigenvalue problem for type B is more complicated since it involves multiple wavenumber pairs coupled to a mean field perturbation. However, it is straightforward to see that a type B eigenfunction in the $\delta$QL problem,

where eigenvalue $\lambda _B$ corresponds to the eigenfunction

with same eigenvalue in the CE2 system.

There are no hybrid eigenfunctions in the $\delta$QL problem with $\delta \varLambda \not \subseteq \varLambda$ as those wavenumber pairs not in $\varLambda$ simply decouple even if they have the same eigenvalue (which is then degenerate). It is worth briefly illustrating this partitioning into type A and type B perturbations in a simple example. Let $q(x,t):=\overline {q}(t)+q^\prime (x,t)$ over the periodic domain $x\in [-{\rm \pi},{\rm \pi} ]$ where $\overline {(\boldsymbol {\cdot })}=({1 }/{2{\rm \pi} })\int ^{{\rm \pi} }_{-{\rm \pi} } (\boldsymbol {\cdot }) \,{\textrm {d}\kern 0.05em x}$ so $\overline {q^\prime }=0$ (and we assume that $q^\prime (-x)=q^\prime (x)$ for simplicity) then consider the QL-like system

where the mean flow is forced by an ‘imposed pressure gradient’ $3/Re^2$, there are the usual dissipation terms and quadratic feedback of the fluctuations onto the mean. There are multiple steady states: we choose to linearise around

*a*–

*c*)\begin{equation} \overline{q}=\frac{2}{Re}, \quad q^\prime_2=\frac{1}{Re}, \quad q^\prime_n=0,\quad n \neq 2, \end{equation}

which leads to the linear perturbation equations

These highlight the difference between ‘new’ wavenumbers ($n \neq 2$) and ‘already-present’ wavenumbers ($n=2$). Using the base flow and restricting to the perturbation to just $n \in \{1,2\}$, then

The associated statistical equations for $\overline {q}$ and $C:= \overline {q^{\prime 2}}$ are

The basic state in this formulation is $\overline {q}=2/Re$ and $C_2=1/Re^2$ (all other $C_n=0$) and the linear ‘statistical’ stability equations (again just for $n \in \{1,2 \}$) are

This leads to the problem

which is similar but different to that in (A12). Ignoring the factor of $1/Re$, the ‘QL’ eigenvalues are $\lambda _A=3/2$ (type A), $\lambda _B=\tfrac {1}{2}(-1+\textrm {i}\sqrt {7})$ and $\lambda _B^*$ (type B) with corresponding eigenvectors (respectively)

*a*–

*c*)\begin{equation} \left[ \begin{array}{@{}c@{}} 0 \\ 1 \\ 0 \end{array} \right], \quad \left[ \begin{array}{@{}c@{}} \lambda_B \\ 0 \\ 1 \end{array} \right] \quad \text{and} \quad \left[ \begin{array}{@{}c@{}} \lambda_B^* \\ 0 \\ 1 \end{array} \right], \end{equation}

whereas the ‘CE2’ eigenvalues are $2\mathrm {Re}(\lambda _A)=3$ (type A), $\lambda _B$ and $\lambda _B^*$ (type B) with eigenvectors (respectively)

*a*–

*c*)\begin{equation} \frac{1}{\lambda_A^2+\lambda_A+2}\left[ \begin{array}{@{}c@{}} -2 \\ \lambda_A^2+\lambda_A+2 \\ -\lambda_A \end{array} \right], \quad \left[ \begin{array}{@{}c@{}} \tfrac{1}{2} \lambda_B \\ 0 \\ 1 \end{array} \right] \quad \text{and} \quad \left[ \begin{array}{@{}c@{}} \tfrac{1}{2} \lambda_B^* \\ 0 \\ 1 \end{array} \right]. \end{equation}

This example shows: (*a*) the partitioning of the perturbations into type A and type B, and that (*b*) CE2 has double the corresponding type A QL eigenvalue (as it is real) and the corresponding CE2 eigenfunction contains other forced components beyond the new wavenumber $n=1$. However, these extra components, which get excited in the QL equations at (next) quadratic order, are passive and do not affect the eigenvalue.

Now we discuss how the $\delta$QL problem is a strict subset of the $\delta$CE2 problem. Things are clearest for the $\delta$CE2 type A case where the eigenvalue problem is solely determined by the statistical stability equation for the second rank cumulant,

Now if $\boldsymbol {e}^{(p)}$ is the $p$th eigenfunction to the OS problem with corresponding eigenvalue $\lambda _p$, it is straightforward to see that $\delta C_{ij}(1,2) = \boldsymbol {e}^{(p)}_i(1) \boldsymbol {e}^{*(p)}_j(2)$ is an eigenfunction of (A21)–(A22) with the corresponding eigenvalue $\lambda _p+\lambda _p^*=2\mathrm {Re}(\lambda _p)$ as already mentioned (just above (A2)). Less clear but still true is the fact that $\delta C_{ij}(1,2) = \boldsymbol {e}^{(p)}_i(1) \boldsymbol {e}^{*(q)}_j(2)$ is also an eigenfunction with eigenvalue $\lambda _p+\lambda _q^*$ which has no counterpart in the OS problem as it does not correspond to one physical flow field (note: first, it has also forced components via (A2) and (A3); and second, that for complex $\lambda _p+\lambda _q^*$, $\delta C_{ij}$ is not Hermitian but the average with its Hermitian conjugate – an eigenfunction with eigenvalue $\lambda _p^*+\lambda _q$ – is). This has the consequence that CE2 has at least the number of unstable directions as the OS problem and probably more since an unstable eigenvalue linearly combined with a weakly stable eigenvalue could produce a new unstable direction. That is, if the OS problem is unstable so is CE2 and *vice versa*. A simple example makes this clear: consider the QL (fluctuation) system $\dot {x}=2x$ and $\dot {y}=-y$ which has one stable and one unstable direction. The equivalent CE2 system where $C_{xx}:=x^2;\, C_{xy}=C_{yx}:=xy$ and $C_{yy:}=y^2$ (and the symmetry $C_{xy}=C_{yx}$ built in),

has two unstable directions and one stable direction. The point is in CE2, it is possible to have a perturbation of form $[0\ 1\ 0]^\textrm {T}$ in the extra unstable direction but this perturbation has no equivalent in the OS problem as it is inconsistent with any one choice of $x$ and $y$. This feature holds for any type A perturbation in CE2 as long as the corresponding velocity is a vector, i.e. there is a non-trivial structure in the cross-stream direction (this was absent by design in the previous example (A6)–(A7)). This is also clear from the dimensions of the respective problems – the QL eigenvalue problem is $O(N_y \times N_y)$ whereas the corresponding CE2 problem is $O(N_y^2 \times N_y^2)$ so it is clear there have to be $O(N_y^2-N_y)$ additional eigenvalues in $\delta$CE2.

Importantly, in type A perturbations, it is clear that any additional instabilities in $\delta$CE2 are dependent on the existence of an instability in the $\delta$QL equations and are also always weaker (i.e. their growths rates are smaller) than this QL instability. By exactly similar reasoning, there is also the possibility of additional unstable type B perturbations in $\delta$CE2. However, here it is less straightforward to deduce their reliance on a QL instability being present or that this dominates their growth rates due to the mean flow perturbation feeding back onto the fluctuation equations. Nevertheless, this seems a plausible assumption: if it is not true, CE2 could exhibit a type B bifurcation when the QL problem does not. Computations to illustrate some of these features have now been done (Nivarti *et al.* Reference Nivarti, Kerswell, Marston and Tobias2022).

In summary, instability in the QL problem is sufficient to conclude the same for the larger CE2 problem whereas concluding stability of CE2 given the stability of the QL problem is true for type A perturbations but remains a conjecture for type B perturbations.

## Appendix B. Solving the eigenvalue problem

Here we explain how we implemented mEOS and EOS calculations as an eigenvalue problem in 2-D channel flow. We expand the fluctuation field in terms of sine and cosine modes,

The nonlinear term in the mean equation then becomes

The EOS equations then take the form

where the scalar variable $\delta G$ is used to impose a zero-flux condition on the mean velocity perturbation $\delta U$: $\int ^1_{-1} \,\delta U \, {\textrm {d}\kern 0.05em y} =0$ and, for each wavenumber $(m,0)$,

The eigenvalue problem is discretised using Chebyshev extreme points with the same resolution as in base flow simulations. The no-slip boundary conditions on the velocity variables are imposed by replacing the appropriate rows of the matrix.

This eigenvalue problem is real and so eigenvalues are either real or come in complex conjugate pairs where a real physical velocity field is formed by addition,