# Supplemental Material: Model for transitional turbulence in a planar shear
flow

Santiago J. Benavides  Dwight Barkley  Mathematics Institute, University of
Warwick, Coventry CV4 7AL, United Kingdom

###### Contents

  1. I Full Model Specification
  2. II Linear Stability Equations
  3. III Details of Model Derivation
    1. III.1 Waleffe flow
    2. III.2 Reynolds Average
    3. III.3 Truncation of Vertical Modes
    4. III.4 Turbulence Closures
      1. III.4.1 Reynolds Stress
      2. III.4.2 Production
      3. III.4.3 Transport and Dissipation
      4. III.4.4 $q_{1}$ and $\alpha$
    5. III.5 Parameter fits
  4. IV Details of Stability Analysis
    1. IV.1 Uniform steady states
    2. IV.2 Linearized equations
    3. IV.3 Normal-mode decomposition
    4. IV.4 Elimination of pressure and full eigenvalue problem
  5. V Angle Selection
    1. V.1 Long-wavelength simplification
    2. V.2 Critical angle at onset of instability
  6. VI Description of Supplemental Movies

##  I Full Model Specification

The model is given by the following equations for the six fields
$u_{0},u_{1},v_{1},w_{0},w_{1},q_{0}$, each depending on streamwise coordinate
$x$, spanwise coordinate $z$, and time $t$,

$\displaystyle\frac{\partial u_{0}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}u_{0}+\frac{1}{2}u_{1}\partial_{x}u_{1}+\frac{1}{2}\beta
v_{1}u_{1}+w_{0}\partial_{z}u_{0}+\frac{1}{2}w_{1}\partial_{z}u_{1}$ (1)
$\displaystyle=-\partial_{x}p_{0}+\nabla^{2}_{H}u_{0}/Re-\alpha_{0}u_{0},$
$\displaystyle\frac{\partial u_{1}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}u_{1}+u_{1}\partial_{x}u_{0}+w_{0}\partial_{z}u_{1}+w_{1}\partial_{z}u_{0}$
$\displaystyle=-\partial_{x}p_{1}+\nabla^{2}_{H}u_{1}/Re-(\alpha+\beta^{2}/Re)(u_{1}-1)-\beta
A(q_{0}),$ $\displaystyle\frac{\partial v_{1}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}v_{1}+w_{0}\partial_{z}v_{1}$
$\displaystyle=-\beta
p_{1}+\partial_{x}A(q_{0})+\nabla^{2}_{H}v_{1}/Re-\beta^{2}v_{1}/Re,$
$\displaystyle\frac{\partial w_{0}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}w_{0}+\frac{1}{2}u_{1}\partial_{x}w_{1}+\frac{1}{2}\beta
v_{1}w_{1}+w_{0}\partial_{z}w_{0}+\frac{1}{2}w_{1}\partial_{z}w_{1}$
$\displaystyle=-\partial_{z}p_{0}+\nabla^{2}_{H}w_{0}/Re-\alpha_{0}w_{0},$
$\displaystyle\frac{\partial w_{1}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}w_{1}+u_{1}\partial_{x}w_{0}+w_{0}\partial_{z}w_{1}+w_{1}\partial_{z}w_{0}$
$\displaystyle=-\partial_{z}p_{1}+\nabla^{2}_{H}w_{1}/Re-(\alpha+\beta^{2}/Re)w_{1},$
$\displaystyle\frac{\partial q_{0}}{\partial t}$
$\displaystyle+u_{0}\partial_{x}q_{0}+\frac{1}{2}u_{1}\partial_{x}q_{1}(u_{1},w_{1},q_{0})+\frac{1}{2}\beta
v_{1}q_{1}(u_{1},w_{1},q_{0})+w_{0}\partial_{z}q_{0}+\frac{1}{2}w_{1}\partial_{z}q_{1}(u_{1},w_{1},q_{0})$
$\displaystyle=\frac{1}{2}A(q_{0})[\beta u_{1}+\partial_{x}v_{1}]-2\alpha
q_{0}-\varepsilon_{0}(q_{0};Re)+\frac{1}{Re}\nabla^{2}_{H}q_{0}+\nabla_{H}\cdot\left(\nu_{T}(q_{0};Re)\nabla_{H}q_{0}\right),$

where $\beta=\pi/2$, and $\nabla_{H}=(\partial_{x},\partial_{z})$. The
pressures $p_{0}$ and $p_{1}$ ensure incompressibility for each set of large-
scale modes,

$\partial_{x}u_{0}+\partial_{z}w_{0}=0,\quad\partial_{x}u_{1}-\beta
v_{1}+\partial_{z}w_{1}=0.$ (2)

The equations are closed with the following,

$\displaystyle A(q_{0})$ $\displaystyle=a\,((q_{0}^{2}+\eta^{2})^{1/2}-\eta),$
$\displaystyle\varepsilon_{0}(q_{0};Re)$ $\displaystyle=\frac{c}{Re}q_{0},$
(3) $\displaystyle q_{1}(u_{1},w_{1},q_{0})$
$\displaystyle=-\frac{u_{1}\partial_{x}q_{0}+w_{1}\partial_{z}q_{0}}{2\alpha+2\kappa+\beta^{2}/Re},$
$\displaystyle\nu_{T}(q_{0};Re)$ $\displaystyle=d\,Re\,q_{0}.$

Model parameters have the values

$a=0.3,\,c=9.65,\,d=0.05,\,\eta=5\times
10^{-3},\,\alpha=0.01,\,\alpha_{0}=0.03,\,\kappa=0.045$ (4)

Note that in the main paper, we lighten notation by dropping the $0$ subscript
from $\varepsilon_{0}$, and write simply $\varepsilon$. Since the model fields
depend only on the horizontal spatial coordinates $(x,z)$, the notation
$\nabla_{H}$ is redundant in stating the model. We use it to clearly
distinguish $\nabla_{H}$ from the three-dimensional gradient
$\nabla=(\partial_{x},\partial_{y},\partial_{z})$ appearing in the derivation
below.

##  II Linear Stability Equations

The uniform turbulent state in the model is a steady state with no $(x,z)$
dependence. These are of the form
$(u_{0},u_{1},v_{1},w_{0},w_{1},q_{0})=(0,u_{ss},0,0,0,q_{ss})$, where
$u_{ss}$ and $q_{ss}$ satisfy the $u_{1}$ and $q_{0}$ equations after dropping
$x,z$, and $t$ derivatives:

$\displaystyle 0$ $\displaystyle=(\alpha+\beta^{2}/Re)(1-u_{ss})-\beta
A(q_{ss}),$ (5a) $\displaystyle 0$ $\displaystyle=\frac{1}{2}A(q_{ss})\beta
u_{ss}-2\alpha q_{ss}-\varepsilon_{0}(q_{ss};Re).$ (5b)

We denote uniform steady states by $(u_{ss},q_{ss})$.

Straightforward linearization about the uniform turbulent state gives the
following eigenvalue problem for the temporal growth rate $\sigma$ of a linear
mode $\hat{u}_{0}\exp(ik_{x}x+ik_{z}z)$, and similarly for
$\hat{u}_{1},\hat{v}_{1},\hat{w}_{0},\hat{w}_{1},\hat{q}_{0}$.

$\displaystyle\sigma\hat{u}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{u}_{1}+\frac{1}{2}\beta
u_{ss}\hat{v}_{1}$ $\displaystyle=$ $\displaystyle\frac{k_{x}^{2}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}+\left(\hat{D}_{2}-\alpha_{0}\right)\hat{u}_{0}$
(6a) $\displaystyle\sigma\hat{u}_{1}+ik_{x}u_{ss}\hat{u}_{0}$ $\displaystyle=$
$\displaystyle{-\frac{ik_{x}\alpha\beta}{k^{2}+\beta^{2}}\,\hat{v}_{1}}+\left(\hat{D}_{3}-\alpha\right)\hat{u}_{1}+\left({\frac{2k_{x}^{2}}{k^{2}+\beta^{2}}}-1\right)\beta\delta
A_{ss}\,\hat{q}_{0}$ (6b) $\displaystyle\sigma\hat{v}_{1}$ $\displaystyle=$
$\displaystyle\left(\hat{D}_{3}-\frac{\alpha\beta^{2}}{k^{2}+\beta^{2}}\right)\hat{v}_{1}+\left(1-\frac{2\beta^{2}}{k^{2}+\beta^{2}}\right)ik_{x}\delta
A_{ss}\hat{q}_{0}$ (6c)
$\displaystyle\sigma\hat{w}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{w}_{1}$
$\displaystyle=$ $\displaystyle\frac{k_{x}k_{z}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}+\left(\hat{D}_{2}-\alpha_{0}\right)\hat{w}_{0}$
(6d) $\displaystyle\sigma\hat{w}_{1}+u_{ss}ik_{x}\hat{w}_{0}$ $\displaystyle=$
$\displaystyle{-\frac{ik_{z}\alpha\beta}{k^{2}+\beta^{2}}\,\hat{v}_{1}}+\left(\hat{D}_{3}-\alpha\right)\hat{w}_{1}+\frac{2k_{x}k_{z}\beta\delta
A_{ss}}{k^{2}+\beta^{2}}\,\hat{q}_{0}$ (6e) $\displaystyle\sigma\hat{q}_{0}$
$\displaystyle=$
$\displaystyle\frac{1}{2}A_{ss}\left(\beta\hat{u}_{1}{+ik_{x}\hat{v}_{1}}\right)+\hat{L}\hat{q}_{0}$
(6f)

where

$k^{2}=k_{x}^{2}+k_{z}^{2},\quad\hat{D}_{2}=-\frac{k^{2}}{Re},\quad\hat{D}_{3}=-\frac{k^{2}+\beta^{2}}{Re},$
(7) $A_{ss}\equiv A(u_{ss},q_{ss}),\quad\delta A_{ss}\equiv\frac{\partial
A}{\partial
q_{0}}(q_{ss}),\quad\delta\varepsilon_{ss}=\frac{d\varepsilon_{0}}{dq_{0}}(q_{ss},Re),$
(8) $\hat{L}=\frac{1}{2}\beta u_{ss}\delta
A_{ss}-2\alpha-\delta\varepsilon_{ss}-\frac{k^{2}}{Re}-\nu_{T}(q_{ss};Re)k^{2}-\frac{u_{ss}^{2}k_{x}^{2}}{2(2\alpha+2\kappa+\beta^{2}/Re)}.$
(9)

See Sec. IV for further details.

For numerical computation of the growth rate $\sigma$ in the main paper, all
terms except those involving $\sigma$ are put onto the right-hand side. The
eigenvalue problem problem becomes $\sigma{\bf x}={\bf M}{\bf x}$, where ${\bf
x}=(\hat{u}_{0},\hat{u}_{1},\hat{v}_{1},\hat{w}_{0},\hat{w}_{1},\hat{q}_{0})^{T}$
and $\bf M$ is a $6\times 6$ matrix. Given parameters $Re$, $k_{x}$, and
$k_{z}$, and $(u_{ss},q_{ss})$ from a numerical solution of (5), the elements
of $\bf M$ are obtained. The eigenvalues of $\bf M$ are found from a standard
eigenvalue solver. The leading eigenvalue, corresponding to maximum growth
rate, is always found to be real.

##  III Details of Model Derivation

In this section, we will walk through the model derivation. The general
approach for the model creation will be the following: first we will present
the Reynolds-Averaged Navier Stokes (RANS) equations, along with the equation
for the Turbulent Kinetic Energy (TKE). These will constitute our main
dynamical variables. We will then introduce our simplification assumptions and
turbulence closures, all of which are motivated (and validated) by
observations from direct numerical simulations.

###  III.1 Waleffe flow

We start from a planar flow confined between stress-free boundaries and driven
by a sinusoidal body force. In the wall-bounded literature this system is
called Waleffe flow (Wf), after Waleffe who used it in his work on the self-
sustaining process [1]. The flow reproduces the same transitional phenomena as
plane Couette flow (pCf) [2, 3], while the stress-free boundary conditions
significantly simplifying analytical and computational treatment.

We use coordinates $(x,y,z)$ aligned with the streamwise, wall-normal, and
spanwise directions, respectively. The components of the velocity
${\bf\tilde{u}}$ are denoted by $(\tilde{u},\tilde{v},\tilde{w})$. We non-
dimensionalize using a length scale $h$ equal to the half gap between the
stress-free walls, and a velocity scale $U$, equal to the maximum laminar
velocity. In non-dimensional coordinates, the walls are at $y=\pm 1$. We use
$\beta\equiv\pi/2$ as a shorthand for the first wall-normal (vertical)
wavenumber.

The momentum equation for Wf is,

$\frac{\partial{\bf\tilde{u}}}{\partial
t}+{\bf\tilde{u}}\cdot\nabla{\bf\tilde{u}}=-\nabla\tilde{p}+\frac{1}{Re}\nabla^{2}{\bf\tilde{u}}-\alpha{\bf\tilde{u}}_{H}+{\bf
f},$ (10)

where $\tilde{p}$ is the pressure, $Re\equiv Uh/\nu$ is the Reynolds number,
$\nu$ is the kinematic viscosity, ${\bf f}$ is the body force, and $\alpha$ is
the drag coefficient. The drag term only acts on the horizontal velocity
components ${\bf\tilde{u}}_{H}\equiv\tilde{u}{\bf e}_{x}+\tilde{w}{\bf e}_{z}$
and is only a significant source of dissipation relative to viscous effects
for $y$-independent modes and large horizontal scales. It represents the
effect of the no-slip boundary conditions in pCf, which acts as an effective
drag on the $y$-averaged horizontal velocities [4, 3]. The body force is
defined to be,

${\bf f}\equiv\left(\alpha+\frac{\beta^{2}}{Re}\right)\sin\left(\beta
y\right){\bf e}_{x}.$ (11)

This body force results in a (linearly stable) laminar solution of the form
${\bf\tilde{u}}_{lam}=\sin(\beta y){\bf e}_{x}$. The momentum equation is
accompanied by an incompressibility condition

$\nabla\cdot{\bf\tilde{u}}=0,$ (12)

periodic boundary conditions in the horizontal $(x,z)$ directions, and stress-
free boundary conditions in the vertical direction

$\partial_{y}\tilde{u}(x,\pm 1,z)=\tilde{v}(x,\pm
1,z)=\partial_{y}\tilde{w}(x,\pm 1,z)=0,$ (13)

which includes a no-penetration condition on the wall-normal velocity.

Throughout this work, we will use direct numerical simulations (DNS) of Wf as
a way to motivate model choices and calibrate parameter values in our model.
We simulated the Wf equations (with periodic boundary conditions in the
horizontal) using the open source pseudo-spectral code Dedalus [5]. A first-
order Crank-Nicolson Adams-Bashforth time stepping scheme was used with a time
step $\Delta t=0.003$ advective time units. A Fourier-spectral method with
$3/2$ dealiasing is used for the horizontal directions, with a resolution of
eight grid points per space unit. A SinCos basis was used in the vertical
direction with the same spatial resolution. Random initial conditions were
used at higher $Re$ and the banded solution was obtained by incrementally
decreasing $Re$ from the resulting uniform state.

###  III.2 Reynolds Average

We apply a Reynolds average to the momentum equation by letting
${\bf\tilde{u}}={\bf u}+{\bf u^{\prime}}$, where ${\bf
u}\equiv\langle{\bf\tilde{u}}\rangle$ is the average of the velocity field and
$\langle{\bf u^{\prime}}\rangle=0$. We do not specify the averaging operation,
although presumably it is either a temporal average or spatial coarse graining
such that scales smaller than $h$, as well as the turbulent space and time
scales, are averaged over, leading to slowly evolving, non-fluctuating large
scale quantities. For this reason, we will refer to ${\bf u}$ as the ‘large-
scale flow’. We also define the Turbulent Kinetic Energy (TKE)
$q\equiv\langle{\bf u^{\prime}}\cdot{\bf u^{\prime}}\rangle/2$. Following
standard procedures [6], we arrive at the equations of motion for the large-
scale flow,

$\frac{\partial{\bf u}}{\partial t}+{\bf u}\cdot\nabla{\bf u}=-\nabla
p+\frac{1}{Re}\nabla^{2}{\bf u}-\alpha{\bf u}_{H}+{\bf
f}+\nabla\cdot{\bf\mathcal{R}},$ (14a) and TKE, $\frac{\partial q}{\partial
t}+{\bf u}\cdot\nabla
q+\nabla\cdot\textbf{T}=\mathcal{P}-\varepsilon+\frac{1}{Re}\nabla^{2}q-2\alpha
q,$ (14b)

where ${\bf u}$ satisfies the same incompresibility constraint and boundary
conditions as ${\bf\tilde{u}}$, equations (12) and (13). The TKE also has
vanishing gradients at the two walls.

A small approximation is already made in equation (14b). The exact drag term
from averaging is $\alpha\langle(u^{\prime})^{2}+(w^{\prime})^{2}\rangle$, and
we have approximated it as $2\alpha
q=\alpha\langle(u^{\prime})^{2}+(v^{\prime})^{2}+(w^{\prime})^{2}\rangle$.
Since
$\langle(v^{\prime})^{2}\rangle\ll\langle(u^{\prime})^{2}+(w^{\prime})^{2}\rangle$,
this approximation has negligible effect.

The averaging has introduced four new fluctuation correlation terms that will
eventually require closure choices. These are the Reynolds stress,

$\mathcal{R}_{ij}=-\langle u^{\prime}_{i}u^{\prime}_{j}\rangle,$ (15)

the production term,

$\mathcal{P}=-\langle u^{\prime}_{i}u^{\prime}_{j}\rangle\partial_{j}u_{i},$
(16)

the pseudo-dissipation rate,

$\varepsilon=\frac{1}{Re}\left|\nabla{\bf u^{\prime}}\right|^{2},$ (17)

and the transport term,

$T_{i}=\frac{1}{2}\langle
u^{\prime}_{i}u^{\prime}_{j}u^{\prime}_{j}\rangle+\frac{1}{2}\langle
u^{\prime}_{i}p^{\prime}\rangle.$ (18)

We use the Einstein summation convention, with $u_{i}$ representing an
indexing of $u,v,w$ and similarly for $u^{\prime}_{i}$. Note that closing the
Reynolds stress term also closes the production term.

###  III.3 Truncation of Vertical Modes

We turn to modeling and continue to use variables ${\bf u}$ and $q$, rather
than introducing separate notation for the model. Our first step is to
introduce a Galerkin truncation in the vertical direction. Given the
sinusoidal body force ${\bf f}$, corresponding sinusoidal laminar flow, and
the stress-free boundary conditions, we use a Fourier basis and project onto
the first two vertical modes (that is, the $y$-independent and $\sin(\beta y)$
or $\cos(\beta y)$ modes). The model large-scale flow and TKE are represented
by

$\displaystyle u(x,y,z,t)$ $\displaystyle=$ $\displaystyle
u_{0}(x,z,t)+u_{1}(x,z,t)\sin(\beta y),$ (19a) $\displaystyle v(x,y,z,t)$
$\displaystyle=$ $\displaystyle v_{1}(x,z,t)\cos(\beta y),$ (19b)
$\displaystyle w(x,y,z,t)$ $\displaystyle=$ $\displaystyle
w_{0}(x,z,t)+w_{1}(x,z,t)\sin(\beta y),$ (19c) $\displaystyle q(x,y,z,t)$
$\displaystyle=$ $\displaystyle q_{0}(x,z,t)+q_{1}(x,z,t)\sin(\beta y),$ (19d)

The no-penetration condition $v(x,\pm 1,z)=0$ implies $v_{0}=0$ and that the
lowest $v$-mode goes as cosine. Incompressibility of the large-scale flow,
$\nabla\cdot{\bf u}=0$, becomes two statements

$\partial_{x}u_{0}+\partial_{z}w_{0}=0,\quad\partial_{x}u_{1}-\beta
v_{1}+\partial_{z}w_{1}=0.$

The pressure field $p$ in (14a) that enforces incompressibility is represented
in the form $p(x,y,z,t)=p_{0}(x,z,t)+p_{1}(x,z,t)\sin(\beta y)$. The fields
$p_{0}$ and $p_{1}$ separately ensure incompressibility for each set of large-
scale modes. Although every fluctuation correlation term has an equivalent
decomposition in vertical modes, $\varepsilon$ is the only one that will
require the consideration of both vertical modes (see the next section below).
Because of this we explicitly define
$\varepsilon=\varepsilon_{0}(x,z,t)+\varepsilon_{1}(x,z,t)\sin(\beta y)$, but
will not expand the remaining fluctuation correlation terms in vertical modes.
Instead, we will write their $y$-dependence directly.

While $v_{1}$ is small, it is essential for the flow dynamics and the retained
modes are the minimal set describing a three-dimensional vector field. For the
TKE field, $q_{0}$ represents the $y$-averaged turbulence and $q_{1}$ allows
for different levels of turbulence in the upper and lower halves of the
geometry. This is important in capturing the so-called ‘overhangs’ at the
edges of turbulent bands where turbulence is not vertically uniform due to the
mean shear. Finally, DNS of Wf confirms that the modes in (14) are the most
energetic modes present in turbulent bands. It is therefore reasonable to base
our model on this minimal set of modes. Substitution of (19) into (14),
followed by a Galerkin projection, will give evolution equations for the
fields, but first we need to consider the turbulence closures.

###  III.4 Turbulence Closures

We will now present the closures. These will be based on fully-resolved
simulations and four-mode simulations of Wf from previous work [2], as well as
simulations performed by the authors for the current study.

####  III.4.1 Reynolds Stress

We start with the Reynolds stress tensor. Only two terms play a significant
role in the dynamics, namely $-\langle u^{\prime}u^{\prime}\rangle$ and
$-\langle u^{\prime}v^{\prime}\rangle$. Fully-resolved simulations reveal only
one dominant vertical mode in each of these Reynolds stress tensor components
– the zeroth vertical mode (invariant in $y$) for $-\langle
u^{\prime}u^{\prime}\rangle$ and the first vertical mode for $-\langle
u^{\prime}v^{\prime}\rangle$. The latter fact can be seen clearly in Figure 4
of Chantry et al. (2016) [2]. It is assumed that the amplitudes of these terms
depend only on the TKE and have the form:

$\displaystyle\mathcal{R}_{12}=\mathcal{R}_{21}=-\langle
u^{\prime}v^{\prime}\rangle$ $\displaystyle=$ $\displaystyle A(q)\cos(\beta
y),$ (20) $\displaystyle\mathcal{R}_{11}=-\langle u^{\prime}u^{\prime}\rangle$
$\displaystyle=$ $\displaystyle-B(q).$ (21)

This results in two Reynolds-stress body forces – one in the streamwise
direction $x$, and one in the wall-normal direction $y$:

$\nabla\cdot{\bf\mathcal{R}}=\left(-\partial_{x}B(q)-A(q)\beta\sin(\beta
y)\right){\bf e}_{x}+\left(\cos(\beta y)\partial_{x}A(q)\right){\bf e}_{y}.$
(22)

We have chosen the signs so that $A(q),B(q)>0$.

We now neglect the $-\langle u^{\prime}u^{\prime}\rangle$ stress, and take
$B(q)=0$ for the remainder. The reason we have included it in the discussion
to this point is because this stress accounts for the largest effect that we
have measured in DNS that we do not include in the present model. This stress
couples the $u_{0}$ field to $q_{0}$ and $q_{1}$. Its associated Reynolds
stress force acts to reduce the amplitude of $u_{0}$, while its associated
production term increase $q_{0}$ (its effect is small compared to other
production terms). We have briefly investigated the role of this stress in the
model and have found that including it does not alter the model dynamics in a
significant way. Therefore, for simplicity we neglect it in the present model.
However, since it is not a negligible term in the $u_{0}$ equation, this has
contributed to the need for an increased drag coefficient $\alpha_{0}$ (see
III.4.4). Future studies may wish to include it. (We find that taking $B(q)$
linear in $q$ is a reasonable approximation.)

Since $A(q)$ is the amplitude of an even mode in $y$, it can only depend on
the even part of the TKE field, that is $q_{0}$ in our model expansion. The
Reynolds-stress closure is specified by a function $A(q_{0})$. We use the
following,

$A(q_{0})=a((q_{0}^{2}+\eta^{2})^{1/2}-\eta),$ (23)

where $a$ and $\eta$ are parameters with $\eta$ small. For $q_{0}\gg\eta$,
$A(q_{0})\approx aq_{0}$, meaning that $\mathcal{R}_{12}$ is simply
proportional to the TKE. From DNS, we find that $A(q_{0})$ scales with a
slightly larger exponent, (approximately $q_{0}^{1.6}$; Supplemental Fig.
1(b)) and does depend somewhat on $Re$, but for simplicity we choose not to
take these factors into account in the simple model. For small $q_{0}$, that
is, $q_{0}\ll\eta$, we have $A(q_{0})\approx aq_{0}^{2}/(2\eta)$, and
$A(q_{0})$ goes to zero quadratically in $q_{0}$. This serves to cut off
turbulent production at small $q_{0}$ (see section III.4.3). Supplemental Fig.
2(c) compares DNS-measured values of the Reynolds stress force across a band
with our closure, as well as with the resulting force across the band from
model simulations (Fig. 2 (a) t = 25350 in the main text).

Supplemental Fig. 1:  Functional dependence of (a) the vertically averaged
pseudo-dissipation rate $\varepsilon_{0}$ on $q_{0}$, and (b) the first mode
amplitude of $\langle u^{\prime}v^{\prime}\rangle$ on $q_{0}$. Taken from a
steady state band found in a DNS of Wf at $Re=140$ (averaged for 1000 time
units).  Supplemental Fig. 2: Comparison of turbulent band profiles and
closures between a DNS of Wf at $Re=140$ (averaged for 1000 time units) and
our model at $Re=75$ (steady state). Across-band profiles of vertical mode
amplitudes for the DNS, (a), and for the model (b). Panel (c) shows the
Reynold stress force for the $u_{1}$ mode equation, measured directly from the
DNS at $Re=140$ (blue), its closure (dashed, orange), and the corresponding
Reynolds stress force for the model (green) at $Re=75$. Panel (d) shows a
similar plot but for the pseudo-dissipation rate.

####  III.4.2 Production

The production $\mathcal{P}$ is just the contraction of the Reynolds stress
tensor $\mathcal{R}$ and large-scale velocity gradient tenor $\nabla{\bf u}$,
followed by projection back to the $q$ equation. Given our simplification of
$\mathcal{R}$, the only terms are
$\mathcal{R}_{12}\partial_{y}u+\mathcal{R}_{21}\partial_{x}v$. Using our
expansions and projecting back to the retained modes gives

$\mathcal{P}=\frac{1}{2}A(q_{0})[\beta u_{1}+\partial_{x}v_{1}].$

It is constant in $y$ and only couples to $q_{0}$.

####  III.4.3 Transport and Dissipation

There are two closures left in the TKE equation, the turbulent pseudo-
dissipation rate and the turbulent transport term. We first focus on the
pseudo-dissipation, where we look at both vertical modes separately, and then
consider the turbulent transport term.

As with $A$, the pseudo-dissipation rate $\varepsilon_{0}$ must only depend on
$q_{0}$, and we further assume that is linear in $q_{0}$. DNS shows a slightly
larger power, closer to $q_{0}^{1.2}$ (Supplemental Fig. 1 (a)), but model
tests confirm that a linear approximation is adequate (Supplemental Fig.2
(d)). DNS results also inform us on the Reynolds number dependence of the
coefficient (found by matching $(u_{ss},q_{ss})$ with DNS as best as possible
for multiple values of $Re$). Our closure for $\varepsilon_{0}$ is,

$\varepsilon_{0}(q_{0};Re)=\frac{c}{Re}q_{0},$ (24)

where $c$ is a parameter. Note that in the main paper, we lighten notation by
dropping the $0$ subscript from $\varepsilon_{0}$, and write simply
$\varepsilon(q_{0};Re)=cq_{0}/Re$.

Having specified the dissipation, we can comment on the functional form of
$A(q_{0})$. The dependence of $A(q_{0})$ on $q_{0}$ gives a quadratic cutoff
in production as $q_{0}\to 0$. Since $\varepsilon_{0}$ is linear in $q_{0}$,
production will necessarily fall below dissipation at small $q_{0}$. This
ensures stability of laminar flow to small perturbations in $q_{0}$.
Physically this rules out “infinitesimal turbulence”, accounting for the fact
that very small fluctuations (deviations from laminar flow), are not turbulent
but rather are small laminar disturbances that we know do not sustain
nonlinearly.

A similarly simple closure exists for $\varepsilon_{1}(q)$, that is
$\varepsilon_{1}(q_{1})=2\kappa q_{1}$, where $\kappa$ does not depend on
$Re$. Although we have measured $\kappa$ from the pseudo-dissipation rate of a
DNS, and not a transport term, this closure can be seen intuitively as the
vertical (turbulent) diffusion of turbulence, modeled as a ‘drag’ term with
drag coefficient $\kappa$ in a two-layer simplification of our system, such
that $\varepsilon_{top}=-\kappa(q_{top}-q_{bottom})$ and
$\varepsilon_{bottom}=-\kappa(q_{bottom}-q_{top})$ [7].

Moving on to the transport term, DNS results and model tests suggest that only
the zeroth vertical mode of the transport term plays a significant role in the
turbulence budget. Apart from this simplification in vertical mode number, we
invoke the gradient-diffusion hypothesis [6, Section 10.3] which states that,

${\bf T}(q_{0};Re)=-\nu_{T}(q_{0};Re)\nabla_{H}q_{0}.$ (25)

We have again used the fact that, since ${\bf T}$ does not depend on $y$, it
must only be a function of $q_{0}$. The turbulent diffusivity is defined in a
standard way, following Pope [6, Section 10.3], as $\nu_{T}\propto
q^{2}_{0}/\varepsilon_{0}$. Since we have already assumed that
$\varepsilon_{0}=cq_{0}/Re$, we can then substitute this into the turbulent
diffusivity arriving at,

$\nu_{T}(q_{0};Re)=d\,Re\,q_{0}.$ (26)

####  III.4.4 $q_{1}$ and $\alpha$

We make two final model approximations. The first concerns $q_{1}$. Projecting
the TKE dynamics (14b) into the $\sin(\beta y)$ mode gives the evolution
equation for $q_{1}$. Ignoring streamwise and spanwise derivatives of $q_{1}$
(which we have verified with DNS), the $q_{1}$ equation simplifies to

$\frac{\partial q_{1}}{\partial
t}+u_{1}\partial_{x}q_{0}+w_{1}\partial_{z}q_{0}=-(2\alpha+\beta^{2}/Re)q_{1}-2\kappa
q_{1}.$ (27)

In an effort to achieve the simplest possible model, we assume that the
timescales of $q_{1}$ are short, resulting in an quasi-static condition for
$q_{1}$. Setting $\partial_{t}q_{1}=0$, we can express $q_{1}$ as a function
of ${\bf u}_{H}$ and $q_{0}$,

$q_{1}(u_{1},w_{1},q_{0})=-\frac{u_{1}\partial_{x}q_{0}+w_{1}\partial_{z}q_{0}}{2\alpha+2\kappa+\beta^{2}/Re}.$
(28)

Thus, $q_{1}$ is not treated as a dynamical field in our model.

The last model approximation we make is to use a larger value of $\alpha$ in
the $u_{0},w_{0}$ equations. In the process of neglecting $\langle
u^{\prime}u^{\prime}\rangle$, which plays a significant role in the budget of
$u_{0}$, we found it necessary to increase the drag for the zero modes of the
large-scale flow. Furthermore, we found that, without this increase, the
large-scale flow has a greater tendency to become unstable and chaotic, making
the banded state unattainable. We therefore introduce a final model parameter
$\alpha_{0}$ to use in the $u_{0},w_{0}$ equations.

###  III.5 Parameter fits

The numerical value of the parameters were first calibrated using measurements
of the Reynolds stresses, turbulent dissipation, and transport term in a time
and along-band average of a single turbulent band found in a Wf simulation run
at $Re=140$ (Supplemental Fig. 2(a)). Due to the nature of our simplifications
and assumptions in our closure choices, we had to further refine these
parameter values in order to ensure simulation results which better matched
those of the DNS (Supplemental Fig. 2(b)).

The first step in choosing parameter values was measuring the Reynolds stress
amplitudes, turbulent dissipation, and turbulent diffusion in a DNS, giving us
a first estimate for what the coefficients in our closures should be. From
there, we tweaked $a$ and $c$ so that the uniform turbulent fixed point in the
homogeneous dynamics (subsection IV.1) would appear at $Re=72$, a value which
we found results in turbulent bands at about $Re=68$. Despite being able to
arbitrarily choose the Reynolds number corresponding to the appearance of
bands in our model, the one we have settled on is about half of that for Wf
($Re\approx 130$). Our attempts at having realistic Wf Reynolds numbers
resulted in the formation of sharp gradients in the large-scale flow which
caused numerical instability, as well as the nonlinear evolution being more
prone to instability and chaos, resulting in a lack of banded states. Although
our DNS budgets do not suggest that a turbulent viscosity plays a significant
role, the lower $Re$ values chosen could be interpreted in such a way.

After calibrating $a$ and $c$, we performed a series of simulations of the
model in a simplified geometry, following that of the ‘tilted’ domain
simulations from Barkley & Tuckerman [8]. That is, simulations were perfomred
in a domain tilted at a fixed angle (24 degrees in our case) to the streamwise
direction, extended in the across-band direction but assumed to be invariant
in the along-band direction, making it a one-dimensional calculation for the
model. We performed a series of simulations to observe the formation of bands,
and to measure the large scale flow and TKE amplitudes, as well as the width
of the bands. From these 1D runs, we modified the rest of the parameters, in
particular the turbulent diffusivity coefficient $d$, so as to best match the
band width found in the DNS.

The final check on the parameter choices was to run a 2D simulation of the
full equations (1) and observe the behavior. Throughout our experimentation
with parameter values, we found that the presence of a linear instability of
the spatially uniform steady state, the formation and stability of bands in
the 1D model, and the angle of the initial bands in the 2D model all to be
extremely robust. However, the two-dimensional dynamics of the full model
proved to be more sensitive to parameter values. In particular, we found
various instabilities of the turbulent bands in the along-band direction, as
well as instances when the large scale flow’s amplitude was too large,
resulting in chaotic motion. This restricted our parameter values, ultimately
resulting in our introduction of $\alpha_{0}$ as well as a forced mismatch
between observed Reynolds number transitions in the DNS and in our model.

##  IV Details of Stability Analysis

In this section we analyze the stability of the the uniform turbulent state in
the model system. The content is standard, straightforward linear stability
analysis, but since there are many steps in the calculations, for completeness
we present a summary of key steps.

###  IV.1 Uniform steady states

The uniform turbulent state in the model is a steady state, or fixed point, of
the equations with no $(x,z)$ dependence. Dropping all horizontal derivatives
from the model equations, it is immediately evident that any uniform steady
state must have $u_{0}=w_{0}=v_{1}=w_{1}=0$, and hence uniform steady states
can have at most $u_{1}$ and $q_{0}$ different from zero. We denote uniform
steady states by $(u_{ss},q_{ss})$, where $u_{ss}$ and $q_{ss}$ satisfy the
$u_{1}$ and $q_{0}$ equations after dropping $x,z$, and $t$ derivatives:

$\displaystyle 0$ $\displaystyle=(\alpha+\beta^{2}/Re)(1-u_{ss})-\beta
A(q_{ss}),$ (29a) $\displaystyle 0$ $\displaystyle=\frac{1}{2}A(q_{ss})\beta
u_{ss}-2\alpha q_{ss}-\varepsilon_{0}(q_{ss};Re).$ (29b)

One sees that laminar flow, $(u_{ss},q_{ss})=(1,0)$, is a solution to these
equations for any positive $Re$. We are primarily concerned with the non-
trivial steady states with $(u_{ss},q_{ss})\neq(1,0)$ corresponding to uniform
turbulence. For $Re>72.4$, a pair of non-trivial steady states exist, with
$u_{ss}<1$ and $q_{ss}>0$, corresponding to uniform turbulence. One of these
states (the ‘upper state’) is stable to spatially uniform perturbations.

###  IV.2 Linearized equations

Linear stability of $(u_{ss},q_{ss})$ is determined by first linearizing the
model equations about this state. Using primes to denote perturbation fields,
we let

$\displaystyle u_{0}=u^{\prime}_{0},w_{0}=w^{\prime}_{0},$ (30) $\displaystyle
u_{1}=u_{ss}+u^{\prime}_{1},v_{1}=v^{\prime}_{1},w_{1}=w^{\prime}_{1}$ (31)
$\displaystyle q_{0}=q_{ss}+q^{\prime}_{0},$ (32)

and substitute these into the model equations and then drop terms quadratic in
the prime variables. Initially, it will be convenient to consider the
linearlized equation in three groups: the vertical mode-0 large-scale-flow
fields, the vertical mode-1 large-scale-flow fields, and the turbulent field.
We define two operators corresponding to the viscous diffusion terms

$\displaystyle D_{2}$ $\displaystyle\equiv\frac{1}{Re}\nabla^{2}_{H}$ (33)
$\displaystyle D_{3}$
$\displaystyle\equiv\frac{1}{Re}\left(\nabla^{2}_{H}-\beta^{2}\right)$ (34)

The vertical mode-0 linear equations are

$\displaystyle\frac{\partial u^{\prime}_{0}}{\partial
t}+\frac{1}{2}u_{ss}\frac{\partial u^{\prime}_{1}}{\partial
x}+\frac{1}{2}\beta u_{ss}v^{\prime}_{1}$ $\displaystyle=$
$\displaystyle-\frac{\partial p^{\prime}_{0}}{\partial
x}+D_{2}u^{\prime}_{0}-\alpha_{0}u^{\prime}_{0}$ (35)
$\displaystyle\frac{\partial w^{\prime}_{0}}{\partial
t}+\frac{1}{2}u_{ss}\frac{\partial w^{\prime}_{1}}{\partial x}$
$\displaystyle=$ $\displaystyle-\frac{\partial p^{\prime}_{0}}{\partial
z}+D_{2}w^{\prime}_{0}-\alpha_{0}w^{\prime}_{0}$ (36) $\frac{\partial
u^{\prime}_{0}}{\partial x}+\frac{\partial w^{\prime}_{0}}{\partial z}=0$ (37)

where $p^{\prime}_{0}$ is the perturbation pressure needed to maintain
incompressibility in the perturbation fields.

The vertical mode-1 linear equations are

$\displaystyle\frac{\partial u^{\prime}_{1}}{\partial t}+u_{ss}\frac{\partial
u^{\prime}_{0}}{\partial x}$ $\displaystyle=$ $\displaystyle-\frac{\partial
p^{\prime}_{1}}{\partial x}+D_{3}u^{\prime}_{1}-\alpha
u^{\prime}_{1}-\beta\delta A_{ss}q^{\prime}_{0}$ (38)
$\displaystyle\frac{\partial v^{\prime}_{1}}{\partial t}$ $\displaystyle=$
$\displaystyle-\beta p^{\prime}_{1}+D_{3}v^{\prime}_{1}+\delta
A_{ss}\frac{\partial q^{\prime}_{0}}{\partial x}$ (39)
$\displaystyle\frac{\partial w^{\prime}_{1}}{\partial t}+u_{ss}\frac{\partial
w^{\prime}_{0}}{\partial x}$ $\displaystyle=$ $\displaystyle-\frac{\partial
p^{\prime}_{1}}{\partial z}+D_{3}w^{\prime}_{1}-\alpha w^{\prime}_{1}$ (40)
$\frac{\partial u^{\prime}_{1}}{\partial x}-\beta
v^{\prime}_{1}+\frac{\partial w^{\prime}_{1}}{\partial z}=0$ (41)

where $p^{\prime}_{1}$ is the perturbation pressure needed to maintain
incompressibility and

$\delta A_{ss}\equiv\frac{\partial A}{\partial q_{0}}(q_{ss})$ (42)

The linearized turbulent field equation is

$\frac{\partial q^{\prime}_{0}}{\partial t}=\frac{1}{2}A_{ss}\left(\beta
u^{\prime}_{1}+\frac{\partial v^{\prime}_{1}}{\partial
x}\right)+Lq^{\prime}_{0}$ (43)

where

$A_{ss}\equiv A(u_{ss},q_{ss})$ (44)

and $L$ is the operator

$\displaystyle L$ $\displaystyle\equiv\frac{1}{2}\beta u_{ss}\delta
A_{ss}-2\alpha-\delta\varepsilon_{ss}$
$\displaystyle+\frac{1}{Re}\nabla^{2}_{H}+\nu_{T}(q_{ss};Re)\nabla^{2}_{H}+\frac{u_{ss}^{2}}{2(2\alpha+2\kappa+\beta^{2}/Re)}\frac{\partial^{2}}{\partial
x^{2}}$ (45)

with

$\delta\varepsilon_{ss}=\frac{d\varepsilon_{0}}{dq_{0}}(q_{ss},Re).$

We will discuss the important operator $L$ below.

###  IV.3 Normal-mode decomposition

We next employ a spectral normal-mode decomposition to transform the
differential operators into algebraic form. Let

$u^{\prime}_{0}=\hat{u}_{0}e^{i(k_{x}x+k_{z}z)}e^{\sigma t}$ (46)

and similarly for the other fields, where $k_{x}$ and $k_{z}$ are the
wavenumbers of the normal modes and $\sigma$ is a temporal eigenvalue to be
determined. In principle $\sigma$ may be complex. Upon substitution we obtain
the following.

The vertical mode-0 linear equations become

$\displaystyle\sigma\hat{u}_{0}+\frac{1}{2}u_{ss}ik_{x}\hat{u}_{1}+\frac{1}{2}\beta
u_{ss}\hat{v}_{1}$ $\displaystyle=$ $\displaystyle-
ik_{x}\hat{p}_{0}+\hat{D}_{2}\hat{u}_{0}-\alpha_{0}\hat{u}_{0}$ (47)
$\displaystyle\sigma\hat{w}_{0}+\frac{1}{2}u_{ss}ik_{x}\hat{w}_{1}$
$\displaystyle=$ $\displaystyle-
ik_{z}\hat{p}_{0}+\hat{D}_{2}\hat{w}_{0}-\alpha_{0}\hat{w}_{0}$ (48)
$ik_{x}\hat{u}_{0}+ik_{z}\hat{w}_{0}=0$ (49)

where

$\hat{D}_{2}=-\frac{k_{x}^{2}+k_{z}^{2}}{Re}=-\frac{k^{2}}{Re}$ (50)

The vertical mode-1 linear equations become

$\displaystyle\sigma\hat{u}_{1}+u_{ss}ik_{x}\hat{u}_{0}$ $\displaystyle=$
$\displaystyle-
ik_{x}\hat{p}_{1}+\hat{D}_{3}\hat{u}_{1}-\alpha\hat{u}_{1}-\beta\delta
A_{ss}\hat{q}_{0}$ (51) $\displaystyle\sigma\hat{v}_{1}$ $\displaystyle=$
$\displaystyle-\beta\hat{p}_{1}+\hat{D}_{3}\hat{v}_{1}+\delta
A_{ss}ik_{x}\hat{q}_{0}$ (52)
$\displaystyle\sigma\hat{w}_{1}+u_{ss}ik_{x}\hat{w}_{0}$ $\displaystyle=$
$\displaystyle-ik_{z}\hat{p}_{1}+\hat{D}_{3}\hat{w}_{1}-\alpha\hat{w}_{1}$
(53) $ik_{x}\hat{u}_{1}-\beta\hat{v}_{1}+ik_{z}\hat{w}_{1}=0$ (54)

where

$\hat{D}_{3}=-\frac{k_{x}^{2}+k_{z}^{2}+\beta^{2}}{Re}=-\frac{k^{2}+\beta^{2}}{Re}$
(55)

The turbulent-field equation becomes

$\sigma\hat{q}_{0}=\frac{1}{2}A_{ss}\left(\beta\hat{u}_{1}+ik_{x}\hat{v}_{1}\right)+\hat{L}\hat{q}_{0}$
(56)

where $\hat{L}$ is obtained from (45) by replacing
$(\partial_{x},\partial_{z})\rightarrow(ik_{x},ik_{z})$.

###  IV.4 Elimination of pressure and full eigenvalue problem

We use incompressibility to obtain expressions for the pressure fields
$\hat{p}_{0}$ and $\hat{p}_{1}$. Multiplying (47) by $ik_{x}$ and (48) by
$ik_{z}$, adding these, and using the incompressibility constraints (49) and
(54), one obtains the following expression for $\hat{p}_{0}$

$\hat{p}_{0}=\frac{ik_{x}\beta u_{ss}}{k^{2}}\,\hat{v}_{1}$ (57)

Similarly, multiplying (51), (52), and (53), by $ik_{x}$, $-\beta$, and
$ik_{z}$, respectively and using incompressibility, one obtains

$\hat{p}_{1}=\frac{2ik_{x}\beta\delta
A_{ss}}{k^{2}+\beta^{2}}\,\hat{q}_{0}+\frac{\alpha\beta}{k^{2}+\beta^{2}}\,\hat{v}_{1}$
(58)

Substituting these expressions for pressure, we obtain the following full
eigenvalue problem for the temporal growth rate $\sigma$

$\displaystyle\sigma\hat{u}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{u}_{1}+\frac{1}{2}\beta
u_{ss}\hat{v}_{1}$ $\displaystyle=$ $\displaystyle\frac{k_{x}^{2}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}+\left(\hat{D}_{2}-\alpha_{0}\right)\hat{u}_{0}$
(59a) $\displaystyle\sigma\hat{u}_{1}+ik_{x}u_{ss}\hat{u}_{0}$
$\displaystyle=$
$\displaystyle{\color[rgb]{1,0,0}-\frac{ik_{x}\alpha\beta}{k^{2}+\beta^{2}}\,\hat{v}_{1}}+\left(\hat{D}_{3}-\alpha\right)\hat{u}_{1}+\left({\color[rgb]{1,0,0}\frac{2k_{x}^{2}}{k^{2}+\beta^{2}}}-1\right)\beta\delta
A_{ss}\,\hat{q}_{0}$ (59b) $\displaystyle\sigma\hat{v}_{1}$ $\displaystyle=$
$\displaystyle\left(\hat{D}_{3}-\frac{\alpha\beta^{2}}{k^{2}+\beta^{2}}\right)\hat{v}_{1}+\left(1-\frac{2\beta^{2}}{k^{2}+\beta^{2}}\right)ik_{x}\delta
A_{ss}\hat{q}_{0}$ (59c)
$\displaystyle\sigma\hat{w}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{w}_{1}$
$\displaystyle=$ $\displaystyle\frac{k_{x}k_{z}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}+\left(\hat{D}_{2}-\alpha_{0}\right)\hat{w}_{0}$
(59d) $\displaystyle\sigma\hat{w}_{1}+u_{ss}ik_{x}\hat{w}_{0}$
$\displaystyle=$
$\displaystyle{\color[rgb]{1,0,0}-\frac{ik_{z}\alpha\beta}{k^{2}+\beta^{2}}\,\hat{v}_{1}}+\left(\hat{D}_{3}-\alpha\right)\hat{w}_{1}{\color[rgb]{1,0,0}+\frac{2k_{x}k_{z}\beta\delta
A_{ss}}{k^{2}+\beta^{2}}\,\hat{q}_{0}}$ (59e) $\displaystyle\sigma\hat{q}_{0}$
$\displaystyle=$
$\displaystyle\frac{1}{2}A_{ss}\left(\beta\hat{u}_{1}{\color[rgb]{1,0,0}+ik_{x}\hat{v}_{1}}\right)+\hat{L}\hat{q}_{0}$
(59f)

We have written these equations such as to indicate the origins of terms.
Those are the left are associated with advection, while those on the left are
associated with pressure gradients, dissipative effects and turbulent
production. The five terms shown in red can be neglected for long-wavelength
(small-$k$) modes as discussed below.

##  V Angle Selection

In this section we consider the selected angle of turbulent bands at onset of
instability of uniform turbulence as $Re$ is decreased. We first consider a
simplification of the linear stability equations for modes with long-
wavelengths relative to the wall spacing. We then show explicitly that in this
long-wavelength case that the selected angle of turbulent bands at onset, the
critical angle, is necessarily less than $45^{\circ}$.

###  V.1 Long-wavelength simplification

Turbulent bands have long wavelengths relative to the wall spacing, which by
choice of length scales is $2h=2$. Considering long-wavelength modes allows
use to simplify the linear stability equations. Letting the wavenumber be
$k=\epsilon\ll 1$, we will drop terms in the eigenvalue problem (59) that are
$O(\epsilon^{2})$.

While $\alpha$ is a small parameter, we assume that $\alpha$ is no smaller
than $\epsilon$. Thus, we make the approximations

$\displaystyle k^{2}+\beta^{2}$ $\displaystyle\rightarrow\beta^{2},$ (60)
$\displaystyle\hat{D}_{2}-\alpha_{0}=-\frac{k^{2}}{Re}-\alpha_{0}$
$\displaystyle\rightarrow-\alpha_{0}$ (61)
$\displaystyle\hat{D}_{3}-\alpha=-\frac{k^{2}+\beta^{2}}{Re}-\alpha$
$\displaystyle\rightarrow-(\alpha+\beta^{2}/{Re})$ (62)

In addition, incompressibility of the mode-1 fields (54) gives

$\frac{ik_{x}}{\beta}\hat{u}_{1}+\frac{ik_{z}}{\beta}\hat{w}_{1}=\hat{v}_{1}.$
(63)

Assuming eigenvector components $\hat{u}_{1}$ and $\hat{w}_{1}$ are $O(1)$, we
have $\hat{v}_{1}=O(\epsilon)$. Thus, terms involving $k_{x}^{2}$,
$k_{x}k_{z}$, $k_{x}\hat{v}_{1}$, and $k_{z}\hat{v}_{1}$ are
$O(\epsilon^{2})$. These term are shown in red in (59a)-(59f). Neglecting
these terms, the eigenvalue equations become:

$\displaystyle\sigma\hat{u}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{u}_{1}+\frac{1}{2}\beta
u_{ss}\hat{v}_{1}$ $\displaystyle=$ $\displaystyle\frac{k_{x}^{2}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}-\alpha_{0}\hat{u}_{0}$ (64a)
$\displaystyle\sigma\hat{u}_{1}+ik_{x}u_{ss}\hat{u}_{0}$ $\displaystyle=$
$\displaystyle-(\alpha+\beta^{2}/{Re})\hat{u}_{1}-\beta\delta
A_{ss}\hat{q}_{0}$ (64b) $\displaystyle\sigma\hat{v}_{1}$ $\displaystyle=$
$\displaystyle-(\alpha+\beta^{2}/{Re})\hat{v}_{1}-ik_{x}\delta
A_{ss}\hat{q}_{0}$ (64c)
$\displaystyle\sigma\hat{w}_{0}+\frac{1}{2}ik_{x}u_{ss}\hat{w}_{1}$
$\displaystyle=$ $\displaystyle\frac{k_{x}k_{z}\beta
u_{ss}}{k^{2}}\,\hat{v}_{1}-\alpha_{0}\hat{w}_{0}$ (64d)
$\displaystyle\sigma\hat{w}_{1}+ik_{x}u_{ss}\hat{w}_{0}$ $\displaystyle=$
$\displaystyle-(\alpha+\beta^{2}/{Re})\hat{w}_{1}$ (64e)
$\displaystyle\sigma\hat{q}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}A_{ss}\beta\hat{u}_{1}+\hat{L}\hat{q}_{0}$ (64f)

This relatively simple eigenvalue problem captures well the onset of
instability of the uniform turbulent fixed points $(u_{ss},q_{ss})$. Table 1
shows a comparison of eigenvalues for the full and simplified systems for
relevant parameters.

full |  $-0.05972$ |  $-0.03587\pm 0.04809i$ |  $-0.02906\pm 0.07571i$ |  0.00069
---|---|---|---|---
simplified | $-0.05889$ | $-0.03505\pm 0.04809i$ | $-0.02841\pm 0.07597i$ | 0.00104
Table 1: Eigenvalues for the full eigenvalue problem (59) and simplified
eigenvalue problem (64) for $Re=82$, $k_{x}=0.1$, $k_{z}=0.24$.

###  V.2 Critical angle at onset of instability

The analysis can continue by focusing on the case of marginal stability, where
$\sigma=0$. No further simplifying assumptions are made and all calculations
below follow exactly from eigenvalue problem (64). Setting $\sigma=0$,
equations (64) become algebraic equations

$\displaystyle 0$ $\displaystyle=\frac{1}{2}ik_{x}u_{ss}\hat{u}_{1}+\beta
u_{ss}\left(\frac{1}{2}-\frac{k_{x}^{2}}{k^{2}}\right)\hat{v}_{1}+\alpha_{0}\hat{u}_{0}$
(65a) $\displaystyle 0$
$\displaystyle=ik_{x}u_{ss}\hat{u}_{0}+(\alpha+\beta^{2}/{Re})\hat{u}_{1}+\beta\delta
A_{ss}\hat{q}_{0}$ (65b) $\displaystyle 0$
$\displaystyle=(\alpha+\beta^{2}/{Re})\hat{v}_{1}+ik_{x}\delta
A_{ss}\hat{q}_{0}$ (65c) $\displaystyle 0$ $\displaystyle=\frac{1}{2}\beta
A_{ss}\hat{u}_{1}+\hat{L}\hat{q}_{0}$ (65d)

We only need to consider these four equations, and not six, because
$\hat{w}_{0}$ and $\hat{w}_{1}$ do not appear in these equations. We are not
interested here in the critical eigenvector, but only in the parameter
relations implied by a zero eigenvalue. This information is contained in (65).

Following straightforward elimination and back substitution, (65) reduces to:

$\displaystyle\left\\{\frac{k_{x}^{2}u_{ss}^{2}}{\left(\alpha+\beta^{2}/Re\right)}\left(\frac{1}{2}-\frac{k_{x}^{2}}{k^{2}}\right)+\left(\frac{1}{2}\frac{k_{x}^{2}u_{ss}^{2}}{\left(\alpha+\beta^{2}/Re\right)}+\alpha_{0}\right)\frac{2\left(\alpha+\beta^{2}/Re\right)\hat{L}}{\beta^{2}A_{ss}\delta
A_{ss}}-\alpha_{0}\right\\}\hat{q}_{0}=0.\quad$

The size of eigenvector component $q_{0}$ is arbitrary, hence the expression
in brackets must vanish whenever there is a zero eigenvalue.

Expressions can be made more compact by introducing the following functions of
$Re$

$c_{1}(Re)=\frac{u_{ss}^{2}}{\left(\alpha+\beta^{2}/Re\right)}>0,\quad
c_{2}(Re)=\frac{2\left(\alpha+\beta^{2}/Re\right)}{\beta^{2}A_{ss}\delta
A_{ss}}>0$ (67)

Note that since $(u_{ss},q_{ss})$ depend on $Re$, so do $A_{ss}$ and $\delta
A_{ss}$. We emphasize that both are positive functions.

Rather than parameterizing eigenmodes modes by wavenumber pairs
$(k_{x},k_{z})$, it is more useful to parameterize modes by $(k,\theta)$,
where

$\tan\theta=\frac{k_{x}}{k_{z}},\qquad k^{2}=k_{x}^{2}+k_{z}^{2}$ (68)

The striped pattern with wavenumbers $(k_{x},k_{z})$ makes angle $\theta$ with
the streamwise direction. Using all the above, the condition for a zero
eigenvalue becomes

$\displaystyle
c_{1}(Re)k^{2}\sin^{2}\theta\left(\frac{1}{2}-\sin^{2}\theta\right)+\left(\frac{c_{1}(Re)}{2}k^{2}\sin^{2}\theta+\alpha_{0}\right)c_{2}(Re)\hat{L}(k,\theta,Re)-\alpha_{0}=0$
(69)

We have explicitly indicated here that $\hat{L}$ is a function of
$(k,\theta,Re)$. It is obtained from (45) by replacing
$(\partial_{x},\partial_{z})\rightarrow(ik_{x},ik_{z})$ to give

$\hat{L}(k,\theta,Re)=\frac{1}{2}\beta u_{ss}\delta
A_{ss}-2\alpha-\delta\varepsilon_{ss}-k^{2}\left(\frac{1}{Re}+\nu_{T}(q_{ss};Re)+c_{3}(Re)\sin^{2}\theta\right),$
(70)

where

$c_{3}(Re)=\frac{u_{ss}^{2}}{2(2\alpha+2\kappa+\beta^{2}/Re)}>0,$ (71)

is a third positive function of $Re$.

Letting $F_{\sigma=0}(k,\theta,Re)$ denote the left-hand side of (69), the
condition for a zero eigenvalue, $\sigma=0$, is written

$\displaystyle F_{\sigma=0}(k,\theta,Re)=0.$ (72)

We are interested in the critical condition where a zero eigenvalue first
appears with decreasing $Re$. This will demarcate the onset of turbulent
bands. At such a point we much have

$\displaystyle F_{\sigma=0}(k_{c},\theta_{c},Re_{c})=0,$ (73a)
$\displaystyle\left.\frac{\partial}{\partial
k}F_{\sigma=0}\right|_{(k_{c},\theta_{c},Re_{c})}=0,$ (73b)
$\displaystyle\left.\frac{\partial}{\partial\theta}F_{\sigma=0}\right|_{(k_{c},\theta_{c},Re_{c})}=0,$
(73c)

where $Re_{c}$ is the critical Reynolds number at which instability first
occurs, and $(k_{c},\theta_{c}$) are the associated critical wavenumber and
critical angle.

(To understand why (73c) must hold, fix $k=k_{c}$ and suppose for that
$\partial_{\theta}F_{\sigma=0}\neq 0$ at the critical point
$(k_{c},\theta_{c},Re_{c})$. Then by the implicit function theorem, there is a
function $\theta=g(Re)$ such that $F_{\sigma=0}(k_{c},g(Re),Re)=0$ for $Re$ on
an open interval containing $Re_{c}$. This implies that there is a zero
eigenvalue both above and below $Re_{c}$. The contradicts the assumption that
$Re_{c}$ is the first value with decreasing $Re$ for which there is a zero
eigenvalue. A similar argument applies fixing $\theta$ and varying $k$.)

In principle, given the model closure parameters, equations (73) could be
evaluated and used to find the critical point for the onset of bands. We will
not pursue this here. However, we consider just the condition on the angle
(73c). Straightforward algebra gives

$\sin^{2}\theta_{c}=\frac{1}{4+k^{2}c_{3}(Re)}\left(1+c_{2}(Re_{c})\hat{L}(k_{c},\theta_{c},Re_{c})-k^{2}\frac{\alpha_{0}c_{3}(Re)}{c_{1}(Re)}\right)$
(74)

To get a bound on the critical angle $\theta_{c}$, let $\hat{L}_{0}$ denote
the terms in $\hat{L}$ that are independent of wavenumber:

$\hat{L}_{0}(Re)=\frac{1}{2}\beta u_{ss}\delta
A_{ss}-2\alpha-\delta\varepsilon_{ss}$ (75)

All the other term in $\hat{L}$ contain a factor $-k^{2}$. Hence,
$\hat{L}<\hat{L}_{0}$. (In fact, the only positive term in $\hat{L}$ is
$\frac{1}{2}\beta u_{ss}\delta A_{ss}$.) Thus,

$\sin^{2}\theta_{c}\leq\frac{1}{4}\left(1+c_{2}(Re_{c})\hat{L}_{0}(Re_{c})\right)$
(76)

Consider now the system’s homogeneous dynamics – the dynamics of spatially
uniform states. In the absence of spatial dependence, we let
$u_{1}(x,z,t)\rightarrow u(t)$ and $q_{0}(x,z,t)\rightarrow q(t)$. The
$(u,q)$-dynamics are governed by the ODEs

$\displaystyle\dot{u}$
$\displaystyle=\left(\alpha+\beta^{2}/Re\right)(1-u)-\beta A(q),$ (77a)
$\displaystyle\dot{q}$ $\displaystyle=\frac{1}{2}\beta uA(q)-2\alpha
q-\varepsilon_{0}(q).$ (77b)

Setting $\dot{u}_{1}=\dot{q}_{0}=0$ gives the steady state equations (29),
whose non-laminar solution is the uniform turbulent fixed point
$(u_{ss},q_{ss})$. The stability of this state to uniform perturbation is
determined by the linearization of (77). For $(u_{ss},q_{ss})$ to be a
linearly stable steady state of (77), the determinant of the $2\times 2$
stability matrix must be positive. This gives

$\displaystyle-\left(\alpha+\beta^{2}/Re\right)\left(\frac{1}{2}\beta
u_{ss}\delta
A_{ss}-2\alpha-\delta\varepsilon_{ss}\right)+\frac{1}{2}\beta^{2}A_{ss}\delta
A_{ss}$ $\displaystyle>0$ (78)

or

$\displaystyle\frac{2\left(\alpha+\beta^{2}/Re\right)\hat{L}_{0}}{\beta^{2}A_{ss}\delta
A_{ss}}=c_{2}(Re)\hat{L}_{0}(Re)<1$ (79)

Using this bound in (76) gives

$\sin^{2}\theta_{c}<\frac{1}{2}$ (80)

Hence, stability of the uniform turbulent state to uniform perturbations
implies that the critical angle for turbulent bands at onset is less than
$45^{\circ}$. While this is a statement about the onset of bands in the model,
the result is independent of the particular form of the closures used and
hence is a very general result. It is possible to improve the bound (80) by
exploiting the specific closures used in the model, but we leave this for a
future publication.

The above result might be the start for an explanation as to why nearly
periodic turbulent-laminar patterns form in planar shear flow, while periodic
arrays of puffs are not observed to spontaneously form in pipe flow [9, 10].
Puffs in pipe flow might be loosely viewed as bands at angle $90^{\circ}$ with
respect to the streamwise direction. The above result rules out linear
instability of turbulence to periodic patterns at $90^{\circ}$. This requires
further study.

##  VI Description of Supplemental Movies

  1. Movie S1

`Lx224Lz100_Re75_localized_initial.mp4`. Model run at $Re=75$ from the
localized initial condition (IC), with $(L_{x},L_{z})=(224,100)$. Demonstrates
growth from tips of band. Corresponds to panel (a) in Figure 2 of the main
text.

  2. Movie S2

`Lx224Lz100_Re75_localized.mp4`. Model run at $Re=75$ from localized IC.
Demonstrates evolution to steady band. Corresponds to panel (a) in Figure 2 of
the main text.

  3. Movie S3

`Lx224Lz100_Re78_localized_initial.mp4`. Model run at $Re=78$ from localized
IC, with $(L_{x},L_{z})=(224,100)$. Demonstrates splitting mechanism.
Corresponds to panel (b) in Figure 2 of the main text.

  4. Movie S4

`Lx224Lz100_Re78_localized.mp4`. Model run at $Re=78$ from localized IC.
Demonstrates evolution to steady bands. Corresponds to panel (b) in Figure 2
of the main text.

  5. Movie S5

`Lx224Lz100_Re78_uniform.mp4`. Model run at $Re=78$ from the unfirom turbulent
fixed point IC ($(u_{ss},q_{ss})$ at $Re=78$), with $(L_{x},L_{z})=(224,100)$.
Demonstrates breakdown of uniform turbulence at $Re<85.1$. Corresponds to
panel (c) in Figure 2 of the main text.

  6. Movie S6

`Lx224Lz100_Re80_localized.mp4`. Model run at $Re=80$ from localized IC, with
$(L_{x},L_{z})=(224,100)$. Demonstrates complex steady state “criss-cross”
pattern.

  7. Movie S7

`Lx224Lz100_Re88_localized.mp4`. Model run at $Re=88$ from localized IC, with
$(L_{x},L_{z})=(224,100)$. Demonstrates evolution to uniform turbulence beyond
$Re=85.1$.

  8. Movie S8

`Lx224Lz200_Re75_localized_gauss.mp4`. Model run at $Re=75$ from gaussian IC,
with $(L_{x},L_{z})=(224,200)$. Demonstrates the robustness of the band angle,
even with different domain aspect ratios.

  9. Movie S9-S11

`Lx448Lz200_Re{75,78,80}_localized.mp4`. Models run at $Re=75,78,80$ from
localized IC, with $(L_{x},L_{z})=(448,200)$. The larger domain runs
demonstrate the complexity of the transients (including different domains with
opposite band angles), as well as the variety of final banded states that can
occur. `Lx448Lz200_Re80_localized.mp4` corresponds to Figure 2(e) from the
main text.

## References

  * Waleffe [1997] F. Waleffe, On a self-sustaining process in shear flows, Phys. Fluids 9, 883 (1997).
  * Chantry  _et al._ [2016] M. Chantry, L. S. Tuckerman, and D. Barkley, Turbulent–laminar patterns in shear flows without walls, J. Fluid Mech. 791, R8 (2016).
  * Tuckerman  _et al._ [2020] L. S. Tuckerman, M. Chantry, and D. Barkley, Patterns in wall-bounded shear flows, Annu. Rev. Fluid Mech. 52, 343 (2020).
  * Suri  _et al._ [2014] B. Suri, J. Tithof, J. Mitchell, Radford, R. O. Grigoriev, and M. F. Schatz, Velocity profile in a two-layer Kolmogorov-like flow, Phys. Fluids 26, 053601 (2014).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res. 2, 023068 (2020).
  * Pope [2000] S. B. Pope,  _Turbulent flows_ (Cambridge university press, 2000).
  * Barkley [2011] D. Barkley, Modeling the transition to turbulence in shear flows, Journal of Physics: Conference Series 318, 032001 (2011).
  * Barkley and Tuckerman [2005] D. Barkley and L. S. Tuckerman, Computational study of turbulent-laminar patterns in Couette flow, Phys. Rev. Lett. 94, 014502 (2005).
  * Moxey and Barkley [2010] D. Moxey and D. Barkley, Distinct large-scale turbulent-laminar states in transitional pipe flow, Proceedings of the National Academy of Sciences 107, 8091 (2010).
  * Avila and Hof [2013] M. Avila and B. Hof, Nature of laminar-turbulence intermittency in shear flows, Phys. Rev. E 87, 063012 (2013).