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.

Figure 1. (a) Schematic of stably stratified PCF with laminar base flow $\boldsymbol {U}(y)=y\boldsymbol {e}_x$ and background density $\bar {\rho }=-y$. The gravity vector $\boldsymbol {g}=-g\boldsymbol {e}_y$ is orthogonal to the channel walls. The blue and red colours represent high-density and low-density fluids, respectively. (b) Block diagram of the feedback interconnection between the linearized dynamics and structured forcing (outlined by the blue dashed line) in (2.4) modelling the nonlinearity.

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

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

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

(2.1c)\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

(2.2a–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 Ioannou1993a), 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

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

(2.3b)\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

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

(2.4b)\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.1c) 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:

(2.5)\begin{equation} \hat{\psi}(y;k_x,k_z,\omega):=\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\psi(x,y,z,t)\exp({-\text{i}(k_x x + k_z z+\omega t)})\,{{\rm d}\kern0.7pt x}\,{\rm d} z\,{\rm d} t,\end{equation}

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

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

(2.6b)\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

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

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

(2.7c)\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.7a) 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

(2.8)\begin{equation} \mathcal{H}^S(y;k_x,k_z,\omega):=\hat{\mathcal{C}}^S \left(\text{i}\omega\mathcal{I}_{3\times 3}-\hat{\mathcal{A}}^S\right)^{{-}1}\hat{\mathcal{B}}^S. \end{equation}

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.4a)–(2.4b) allows us to perform the same spatio-temporal Fourier transform on the model of the nonlinearity, which can be decomposed as

(2.9)\begin{equation} \begin{bmatrix} \hat{f}_{x,\xi}\\ \hat{f}_{y,\xi}\\ \hat{f}_{z,\xi}\\ \hat{f}_{\rho,\xi} \end{bmatrix}=\text{diag} \left(-\hat{\boldsymbol{u}}_{\xi}^{\text{T}},-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} \right) \text{diag}\left(\widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}, \widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}\right) \begin{bmatrix} \hat{u}\\ \hat{v}\\ \hat{w}\\ \hat{\rho} \end{bmatrix}. \end{equation}

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.

Figure 2. Block diagram showing the feedback interconnection of the SIOA framework applied to stratified PCF. (a) Redraws figure 1(b), where the blocks inside of ($\text {- -}$, blue) lines represent the modelled forcing in (2.9) corresponding to the bottom block in figure 1(b) also inside of ($\text {- -}$, blue). (b) Redraws (a) after discretization with the top block corresponding to the combination of the two top blocks in (a) and the bottom block corresponding to the bottom block of (a).

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

(2.10)\begin{equation} \mathcal{H}^S_{\boldsymbol{\nabla}}(y;k_x,k_z,\omega):=\text{diag} \left(\widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}, \widehat{\boldsymbol{\nabla}},\widehat{\boldsymbol{\nabla}}\right) \mathcal{H}^S(y;k_x,k_z, \omega). \end{equation}

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

(2.11)\begin{equation} \hat{\boldsymbol{u}}^S_{\varXi}:=\text{diag}\left(-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} ,-\hat{\boldsymbol{u}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{u}}_{\xi}^{\text{T}} \right). \end{equation}

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

(2.12)\begin{equation} \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}} \left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right]:= \frac{1}{\min\{\bar{\sigma}[\hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}]: \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}\in \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}, \det[\boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega) \hat{\boldsymbol{\mathsf{u}}}^S_{\varXi}]=0\}}.\end{equation}

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.

(2.13)\begin{equation} \hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}:=\left\{{\rm {diag}} \left(-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},-\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}}, -\hat{\boldsymbol{\mathsf{u}}}_{\xi}^{\text{T}},-\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{equation}

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

(2.14)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z):=\underset{\omega \in \mathbb{R}}{\text{sup}} \mu_{\hat{\boldsymbol{\mathsf{U}}}^S_{\varXi}}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right], \end{equation}

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

(2.15)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\infty}(k_x,k_z):=\underset{\omega \in \mathbb{R}}{\text{sup}} \bar{\sigma}\left[\boldsymbol{\mathsf{H}}^S_{\boldsymbol{\nabla}}(k_x,k_z,\omega)\right].\end{equation}

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 6a), 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.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

(3.1)\begin{equation} \|\mathcal{H}^S\|_{\infty}(k_x,k_z):= \underset{\omega \in \mathbb{R}}{\text{sup}} \bar{\sigma}\left[\boldsymbol{\mathsf{H}}^S(k_x,k_z,\omega)\right].\end{equation}

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 2b) 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.

Figure 3. Plots of (a) $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$, (b) $\log _{10}[\|\mathcal {H}^S\|_{\infty }(k_x,k_z)]$ and (c) $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\infty }(k_x,k_z)]$ for stratified PCF at $Re=865$, $Ri_b=0.02$ and $Pr=0.7$. Here, the symbols (${\triangle }$) are characteristic wavelengths ($\lambda _x=32{\rm \pi}$, $\lambda _z=16{\rm \pi}$) corresponding to the oblique turbulent band observed in DNS in the same flow regime (Deusebio et al. Reference Deusebio, Caulfield and Taylor2015; Taylor et al. Reference Taylor, Deusebio, Caulfield and Kerswell2016). The lines (–) are $\lambda _z=\lambda _x\tan (27^\circ )$ indicating the $27^\circ$ angle of the oblique turbulent bands.

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 4b). 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$).

Figure 4. Plots of $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$ at $Pr=0.7$ and: (a) $Re=4250$ and $Ri_b=0.02$; (b) $Re=865$ and $Ri_b=0.2541$ and (c) $Re=52\,630$ and $Ri_b=0.15$.

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

(3.2)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}^M:=\underset{k_z,k_x}{\text{max}} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z),\end{equation}

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

(3.3)\begin{equation} \|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}^{sc}:= \underset{k_z,k_x=10^{{-}6}}{\text{max}}\|\mathcal{H}^S_{\boldsymbol{\nabla}}\|_{\mu}(k_x,k_z). \end{equation}

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

Figure 5. The dependence on $Ri_b$ and $Re$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ (–, black) and $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ ($\text {- -}$, red) at $Pr=0.7$. Each black marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{M}$ indicates the associated value of (a) $Re$ and (b) $Ri_b$ defined in the legend. Each red marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ indicates the same parameter value for (a) $Re$ and (b) $Ri_b$ as the corresponding black one.

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

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

(3.4c,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 Caulfield2017a,Reference Zhou, Taylor, Caulfield and Lindenb; 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 Caulfield2017a,Reference Zhou, Taylor, Caulfield and Lindenb). 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 Caulfield2017a,Reference Zhou, Taylor, Caulfield and Lindenb).

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

Figure 6. The dependence on $PrRi_b$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ at $Re=4250$.

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.

Figure 7. The dependence on $Ri_b$ and $Pr$ of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ (–, black) and $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ ($\text {- -}$, red) at $Re=4250$. Each black marker on the lines of $\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }^M$ indicates the associated value of (a) $Pr$ and (b) $Ri_b$ indicated in the legend. Each red marker on the lines of $\|\mathcal {H}^{S}_{\boldsymbol {\nabla }}\|_{\mu }^{sc}$ indicates the same parameter as the corresponding black marker.

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

Figure 8. Plots of $\log _{10}[\|\mathcal {H}^S_{\boldsymbol {\nabla }}\|_{\mu }(k_x,k_z)]$ at $Re=4250$, $Ri_b=0.02$ with three different Prandtl numbers at (a) $Pr=10^{-4}$, (b) $Pr=7$ and (c) $Pr=70$.

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 4c) 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

(3.5)\begin{equation} \mathcal{H}^S_{ij}=\hat{\mathcal{C}}^S_{i}\left(\text{i}\omega \mathcal{I}_{3\times 3}-\hat{\mathcal{A}}^S\right)^{{-}1}\hat{\mathcal{B}}^S_{j},\end{equation}

with

(3.6a,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}

(3.6c,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}

(3.6e,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}

(3.6g,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.4a)–(2.4b) 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

(3.7)\begin{equation} \mathcal{H}_{\boldsymbol{\nabla} ij}^S=\widehat{\boldsymbol{\nabla}}\mathcal{H}_{ij}^S. \end{equation}

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

Figure 9. Componentwise values of $\log _{10}[\|\mathcal {H}^S_{ij}\|_{\infty }](k_x,k_z)$ at $Re=865$, $Ri_b=0.02$, $Pr=0.7$ ($ij=ux,vy,wz,\rho \rho$). Note the range of the colour bar for each panel is modified based on the maximum value of the particular quantity: (a) $\|\mathcal {H}^S_{u x}\|_{\infty }$; (b) $\|\mathcal {H}^S_{v y}\|_{\infty }$; (c) $\|\mathcal {H}^S_{w z}\|_{\infty }$; (d) $\|\mathcal {H}^S_{\rho \rho }\|_{\infty }$.

Figure 10. Componentwise values of $\log _{10}[\|\mathcal {H}^S_{ij}\|_{\infty }](k_x,k_z)$ at $Re=4250$, $Ri_b=0.02$, $Pr=70$ ($ij=ux,vy,wz$, and $\rho \rho$). Note the range of the colour bar for each panel is modified based on the maximum value of the particular quantity: (a) $\|\mathcal {H}^S_{u x}\|_{\infty }$; (b) $\|\mathcal {H}^S_{v y}\|_{\infty }$; (c) $\|\mathcal {H}^S_{wz}\|_{\infty }$; (d) $\|\mathcal {H}^S_{\rho \rho }\|_{\infty }$.