2. Direct-adjoint-looping (DAL) method

We consider stably stratified Boussinesq PCF with a linear equation of state in which fluid flows between two parallel horizontal plates moving in opposite directions with relative speed $2{\rm\Delta}U$ , separated by a distance $2H$ , across which a constant density difference $2{\rm\Delta}{\it\rho}$ is maintained from a reference density ${\it\rho}_{r}\gg {\rm\Delta}{\it\rho}$ . Non-dimensionalising with respect to ${\rm\Delta}U$ , $H$ and ${\rm\Delta}{\it\rho}$ , and decomposing the total velocity and density fields as $(\boldsymbol{u}_{\mathit{tot}},{\it\rho}_{\mathit{tot}})=(\boldsymbol{U},\bar{{\it\rho}})+(\boldsymbol{u},{\it\rho})$ , where $(\boldsymbol{U},\bar{{\it\rho}})$ is the laminar state, we obtain

(2.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\partial _{t}\boldsymbol{u}+(\boldsymbol{u}+\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\boldsymbol{u}+\boldsymbol{U})=-\boldsymbol{{\rm\nabla}}p-\mathit{Ri}_{B}{\it\rho}\hat{\boldsymbol{y}}+\mathit{Re}^{-1}{\rm\nabla}^{2}\boldsymbol{u},\quad (\boldsymbol{u},{\it\rho})(y=\pm 1)=(\mathbf{0},0)\\ \partial _{t}{\it\rho}+(\boldsymbol{u}+\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}({\it\rho}+\bar{{\it\rho}})=(\mathit{Re}\mathit{Pr})^{-1}{\rm\nabla}^{2}{\it\rho},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \boldsymbol{U}=y\hat{\boldsymbol{x}},\quad \bar{{\it\rho}}=-y\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

with the three non-dimensional parameters: Reynolds number $\mathit{Re}$ ; Prandtl number $\mathit{Pr}$ ; and bulk Richardson number $\mathit{Ri}_{B}$ defined as

(2.2a−c ) $$\begin{eqnarray}\mathit{Re}=\frac{{\rm\Delta}UH}{{\it\nu}},\quad \mathit{Pr}=\frac{{\it\nu}}{{\it\kappa}},\quad \mathit{Ri}_{B}=\frac{g{\rm\Delta}{\it\rho}H}{{\it\rho}_{r}{\rm\Delta}U^{2}},\end{eqnarray}$$

where ${\it\nu}$ is the kinematic viscosity and ${\it\kappa}$ is the diffusivity of density. For the calculations presented here, following Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) we choose $\mathit{Re}=1000$ , set $\mathit{Pr}=1$ , and vary $\mathit{Ri}_{B}$ .

We define a nonlinear optimal perturbation as an initial condition $(\boldsymbol{u},{\it\rho})(t=0)=(\boldsymbol{u}_{0},{\it\rho}_{0})$ of given initial total energy (kinetic energy plus potential energy) $E_{0}$ , defined as

(2.3) $$\begin{eqnarray}E_{0}=\langle \boldsymbol{u}_{0}\boldsymbol{\cdot }\boldsymbol{u}_{0}+\mathit{Ri}_{B}{\it\rho}_{0}^{2}\rangle /2,\quad \text{where }\langle \boldsymbol{a},\boldsymbol{b}\rangle \equiv \frac{1}{V}\int _{V}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{b}\,\text{d}V,\end{eqnarray}$$

that maximises, over a given time horizon $T$ , a given quantity of interest ${\mathcal{J}}(\boldsymbol{u},{\it\rho},T)$ . Here, $V$ is the volume of the computational domain. In order to find minimal seeds for turbulence, which are initial conditions whose trajectories eventually transition to turbulence, we typically choose an appropriately large large value $T=300$ , and we consider a functional ${\mathcal{J}}$ that takes heightened values in the turbulent state. With the presence of a density field, it is possible to find large amplitude waves that have an instantaneously large total energy, but are not a turbulent state. Therefore, rather than choosing the total energy density at $T$ as ${\mathcal{J}}$ , as was done by Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) for unstratified minimal seeds, we choose the total time-averaged dissipation of total energy density, the stratified generalisation of the objective functional chosen by Monokrousos et al. (Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011), for an equivalent minimal seed calculation.

We thus write

(2.4) $$\begin{eqnarray}{\mathcal{J}}=(T\mathit{Re})^{-1}[\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\mathit{Ri}_{B}\mathit{Pr}^{-1}\boldsymbol{{\rm\nabla}}{\it\rho}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}],\end{eqnarray}$$

where $[\boldsymbol{a},\boldsymbol{b}]\equiv \int _{0}^{T}\langle \boldsymbol{a},\boldsymbol{b}\rangle \,\text{d}t$ , and the maximisation of ${\mathcal{J}}$ can be conducted by taking variations of the augmented and constrained functional  ${\mathcal{L}}$ :

(2.5) $$\begin{eqnarray}\displaystyle {\mathcal{L}} & = & \displaystyle {\mathcal{J}}(\boldsymbol{u},{\it\rho},T)-[\partial _{t}\boldsymbol{u}+N(\boldsymbol{u})+\boldsymbol{{\rm\nabla}}p+\mathit{Ri}_{B}{\it\rho}\hat{\boldsymbol{y}}-\mathit{Re}^{-1}{\rm\nabla}^{2}\boldsymbol{u},\boldsymbol{v}]-[\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u},q]\nonumber\\ \displaystyle & & \displaystyle -\,[\partial _{t}{\it\rho}+(\boldsymbol{u}+\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\bar{{\it\rho}}+{\it\rho})-(\mathit{Re}\mathit{Pr})^{-1}{\rm\nabla}^{2}{\it\rho},{\it\eta}]+\langle \boldsymbol{u}_{0}-\boldsymbol{u}(0),\boldsymbol{v}_{0}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\langle {\it\rho}_{0}-{\it\rho}(0),{\it\eta}_{0}\rangle -(\langle \boldsymbol{u}_{0},\boldsymbol{u}_{0}\rangle /2+\langle \mathit{Ri}_{B}{\it\rho}_{0}^{2}\rangle /2-E_{0})c,\end{eqnarray}$$

where $N(\boldsymbol{u})=(\boldsymbol{u}+\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(\boldsymbol{u}+\boldsymbol{U})$ . The Lagrange multipliers $\boldsymbol{v}$ , $q$ and ${\it\eta}$ are termed the adjoint velocity, pressure and density and together enforce the Boussinesq Navier–Stokes equations (2.1) on $\boldsymbol{u}$ and ${\it\rho}$ . The initial conditions and initial energy are enforced by $\boldsymbol{v}_{0}$ , ${\it\eta}_{0}$ and $c$ . Taking variations with respect to $\boldsymbol{u}$ , ${\it\rho}$ , $\boldsymbol{u}_{0}$ and ${\it\rho}_{0}$ yields the following system of ‘adjoint’ or ‘dual’ equations that must also be satisfied at all times and points in space by a nonlinear optimal perturbation:

(2.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\partial _{t}\boldsymbol{v}+N^{\dagger }(\boldsymbol{v},\boldsymbol{u})+\mathit{Re}^{-1}{\rm\nabla}^{2}\boldsymbol{v}+\boldsymbol{{\rm\nabla}}q-{\it\eta}\boldsymbol{{\rm\nabla}}(\bar{{\it\rho}}+{\it\rho})-(\mathit{Re}T)^{-1}{\rm\nabla}^{2}\boldsymbol{u}=\mathbf{0},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0,\\ \partial _{t}{\it\eta}+(\boldsymbol{U}+\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\eta}+(\mathit{Re}\mathit{Pr})^{-1}{\rm\nabla}^{2}{\it\eta}-\mathit{Ri}_{B}\,\hat{\boldsymbol{y}}\boldsymbol{\cdot }\boldsymbol{v}-\mathit{Ri}_{B}(\mathit{Re}\mathit{Pr}T)^{-1}{\rm\nabla}^{2}{\it\rho}=0,\end{array}\right\}\vphantom{\mathop{\sum }_{\sum }}\end{eqnarray}$$

(2.7a,b ) $$\begin{eqnarray}\boldsymbol{v}(T)=\mathbf{0},\quad {\it\eta}(T)=0,\end{eqnarray}$$

(2.8a−c ) $$\begin{eqnarray}\boldsymbol{v}(0)-c\boldsymbol{u}_{0}=\mathbf{0},\quad {\it\eta}(0)-c{\it\rho}_{0}\mathit{Ri}_{B}=0,\quad \langle \boldsymbol{u}_{0},\boldsymbol{u}_{0}\rangle +\langle \mathit{Ri}_{B}{\it\rho}_{0}^{2}\rangle -2E_{0}=0,\end{eqnarray}$$

where $N^{\dagger }(v_{i},\boldsymbol{u})=\partial _{j}((U_{j}+u_{j})v_{i})-v_{j}\partial _{i}(U_{j}+u_{j})$ .

The first step of the DAL method to find a nonlinear optimal perturbation is to choose an initial condition guess $(\boldsymbol{u}_{0},{\it\rho}_{0})$ of energy density $E_{0}$ . This initial condition is then integrated forwards in time to $t=T$ using the direct Boussinesq Navier–Stokes equations (2.1). These flow fields are stored, and used to integrate backwards in time the adjoint fields $\boldsymbol{v}$ and ${\it\eta}$ from the null ‘end’ conditions (2.7) using the adjoint equations (2.6), which actually depend on the direct fields $(\boldsymbol{u},{\it\rho})$ . Once values for the adjoint variables are obtained at $t=0$ , we have compatibility conditions (2.8) relating $\boldsymbol{v}(0)$ to $\boldsymbol{u}_{0}$ and ${\it\eta}(0)$ to ${\it\rho}_{0}$ that must be satisfied by a nonlinear optimal perturbation. If not satisfied, these compatibility conditions yield gradient information for the objective functional ${\mathcal{J}}$ with respect to changes in the initial condition $(\boldsymbol{u}_{0},{\it\rho}_{0})$ , allowing a refined guess for the nonlinear optimal perturbation to be made. This sequence of steps is continued until convergence.

For large $T$ , any initial condition that eventually becomes turbulent will be a turbulent seed using this method. The minimal seed, however, has the special property of having the lowest critical initial energy density $E_{0}=E_{c}$ of all such seeds. We first identify a perturbation that causes turbulence, with initial energy density $E_{0}=E_{{\it\alpha}}\gg E_{c}$ . This initial condition is then used as an initial guess for the DAL method at a smaller initial energy density $E_{0}=E_{{\it\beta}}<E_{{\it\alpha}}$ , with uniformly rescaled energy. The DAL method then finds a more efficient route to turbulence at $E_{0}=E_{{\it\beta}}$ . This more efficient initial condition is then uniformly rescaled to have energy density $E_{0}=E_{{\it\gamma}}<E_{{\it\beta}}$ , and the process is repeated. This ‘laddering down’ is continued until $E_{0}=E_{c}$ , at which point any further reduction in $E_{0}$ cannot find an initial condition leading to turbulence.

3. Stratified minimal seeds for turbulence

We investigate two geometries, namely the ‘narrow’ geometry ‘N’ investigated by Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) in the unstratified case, $4.35{\rm\pi}\times 2\times 1.05{\rm\pi}$ , one of the geometries considered by Butler & Farrell (Reference Butler and Farrell1992) for linear, unstratified optimal perturbations, in which the unstratified minimal seeds are localised in the streamwise direction but essentially fill the spanwise width, and $4.35{\rm\pi}\times 2\times 2.10{\rm\pi}$ , a ‘wide’ geometry ‘W’, which is twice as wide, and allows for localisation in the spanwise direction also. We used a modified version of the parallelised CFD solver Diablo (Taylor Reference Taylor2008) to solve the forward and adjoint equations, which uses Fourier modes in the streamwise $x$ and spanwise $z$ directions, and finite differences in the wall-normal $y$ direction, and a combined implicit–explicit Runge–Kutta–Wray Crank–Nicholson time integration scheme. The resolution for geometry N was $128\times 256\times 32$ and for geometry W was $128\times 256\times 64$ .

Using the laddering down approach described above, we have converged to the minimal seed for $\mathit{Ri}_{B}=0$ (confirming quantitatively the unstratified case of Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) using a different code, and objective functional ${\mathcal{J}}$ ), $\mathit{Ri}_{B}=10^{-4}$ , $10^{-3}$ , $3\times 10^{-3}$ and $10^{-2}$ in geometry N and $\mathit{Ri}_{B}=0$ , $3\times 10^{-3}$ and $10^{-2}$ in geometry W, using the fixed value $T=300$ for all bulk Richardson numbers except for the largest, $\mathit{Ri}_{B}=10^{-2}$ , for which we use $T=400$ . The respective values of $E_{c}(\mathit{Ri}_{B})$ are shown in table 1. We immediately see that $E_{c}$ is an increasing function of $\mathit{Ri}_{B}$ , as expected, since a stable stratification inhibits vertical motions, and so a transition process involving vertical motions should be expected to require a larger energy input. Interestingly, $E_{c}$ in geometry W is approximately 40 % of $E_{c}$ in geometry N for the same $Ri_{B}$ . Since $E_{c}$ , as defined in (2.3), is an energy density, and the volume of geometry W is twice that of geometry N, this suggests that the minimal seeds in geometry W are spanwise localised, and that the narrow geometry N actually requires higher maximum amplitudes of perturbation in the minimal seed due to the enforced spanwise periodicity.

Table 1. Values of the critical energy density $E_{c}$ , the energy of the minimal seed, for various bulk Richardson numbers $\mathit{Ri}_{B}$ in the two geometries N and W. The upper bound corresponds to the flow evolutions shown in subsequent figures, and in supplementary movies 1 and 2 available at http://dx.doi.org/10.1017/jfm.2015.596. The lower bound corresponds to an $E_{0}$ at which a turbulent state cannot be attained.

3.1. Minimal seeds in narrow geometry N

Figure 2 shows the time evolution from the minimal seed initial conditions for geometry N of the total energy density $E(t)$ defined as

(3.1) $$\begin{eqnarray}E(t)={\textstyle \frac{1}{2}}\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}+\mathit{Ri}_{B}{\it\rho}^{2}\rangle =K(t)+P(t),\end{eqnarray}$$

where $K(t)$ is the kinetic energy, $P(t)$ is the potential energy, and angled brackets denote volume averaging as defined in (2.3). Streamwise velocity $u=\pm 0.6\max (u)$ isosurfaces are plotted in figures 3–5 at times $t=0$ and 25, $t=70$ and 150, and $t=210$ and 280 respectively, and videos of the flow evolution are available as supplementary material.

Figure 2. Time variation of energy density $E(t)$ as defined in (3.1) for the minimal seed trajectories in geometry N for $\mathit{Ri}_{B}=0$ , black; $10^{-4}$ , red; $10^{-3}$ , blue; $3\times 10^{-3}$ , green; and $10^{-2}$ , purple.

Figure 3. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=0$ (a–e) and $t=25$ (f–j) for the minimal seed trajectories in geometry N for $\mathit{Ri}_{B}=0$ , $10^{-4}$ , $10^{-3}$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

Figure 4. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=70$ (a–e) and $t=150$ (f–j) for the minimal seed trajectories in geometry N for $\mathit{Ri}_{B}=0$ , $10^{-4}$ , $10^{-3}$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

Figure 5. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=210$ (a–e) and $t=280$ (f–j) for the minimal seed trajectories in geometry N for $\mathit{Ri}_{B}=0$ , $10^{-4}$ , $10^{-3}$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

The unstratified minimal seed, as shown in the leftmost column, consists of a localised patch of flow structures aligned against the mean shear, unwrapping via the Orr mechanism (see Orr Reference Orr1907) into an array of streamwise aligned structures with a distinct oblique component, shown in figure 3(f). These oblique structures then transfer energy into streamwise-independent streaks through the oblique wave mechanism in which structures with wavenumbers $(k_{x},k_{z})=(0,\pm a)$ interact nonlinearly to move energy into the wavenumber $(k_{x},k_{z})=(0,0)$ . These streaks are then able to ‘self-sustain’ for a long period of the flow’s evolution, as shown in the leftmost column of figure 4, by utilising the lift-up mechanism (described by Landahl Reference Landahl1980) and their own instability to offset viscous decay, (see Waleffe Reference Waleffe1997), before eventually being of large enough amplitude to transition to a high-energy oblique structure, as shown in figure 5(a,f), which is visited only transiently, leading to a break down to small-scale turbulence due to an instability reminiscent of the Kelvin–Helmholtz instability. The energetics of these sequential growth mechanisms in unstratified minimal seed trajectories were discussed by Pringle, Willis & Kerswell (Reference Pringle, Willis and Kerswell2012) and Duguet et al. (Reference Duguet, Monokrousos, Brandt and Henningson2013).

This transition mechanism has been interpreted by Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) as an initial concentration of energy in a small region of the flow, allowing a minimisation of the total energy input whilst maximising the local energy input. From this point, the most efficient way to transition to turbulence is by exploiting the lift-up mechanism on a high-energy flow structure and thus increasing the total energy of the flow to a point beyond which it is sufficiently unstable to become turbulent. As shown in figure 2, there is a sustained period of near-constant energy density $E(t)$ , and during this time the flow is only very slowly changing, indicating that the flow trajectory is near a coherent state. The isosurfaces plotted in figure 4(a,f) show that this state is very reminiscent of an SSP/VWI state. As $\mathit{Ri}_{B}$ is increased to $\mathit{Ri}_{B}=10^{-4}$ and $10^{-3}$ , the minimal seed, its trajectory, and the coherent state to which it approaches, remain very similar to the unstratified case, as can be seen from the flow isosurfaces shown in figures 3(b,g,c,h)–5(b,g,c,h), although there are the beginnings of an oblique structure appearing in the coherent states to which these flows evolve.

As $\mathit{Ri}_{B}$ is increased further, the behaviour changes qualitatively. The minimal seed for $\mathit{Ri}_{B}=3\times 10^{-3}$ (as shown in figure 3 d) still consists of a localised patch of flow structures aligned against the mean shear, which unwrap via the Orr mechanism into the same array of streamwise aligned structures with a distinct oblique component, followed by the oblique wave mechanism. However, the newly created streamwise-independent streaks are visited only transiently, no longer able to be sustained, and the flow evolves into a new, fully three-dimensional coherent state, which is shown in figure 4(d,i). This new stratified coherent state is different from the unstratified SSP/VWI state visited by the unstratified minimal seed trajectory, with the stratified coherent state being clearly more three-dimensional and very reminiscent of the stratified edge states recently reported by Olvera & Kerswell (Reference Olvera and Kerswell2014).

Figure 2 shows that in the cases $\mathit{Ri}_{B}=10^{-4}$ , $10^{-3}$ and $3\times 10^{-3}$ , transition to turbulence is apparently faster than the unstratified case, with the same accuracy in the estimated value of $E_{c}$ , and the same time target time $T=300$ . However, the time taken to transition to turbulence is a function of how close our numerical estimates for the minimal seed initial conditions are to the edge manifold. In the limit of successively better approximations to the ‘ideal’ minimal seed, which lies exactly on the edge manifold, $t_{\mathit{transition}}\rightarrow \infty$ . Therefore, since the value of $E_{c}$ is only bracketed to a certain accuracy, the time to transition for the minimal seeds presented here is most likely to be affected by the actual value of $E_{c}$ relative to the energy bracket for $E_{c}$ found here.

The minimal seed for $\mathit{Ri}_{B}=10^{-2}$ is once again of a qualitatively different character. Although it still consists of a localised patch of flow structures that unwrap via the Orr mechanism, there appears to be no quasi-steady flow structure into which it evolves. The flow is now chaotic with a weak oscillation, with no single structure dominating the flow evolution. Indeed, due to the chaotic nature of the dynamics on the edge manifold that this trajectory follows, the DAL method struggled to identify turbulence solutions without extending the optimisation time interval to $T=400$ . Even after extending the optimisation time interval, the DAL method required up to ten times as many iterations to identify initial conditions that transition to turbulence for the flow with $\mathit{Ri}_{B}=10^{-2}$ compared to the number of iterations required for flows with smaller bulk Richardson numbers, apparently because of the chaotic nature of the edge manifold. As noted in the figure captions, videos of the evolution of the minimal seeds for the different values of $\mathit{Ri}_{B}$ are available as supplementary material.

3.2. ‘Wide’ geometry W

As already noted, the minimal seeds in the narrow geometry N are streamwise localised, but fill much of the spanwise extent of the computational domain. To investigate to some extent the sensitivity of the identified minimal seeds to the flow geometry, we also calculate minimal seeds in geometry W, which is twice as wide in the spanwise direction. Figure 6 shows the time evolution from the minimal seed initial conditions of the total energy density $E(t)$ as defined in (3.1) for the wide geometry W. Streamwise velocity $u=\pm 0.6\max (u)$ isosurfaces are plotted in figures 7–9 at times $t=0$ and 25, $t=70$ and 150, and $t=210$ and 280 respectively, and videos of the flow evolution are available as supplementary material.

Figure 6. Time variation of energy density $E(t)$ as defined in (3.1) for the minimal seed trajectories in geometry W for $\mathit{Ri}_{B}=0$ , black; $3\times 10^{-3}$ , green; and $10^{-2}$ , purple.

Figure 7. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=0$ (a–c) and $t=25$ (d–f) for the minimal seed trajectories in geometry W for $\mathit{Ri}_{B}=0$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

Figure 8. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=70$ (a–c) and $t=150$  (d–f) for the minimal seed trajectories in geometry W for $\mathit{Ri}_{B}=0$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

Figure 9. Isosurfaces of streamwise perturbation velocity $u=\pm 0.6\max (u)$ at $t=210$ (a–c) and $t=280$ (d–f) for the minimal seed trajectories in geometry W for $\mathit{Ri}_{B}=0$ , $3\times 10^{-3}$ and $10^{-2}$ , from left to right. Videos of the flow evolution are available as supplementary material.

The $\mathit{Ri}_{B}=0$ unstratified minimal seed trajectory in this wider geometry shares the same characteristic evolution as the minimal seed trajectory in the narrower geometry, but with the addition of spanwise localisation, that is a streamwise and spanwise localised patch of flow structures aligned against the mean shear which unwrap via the Orr mechanism into a streamwise aligned structure with a distinct oblique component, as shown in figure 7(a,d), before utilising the oblique wave mechanism to produce a spanwise isolated pair of streaks, reminiscent of the non-localised structure seen in geometry N. As is apparent in the tabulated values of $E_{c}$ listed in table 1, these structures are slightly less energetic than the minimal seeds identified in the narrow geometry N, in that, as already noted, the critical values of the energy density $E_{c}$ (i.e. the energy divided by the volume of the computational domain) in geometry W are approximately 40 % of the equivalent values determined in geometry N.

These streaks survive in the flow for an extended period of time, as shown in figure 8(a,d), and are another realisation of an SSP/VWI coherent structure. These streaks eventually transition to a high-energy spanwise isolated oblique structure, as shown in figure 9(a,d), which is visited only transiently, before breaking down to small-scale turbulence. This sequence of events, as well as the spanwise localisation, are both consistent with those reported by Monokrousos et al. (Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011) and verified using a different objective functional by Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) in the domain $4{\rm\pi}\times 2\times 2{\rm\pi}$ at the larger Reynolds number $\mathit{Re}=1500$ .

The minimal seed trajectory for $\mathit{Ri}_{B}=3\times 10^{-3}$ (shown in figures 7 b,e–9 b,e) again shares the same characteristic evolution as the equivalent minimal seed trajectory in the narrow flow geometry N. The initial condition consists of the same spanwise and streamwise localised patch of flow structures aligned against the mean shear as in the unstratified flow, which again unwrap via the Orr mechanism into a streamwise aligned structure with a distinct oblique component, as shown in the figure 7(e). The oblique wave mechanism transfers energy into a pair of spanwise isolated streaks that are visited only transiently, and the flow evolves onto a long-lived spanwise localised three-dimensional coherent state which has the same oblique characteristics as the one found in the narrow geometry N. This structure is eventually no longer able to be maintained, and transition to turbulence occurs.

Furthermore, the minimal seed trajectory for $\mathit{Ri}_{B}=10^{-2}$ in geometry W is also a spanwise localised version of the equivalent trajectory in the narrower geometry N, as can be seen by comparison of figures 3(e,j)–5(e,j) and 7(c,f)–9(c,f). The initial condition in geometry W consists of a streamwise and spanwise localised patch of flow structures that unwrap via the Orr mechanism into a flow which quickly becomes chaotic. The flow continues in much the same way as its small domain version, until eventually transitioning to turbulence. Once again, it appeared computationally more difficult to identify the minimal seed for the flow with $Ri_{B}=10^{-2}$ than for the flows with smaller bulk Richardson numbers, requiring two or three times as many iterations, although this convergence was still faster than for the equivalent flow with the same $Ri_{B}=10^{-2}$ in the narrower geometry N.

4. Rolls, streaks and waves

To analyse the effect that stratification has on the SSP/VWI states, we decompose the full perturbation velocity and density fields into roll, streak and wave components. First, we define $\boldsymbol{{\mathcal{U}}}=\langle \boldsymbol{u}\rangle _{x}$ and ${\it\Theta}=\langle {\it\rho}\rangle _{x}$ , where $\langle \boldsymbol{a}\rangle _{x}=(1/L_{x})\int _{0}^{L_{x}}\boldsymbol{a}\,\text{d}x$ , and decompose $\boldsymbol{u}=\boldsymbol{{\mathcal{U}}}+\hat{\boldsymbol{u}}=\boldsymbol{{\mathcal{U}}}_{r}+\boldsymbol{{\mathcal{U}}}_{s}+\hat{\boldsymbol{u}}$ and ${\it\rho}={\it\Theta}+\hat{{\it\rho}}$ , so that

(4.1) $$\begin{eqnarray}(\boldsymbol{u},{\it\rho})=(\boldsymbol{{\mathcal{U}}}_{r}+\boldsymbol{{\mathcal{U}}}_{s}+\hat{\boldsymbol{u}},{\it\rho})=(0,{\mathcal{V}},{\mathcal{W}},0)_{\mathit{roll}}+({\mathcal{U}},0,0,{\it\Theta})_{\mathit{streak}}+(\hat{u} ,\hat{v},{\hat{w}},\hat{{\it\rho}})_{wave},\end{eqnarray}$$

where the $r$ subscript denotes the streamwise-independent wall-normal and spanwise roll velocity and the $s$ subscript denotes the streamwise-independent streamwise streak velocity. As originally argued independently by Hall & Smith (Reference Hall and Smith1991) and Waleffe (Reference Waleffe1997), at high Reynolds number, if there are rolls in the flow of typical amplitude $O({\it\epsilon})$ with ${\it\epsilon}\ll 1$ , then their decay rate due to viscosity is $O(\mathit{Re}^{-1})$ . During the $O(\mathit{Re})$ time in which they survive in the flow, they can advect streamwise velocity through the $O(1)$ shear of PCF a distance $O({\it\epsilon}\mathit{Re})$ and so produce $O({\it\epsilon}\mathit{Re})$ streaks in the flow. If the amplitude of these streaks is sufficiently large, they can undergo an instability which creates a wave field. The nonlinear self-interaction of this wave field then puts energy back into the rolls, and provided that this input of energy is sufficient to balance the viscous decay of the rolls, a ‘self-sustaining process’ associated with this ‘vortex–wave interaction’ is possible. Supposing that the streaks need to be $O(1)$ to become unstable requires ${\it\epsilon}=\mathit{Re}^{-1}$ and so a quadratic self-interaction of the wave field of $O(\mathit{Re}^{-1})$ is sufficient. Hall & Sherwin (Reference Hall and Sherwin2010) showed that there is a numerically exact solution in unstratified PCF for asymptotically large $\mathit{Re}$ such that the rolls have ${\mathcal{V}},{\mathcal{W}}=O(\mathit{Re}^{-1})$ and the streaks have ${\mathcal{U}}=O(1)$ throughout the domain. The waves inject energy into the rolls primarily in an $O(\mathit{Re}^{-1/3})$ critical layer where the background PCF plus the streak flow has the same velocity as the wave velocity, within which $(\hat{u} ,\hat{v},{\hat{w}})=O(\mathit{Re}^{-5/6},\mathit{Re}^{-7/6},\mathit{Re}^{-5/6})$ , giving an integrated Reynolds stress contribution over the critical layer of $O(\mathit{Re}^{-1})$ . The waves essentially provide an $O(\mathit{Re}^{-1})$ jump across the critical layer in the component of the roll velocity normal to the critical layer through their Reynolds stresses. Such scalings have been used successfully to provide initial guesses for a search for numerically exact solutions of the Navier–Stokes equations (see Waleffe Reference Waleffe2001; Wedin & Kerswell Reference Wedin and Kerswell2004).

To see how this scaling is affected by the addition of a stable stratification, we examine the energetics for the total energy density of the rolls and of the streaks. We obtain

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}K_{r}(t)=\frac{1}{2}\frac{\text{d}}{\text{d}t}\langle |\boldsymbol{{\mathcal{U}}}_{r}|^{2}\rangle _{y,z}=\langle -\mathit{Ri}_{B}{\mathcal{V}}{\it\Theta}-\boldsymbol{{\mathcal{U}}}_{r}\boldsymbol{\cdot }\langle \hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{u}}\rangle _{x}-\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}\boldsymbol{{\mathcal{U}}}_{r}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{{\mathcal{U}}}_{r}\rangle _{y,z}, & \displaystyle\end{eqnarray}$$

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}K_{s}(t)=\frac{1}{2}\frac{\text{d}}{\text{d}t}\langle |\boldsymbol{{\mathcal{U}}}_{s}|^{2}\rangle _{y,z}=\langle -{\mathcal{V}}{\mathcal{U}}-\boldsymbol{{\mathcal{U}}}_{s}\boldsymbol{\cdot }\langle \hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{\boldsymbol{u}}\rangle _{x}-\mathit{Re}^{-1}\boldsymbol{{\rm\nabla}}\boldsymbol{{\mathcal{U}}}_{s}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{{\mathcal{U}}}_{s}\rangle _{y,z}, & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}t}P_{s}(t)=\frac{1}{2}\frac{\text{d}}{\text{d}t}\langle \mathit{Ri}_{B}{\it\Theta}^{2}\rangle _{y,z}=\langle \mathit{Ri}_{B}{\mathcal{V}}{\it\Theta}-\mathit{Ri}_{B}{\it\Theta}\langle \hat{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{{\it\rho}}\rangle _{x}-\mathit{Ri}_{B}(\mathit{Re}\mathit{Pr})^{-1}\boldsymbol{{\rm\nabla}}{\it\Theta}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\Theta}\rangle _{y,z}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\langle \boldsymbol{a}\rangle _{y,z}=(1/2L_{z})\int _{-1}^{1}\int _{0}^{L_{z}}\boldsymbol{a}\,\text{d}z\,\text{d}y$ , $K_{r}(t)$ is the roll kinetic energy density, $K_{s}(t)$ is the streak kinetic energy density, and $P_{s}(t)$ is the streak potential energy density. Note that $K_{r}(t)+K_{s}(t)+K_{w}(t)=K(t)$ , where $K_{w}(t)$ is the wave kinetic energy density

(4.5) $$\begin{eqnarray}K_{w}(t)={\textstyle \frac{1}{2}}\langle |\hat{\boldsymbol{u}}|^{2}\rangle .\end{eqnarray}$$

Figure 10. Time dependence of normalised streak kinetic energy density $K_{s}/K$ as defined in (4.3) (plotted with a blue line) and normalised wave kinetic energy density $K_{w}/K$ as defined in (4.5) (plotted with a red line) for the minimal seed trajectories in geometry N for $\mathit{Ri}_{B}=0$ (a), $10^{-4}$ (b), $10^{-3}$ (c), $3\times 10^{-3}$ (d), $10^{-2}$ (e).

Taking ${\it\epsilon}=\mathit{Re}^{-1}$ , the correct scaling for the viscous decay time is obtained from (4.2). We also see from (4.3) that when ${\mathcal{V}}=O(\mathit{Re}^{-1})$ is sustained over a period of $O(\mathit{Re})$ , we obtain streaks ${\mathcal{U}}=O(1)$ . We can also obtain the scaling of Hall & Sherwin (Reference Hall and Sherwin2010) for the wave velocities, and their Reynolds stress contribution, under the assumption that they act over an $O(\mathit{Re}^{-1/3})$ critical layer, and balance the viscous dissipation there.

For stratified PCF with $\mathit{Ri}_{B}\neq 0$ there is a new buoyancy flux term $-\mathit{Ri}_{B}{\mathcal{V}}{\it\Theta}$ entering the energetics of the rolls. Comparing the first term in the right-hand side of (4.3) to that of (4.4), which are production terms for the streak kinetic energy density and streak potential energy density respectively, we see that ${\it\Theta}=O(1)$ everywhere for a passive scalar field placed in a self-sustained process. Thus, (4.2) shows that a buoyancy flux associated with the streak flow effective across the whole domain provides a contribution to the kinetic energy density $K_{r}$ of the rolls of $O(\mathit{Ri}_{B}\mathit{Re}^{-1})$ , and is able to disrupt fully the flux of wave Reynolds stress input of $O(\mathit{Re}^{-3})$ over the $O(\mathit{Re}^{-1/3})$ critical layer when $\mathit{Ri}_{B}=O(\mathit{Re}^{-2})$ . This simple minded high Reynolds number scaling argument, which is domain size independent, is not inconsistent with the moderate Reynolds number minimal seed calculations presented above for which the unstratified coherent state visited by the unstratified minimal seed trajectory is no longer a viable solution for sufficiently large $\mathit{Ri}_{B}\sim 3\times 10^{-3}$ , and the new coherent states have a modified roll structure that is no longer streamwise-independent.

Figures 10 and 11 show the time evolution of the streak kinetic energy density $K_{s}(t)$ , as defined in (4.3) and plotted with a blue line, and the wave kinetic energy density $K_{w}(t)$ , as defined in (4.5) and plotted with a red line, normalised by the total kinetic energy $K(t)$ for each of the minimal seed trajectories in geometries N and W respectively. It is clear that for $\mathit{Ri}_{B}=0$ , $10^{-4}$ and $10^{-3}$ , there is a balance (shown by the approximate plateaux in the streak and wave energy components) between these energies (and also the residual roll kinetic density density $K_{r}(t)$ , which is not shown as it is, unsurprisingly, appreciably smaller in magnitude) for a large period of the flow evolution. For $\mathit{Ri}_{B}=3\times 10^{-3}$ , this balance has been significantly disrupted, and is not maintained purely by the velocity fields. There also is a small amplitude oscillation, and for $\mathit{Ri}_{B}=10^{-2}$ the above-mentioned weakly oscillatory nature of the flow is clear. The unstratified self-sustaining process is completely disrupted for $\mathit{Ri}_{B}\gtrsim 3\times 10^{-3}$ , as expected.

Figure 11. Time dependence of normalised streak kinetic energy density $K_{s}/K$ as defined in (4.3) (plotted with a blue line) and normalised wave kinetic energy density $K_{w}/K$ as defined in (4.5) (plotted with a red line) for the minimal seed trajectories in geometry W for $\mathit{Ri}_{B}=0$ (a), $3\times 10^{-3}$ (b) and $10^{-2}$ (c).