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

(2.6)\begin{gather} -{\rm i}m \alpha c\, \tilde{u} = \frac{1}{Re}(\tilde{u}_{yy}-m^2\alpha^2 \tilde{u})-{\rm i}m \alpha \tilde{p} -{\rm i}m \alpha U\tilde{u}-U_y \tilde{v}, \end{gather}

(2.7)\begin{gather}-{\rm i}m \alpha c \,\tilde{v} = \frac{1}{Re}(\tilde{v}_{yy}-m^2 \alpha^2\tilde{v})-\tilde{p}_y -{\rm i}m \alpha U\tilde{v}, \end{gather}

(2.8)\begin{gather}0 = {\rm i}m \alpha \, \tilde{u}+\tilde{v}_y , \end{gather}

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

(2.9)\begin{equation} (U-c)(\partial_y^2-m^2 \alpha^2)\psi-U_{yy} \psi=\frac{1}{{\rm i} m \alpha Re} (\partial_y^2-m^2 \alpha^2)^2 \psi. \end{equation}

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),

(2.10)\begin{align} \frac{1}{Re} [(1+E(y))U_y]_y&:=\frac{1}{Re}U_{yy}- (\overline{\tilde{u} \tilde{v}} )_y={-}G \Rightarrow E(y) \nonumber\\ &:={-}Re \frac{ \overline{ \tilde{u} \tilde{v}} }{U_y} ={-}\frac{Re}{U_y} \int^y_0 G(\bar{y})\,{\rm d}\bar{y}-1. \end{align}

Then, the eigenproblem can be modified to

(2.11)\begin{gather} -{\rm i}m \alpha c\, \tilde{u} = \frac{1}{Re} \left( \left[(1+E)\tilde{u}_{y} \right]_y-m^2\alpha^2 (1+E) \tilde{u} \right)-{\rm i}m \alpha \tilde{p} -{\rm i}m \alpha U\tilde{u}-U_y \tilde{v}, \end{gather}

(2.12)\begin{gather}-{\rm i}m \alpha c \,\tilde{v} = \frac{1}{Re}\left( \left[(1+E)\tilde{v}_{y} \right]_y-m^2\alpha^2 (1+E) \tilde{v} \right)-\tilde{p}_y -{\rm i}m \alpha U\tilde{v}, \end{gather}

(2.13)\begin{gather}0 = {\rm i}m \alpha \, \tilde{u}+\tilde{v}_y. \end{gather}

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

(2.14)\begin{equation} \boldsymbol{C}(y_1,y_2,t):= \overline{\boldsymbol{\tilde{u}}(x,y_1,z,t) \otimes \boldsymbol{\tilde{u}}(x,y_2,z,t)} =\left(\begin{array}{@{}ccc@{}} C_{11} & C_{12} & C_{13} \\ C_{21} & C_{22} & C_{23} \\ C_{31} & C_{32} & C_{33} \end{array} \right), \end{equation}

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

(2.15)\begin{equation} C_{ijk}^{(3)}(1,2,3):=\overline{[\boldsymbol{\tilde{u}}(x,y_1,z,t)]_i [\boldsymbol{\tilde{u}}(x,y_2,z,t)]_j [\boldsymbol{\tilde{u}}(x,y_3,z,t)]_k}. \end{equation}

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.

(2.16)\begin{align} C_{ijkl}^{(4)}(1,2,3,4)&:= \overline{ [\boldsymbol{\tilde{u}}(x,y_1,z,t)]_i [\boldsymbol{\tilde{u}}(x,y_2,z,t)]_j [\boldsymbol{\tilde{u}}(x,y_3,z,t)]_k [\boldsymbol{\tilde{u}}(x,y_4,z,t)]_l } \end{align}

(2.17)\begin{align} & \quad -C_{ij}(1,2)C_{kl}(3,4) -C_{ik}(1,3)C_{jl}(2,4)-C_{il}(1,4)C_{jk}(2,3). \end{align}

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

(2.18)\begin{equation} \boldsymbol{\tilde{u}}(x,y,z,t):= \sum_m \sum_n \boldsymbol{\tilde{u}}^{mn}(y,t) \exp({{\rm i}m \alpha x+{\rm i}n \beta z}), \end{equation}

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

(2.19)\begin{align} C_{ij}(1,2) &= \sum_m \sum_n \left\{ C_{ij}^{mn}(1,2):= [\boldsymbol{\tilde{u}}^{mn} (y_1,t)]_i [\boldsymbol{\tilde{u}}^{{-}m-n} (y_2,t)]_j \right\}, \end{align}

(2.20)\begin{align}C^{(3)}_{ijk}(1,2,3) &= \sum_m \sum_n \left\{ C_{ijk}^{(3)mn}(1,2,3) :=\sum_p \sum_q [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i [\boldsymbol{\tilde{u}}^{pq}(y_2,t)]_j\right. \nonumber\\ &\quad \left. \times [\boldsymbol{\tilde{u}}^{-(m+p)-(n+q)}(y_3,t)]_k \vphantom{\sum_p}\right\} \end{align}

(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,

(2.21)\begin{align} &\partial_t C_{ij}^{mn}(1,2) = [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i \partial_t [\boldsymbol{\tilde{u}}^{{-}m-n}(y_2,t)]_j +\partial_t [\boldsymbol{\tilde{u}}^{mn}(y_1,t)]_i [\boldsymbol{\tilde{u}}^{{-}m-n}(y_2,t)]_j \nonumber\\ &\quad = \frac{1}{Re} \left(\partial^2_1+\partial^2_2 -2m^2\alpha^2-2n^2\beta^2\right) C_{ij}^{mn}(1,2) \nonumber\\ & \qquad -\left[\begin{array}{@{}c@{}} -{\rm i}m\alpha \\ \partial_2 \\ -{\rm i}n\beta \end{array}\right]_j C_{i4}^{mn}(1,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_i C^{mn}_{4j}(1,2) \nonumber\\ &\qquad +{\rm i}m\alpha [U(2)-U(1)]C_{ij}^{mn}(1,2) -U_y(2)C_{i2}^{mn}(1,2) \delta_{1j} \nonumber\\ &\qquad -U_y(1)C_{2j}^{mn}(1,2) \delta_{i1} -\left[\begin{array}{@{}c@{}} -{\rm i}m\alpha \\ \partial_2 \\ -{\rm i}n\beta \end{array}\right]_k C_{ijk}^{(3)mn}(1,2,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_k C_{ikj}^{(3)-m-n}(2,1,1), \end{align}

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

(2.22)\begin{equation} \left[\begin{array}{@{}c@{}} {\rm i}m \alpha\\ \partial_1 \\ {\rm i}n \beta \end{array}\right]_i C_{ij}^{mn}(1,2)=0 \quad j \in \{1,2,3\} \end{equation}

along with the equation for the mean equation (2.3)

(2.23)\begin{equation} U_t -\frac{1}{Re}U_{yy} = G - \sum_m \sum_n \partial_y C_{12}^{mn}(y,y,t) \end{equation}

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),

(2.24)\begin{gather} U_t = \frac{1}{Re}U_{yy}+G - ( \overline{\tilde{u} \tilde{v}} )_y, \end{gather}

(2.25)\begin{gather}\boldsymbol{\tilde{u}}_t = \frac{1}{Re} \nabla^2 \boldsymbol{\tilde{u}} -\boldsymbol{\nabla} \tilde{p} -{U} \boldsymbol{\tilde{u}}_x - \tilde{v} U_y \boldsymbol{\hat{x}}, \end{gather}

(2.26)\begin{gather}0 =\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tilde{u}} . \end{gather}

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

(2.27)\begin{align} \partial_t \,\delta C_{ij}^{mn}(1,2) &= \frac{1}{Re} \left(\partial^2_1+\partial^2_2 -2m^2\alpha^2-2n^2\beta^2\right) \delta C_{ij}^{mn}(1,2) \nonumber\\ &\qquad -\left[\begin{array}{@{}c@{}} -{\rm i}m \alpha\\ \partial_2 \\ -{\rm i}n \beta \end{array}\right]_j \delta C_{i4}^{mn}(1,2) -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_1 \\ {\rm i}n\beta \end{array}\right]_i \delta C^{mn}_{4j}(1,2) \nonumber\\ &\qquad +{\rm i}m \alpha [U(2)-U(1)] \delta C_{ij}^{mn}(1,2) -U_y(2) \delta C_{i2}^{mn}(1,2) \delta_{1j} \nonumber\\ &\qquad -U_y(1) \delta C_{2j}^{mn}(1,2) \delta_{i1} \nonumber\\ &\qquad +{\rm i}m \alpha [\delta U(2)- \delta U(1)] C_{ij}^{mn}(1,2) -\delta U_y(2) C_{i2}^{mn}(1,2) \delta_{1j}\nonumber\\ &\qquad -\delta U_y(1) C_{2j}^{mn}(1,2) \delta_{i1}, \end{align}

(2.28)\begin{align} 0 &= \left[\begin{array}{@{}c@{}} {\rm i}m\alpha\\ \partial_1 \\ {\rm i}n \beta \end{array}\right]_i \delta C_{ij}^{mn}(1,2), \end{align}

(2.29)\begin{align} \partial_t \delta U &= \frac{1}{Re}\delta U_{yy} +\delta G - \sum_m \sum_n \partial_y \delta C_{12}^{mn}(y,y,t). \end{align}

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,

(2.30)\begin{align} &\|\boldsymbol{C}^{mn}-\boldsymbol{\tilde{u}}^{mn}_0 \otimes \boldsymbol{\tilde{u}}^{*mn}_0\|^2_F := \sum_{i,j=1}^3 \sum_{p,q=1}^N (C^{mn}_{ij}(p,q)-[\boldsymbol{\tilde{u}}^{mn}_0(y_p)]_i[\boldsymbol{\tilde{u}}^{*mn}_0(y_q)]_j) \nonumber\\ &\qquad \times (C^{*mn}_{ij}(p,q)-[\boldsymbol{\tilde{u}}^{*mn}_0(y_p)]_i[\boldsymbol{\tilde{u}}^{mn}_0(y_q)]_j) \nonumber\\ &\quad =\left\|\left(\begin{array}{@{}ccc@{}} C^{mn}_{11} & C^{mn}_{12} & C^{mn}_{13} \\ C^{mn}_{21} & C^{mn}_{22} & C^{mn}_{23} \\ C^{mn}_{31} & C^{mn}_{32} & C^{mn}_{33} \end{array} \right) - \left[\begin{array}{@{}c@{}} \tilde{u}^{mn}_{(1:N)}\\ \tilde{v}^{mn}_{(1:N)} \\ \tilde{w}^{mn}_{(1:N)} \end{array} \right] \left[\begin{array}{@{}ccc@{}} \tilde{u}^{*mn}_{(1:N)} & \tilde{v}^{*mn}_{(1:N)} & \tilde{w}^{*mn}_{(1:N)} \end{array} \right] \right\|^2_F, \end{align}

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

(2.31)\begin{gather} \partial_t \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} = {\mathbb{L}}^{mn}(U) \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} +{\mathbb{J}}^{mn}(\delta U) \boldsymbol{\tilde{u}}_0^{mn} \quad \forall (m,n) \neq (0,0), \end{gather}

(2.32)\begin{gather}\partial_t \delta U = \frac{1}{Re}\delta U_{yy} +\delta G - \sum_m \sum_n \partial_y \left\{ \tilde{u}_0^{mn} \delta \tilde{v}^{{-}m-n} +\delta \tilde{u}^{mn} \tilde{v}_0^{{-}m-n} \right\}, \end{gather}

where

(2.33)\begin{gather} {\mathbb{L}}^{mn}(U)\, \boldsymbol{\tilde{u}}^{mn} := \frac{1}{Re} [\partial_y^2-m^2\alpha^2 -n^2 \beta^2]\boldsymbol{\tilde{u}}^{mn} -\left[\begin{array}{@{}c@{}} {\rm i}m\alpha \\ \partial_y \\ {\rm i}n\beta \end{array} \right] \tilde{p}^{mn} -{\rm i}m \alpha{U} \boldsymbol{\tilde{u}}^{mn} - \tilde{v}^{mn} U_y \boldsymbol{\hat{x}}, \end{gather}

(2.34)\begin{gather}{\mathbb{J}}^{mn}(\delta U) :=\lim_{\epsilon \rightarrow 0} \, \frac{1}{\epsilon}[{\mathbb{L}}^{mn}(U+ \epsilon\delta U)-{\mathbb{L}}^{mn}(U)]={\mathbb{L}}^{mn}(U+\delta U)-{\mathbb{L}}^{mn}(U) \end{gather}

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

(2.35)\begin{align} &(\delta U(y,t), \boldsymbol{\delta} \boldsymbol{\tilde{u}}(x,y,z,t)) \nonumber\\ &\quad =\left( \delta \hat{U}(y) {\rm e}^{\lambda t}, \,\sum_m \sum_n \left[ \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn}(y,t):=\boldsymbol{\delta} \widehat{\boldsymbol{u}}^{mn}(y) {\rm e}^{\lambda t} \right]\exp({{\rm i}(m \alpha x+ n \beta z)}) \right). \end{align}

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,

(2.36)\begin{equation} (\delta \hat{U}, \,\boldsymbol{\delta} \widehat{\boldsymbol{C}}^{mn}) \,{\rm e}^{\lambda t}. \end{equation}

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

(2.37)\begin{equation} (\delta \hat{U}, \,\boldsymbol{\delta} \boldsymbol{\hat{u}}^{mn},\,\boldsymbol{\delta} \boldsymbol{\hat{u}}^{{-}m-n}) \,{\rm e}^{\lambda t} \end{equation}

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

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

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

(3.5)\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 \alpha U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{mn} - \delta \tilde{v}^{mn}U_y\boldsymbol{\hat{x}} -{\rm i}m \alpha {\delta U} \boldsymbol{\tilde{u}}_0^{mn} -\tilde{v}_0^{mn} \delta U_y \boldsymbol{\hat{x}} \nonumber\\ -[\boldsymbol{\delta} \boldsymbol{\tilde{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0+\boldsymbol{\tilde{u}}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}]^{mn}, \nonumber\\ {\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}

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

Figure 3. Fluctuation energy as a ratio to the total flow energy for the four 2-D channel test states at $Re=36\,300$: symmetric S (light blue), asymmetric A (pink), F1 (blue) and F2 (red). (a) Decomposition by streamwise wavenumber $k$. (b) Temporal fluctuations over a typical time-averaging window.

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.

Figure 4. Typical eigenvalue convergence with respect to the total number of the streamwise wavenumbers $N_x$ included in the EOS model, shown for $k=2$ (a) and $k=3$ (b) for the F2 state.

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 5a) and asymmetric (figure 5b) 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.

Figure 5. Comparison of the eigenvalues $\sigma =-\textrm {i}m\alpha (c_r+\textrm {i}\,c_i)$ – so growth rate is $m \alpha c_i$ – for OS ($\bullet$), OS with eddy viscosity ($\square$), EOS ($\Diamond$) and mEOS ($\bigcirc$) analyses. Only the leading eigenvalues for $m=1$ (blue), $m=2$ (red) and $m=3$ (green) are shown. (a) Symmetric S, (b) asymmetric A, (c) F1, (d) F2 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.

Figure 6. Power spectra for the EOS eigenvectors for the leading eigenvalues. (a) Symmetric S state ($\times$) and asymmetric A state ($\bullet$, red) with $m=1$; (b) states F1 ($\blacksquare$, cyan) and F2 ($+$, blue) with $m=2$. The mean velocity perturbation ($m=0$) is at least an order of magnitude smaller than the leading wavenumber contribution to the power spectrum in all 4 cases.

The leading eigenvalue for the body-forced states states F1 (figure 5c) and F2 (figure 5d) 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

(4.4)\begin{equation} A \boldsymbol{y_R} = \sigma B \boldsymbol{y_R} \end{equation}

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

(4.5)\begin{equation} \left(A + {\delta A_A} +\delta A_S \right)\left(\boldsymbol{y_R} + \boldsymbol{\delta y_R } \right) = \left( \sigma +\delta \sigma \right) B \left( \boldsymbol{y_R} +\boldsymbol{\delta y_R} \right), \end{equation}

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

(4.6)\begin{equation} {\delta \sigma} = {\delta \sigma_A} + {\delta \sigma_S}:= \frac{\boldsymbol{y_L}^{{\dagger}}{\delta A_A\boldsymbol{y_R}}}{\boldsymbol{y_L}^{{\dagger}}B\boldsymbol{y_R}}+\frac{\boldsymbol{y_L}^{{\dagger}}{\delta A_S\boldsymbol{y_R}}}{\boldsymbol{y_L}^{{\dagger}}B\boldsymbol{y_R}}. \end{equation}

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.

Figure 7. Perturbation analysis results for F2 state. (a–c) Streamwise component of the unperturbed left eigenvector $\delta u_{sL}$, streamwise component of the perturbed advection term $-\textrm {i}m\alpha u_{0s}^{m0} \delta U_{sR}$ and streamwise component of the perturbed shear term $-v_{0s}^{m0} \partial _y \delta U_{sR}$. (d–f) Wall-normal component of the unperturbed left eigenvector $\delta v_{sL}$, wall-normal component of the perturbed advection term $-\textrm {i}m\alpha v_{0s}^{m0} \delta U_R$ and mean profile component of the unperturbed right eigenvector $-\delta U_{sR}$. Solid (dashed) lines indicate real (imaginary) parts of the vectors. Subscript $(\boldsymbol {\cdot })_s$ refers to the sine components of the eigenfunctions as explained in Appendix B.

Figure 8. Perturbation analysis results for S state. (a–c) Streamwise component of the unperturbed left eigenvector $\delta u_{cL}$, streamwise component of the perturbed advection term $-\textrm {i}m\alpha u_{0c}^{m0} \delta U_{cR}$ and streamwise component of the perturbed shear term $-v_{0c}^{m0} \partial _y \delta U_{cR}$. (d–f) Wall-normal component of the unperturbed left eigenvector $\delta v_{cL}$, wall-normal component of the perturbed advection term $-\textrm {i}m\alpha v_{0c}^{m0} \delta U_R$ and mean profile component of the unperturbed right eigenvector $-\delta U_{cR}$. Solid (dashed) lines indicate real (imaginary) parts of the vectors. Subscript $(\boldsymbol {\cdot })_c$ refers to the cosine components of the eigenfunctions as explained in Appendix B.

Table 1. Perturbation analysis results for the leading eigenvalue in the F2 state (left) and the symmetric S state (right). The first-order change in the eigenvalue caused by the advection ($\sigma_A$) and shear ($\sigma_S$) terms in the mEOS equations. Unperturbed (standard OS) eigenvalues are given by $\sigma _{OS}$ and mEOS eigenvalues are given by $\sigma _{mEOS}$. The table shows that F2 changes to being stable when using mEOS analysis, whereas S remains unstable.

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

Figure 9. Leading singular values $\sigma _i$ of the correlation matrices $\boldsymbol {C}^{m0}$ time averaged over 1000 time units. Shown for $m=1$ (blue), $m=2$ (red) and $m=3$ (green). (a) Symmetric S, (b) asymmetric A, (c) F1, (d) F2 states. Note the gap after the leading singular value for $m=2$ (red) across all states confirming rank 1 approximation of the correlation matrix for the cases when mEOS and EOS analyses show stabilisation.

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,

(5.1)\begin{equation} \boldsymbol{\tilde{u}}^{mn}_0=\sum_{i=1}^{N} \frac{\sigma_i}{\sum_{i=1}^{N} \sigma_i} \,\sqrt{\sigma_i} \boldsymbol{\hat{e}_i}, \end{equation}

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

(5.2)\begin{gather} \partial_t\delta U=\frac{1}{Re}\partial^2_y \delta U+\delta G -[\boldsymbol{\tilde{u}}_0(t).\boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}+ \boldsymbol{\delta}\boldsymbol{\tilde{u}}.\boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0(t)]^{00}, \end{gather}

(5.3)\begin{align}\partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{m0} &= \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0}-\boldsymbol{\nabla} \delta \tilde{p}^{m0} - {\rm i}m U \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0} - \delta \tilde{v}^{m0}U_y\boldsymbol{\hat{x}}\nonumber\\ &\quad -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{m0} -\tilde{v}_0^{m0} \delta U_y \boldsymbol{\hat{x}}, \quad m>0, \end{align}

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.

Figure 10. Results of various numerical experiments: time evolution of the perturbation energy $E_{pert}$ as a ratio to the base flow energy $E_{base}$ averaged over multiple DNS runs. Perturbations were time stepped using linearised QL (blue), LNS (orange) or full nonlinear Navier–Stokes (green) equations. For each case, base flow was time dependent and generated using full Navier–Stokes DNS. Results shown for the symmetric S (a) and body-force-driven BF2 (b) states. Black dashed lines are provided as guides to emphasize the zero gradient of the saturated states and to indicate the growth rate ${\rm Re}(\sigma )$ of the unstable states. For comparison, the real part of the leading mEOS eigenvalue ${\rm Re}(\sigma _{mEOS})$ is given for each of the states.

Figure 11. Comparison of the leading eigenvectors from the EOS (a), mEOS $m=1$ (b) analyses and an instantaneous snapshot of the velocity perturbation field in the QL experiment at $t=100$ (c) for the symmetric (S) state. Only the streamwise velocity component is shown. There is a good correspondence in the leading eigenvector across EOS and mEOS analyses and the QL experiment using a time-dependent base flow.

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,

(5.4)\begin{gather} \partial_t{\boldsymbol{\delta} \boldsymbol{\tilde{u}}}^{m0} = \frac{1}{Re} \nabla^2 \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0}-\boldsymbol{\nabla} \delta \tilde{p}^{m0} - {\rm i}m U(t) \boldsymbol{\delta} \boldsymbol{\tilde{u}}^{m0} - \delta \tilde{v}^{m0}U_y(t)\boldsymbol{\hat{x}} -{\rm i}m{\delta U} \boldsymbol{\tilde{u}}_0^{m0}(t)\nonumber\\ -\tilde{v}_0^{m0}(t) \delta U_y \boldsymbol{\hat{x}} -[\boldsymbol{\tilde{u}}_0(t).\boldsymbol{\nabla} \boldsymbol{\delta} \boldsymbol{\tilde{u}}+ \boldsymbol{\delta}\boldsymbol{\tilde{u}}.\boldsymbol{\nabla} \boldsymbol{\tilde{u}}_0(t)]^{m0}. \end{gather}

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 10a). 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.

Figure 12. Comparison of the leading eigenvectors from the mEOS $m=1$ (a) and mEOS $m=3$ (b) analyses with the instantaneous snapshots of the perturbation fields obtained from the QL experiment at $t=100$ (c) and linear NS experiment at $t=100$ (d) for the F2 state. The structures of the leading mEOS eigenvectors seem unrelated to the growing perturbations in the QL or LNS regimes.

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.

Figure 13. Time evolution of the relative difference between correlation matrices of the full flow $C^{m0}$ and base flow $C^{m0}_0$ shown for $m=1$ (blue), $m=2$ (red) and $m=3$ (green). The correlation matrices were time averaged over a rolling time window to display the statistical behaviour. Shown for the full nonlinear NS experiment where the difference slowly approaches $0$ signifying the recovery of the same statistical state.

Figure 14. Snapshots of the streamwise velocity component of the perturbation field taken from the full nonlinear Navier–Stokes experiment with the F2 state. At time $t = 10$ the unstable structure is still growing until time $t = 20$ when it reaches its maximum size and starts to spread out through the channel, at time $t = 30$ the perturbation saturates and shows no qualitative change at a much later time $t = 200$. Results are shown for (a) $t = 10$, (b) $t = 20$, (c) $t = 30$, (d) $t = 200$.