## 1. Introduction

Statically stable density stratification in wall-bounded shear flows plays an important role in many industrial and environmental applications, e.g. in cooling equipment (Zonta & Soldati Reference Zonta and Soldati2018), and the turbulent boundary layers governing atmospheric and oceanic flows (Pedlosky Reference Pedlosky2013; Vallis Reference Vallis2017). In the atmospheric boundary layer, stable stratification arising from strong ground cooling effects is of particular importance at night (Nieuwstadt Reference Nieuwstadt1984; Mahrt Reference Mahrt1999, Reference Mahrt2014) and near the polar region (Grachev *et al.* Reference Grachev, Fairall, Persson, Andreas and Guest2005). At the ocean floor, stable density stratification is known to influence the boundary layer thickness (Weatherly & Martin Reference Weatherly and Martin1978; Lien & Sanford Reference Lien and Sanford2004).

(Stably) stratified plane Couette flow (PCF) is a canonical model for stratified wall-bounded shear flow. When the density as well as the velocity is maintained at different values at the two horizontal boundary planes, with gravity acting vertically, stratified PCF has the added benefit (as defined more precisely below) that a natural bulk Richardson number $Ri_b$ can be defined, thereby enabling the relative significance of the imposed stratification and shear to be captured. Furthermore, unstratified PCF has no linear instability for any Reynolds number ($Re$, again defined more precisely below) (Romanov Reference Romanov1973), and yet is observed to transition at Reynolds numbers as low as $Re=360\pm 10$ (Tillmark & Alfredsson Reference Tillmark and Alfredsson1992). Stratified PCF is a convenient model flow for investigating the effect of stable stratification on the transition dynamics (Deusebio, Caulfield & Taylor Reference Deusebio, Caulfield and Taylor2015).

Stable stratification provides a restoring buoyancy force inhibiting vertical motion (Turner Reference Turner1979; Davidson Reference Davidson2013). Thus, transition to turbulence in stably stratified PCF typically occurs at a higher Reynolds number than unstratified PCF; see e.g. Deusebio *et al.* (Reference Deusebio, Caulfield and Taylor2015); Eaves & Caulfield (Reference Eaves and Caulfield2015); Deguchi (Reference Deguchi2017) and Olvera & Kerswell (Reference Olvera and Kerswell2017). In the transitional regime, both stratified PCF and unstratified PCF exhibit spatial intermittency; i.e. the coexistence of laminar and turbulent regions. In the relatively low-$Re$ low-$Ri_b$ intermittent regime, the spatial intermittency in stratified PCF is characterized by oblique turbulent bands (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016) at least qualitatively similar to those seen in unstratified PCF (Prigent *et al.* Reference Prigent, Grégoire, Chaté and Dauchot2003; Duguet, Schlatter & Henningson Reference Duguet, Schlatter and Henningson2010) with a very large channel size (${\sim }O(100)$ times the channel half-height). In the high-$Re$ high-$Ri_b$ intermittent regime, flow structures are instead characterized by turbulent and laminar layers over the vertical direction due to the strong effect of buoyancy (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015). This spatial intermittency directly imposes challenges for the computation of averaged measurements of flow behaviour (such as the efficiency of mixing or the dissipation rate), and thus an understanding of the underlying mechanisms is important for the parameterization of turbulence properties, in particular the irreversible mixing in stratified flows (Caulfield Reference Caulfield2020, Reference Caulfield2021).

The existence of a unique critical Richardson number that separates flow into laminar and turbulent regimes is questionable, to put it mildly (Andreas Reference Andreas2002; Galperin, Sukoriansky & Anderson Reference Galperin, Sukoriansky and Anderson2007). A threshold value close to $1/4$ is supported by some field measurements (Kundu & Beardsley Reference Kundu and Beardsley1991) and experiments (Rohr *et al.* Reference Rohr, Itsweire, Helland and Van Atta1988), although other field measurements reported sustained turbulence in flows with Richardson numbers ${\simeq }1$ (Lyons, Panofsky & Wollaston Reference Lyons, Panofsky and Wollaston1964). More recently, increasing evidence has been found that vertically sheared stably stratified flow appears to self-organize to maintain an appropriately defined Richardson number near $1/4$, both in field observations (Smyth & Moum Reference Smyth and Moum2013; Smyth, Nash & Moum Reference Smyth, Nash and Moum2019) and in simulations (Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018). This threshold value of $1/4$ also appears in the classical ‘Miles–Howard’ theorem (Howard Reference Howard1961; Miles Reference Miles1961), which provides a necessary condition for linear instability in inviscid, non-diffusive steady parallel flow. In particular, it states that instability requires that the local or gradient Richardson number be less than $1/4$ somewhere. Therefore, it is of interest to consider stratified PCF as a well-controlled sheared stratified flow to investigate whether some kind of self-organized criticality and/or marginal stability naturally emerges in a viscous and diffusive flow.

The Prandtl number ($Pr=\nu /\kappa$, where $\nu$ is the kinematic viscosity and $\kappa$ is the diffusivity of the density field) plays a perhaps unsurprisingly important role in determining flow structures. For example, for sufficiently small $Pr$, flows with the same value of the product $Pr Ri_b$ develop the same averaged vertical density profile (Langham, Eaves & Kerswell Reference Langham, Eaves and Kerswell2020). This observation that the product $PrRi_b$ determines flow behaviour in the low Prandtl number limit is widely observed in stratified shear flows; see e.g. Lignieres (Reference Lignieres1999); Garaud, Gallet & Bischoff (Reference Garaud, Gallet and Bischoff2015); Garaud, Gagnier & Verhoeven (Reference Garaud, Gagnier and Verhoeven2017) and Garaud (Reference Garaud2021). Conversely, in the high Prandtl number regime, studies of exact coherent structures in stratified PCF (Langham *et al.* Reference Langham, Eaves and Kerswell2020) show that a nearly uniform density region forms near the channel centre, and the influence of bulk Richardson number on the averaged properties of the flow are significantly reduced. Moreover, Taylor & Zhou (Reference Taylor and Zhou2017) proposed a multi-parameter criterion for the formation of a ‘staircase’ in the density distribution (i.e. a distribution with relatively deep ‘layers’ of nearly uniform density separated by relatively thin interfaces of enhanced density gradient) which suggests that this staircase formation is actually favoured for a large Prandtl number (Taylor & Zhou Reference Taylor and Zhou2017). The sharpness of the density interfaces also appears to increase as the Prandtl number increases (Zhou *et al.* Reference Zhou, Taylor, Caulfield and Linden2017*b*). In addition, increasing the Prandtl number has a larger influence on the mean density profiles than on the mean velocity profiles (Zhou, Taylor & Caulfield Reference Zhou, Taylor and Caulfield2017*a*).

The oblique turbulent bands observed in the intermittent regime of stratified PCF (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016) require a very large channel size to accommodate them fully, which poses challenges for both simulations and experiments. The three different flow parameters of interest, $Re$, $Ri_b$ and $Pr$ also lead to computational challenges in exploring the full range of flow regimes. To overcome these challenges to direct numerical simulation (DNS), we use an input–output (resolvent) analysis based approach. Such methods, built upon the spatio-temporal frequency response, have been widely employed in unstratified wall-bounded shear flows (Farrell & Ioannou Reference Farrell and Ioannou1993*a*; Bamieh & Dahleh Reference Bamieh and Dahleh2001; Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010; McKeon Reference McKeon2017). This analysis framework has the advantage of computational tractability and is not subject to finite channel effects. Related analysis has shown promise in studying stratified flows, including inviscid stratified shear flow with constant shear (Farrell & Ioannou Reference Farrell and Ioannou1993*b*), stratified PCF (Jose *et al.* Reference Jose, Roy, Bale and Govindarajan2015, Reference Jose, Roy, Bale, Iyer and Govindarajan2018) and stratified turbulent channel flow (Ahmed *et al.* Reference Ahmed, Bae, Thompson and McKeon2021).

In this work, we extend the structured input–output analysis (SIOA) originally developed for unstratified PCF (Liu & Gayme Reference Liu and Gayme2021) to stratified PCF. Prior application of the SIOA approach to unstratified transitional wall-bounded shear flows (Liu & Gayme Reference Liu and Gayme2021) demonstrated that including the componentwise structure of the nonlinearity uncovers a wider range of known key flow features identified through nonlinear analysis, experiments and DNS, but not captured through traditional (unstructured) input–output approaches. Here, SIOA for stratified PCF includes the effect of nonlinearity in the momentum and density equations (under the Boussinesq approximation) within a computationally tractable linear framework through a feedback interconnection between the linearized dynamics and a structured forcing that is explicitly constrained to preserve the componentwise structure of the nonlinearity. The structured singular value (Doyle Reference Doyle1982; Safonov Reference Safonov1982) of the spatio-temporal frequency response associated with this feedback interconnection can then be calculated at each streamwise and spanwise length scale. This structured singular value can be interpreted as the flow structure that shows the largest input–output gain (amplification) given the structured feedback interconnection.

Here, we apply the SIOA to characterize highly amplified flow structures in the intermittent regime of stratified PCF and investigate the behaviour of the flow across a range of $Re$, $Ri_b$ and $Pr$. Our aims are twofold. First, we wish to investigate whether the structures predicted by the SIOA can be quantitatively identified with fully nonlinear structures that have been observed in previously reported DNS of stratified PCF with specific values of the control parameters $Re$, $Ri_b$ and $Pr$. Second, we wish to explore the dependence on the control parameters of predictions from the SIOA in parameter regimes which are not (as yet) accessible to DNS. More specifically, to address our first aim, we examine how $Re$ and $Ri_b$ affect flow structures with Prandtl number set at $Pr=0.7$, i.e. the value for air. We demonstrate that SIOA does indeed predict the characteristic wavelengths and angle of the oblique turbulent bands observed in very large channel size DNS of the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF at the same $Pr$ (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016). We further show that, in the high-$Re$ high-$Ri_b$ intermittent regime, the SIOA identifies quasi-horizontal flow structures resembling turbulent–laminar layers (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015).

Having achieved our first aim, and demonstrated the usefulness of the SIOA for identifying realistic nonlinear flow structures, we then turn our attention to our second aim. We demonstrate that increasing the bulk Richardson number reduces the amplification of streamwise-varying flow structures. These results show that the classical Miles–Howard stability criterion ($Ri_b\leq 1/4$) appears (perhaps fortuitously) to be associated with a change in the most amplified flow structures, which is robust for a wide range of $Re$ and valid at $Pr\approx 1$.

We then examine flow behaviour at different $Ri_b$ and $Pr$. For flows with $Pr\ll 1$, a larger bulk Richardson number is required to reduce the amplification of streamwise-varying flow structures to the same level as streamwise-independent ones compared with $Pr\approx 1$. The largest amplification also is predicted to occur at the same value of the product $Pr Ri_b$ consistent with the observation of the averaged density profile only varying with the product $Pr Ri_b$ in the $Pr\ll 1$ regime (Langham *et al.* Reference Langham, Eaves and Kerswell2020). For flows with $Pr\gg 1$, the SIOA identifies another quasi-horizontal flow structure independent of $Ri_b$. By decomposing input–output pathways into separate components associated with velocity and density fluctuations, we show that these quasi-horizontal flow structures at $Pr\gg 1$ are primarily associated with fluctuations in the density field. We further highlight the importance of this density-associated flow structure at $Pr\gg 1$ by constructing an analytical scaling argument for the input–output amplification in terms of $Re$ and $Pr$ under the assumptions of unstratified flow (with $Ri_b=0$) and streamwise invariance. The above observations using SIOA distinguish two types of quasi-horizontal flow structures, one associated with the high-$Re$ high-$Ri_b$ regime and the other one associated with density perturbations that emerges in the high $Pr$ regime.

To achieve our twin aims, and to demonstrate the above summarized results, the remainder of this paper is organized as follows. Section 2 describes the flow configuration of stratified PCF and then develops the SIOA for this flow. Section 3 analyses the results obtained from SIOA focusing on the wall-parallel length scale of flow structures. In § 4, we develop analytical scaling arguments with respect to $Re$ and $Pr$ to investigate behaviour for flows in the high $Pr$ limit. Finally we draw conclusions and suggest some avenues of future work in § 5.

## 2. Structured input–output response of stratified flow

### 2.1. Governing equations

We consider stably stratified PCF between two infinite parallel plates and employ $x$, $y$ and $z$ to denote the streamwise, wall-normal (or vertical) and spanwise directions. The corresponding (assumed incompressible) velocity components are denoted as $u$, $v$ and $w$. The coordinate frames and configurations for this stratified PCF are shown in figure 1. We express the velocity field as a vector $\boldsymbol {u}_{{tot}}=[u_{{tot}} \enspace v_{{tot}} \enspace w_{{tot}}]^\text {T}$ with ${}^\text {T}$ indicating the transpose. We then decompose the velocity field into the sum of a laminar linearly varying base flow $U(y)=y$ and fluctuations about the base flow $\boldsymbol {u}$; i.e. $\boldsymbol {u}_{{tot}}=U(y)\boldsymbol {e}_x+\boldsymbol {u}$ with $\boldsymbol {e}_x$ denoting the $x$-direction (streamwise) unit vector. The pressure field is similarly decomposed as $p_{{tot}}=P+p$. We decompose the density $\rho _{{tot}}$ as the sum of a reference density $\rho _r$, a base, again linearly varying, density $\bar {\rho }=-y$ and a density fluctuation $\rho$; i.e. $\rho _{{tot}}=\rho _r+\bar {\rho }+\rho$. We use $\rho _0$ to denote half of the density difference between the top and bottom walls, which is assumed to be much smaller than the reference density $\rho _0\ll \rho _r$ so that the Boussinesq approximation can be used.

The dynamics of the fluctuations $\boldsymbol {u}$, $p$ and $\rho$ is hence governed by the Navier–Stokes equations for an incompressible velocity field under the Boussinesq approximation and an advection–diffusion equation for the density

*a*)\begin{gather} \partial_{t} \boldsymbol{u} + U\partial_x \boldsymbol{u} + v\frac{{\rm d} U }{{\rm d}y}\boldsymbol{e}_x+Ri_b\rho\boldsymbol{e}_y +\boldsymbol{\nabla} p -\frac{1}{Re}{\nabla}^2 \boldsymbol{u} ={-} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}, \end{gather}

*b*)\begin{gather}\partial_t \rho+U\partial_x \rho+v\frac{{\rm d}\bar{\rho}}{{\rm d}y}-\frac{1}{RePr}\nabla^2 \rho ={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho, \end{gather}

*c*)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0. \end{gather}

Here, the spatial variables are normalized by the channel half-height $h$, and the velocity is normalized by half of the velocity difference between the top and bottom walls $U_w$, where $\pm U_w$ is the velocity at the channel walls. Time and pressure are normalized by $h/U_w$ and $\rho _r U_w^2$, respectively. The base density field $\bar {\rho }(y)$ and the density fluctuations $\rho$ are normalized by $\rho _0$. Under this normalization, the base density profile $\bar {\rho }=-y$ is balanced by a hydrostatic pressure $P=Ri_b y^2/2$.

The non-dimensional control parameters are the Reynolds number $Re$, the Prandtl number $Pr$ and the bulk Richardson number $Ri_b$, naturally defined as

*a*–

*c*)\begin{equation} Re:=\frac{U_wh}{\nu},\quad Pr:=\frac{\nu}{\kappa},\quad Ri_b:=\frac{g \rho_0 h}{\rho_r U_w^2}, \end{equation}

where $\nu$ is the kinematic viscosity, $\kappa$ is the molecular diffusivity of the density scalar and $g$ is the magnitude of gravity. The gravity vector is in the direction orthogonal to the wall $\boldsymbol {g}=-g\boldsymbol {e}_y$ with $\boldsymbol {e}_y$ denoting the $y$-direction (wall-normal, or vertical) unit vector. In equation set (2.1), $\boldsymbol {\nabla }:=[\partial _x \enspace \partial _y \enspace \partial _z]^\text {T}$ represents the gradient operator, and $\nabla ^2:=\partial _{x}^2+\partial _{y}^2+\partial _{z}^2$ represents the Laplacian operator. We impose no-slip boundary conditions at the wall $\boldsymbol {u}(y=\pm 1)=\boldsymbol {0}$ and Dirichlet boundary conditions for density fluctuations $\rho (y=\pm 1)=0$ that can be maintained by e.g. constant temperatures at the wall with a linear equation of state (with the hotter plate at the top).

A large body of linear analysis techniques views the nonlinear terms as a forcing, which enables these terms to be represented as an unmodelled effect (which can be thought of as some type of ‘uncertainty’ in the equations). There are a wide range of such models, but a common approach is a delta-correlated or coloured stochastic forcing that captures a wide range of the unmodelled effects, see e.g. the discussion in Farrell & Ioannou (Reference Farrell and Ioannou1993*a*), Bamieh & Dahleh (Reference Bamieh and Dahleh2001), Jovanović & Bamieh (Reference Jovanović and Bamieh2005), McKeon & Sharma (Reference McKeon and Sharma2010), McKeon (Reference McKeon2017) and Zare, Jovanović & Georgiou (Reference Zare, Jovanović and Georgiou2017). Here, we similarly write the nonlinear terms as the forcing

*a*)\begin{gather} \boldsymbol{f}_{\boldsymbol{u}}:={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}=\begin{bmatrix} -\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}u & -\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}v & -\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}w\end{bmatrix}^\text{T}=:\begin{bmatrix} f_x & f_y & f_z\end{bmatrix}^\text{T}. \end{gather}

*b*)\begin{gather}f_{\rho}:={-}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho, \end{gather}

which turns (2.1) into a set of linear evolution equations subject to the forcing terms $\boldsymbol {f}_{\boldsymbol {u}}$ and $f_{\rho }$.

We now construct a model of the nonlinearity, where the velocity field $-\boldsymbol {u}$ in (2.3) associated with the forcing components can be viewed as the gain operator of an input–output system in which the velocity and density gradients $\boldsymbol {\nabla } u$, $\boldsymbol {\nabla } v$, $\boldsymbol {\nabla } w$ and $\boldsymbol {\nabla }\rho$ act as the respective inputs and the forcing components $f_x$, $f_y$, $f_z$ and $f_{\rho }$ act as the respective outputs. This input–output models of the nonlinear components in the momentum and density equations (2.3), are respectively given by

*a*)\begin{gather} \boldsymbol{f}_{\boldsymbol{u},\xi}:={-}\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}=\begin{bmatrix} -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}u & -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}v & -\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}w \end{bmatrix}^\text{T}=:\begin{bmatrix} f_{x,\xi} & f_{y,\xi} & f_{z,\xi}\end{bmatrix}^\text{T}, \end{gather}

*b*)\begin{gather}f_{\rho,\xi}:={-}\boldsymbol{u}_{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}\rho. \end{gather}

Here, $-\boldsymbol {u}_{\xi }$ in (2.4) maps the corresponding velocity and density gradients into each component of the modelled forcing driving the linearized dynamics. This forcing in (2.4) is referred to as structured forcing because it preserves the componentwise structure of the nonlinear terms in (2.3). Figure 1(*b*) shows a block diagram of the feedback interconnection between the linearized dynamics and this forcing whose block-diagonal structure mirrors the nonlinear interactions in the Navier–Stokes equations, i.e. the forcing does not include terms such as $-\boldsymbol {u}\boldsymbol {\cdot}\boldsymbol {\nabla } v$ in the forcing $f_{x,\xi }$ since this term does not appear in the Navier–Stokes equations.

Although the nonlinearity in (2.3) can be written in many different ways, the current formulation leads to a straightforward and unified formulation for structured forcing in each momentum and density equation in (2.4). We next exploit this form of the equations to construct an input–output map using the structured singular value formalism (Packard & Doyle Reference Packard and Doyle1993; Zhou, Doyle & Glover Reference Zhou, Doyle and Glover1996). This map will enable us to analyse the fluctuations which are prominent in the intermittent regime.

### 2.2. Structured input–output response

We need to define the spatio-temporal frequency response $\mathcal {H}_\nabla ^S(y;k_x,k_z,\omega )$ of stratified PCF that will form the basis of the structured input–output response. We use the superscript $^S$ to distinguish this operator from its counterpart for unstratified wall-bounded shear flow (Liu & Gayme Reference Liu and Gayme2021). We employ the standard transformation to express the velocity field dynamics in (2.1) in terms of the wall-normal velocity $v$ and the wall-normal vorticity $\omega _y:=\partial _z u-\partial _x w$ (Schmid & Henningson Reference Schmid and Henningson2012). This transformation enforces the incompressibility constraint in (2.1*c*) and eliminates the pressure by construction. We exploit shift invariance in the $(x,z)$ spatial directions and assume shift invariance in time $t$, which allows us to perform the following triple Fourier transform:

where $\text {i}=\sqrt {-1}$, $\omega$ is the temporal frequency and $k_x = 2{\rm \pi} /\lambda _x$ and $k_z = 2{\rm \pi} /\lambda _z$ are the $x$ and $z$ wavenumbers, respectively. The sign of the temporal frequency $\omega$ in (2.5) is chosen for ease of employing control-oriented toolboxes in our computations.

The resulting equations describing the transformed linearized equations subject to the forcing $\left [\begin {smallmatrix}\boldsymbol {f}_{\boldsymbol {u}, \xi }\\ f_{\rho,\xi }\end {smallmatrix}\right ]$ are given by

*a*)\begin{gather} \text{i}\omega \begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}=\hat{\mathcal{A}}^S\begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}+\hat{\mathcal{B}}^S\begin{bmatrix}\hat{f}_{x,\xi}\\ \hat{f}_{y,\xi}\\ \hat{f}_{z,\xi}\\ \hat{f}_{\rho,\xi}\end{bmatrix}, \end{gather}

*b*)\begin{gather}\begin{bmatrix}\hat{u}\\ \hat{v}\\ \hat{w}\\ \hat{\rho}\end{bmatrix}=\hat{\mathcal{C}}^S\begin{bmatrix} \hat{v}\\ \hat{\omega}_y\\ \hat{\rho} \end{bmatrix}. \end{gather}

The operators in equation set (2.6) are defined as

*a*)\begin{gather} \hat{\mathcal{A}}^S(k_x,k_z):= \hat{\mathcal{M}}^{{-}1} \begin{bmatrix} -\text{i}k_xU{\widehat{\nabla}}^2+\text{i}k_xU''+\dfrac{\widehat{{\nabla}}^4}{Re} & 0 & Ri_b(k_x^2+k_z^2)\\ -\text{i}k_zU' & -\text{i}k_x U+\dfrac{\widehat{{\nabla}}^2}{Re} & 0\\ -\bar{\rho}' & 0 & -\text{i}k_x U+\dfrac{\widehat{\nabla}^2}{Re Pr} \end{bmatrix}, \end{gather}

*b*)\begin{gather}\mathcal{\hat{B}}^S(k_x,k_z):=\hat{\mathcal{M}}^{{-}1} \begin{bmatrix} -\text{i}k_x\partial_y & -(k_x^2+k_z^2) & -\text{i}k_z \partial_y & 0\\ \text{i}k_z & 0 & -\text{i}k_x & 0\\ 0 & 0 & 0 & \mathcal{I} \end{bmatrix},\quad \hat{\mathcal{M}}:=\begin{bmatrix} \widehat{{\nabla}}^2 & 0 & 0\\ 0 & \mathcal{I} & 0\\ 0 & 0 & \mathcal{I} \end{bmatrix}, \end{gather}

*c*)\begin{gather}\mathcal{\hat{C}}^S(k_x,k_z):= \frac{1}{k_x^2+k_z^2}\begin{bmatrix} \text{i}k_x\partial_y & -\text{i}k_z & 0 \\ k_x^2+k_z^2 & 0 & 0\\ \text{i}k_z \partial_y & \text{i}k_x & 0 \\ 0 & 0 & k_x^2+k_z^2 \end{bmatrix}, \end{gather}

where $U':={\rm d} U(y)/{{\rm d}y}$, $U'':={\rm d}^2U(y)/{{\rm d}y}^2$, $\bar {\rho }':={\rm d}\bar {\rho }(y)/{{\rm d}y}$, $\widehat {{\nabla }}^2:=\partial _{yy}-k_x^2-k_z^2$, $\widehat {\nabla }^4:=\partial ^{(4)}_{y}-2(k_x^2+k_z^2)\partial _{yy}+(k_x^2+k_z^2)^2$ and $\mathcal {I}$ is the identity operator. The equation associated with the $\hat {\mathcal {A}}^S$ operator in (2.7*a*) can also be obtained by generalizing the classical Taylor–Goldstein equation (Goldstein Reference Goldstein1931; Taylor Reference Taylor1931; Smyth & Carpenter Reference Smyth and Carpenter2019) to include viscosity, density diffusivity and coupling with wall-normal vorticity $\hat {\omega }_y$ with $k_z\neq 0$. The boundary conditions associated with (2.6) are $\hat {v}(y=\pm 1)= ({\partial \hat {v}}/{\partial y})(y=\pm 1)=\hat {\omega }_y(y=\pm 1)=\hat {\rho }(y=\pm 1)=0$.

The spatio-temporal frequency response $\mathcal {H}^S$ of the system in (2.6), which maps the input forcing to the velocity and density fields at the same spatial-temporal wavenumber–frequency triplet; i.e. $\left [\begin {smallmatrix}\boldsymbol {\hat {u}}(y;k_x,k_z,\omega )\\ {\hat {\rho }}(y;k_x,k_z,\omega ) \end {smallmatrix}\right ]=\mathcal {H}^S(y;k_x,k_z,\omega )\left [\begin {smallmatrix} \boldsymbol {\hat {f}}_{\boldsymbol {u}, \xi }(y;k_x,k_z,\omega )\\ \hat {f}_{\rho,\xi }(y;k_x,k_z,\omega )\end {smallmatrix}\right ]$, is given by

Here, $\mathcal {I}_{3\times 3}:=\text {diag}(\mathcal {I},\mathcal {I},\mathcal {I})$, where $\text {diag}(\,\cdot\,)$ indicates a block-diagonal operation.

The linear form of (2.4*a*)–(2.4*b*) allows us to perform the same spatio-temporal Fourier transform on the model of the nonlinearity, which can be decomposed as

A block diagram illustrating this decomposition of the modelled nonlinearity is shown inside the blue dashed line ($\text {- -}$, blue) in figure 2(*a*). This block-diagonal structure constrains the modelled nonlinear interactions, i.e. provides structured forcing.

In order to isolate the gain operator $-\boldsymbol {u}_{\xi }$, we combine the linear gradient operator with the spatio-temporal frequency response of the linearized system (2.8). The resulting modified frequency response operator with outputs that are the vectorized gradients of the velocity and density components is defined as

The resulting system model can be redrawn as a feedback interconnection between this linear operator and the structured uncertainty

Here, the structure is introduced in terms of the diagonal form of $\hat {\boldsymbol {u}}_{\varXi }^S$ that enforces the componentwise structure of the nonlinearity in the forcing model defined in (2.4). Figure 2(*b*) illustrates this feedback interconnection between the modified spatio-temporal frequency response and the structured uncertainty, where $\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}$ and $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$, respectively, represent the spatial discretizations (numerical approximations) of $\mathcal {H}^S_{\boldsymbol {\nabla }}$ in (2.10) and $\boldsymbol {\hat {u}}^S_{\varXi }$ in (2.11).

We are interested in characterizing the horizontal length scales of the most amplified flow structures under this structured forcing. This amplification can be quantified in terms of the structured singular value of the modified frequency response operator $\mathcal {H}^S_{\boldsymbol {\nabla }}$; see e.g. Packard & Doyle (Reference Packard and Doyle1993, definition 3.1) and Zhou *et al.* (Reference Zhou, Doyle and Glover1996, definition 11.1), which is defined as follows.

Definition 2.1 Given wavenumber and frequency pair $(k_x,k_z,\omega )$, the structured singular value $\mu _{\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}(k_x,k_z,\omega )]$ is defined as

If no $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }\in \hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ makes $\boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$ singular, then $\mu _{\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}]:=0$.

Here, $\bar {\sigma }[\,\cdot\,]$ is the largest singular value, $\det [\,\cdot\,]$ is the determinant of the argument, and $\boldsymbol{\mathsf{I}}$ is the identity matrix. The subscript of $\mu$ in (2.12) is a set $\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ containing all uncertainties having the same block-diagonal structure as $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$; i.e.

where $N_y$ denotes the number of grid points in $y$.

Definition 2.1 suggests that the inverse of the structured singular value $1/\mu$ is the minimal norm of the perturbation $\hat {\boldsymbol{\mathsf{u}}}^S_{\varXi }$ that destabilizes the feedback interconnection in figure 2(*b*) in the input–output sense defined by the small gain theorem, see Liu & Gayme (Reference Liu and Gayme2021, proposition 2.2) and Zhou *et al.* (Reference Zhou, Doyle and Glover1996, theorem 11.8). This interpretation suggests that the flow field is more sensitive to perturbations with the flow structures associated with a larger amplification measured by the structured singular value $\mu$. A similar notion of destabilizing perturbation was also employed to interpret the largest (unstructured) singular value (Trefethen *et al.* Reference Trefethen, Trefethen, Reddy and Driscoll1993, p. 581).

Here, the form of structured uncertainty in (2.13) allows the full degrees of freedom for the complex matrix $-\hat {\boldsymbol{\mathsf{u}}}^{\text {T}}_{\xi }\in \mathbb {C}^{N_y\times 3N_y}$ for ease of computation. While $\boldsymbol {u}_\xi$ is not constrained to be incompressible, the incompressibility of $\boldsymbol {u}$ and the role of pressure are accounted for within the current $v$-$\omega _y$ formulation. Further refinement to better represent the physics and uncover the form of $\boldsymbol {u}_\xi$ requires an extension of both the analysis method and computational tools. These extensions are beyond the scope of the current work.

We then define the structured response following Liu & Gayme (Reference Liu and Gayme2021) as

where $\text {sup}$ represents a supremum (least upper bound) operation. Here, we abuse the notation and terminology by writing $\|\,\cdot\,\|_{\mu }$ (Packard & Doyle Reference Packard and Doyle1993), although this quantity is not a proper norm. We employ this notation in analogy with the corresponding unstructured response of the feedback interconnection, which is given by

This quantity is the unstructured counterpart of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$, which is obtained by replacing the structured uncertainty set $\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }$ with the set of full complex matrices $\mathbb {C}^{4N_y\times 12N_y}$. In both cases, a larger value indicates that the corresponding flow structures (associated with a particular $k_x$ and $k_z$ pair) have larger amplification under either structured or unstructured feedback forcing. For example, a larger value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)$ indicates that the corresponding flow structures (associated with a particular $k_x$ and $k_z$ pair) have larger amplification under structured feedback in figure 2(*b*).

### 2.3. Numerical method

We employ the Chebyshev differential matrix (Trefethen Reference Trefethen2000; Weideman & Reddy Reference Weideman and Reddy2000) to discretize the operators in equation set (2.7). Our code is validated through comparison with the unstratified PCF and Poiseuille flow results in Jovanović (Reference Jovanović2004), Jovanović & Bamieh (Reference Jovanović and Bamieh2005) and Schmid (Reference Schmid2007). The implementation of stratification is validated by reproducing the maximum growth rate of the linear normal mode in a layered stratified PCF determined by Eaves & Caulfield (Reference Eaves and Caulfield2017, figures 3 and 6*a*), as well as the linear stability predictions for the unstable stratification configuration in Olvera & Kerswell (Reference Olvera and Kerswell2017, figure 1 and appendix B). We use $N_y=60$ collocation points not including the boundary points over the wall-normal extent, as well as $48$ and $36$ logarithmically spaced streamwise and spanwise wavenumbers in the respective spectral ranges $k_x \in [10^{-4},10^{0.48}]$ and $k_z \in [10^{-2},10^{1.2}]$, unless otherwise mentioned. To verify that this resolution is sufficient to achieve grid convergence we recomputed selected results with 1.5 times the number of collocation points in the wall-normal direction and verified that the results did not change. The quantity $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (2.14) for each wavenumber pair $(k_x,k_z)$ is computed using the mussv command in the Robust Control Toolbox (Balas *et al.* Reference Balas, Chiang, Packard and Safonov2005) of MATLAB. The arguments of mussv employed here include the state-space model of $\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla }}$ that samples the frequency domain adaptively. The BlockStructure argument comprises four full $N_y\times 3N_y$ complex matrices, and we use the ‘Uf’ algorithm option.

## 3. Structured spatio-temporal frequency response of stratified flow

In this section, we use the SIOA approach described in § 2.2 to characterize the flow structures that are most amplified in stably stratified PCF.

### 3.1. Low-$Re$ low-$Ri_b$ vs high-$Re$ high-$Ri_b$ intermittency

In this subsection, we analyse flow structures that are prominent in either the low-$Re$ low-$Ri_b$ or the high-$Re$ high-$Ri_b$ intermittent regime (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015). Here, we keep the Prandtl number fixed at $Pr=0.7$. This value corresponds to thermally stratified air and is the same value studied by Deusebio *et al.* (Reference Deusebio, Caulfield and Taylor2015). We first consider a flow with $Re=865$, $Ri_b=0.02$ and $Pr=0.7$, where oblique turbulent bands have been observed (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016). In order to evaluate the relative effect of the feedback interconnection and the imposed structure, we also compute $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }(k_x,k_z)$ defined in (2.15) and

Here, $\boldsymbol{\mathsf{H}}^S$ is the discretization of spatio-temporal frequency response operator $\mathcal {H}^S$ in (2.8), i.e. the spatio-temporal frequency operator governing the linearized dynamics without the feedback interconnection. The values of $\|\mathcal {H}\|_{\infty }$ for unstratified plane Couette and plane Poiseuille flows were previously analysed in Jovanović (Reference Jovanović2004), Schmid (Reference Schmid2007) and Illingworth (Reference Illingworth2020). The quantity in (3.1) describes the maximum singular value of the frequency response operator $\mathcal {H}^S$, which represents the maximal gain of $\mathcal {H}^S$ over all temporal frequencies; i.e. the worst-case amplification over harmonic inputs.

Figure 3 shows $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (*a*) alongside (*b*) $\|\mathcal {H}^S\|_{\infty }$ and (*c*) $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$. We indicate the characteristic wavelength pair $\lambda _x=32{\rm \pi}$, $\lambda _z=16{\rm \pi}$ corresponding to the oblique turbulent bands observed in DNS under the same flow regime (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016, figure 2*b*) in these panels using the symbol ($\triangle$, black). These structures are observed to have a characteristic inclination angle (measured from the streamwise direction in $x$–$z$ plane) of $\theta :=\tan ^{-1}(\lambda _z/\lambda _x)\approx 27^\circ$, which is indicated in all panels by the black solid line (–) that plots $\lambda _z=\lambda _x\tan (27^\circ )$. While there is some footprint of these structures and this angle in all three panels, the correspondence with the peak amplitude is most prominent in (*a*). In fact, the peak value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ in (*a*) occurs at streamwise and spanwise wavenumbers associated with the characteristic wavelengths and angle of the oblique turbulent bands reported in DNS (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016), and the line representing the angle of the oblique turbulent bands crosses through the centre of the narrow roughly elliptical peak region whose principal axis coincides with this angle. The results in figure 3(*a*) suggest that the SIOA captures both the wavelengths and angle of the oblique turbulent bands in the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF. This analysis suggests that these oblique turbulent bands arise in the intermittent regime of stratified PCF due to their large amplification, or equivalently their sensitivity to disturbances.

The traditional input–output analysis results, $\|\mathcal {H}^S\|_\infty$ in (*b*), provide a noticeable improvement compared with growth rate analysis (as presented in more detail in Appendix A) and are also able to identify the preferred wavenumber pair in this intermittent regime. However, this analysis suggests larger amplification of the streamwise elongated modes. Moreover, the inclusion of an unstructured feedback loop quantified through $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ in (*c*) correctly orders the relative amplification between the oblique turbulent bands and streamwise elongated structures ($k_x\approx 0$). The differences between $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ and $\|\mathcal {H}^S\|_{\infty }$ are likely associated with the additional $\widehat {\boldsymbol {\nabla }}$ operator in defining $\mathcal {H}_{\boldsymbol {\nabla }}$ in (2.10), which emphasizes flow structures with a larger horizontal wavenumber. The difference between $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }$ is associated with the structured feedback interconnection that constrains the permissible feedback pathway, which weakens the amplification associated with the lift-up mechanism; see similar discussion on unstratified PCF (Liu & Gayme Reference Liu and Gayme2021, § 3.3). A comparison of the results in figure 3 suggests that it is the imposition of the componentwise structure from the nonlinear terms in (2.3) that further improves the prediction of the oblique turbulent bands.

We now consider the high-$Re$ high-$Ri_b$ intermittent regime, which was shown to be qualitatively different in behaviour from the low-$Re$ low-$Ri_b$ intermittent regime (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015). We first isolate the effect of increasing either $Re$ or $Ri_b$. Figure 4(*a*) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for a flow with $Ri_b=0.02$ and $Re=4250$. The larger colour bar range vs figure 3 highlights the expected higher magnitudes vs those for a flow with a lower Reynolds number ($Re=865$). We can see that the wavenumber pair of the peak region extends towards smaller values (larger wavelengths) than those associated with the oblique turbulent bands that were in the peak region in figure 3(*a*). Figure 4(*b*) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for a higher bulk Richardson number $Ri_b=0.2541$ and the same $Re$ and $Pr$ values as figure 3(*a*). Here, the amplification associated with the streamwise-varying flow structures such as the oblique turbulent bands observed in figure 3(*a*) is reduced and quasi-horizontal flow structures $(k_x\approx 0, k_z\approx 0)$ show a similar level of amplification (see the bottom left corner in figure 4*b*). Note that this flow structure associated with $k_x\approx 0$, $k_z\approx 0$ is referred to as quasi-horizontal to distinguish it from a horizontally uniform mode ($k_x=0$, $k_z=0$).

Armed with these insights, we consider the combined high-$Re$ high-$Ri_b$ intermittent regime ($Re=52\,630$ and $Ri_b=0.15$); these values correspond to the results shown in figure 7 of Deusebio *et al.* (Reference Deusebio, Caulfield and Taylor2015). Figure 4(*c*) presents $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ for these parameter values with an increased wall-normal grid with $N_y=90$. Here, the amplification of the oblique turbulent band is of a similar order to that of flow structures with a wide range of wavenumber pairs ranging from $k_x\lesssim 10^{-2}$ and $k_z\lesssim 1$ down to $k_x\approx 0$ and $k_z\approx 0$. These latter flow structures resemble the quasi-horizontal flow structures that have a horizontal length scale much larger than their vertical length scales, which are limited by the channel height and therefore restricted to non-dimensional scales of the order of $2$. The response in this regime, therefore, shows a large qualitative difference from that in the low-$Re$ low-$Ri_b$ ($Re=865$ and $Ri_b=0.02$) intermittent regime shown in figure 3(*a*). This qualitative difference mirrors the different features in intermittent regimes described by Deusebio *et al.* (Reference Deusebio, Caulfield and Taylor2015), where oblique turbulent bands are prevalent in the low-$Re$ low-$Ri_b$ intermittent regime, but the high-$Re$ high-$Ri_b$ intermittent regime is characterized by turbulent–laminar layers indicating a large horizontal length scale.

In figure 4(*c*), we also observe that the quasi-horizontal flow structures have a streamwise wavelength much larger than their spanwise wavelength ($\lambda _x\gg \lambda _z$), which is also consistent with the observation in Deusebio *et al.* (Reference Deusebio, Caulfield and Taylor2015) that the turbulent and laminar layers in the high-$Re$ high-$Ri_b$ intermittent regime are homogeneous in the streamwise direction. This behaviour can be understood through an order of magnitude estimation of the terms in the continuity equation. We assume highly anisotropic flow with $v\approx 0$ under strong stratification, which simplifies the continuity equation to ${\partial u}/{\partial x}+{\partial w}/{\partial z}=0$. We further assume that the restoring buoyancy force due to stratification does not have a preference between the streamwise or spanwise direction and, therefore, we also assume ${\partial u}/{\partial x}$ and ${\partial w}/{\partial z}$ are the same order of magnitude. However, in the current stratified PCF configuration, streamwise velocity is associated with a characteristic velocity much larger than its spanwise counterpart due to the base flow velocity. As a result, streamwise variation is reduced much faster than spanwise variation (i.e. $k_x\ll k_z$) to keep ${\partial u}/{\partial x}$ and ${\partial w}/{\partial z}$ the same order of magnitude.

### 3.2. Case $Ri_b>1/4$: a change in the most amplified flow structures

The Miles–Howard theorem (Howard Reference Howard1961; Miles Reference Miles1961) implies that the laminar base flow would be linearly stable in the limits where $\nu$ and $\kappa$ both are zero if $Ri_b > 1/4$. Although the theorem is not applicable to unsteady flows with finite $Re$ and $Pr$, it has also been observed that a ‘marginal’ or ‘critical’ Richardson number near this threshold value appears to emerge naturally in simulations (Salehipour *et al.* Reference Salehipour, Peltier and Caulfield2018) and field measurement (Smyth *et al.* Reference Smyth, Nash and Moum2019). In the previous section, we noted that increasing $Ri_b=0.02$ to $Ri_b=0.2541$, for a fixed Reynolds number of $Re=865$, reduces the overall response and changes the types of flow structures, $(k_x, k_z)$ wavenumber pairs, that exhibit the largest response, see figures 3(*a*) and 4(*b*). In this subsection, we further investigate whether this apparently marginal threshold $Ri_b \simeq 1/4$ is associated with a change in flow structures and whether this behaviour is independent of $Re$. As in the previous subsection, the Prandtl number is fixed at $Pr=0.7$.

Here, we aggregate results varying over a range of $(k_x,k_z)$ wavenumber pairs in terms of the maximum value

over the wavenumber domain $k_x\in [10^{-6},10^{0.48}]$ and $k_z\in [10^{-6}, 10^{0.48}]$. Lowering the minimum value of the wavenumber ranges vs those considered in the previous subsection is motivated by the observation in figure 4 that both the $k_x$ and $k_z$ values corresponding to the peak region decrease with increasing Reynolds and Richardson numbers. To separate streamwise-varying and streamwise-independent flow structures, we similarly evaluate

This quantity restricts the streamwise wavenumber to $k_x=10^{-6}$ to approximate the streamwise constant modes and computes the maximum value over $k_z\in [10^{-6}, 10^{0.48}]$, where we use an upper bound of $10^{0.48}$ rather than the larger value of $10^{1.2}$ to save computation time. This change in the upper bound was not found to affect the results since the $k_z$ associated with the maximum value was consistently found to be below this upper bound. The restriction in streamwise wavelengths to $k_x=10^{-6}$ naturally includes the quasi-horizontal flow structures prevalent in the high-$Re$ high-$Ri_b$ regime ($k_x\approx 0$, $k_z\approx 0$) as an extreme case, but excludes streamwise-varying flow structures such as the oblique turbulent bands observed in the low-$Re$ low-$Ri_b$ regime discussed in § 3.1.

Figure 5 shows the variation of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ (solid lines) and $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ (dashed lines) with bulk Richardson number $Ri_b\in [0,6]$ and Reynolds number $Re\in [865,15\,000]$. The quantities $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{M}$ including streamwise-varying flow structures are very close to $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ when $Ri_b \gtrsim 1/4$ for the full range of Reynolds numbers $Re\in [865,15\,000]$ in figure 5(*a*). This phenomenon is also reflected in figure 5(*b*), where for flows with $Ri_b=0.24$ and $Ri_b=0.75$, the curves for $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{M}$ and $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ largely overlap. These trends suggest that the inviscid marginal stability value $Ri_b=1/4$ predicted by the Miles–Howard theorem (Howard Reference Howard1961; Miles Reference Miles1961) for laminar flow is apparently associated with a change in the most amplified flow structure in stratified PCF at finite $Re$ and $Pr \sim 1$.

The plots in figure 5 show that the largest amplification of the streamwise-invariant modes represented by $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}$ (dashed lines) do not appear to be influenced by $Ri_b$, as shown by the horizontal dashed lines in figure 5(*a*) and the overlapping dashed lines in figure 5(*b*). We further explore this $Ri_b$ independence for streamwise constant flow structures by considering the limit of horizontal invariance $\partial _x(\,\cdot\,)=0$ and $\partial _z(\,\cdot\,)=0$ ($k_x=0$ and $k_z=0$), which directly results in $\partial _y v=0$ due to the continuity equation. The boundary condition $v(y=\pm 1)=0$ then directly results in $v=0$. Using these assumptions, the advection terms vanish (i.e. $U\partial _x(\,\cdot\,)=0$, and $\boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla }(\,\cdot\,)=0$) in each momentum and density equation. The terms associated with the background shear and density gradient also vanish due to zero wall-normal velocity; i.e. $vU'=v\bar {\rho }'=0$.

These observations lead to a simplification of the momentum and density equations in (2.1) to

*a*,

*b*)\begin{gather} \partial_t u=\frac{1}{Re}\partial_{yy} u,\quad \partial_y p={-}Ri_b \rho, \end{gather}

*c*,

*d*)\begin{gather}\partial_t w=\frac{1}{Re}\partial_{yy} w,\quad \partial_t \rho=\frac{1}{RePr}\partial_{yy} \rho. \end{gather}

Here, we can see that horizontal momentum and density field equations are all reduced to the diffusion equation, and the wall-normal momentum equation is reduced to a balance between the buoyancy force and the vertical pressure gradient i.e. to a hydrostatic balance. This balance suggests that the only dependence on $Ri_b$ can be absorbed into the pressure by rescaling pressure and thus does not influence the quasi-horizontal flow structures. The results presented in figures 3(*a*) and 4 suggest that flow structures with $k_x=10^{-4}$ lead to the same structured response $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_\mu$ at a wide range of spanwise wavenumbers $k_z\leq 1$, and this value is consistent with that of quasi-horizontal flow structures $(k_x\approx 0,k_z\approx 0)$ that are independent of $Ri_b$.

This analysis provides evidence that quasi-horizontal flow structures are associated with amplification which is independent of $Ri_b$. Instead, high $Ri_b$ (i.e. strong stratification) will reduce the amplification of other horizontally varying flow structures such as the oblique turbulent bands observed in the low-$Re$ low-$Ri_b$ intermittent regime of stably stratified PCF. Furthermore, it appears in figure 5(*b*) that quasi-horizontal flow structures are also increasingly amplified as $Re$ increases with scaling law $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}\sim Re$. This behaviour indicates that quasi-horizontal flow structures may prefer a high-$Re$ high-$Ri_b$ regime. In § 4, we further explore this scaling law $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^{sc}\sim Re$ by developing analytical scaling arguments for $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in unstratified and streamwise-invariant flow.

### 3.3. Effects of low and high $Pr$

The Prandtl number is known to play an important role in the types of flow structures characterizing stratified PCF (Taylor & Zhou Reference Taylor and Zhou2017; Zhou *et al.* Reference Zhou, Taylor and Caulfield2017*a*,Reference Zhou, Taylor, Caulfield and Linden*b*; Langham *et al.* Reference Langham, Eaves and Kerswell2020). The Prandtl number also varies over a wide range in different applications. For example, $Pr\ll 1$ is relevant for flow in the stellar interior; see e.g. Garaud (Reference Garaud2021), while a Prandtl number of $Pr=7$ corresponds to thermally stratified water. The Schmidt number (the analogous parameter to the Prandtl number for compositionally induced density variations) for salt-stratified water is significantly larger $Sc \simeq 700$. Moreover, the Prandtl number is obviously not well defined under the inviscid and non-diffusive assumptions of the Miles–Howard theorem. In this subsection, we explore the effect of low or high Prandtl number on flow structures. Here, we keep the Reynolds number fixed at $Re=4250$ following Zhou *et al.* (Reference Zhou, Taylor and Caulfield2017*a*,Reference Zhou, Taylor, Caulfield and Linden*b*). In order to resolve fully the additional scales introduced at high $Pr$, we increase the number of wall-normal grid points to $N_y=120$ at $Pr=70$, which is chosen as a more computationally accessible ‘large’ value, as previously considered by Zhou *et al.* (Reference Zhou, Taylor and Caulfield2017*a*,Reference Zhou, Taylor, Caulfield and Linden*b*).

We first investigate the effect of low $Pr$. The cross-channel density profiles of exact coherent structures in stratified PCF were shown to match in flows with the same $PrRi_b$ at $Pr\in [10^{-4},10^{-2}]$ (Langham *et al.* Reference Langham, Eaves and Kerswell2020, figure 3). This combined measure $PrRi_b$ has been proposed as the natural control parameter for stably stratified shear flows at the low Prandtl number limit $Pr\ll 1$ (Lignieres Reference Lignieres1999; Garaud *et al.* Reference Garaud, Gallet and Bischoff2015). In order to further explore this dependence, we plot $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ as a function of $PrRi_b$ for Prandtl numbers in the range $Pr\in [10^{-4}, 7]$ in figure 6. Here, the results $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }^M$ for $Pr\in [10^{-4},10^{-1}]$ (magenta dashed lines) again show a natural matching dependence on $PrRi_b$. This behaviour breaks down for flows with $Pr\geq 1$, as shown in figure 6. A similar end to the region of matched dependence on $PrRi_b$ alone is observed in studies using exact coherent structures, where the density profile at $Pr= 0.1$ deviates from the matching profiles for flows with $Pr\in [10^{-4}, 10^{-2}]$ yet fixed $Pr Ri_b$ (Langham *et al.* Reference Langham, Eaves and Kerswell2020, figure 3).

In figure 7, we plot $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ as a function of $Ri_b$ and $Pr$ over the respective ranges $Ri_b\in [0,6]$ and $Pr\in [10^{-4},70]$. Figure 7(*a*) shows that the marginal stability value $Ri_b=1/4$ is not associated with any significant changes in the types of flow features undergoing the largest amplification for flows with $Pr=0.01$ ($\triangle$). Figure 7(*b*) further suggests that flows with smaller $Pr$ require a larger $Ri_b$ to reduce the amplification of streamwise-varying flow structures to the same level as streamwise constant structures. This behaviour is consistent with the observation that the exact coherent structures in the low $Pr$ limit require a larger $Ri_b$ to be affected by stratification in PCF (Langham *et al.* Reference Langham, Eaves and Kerswell2020). Figure 7 further shows that, for flows with high $Pr$, the quantities $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ and $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ are the same over a wide range of $Ri_b\in [0,6]$. In particular, the horizontal lines plotted in figure 7(*a*) for flows with different $Ri_b$ collapse to one line in the high $Pr\gg 1$ regime shown in figure 7(*b*). This observation is also consistent with Langham *et al.* (Reference Langham, Eaves and Kerswell2020), who noted that, in the high Prandtl number limit $Pr\gg 1$, the effect of increasing $Ri_b$ is mitigated.

In order to investigate in isolation the effect of Prandtl number on the amplification of each wavenumber pair $(k_x,k_z)$, in figure 8 we plot $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)$ for flows with (*a*) $Pr=10^{-4}$, (*b*) $Pr=7$ and (*c*) $Pr=70$. The peak region in figure 8(*a*) at $Pr=10^{-4}$ resembles the shape of the peak region for unstratified PCF (Liu & Gayme Reference Liu and Gayme2021, figure 4*a*). We have also computed the results for unstratified PCF at the same Reynolds number $Re=4250$ (not shown here), and we find almost the same results as shown in figure 8(*a*). This similarity suggests that, for the same bulk Richardson number, a lower Prandtl number will result in a weakening of the stabilizing effect of stratification. Comparing $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in figure 4(*a*) at $Pr=0.7$ with the same quantity at $Pr=7$ and $Pr=70$, respectively shown in figures 8(*b*) and 8(*c*), we notice that the amplification associated with the wavenumber pair $k_x=10^{-4}$ and $k_z=10^{-2}$ increases with $Pr$. More specifically, the value of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }$ associated with $k_x=10^{-4}$ and $k_z=10^{-2}$ becomes comparable to the values associated with the $k_x\approx 10^{-2}$ and $k_z\approx 10^{-1}$ at $Pr=7$ as shown in 8(*b*). The wavenumber pair $k_x=10^{-4}$ and $k_z=10^{-2}$ is associated with the largest magnitude over the $(k_x,k_z)$ contour region at $Pr=70$, as shown in figure 8(*c*).

The quasi-horizontal flow structures ($k_x\approx 0$, $k_z\approx 0$) observed in flows at high $Pr$ have different features from those previously observed in the high-$Re$, high-$Ri_b$ regime (e.g. results for flows with $Re=52\,630$ and $Ri_b=0.15$ shown in figure 4*c*) and described in § 3.1. This indicates that a new type of quasi-horizontal flow structure appears in flows with sufficiently high $Pr$. The appearance of this flow structure at a high $Pr$ suggests that this flow structure is associated with fluctuations in the density field. This can be further explored by isolating the input–output pathway for each component of the momentum and density equations, i.e. inputs $f_x, f_y,\ f_z,\ f_{\rho }$ in (2.3) to outputs $u$, $v$, $w$ and $\rho$. These input–output pathways can be studied through the definition of operators $\mathcal {H}^S_{ij}$, where $j$ defines the forcing input component ($j=x,y,z,\rho$) and $i=u,v,w,\rho$ describes each velocity or density output component

with

*a*,

*b*)\begin{gather} \mathcal{\hat{B}}^S_x:=\mathcal{\hat{B}}^S\begin{bmatrix} \mathcal{I} & 0 & 0 & 0 \end{bmatrix}^{\text{T}},\quad \mathcal{\hat{B}}^S_y:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & \mathcal{I} & 0 & 0 \end{bmatrix}^{\text{T}}, \end{gather}

*c*,

*d*)\begin{gather}\mathcal{\hat{B}}^S_z:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & 0 & \mathcal{I} & 0 \end{bmatrix}^{\text{T}},\quad \mathcal{\hat{B}}^S_\rho:=\mathcal{\hat{B}}^S\begin{bmatrix} 0 & 0 & 0 & \mathcal{I} \end{bmatrix}^{\text{T}}, \end{gather}

*e*,

*f*)\begin{gather}\mathcal{\hat{C}}^S_u:=\begin{bmatrix}\mathcal{I} & 0 & 0 & 0 \end{bmatrix} \mathcal{\hat{C}}^S ,\quad \mathcal{\hat{C}}^S_v:=\begin{bmatrix}0 & \mathcal{I} & 0 & 0 \end{bmatrix} \mathcal{\hat{C}}^S, \end{gather}

*g*,

*h*)\begin{gather}\mathcal{\hat{C}}^S_w:=\begin{bmatrix}0 & 0 & \mathcal{I} & 0 \end{bmatrix} \mathcal{\hat{C}}^S,\quad \mathcal{\hat{C}}^S_\rho:=\begin{bmatrix}0 & 0 & 0 & \mathcal{I} \end{bmatrix} \mathcal{\hat{C}}^S. \end{gather}

Figures 9 and 10 present quantities $\|\mathcal {H}^S_{ij}\|_{\infty }$ $ij=ux,vy,wz$ and $\rho \rho$ for the control parameters ($Re=865, Ri_b=0.02, Pr=0.7)$ and $(Re=4250,Ri_b=0.02, Pr=70)$ respectively. The combined effect of these four panels associated with the input and output in the same component resembles the shape of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ at the same flow regime in figures 3(*a*) and 8(*c*). This correspondence is because the structured feedback interconnections in (2.4*a*)–(2.4*b*) constrain the permissible feedback pathways. In figure 10, (*d*) $\|\mathcal {H}_{\rho \rho }^S\|_{\infty }$ (associated with input $f_\rho$ and output $\rho$) is significantly larger than the responses shown in the other panels for a flow with $Pr=70$, suggesting the strong role of density in the amplification for this parameter range. We can further isolate each component of the frequency response operator $\mathcal {H}_{\boldsymbol {\nabla }}^S$ by defining

The values $\|\mathcal {H}_{\boldsymbol {\nabla } ij}^S\|_{\infty }$, not shown here for brevity, show qualitatively similar behaviour to $\|\mathcal {H}^S_{ij}\|_{\infty }$ as plotted in figures 9 and 10. This componentwise analysis demonstrates that the quasi-horizontal flow structures appearing at high $Pr$ are associated with density fluctuations. The appearance of this type of quasi-horizontal flow structure associated with density fluctuations in a high $Pr$ regime is qualitatively consistent with previous observations that sharp density gradients or even density ‘staircases’ can be observed when $Pr$ is increased (Taylor & Zhou Reference Taylor and Zhou2017; Zhou *et al.* Reference Zhou, Taylor, Caulfield and Linden2017*b*).

## 4. Scaling for density-associated flow structures when $Pr\gg 1$

The previous subsection reveals the appearance of quasi-horizontal flow structures associated with density fluctuations in the $Pr\gg 1$ limit. In this section, we construct an analytical scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in terms of $Re$ and $Pr$ to provide further evidence that such flow structures prefer the $Pr\gg 1$ regime. The analytical scaling in terms of $Re$ and $Pr$ can also further provide insight into high $Re$ and $Pr$ flow regimes beyond the current computation range achievable through DNS.

We assume streamwise-invariance ($k_x=0$) and unstratified flow ($Ri_b=0$) to facilitate analytical derivation. The importance of streamwise-invariant flow structures is suggested by the quasi-horizontal flow structures ($k_x\approx 0$ and $k_z\approx 0$), which are nearly streamwise constant. The independence with respect to variations in $Ri_b$ of the amplification of streamwise-invariant flow structures $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|^{sc}_{\mu }$ shown in figures 5 and 7 and the analysis of Langham *et al.* (Reference Langham, Eaves and Kerswell2020) suggest that $Ri_b=0$ (i.e. density fluctuations can be treated as a passive scalar) is a reasonable regime to consider to obtain further insight. The analytically derived $Re$ and $Pr$ scalings of each component of $\mathcal {H}_{ij}^S$ in (3.5) and $\mathcal {H}_{\boldsymbol {\nabla } ij}^S$ in (3.7) are presented in theorems 4.1(*a*) and (*b*), respectively.

Theorem 4.1 Consider streamwise-invariant ($k_x=0$) unstratified ($Ri_b=0$) plane Couette flow with a passive scalar ‘density’ field.

(*a*) Each component of $\|\mathcal {H}^S_{ij}\|_{\infty }$ ($i=u,v,w,\rho$ and $j=x,y,z,\rho$) scales as

where functions $h^S_{ij}(k_z)$ are independent of $Re$ and $Pr$.

(*b*) Each component of $\|\mathcal {H}^S_{\boldsymbol {\nabla } ij}\|_{\infty }$ ($i=u,v,w,\rho$ and $j=x,y,z,\rho$) scales as

where functions $h^S_{\boldsymbol {\nabla } ij}(k_z)$ are independent of $Re$ and $Pr$.

The first three columns and three rows presented in (4.1) are the same as those derived in Jovanović (Reference Jovanović2004, theorem 11) for unstratified wall-bounded shear flows with no passive scalar field. The details of the proof are presented in Appendix B. These results demonstrate that $Pr$ only contributes to the scaling associated with the density field (here of course assumed to be a passive scalar); i.e. the bottom rows of (4.1) and (4.2) corresponding to the density output. We also note that the rightmost columns of (4.1) and (4.2) show that the forcing in the density equation $f_{\rho }$ does not influence the output corresponding to velocity components $u$, $v$ and $w$, which is consistent with the assumption that $Ri_b=0$, in that the density perturbation behaves as a passive scalar.

The effect of imposing a componentwise structure for the nonlinearity within the feedback is analogous to the effect seen in unstratified PCF (Liu & Gayme Reference Liu and Gayme2021, § 3.3). The imposed correlation between each component of the modelled forcing $f_{x,\xi }$, $f_{y,\xi }$, $f_{z,\xi }$, $f_{\rho,\xi }$, and the respective velocity and density components $u$, $v$, $w$, $\rho$ constrain the influence of the forcing to its associated component of the velocity or density field. Thus, the overall scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ is related to the worst-case scaling of the diagonal terms in (4.2) in theorem 4.1. The concept is formalized in theorem 4.2, and we relegate the details of the proof to Appendix B.

Theorem 4.2 Given a wavenumber pair $(k_x,k_z)$

We can combine results in theorems 4.1(*b*) and 4.2 to obtain the scaling of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ in corollary 4.3

Corollary 4.3 Consider streamwise-invariant ($k_x=0$) unstratified ($Ri_b=0$) PCF with a passive scalar ‘density’ field

where functions $h^S_{\boldsymbol {\nabla } ij}(k_z)$ with $ij=ux,vy,wz,\rho \rho$ are independent of $Re$ and $Pr$.

Although corollary 4.3 provides a lower bound on $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$, the numerical results suggest that $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$ follows the same $Re$ and $Pr$ scaling as the right-hand side of (4.4) in corollary 4.3. For example, corollary 4.3 suggests that the lower bound of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)$ will scale as ${\sim }Re$ at a fixed $Pr$, which is consistent with the red dashed lines of figure 5(*b*). At a fixed $Re$, corollary 4.3 also suggests that $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)\sim Pr$ in the limit $Pr\gg 1$, but $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }(0,k_z)$ will become independent of $Pr$ in the limit $Pr\ll 1$. This is also consistent with the numerical results shown in the red dashed lines of figure 7(*b*) that suggest $\|\mathcal {H}_{\boldsymbol {\nabla }}^{S}\|^{sc}_{\mu }\sim Pr$ when $Pr\gg 1$ and independently of $Pr$ when $Pr\ll 1$. For $Pr\gg 1$, theorem 4.2 and corollary 4.3 further suggest that the component $\|\mathcal {H}^{S}_{\boldsymbol {\nabla } \rho \rho }\|_{\infty }$ associated with the density will dominate the overall behaviour of $\|\mathcal {H}_{\boldsymbol {\nabla }}^S\|_{\mu }$, which is consistent with the large amplification of quasi-horizontal flow structures associated with density fluctuations, i.e. $\|\mathcal {H}^S_{\rho \rho }\|_\infty$ shown in figure 10(*p*). Corollary 4.3 further supports the notion that the flow structures associated with density fluctuations prefer the flow regime with $Pr\gg 1$ under the assumptions of streamwise-invariant ($k_x=0$) and unstratified ($Ri_b=0$) flow.

## 5. Conclusions and future work

In this paper, we have extended the SIOA originally developed for unstratified wall-bounded shear flows (Liu & Gayme Reference Liu and Gayme2021) to stratified PCF. We first apply SIOA to characterize highly amplified flow structures in the intermittent regimes of stratified PCF. We examine how variations in $Re$ and $Ri_b$ affect flow structures with $Pr=0.7$. SIOA predicts the characteristic wavelengths and angle of the oblique turbulent bands observed in very large channel size DNS of the low-$Re$ low-$Ri_b$ intermittent regime of stratified PCF (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015; Taylor *et al.* Reference Taylor, Deusebio, Caulfield and Kerswell2016). In the high-$Re$ high-$Ri_b$ intermittent regime, SIOA identifies quasi-horizontal flow structures resembling turbulent–laminar layers (Deusebio *et al.* Reference Deusebio, Caulfield and Taylor2015).

Having validated the ability of the SIOA approach to predict important structures in the intermittent regime, we next investigate the behaviour of the flow across a range of important control parameters $Re$, $Ri_b$ and $Pr$. Increasing $Ri_b$ is shown to reduce the amplification of streamwise-varying flow structures. The results indicate that the classical marginally stable $Ri_b= 1/4$ for the laminar base flow appears to be associated with a change in the most amplified flow structures, an observation which is robust for a wide range of $Re$ and valid at $Pr\approx 1$.

We then examine flow behaviour at different $Ri_b$ and $Pr$. For flows with $Pr\ll 1$, a larger value of $Ri_b$ is required to reduce the amplification of streamwise-varying flow structures to the same level as streamwise-invariant ones compared with flows with $Pr\approx 1$. The largest amplification also occurs at the same value of $Pr Ri_b$, consistent with the observation of matching averaged density profile for flows with the same value of $PrRi_b$ in the $Pr\ll 1$ regime (Langham *et al.* Reference Langham, Eaves and Kerswell2020). For flows with $Pr\gg 1$, the SIOA identifies another quasi-horizontal flow structure that is independent of $Ri_b$. By decomposing input–output pathways into each velocity and density component, we show that these quasi-horizontal flow structures for flows with $Pr\gg 1$ are associated with density fluctuations. The importance of this density-associated flow structure for flows with $Pr\gg 1$ is further highlighted through a derived analytical scaling of amplification with respect to $Re$ and $Pr$ under the assumptions that the flow is streamwise invariant ($k_x=0$) and unstratified (i.e. $Ri_b=0$ and the density behaves as a passive scalar). The above observations using SIOA distinguish two types of quasi-horizontal flow structures, one emerging in the high-$Re$ high-$Ri_b$ regime and the other one (associated with density fluctuations) emerging in the high $Pr$ regime.

The results here suggest the promise of this computationally tractable approach in identifying horizontal length scales of prominent flow structures in stratified wall-bounded shear flows and opens up many directions for future work. For example, this framework may be extended to other stratified wall-bounded shear flows such as stratified channel flow (Garcia-Villalba & del Alamo Reference Garcia-Villalba and del Alamo2011), stratified open channels (Flores & Riley Reference Flores and Riley2011; Brethouwer, Duguet & Schlatter Reference Brethouwer, Duguet and Schlatter2012; Donda *et al.* Reference Donda, Van Hooijdonk, Moene, Jonker, van Heijst, Clercx and van de Wiel2015; He & Basu Reference He and Basu2015, Reference He and Basu2016) and the stratified Ekman layer (Deusebio *et al.* Reference Deusebio, Brethouwer, Schlatter and Lindborg2014), where intermittent regimes of flow dynamics were also observed. This framework may also be extended to configurations where the background density gradient and velocity gradient are orthogonal, e.g. spanwise stratified PCF (Facchini *et al.* Reference Facchini, Favier, Le Gal, Wang and Le Bars2018; Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2019) and spanwise stratified plane Poiseuille flow (Le Gal *et al.* Reference Le Gal, Harlander, Borcia, Le Dizès, Chen and Favier2021).

## Funding

C.L. and D.F.G. gratefully acknowledge support from the US National Science Foundation (NSF) through grant number CBET 1652244. C.L. also greatly appreciates the support from the Chinese Scholarship Council.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Growth rate analysis

Here, we present the growth rate of the dynamics in (2.6) computed as

where $\text {eig}(\,\cdot\,)$ is the eigenvalue of the argument, ${\rm Re}[\,\cdot\,]$ represents the real part, $\text {max}\{\,\cdot \,\}$ is the maximum value of the argument and $\hat {\boldsymbol{\mathsf{A}}}^S$ is the discretization of operator $\hat {\mathcal {A}}^S$. Figure 11 shows the growth rate $R(\hat {\mathcal {A}}^S)(k_x,k_z)$ in (A1). Here, we observe that this modal growth rate analysis $R(\hat {\mathcal {A}}^S)(k_x,k_z)$ cannot distinguish a preferential structure size over a wide range of wavenumbers $k_x\lesssim 1$ and $k_z\lesssim 10$, and there is no identified instability consistent with Davey & Reid (Reference Davey and Reid1977).

## Appendix B. Proofs of theorems 4.1–4.2

#### B.1. Proof of theorem 4.1

Proof. The proof of theorem 4.1 naturally follows the procedure in unstratified flow (Jovanović Reference Jovanović2004; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Jovanović Reference Jovanović2021) and is outlined as a block diagram in figure 12. Under the assumption of streamwise invariance ($k_x=0$) and taking the passive scalar limit ($Ri_b=0$) for stratified PCF in theorem 4.1, the operators $\hat {\mathcal {A}}^S$, $\hat {\mathcal {B}}^S$ and $\hat {\mathcal {C}}^S$ can be simplified and their non-zero elements can be defined as

*a*)$$\begin{gather} \hat{\mathcal{A}}^S(k_x,k_z)=\begin{bmatrix} \dfrac{\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^4}{Re} & 0 & 0\\ -\text{i}k_zU' & \dfrac{\widehat{{\nabla}}^2}{Re} & 0\\ -\bar{\rho}' & 0 & \dfrac{\widehat{\nabla}^2}{Re Pr} \end{bmatrix}, \end{gather}$$

*b*)$$\begin{gather}\mathcal{\hat{B}}^S(k_x,k_z)= \begin{bmatrix} 0 & -k_z^2\widehat{{\nabla}}^{{-}2} & -\text{i}k_z\widehat{{\nabla}}^{{-}2} \partial_y & 0\\ \text{i}k_z & 0 & 0 & 0\\ 0 & 0 & 0 & \mathcal{I} \end{bmatrix}=:\begin{bmatrix} 0 & \hat{\mathcal{B}}^S_{y,1} & \hat{\mathcal{B}}^S_{z,1} & 0\\ \mathcal{B}^S_{x,2} & 0 & 0 & 0\\ 0 & 0 & 0 & \hat{\mathcal{B}}^S_{\rho,3} \end{bmatrix} \end{gather}$$

*c*)$$\begin{gather}\mathcal{\hat{C}}^S(k_x,k_z)=\begin{bmatrix} 0 & -\text{i}/k_z & 0 \\ \mathcal{I} & 0 & 0\\ \text{i} \partial_y/k_z & 0 & 0 \\ 0 & 0 & \mathcal{I} \end{bmatrix}=:\begin{bmatrix} 0 & \mathcal{\hat{C}}^S_{u,2} & 0 \\ \mathcal{\hat{C}}^S_{v,1} & 0 & 0\\ \mathcal{\hat{C}}^S_{w,1} & 0 & 0 \\ 0 & 0 & \mathcal{\hat{C}}^S_{\rho,3} \end{bmatrix} \end{gather}$$

We employ a matrix inverse formula for the lower triangle block matrix

to compute $(\text {i}\omega \mathcal {I}_{3\times 3}-\hat {\mathcal {A}}^S)^{-1}$. Then, we employ a change of variable $\varOmega _1=\omega Re$ and $\varOmega _2=\omega RePr$ to obtain componentwise frequency response operators $\mathcal {H}^S_{ij}$ with $i=u,v,w,\rho$ and $j=x,y,z,\rho$ as

*a*)\begin{gather} \mathcal{H}^S_{ux}=\hat{\mathcal{C}}^S_{u,2}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}\hat{\mathcal{B}}^S_{x,2}, \end{gather}

*b*)\begin{gather}\mathcal{H}^S_{uy}=\hat{\mathcal{C}}^S_{u,2}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}(-\text{i}k_zU')Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{y,1}, \end{gather}

*c*)\begin{gather}\mathcal{H}^S_{uz}=\hat{\mathcal{C}}^S_{u,2}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}(-\text{i}k_zU')Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{z,1}, \end{gather}

*d*)\begin{gather}\mathcal{H}^S_{u\rho}=0, \end{gather}

*e*)\begin{gather}\mathcal{H}^S_{vx}=0, \end{gather}

*f*)\begin{gather}\mathcal{H}^S_{vy}=\hat{\mathcal{C}}^S_{v,1}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{y,1}, \end{gather}

*g*)\begin{gather}\mathcal{H}^S_{vz}=\hat{\mathcal{C}}^S_{v,1}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{z,1}, \end{gather}

*h*)\begin{gather}\mathcal{H}^S_{v\rho}=0, \end{gather}

*i*)\begin{gather}\mathcal{H}^S_{wx}=0, \end{gather}

*j*)\begin{gather}\mathcal{H}^S_{wy}=\hat{\mathcal{C}}^S_{w,1}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{y,1}, \end{gather}

*k*)\begin{gather}\mathcal{H}^S_{wz}=\hat{\mathcal{C}}^S_{w,1}Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{z,1}, \end{gather}

*l*)\begin{gather}\mathcal{H}^S_{w\rho}=0, \end{gather}

*m*)\begin{gather}\mathcal{H}^S_{\rho x}=0, \end{gather}

*n*)\begin{gather}\mathcal{H}^S_{\rho y}=\hat{\mathcal{C}}^S_{\rho,3}RePr\left(\text{i}\varOmega_2 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}(-\bar{\rho}')Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{y,1}, \end{gather}

*o*)\begin{gather}\mathcal{H}^S_{\rho z}=\hat{\mathcal{C}}^S_{\rho,3}RePr\left(\text{i}\varOmega_2 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}(-\bar{\rho}')Re\left(\text{i}\varOmega_1 \mathcal{I}-\widehat{{\nabla}}^{{-}2}\widehat{{\nabla}}^{4}\right)^{{-}1} \hat{\mathcal{B}}^S_{z,1}, \end{gather}

*p*)\begin{gather}\mathcal{H}^S_{\rho \rho}=\hat{\mathcal{C}}^S_{\rho,3}RePr\left(\text{i}\varOmega_2 \mathcal{I}-\widehat{{\nabla}}^2\right)^{{-}1}\hat{\mathcal{B}}^S_{\rho,3}. \end{gather}

Taking the operation that $\|\,\cdot\,\|_{\infty }=\underset {\omega \in \mathbb {R}}{\text {sup}}\bar {\sigma }[\,\cdot\,]=\underset {\varOmega _1\in \mathbb {R}}{\text {sup}}\bar {\sigma }[\,\cdot\,]=\underset {\varOmega _2\in \mathbb {R}}{\text {sup}}\bar {\sigma }[\,\cdot\,]$, we obtain the scaling relation in theorem 4.1(*a*).

Using the relation that $\mathcal {H}^S_{\boldsymbol {\nabla } ij}=\widehat {\boldsymbol {\nabla }}\mathcal {H}^S_{ij}$ in (3.7) with $i=u,v,w,\rho$, and $j=x,y,z,\rho$, and similarly employing the definition of $\|\,\cdot\,\|_{\infty }$, we obtain the scaling relation in theorem 4.1(*b*).

In figure 12, the block $-\text {i}k_zU'$ inside the dashed line ($\text {- -}$, red) contributes to the relatively large scalings of $\|\mathcal {H}^S_{uy}\|_{\infty }\sim Re^2$, $\|\mathcal {H}^S_{uz}\|_{\infty }\sim Re^2$ at high $Re$ in (4.1) of theorem 4.1(*a*), which has been attributed to the lift-up mechanism; see discussion in Jovanović (Reference Jovanović2021). Similarly, the block $-\bar {\rho }'$ outlined by ($\text {-}\cdot \text {-}$, blue) contributes to the relatively large scalings of $\|\mathcal {H}^S_{\rho y}\|_{\infty }\sim Re^2Pr$, and $\|\mathcal {H}^S_{\rho z}\|_{\infty }\sim Re^2Pr$ at high $Re$ or $Pr$. This similarity between streamwise streaks and density streaks is consistent with the observation that passive scalar streaks can be generated by the same lift-up mechanism as the streamwise streaks (Chernyshenko & Baig Reference Chernyshenko and Baig2005).

#### B.2. Proof of theorem 4.2

Proof. We define the set of uncertainties

*a*)\begin{gather} \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,ux}:=\left\{\text{diag}\left(-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}}\right):-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}\right\}, \end{gather}

*b*)\begin{gather}\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,vy}:=\left\{\text{diag}\left(\boldsymbol{\mathsf{0}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}}\right):-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}\right\}, \end{gather}

*c*)\begin{gather}\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,wz}:=\left\{\text{diag}\left(\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},\boldsymbol{\mathsf{0}}\right):-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}\right\}, \end{gather}

*d*)\begin{gather}\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,\rho\rho}:=\left\{\text{diag}\left(\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{0}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}}\right):-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}\right\}. \end{gather}

Here, $\boldsymbol{\mathsf{0}}\in \mathbb {C}^{N_y\times 3N_y}$ is a zero matrix with the size $N_y\times 3N_y$. Then, using the definition of the structured singular value in definition 2.1, we have

*a*)\begin{align} & \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,ux}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]\nonumber\\ &\quad =\frac{1}{\text{min}\{\bar{\sigma}[\hat{\boldsymbol{\mathsf{u}}}^S_{\varXi,ux}]: \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi,ux}\in \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,ux}, \text{det}[\boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega) \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi,ux}]=0\}} \end{align}

*b*)\begin{align} &\quad =\frac{1}{\text{min}\{\bar{\sigma}[-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}]: -\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi}\in \mathbb{C}^{N_y\times 3N_y}, \text{det}[\boldsymbol{\mathsf{I}}_{3N_y}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla} ux}(k_x,k_z,\omega) (-\hat{\boldsymbol{\mathsf{u}}}^{\text{T}}_{\xi})]=0\}} \end{align}

*c*)\begin{align} &\quad =\bar{\sigma}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla} ux}(k_x,k_z,\omega)]. \end{align}

Here, equality (B5*a*) is obtained by substituting the uncertainty set in (B4*a*) into definition 2.1. The equality (B5*b*) is obtained by performing a block-diagonal partition of terms inside of $\text {det}[\,\cdot\,]$ and employing zeros in the uncertainty set in equation (B4*a*). Here, $\boldsymbol{\mathsf{H}}^S_{\boldsymbol {\nabla } ux}$ is the discretization of $\mathcal {H}^S_{\boldsymbol {\nabla } ux}$ and $\boldsymbol{\mathsf{I}}_{3N_y}\in \mathbb {C}^{3N_y\times 3N_y}$ in (B5*b*) is an identity matrix with matching size $(3N_y\times 3N_y)$, where we use the subscripts to distinguish it from $\boldsymbol{\mathsf{I}} \in \mathbb {C}^{12N_y\times 12N_y}$ in (B5*a*). The equality (B5*c*) uses the definition of the unstructured singular value; see e.g. Zhou *et al.* (Reference Zhou, Doyle and Glover1996, (11.1)).

Similarly, we have

*a*)\begin{gather} \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,vy}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]=\bar{\sigma}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla} vy}(k_x,k_z,\omega)], \end{gather}

*b*)\begin{gather}\mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,wz}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]=\bar{\sigma}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla} wz}(k_x,k_z,\omega)], \end{gather}

*c*)\begin{gather}\mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi,\rho\rho}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]=\bar{\sigma}[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla} \rho \rho}(k_x,k_z,\omega)]. \end{gather}

Using the fact that $\hat {\boldsymbol{\mathsf{U}}}^S_{\varXi }\supseteq \hat {\boldsymbol{\mathsf{U}}}^S_{\varXi,ij}$ with $ij=ux,vy,wz,\rho \rho$ and equalities in (B5)–(B6), we have

Applying the supreme operation $\underset {\omega \in \mathbb {R}}{\text {sup}}[\,\cdot\,]$ on (B7) and using the definitions of $\|\,\cdot\,\|_{\mu }$ and $\|\,\cdot\,\|_{\infty }$ we have

*a*,

*b*)