# Near-cancellation of up- and down-gradient momentum transport in forced
magnetized shear-flow turbulence

B. Tripathi1<EMAIL_ADDRESS>A.E. Fraser2 P.W. Terry1 E.G. Zweibel1 M.J.
Pueschel3,4 1University of Wisconsin-Madison, Madison, Wisconsin 53706, U.S.A.
2University of California, Santa Cruz, Santa Cruz, California 95064, U.S.A.
3Dutch Institute for Fundamental Energy Research, 5612 AJ Eindhoven, The
Netherlands
4Eindhoven University of Technology, 5600 MB Eindhoven, The Netherlands

###### Abstract

Visco-resistive magnetohydrodynamic turbulence, driven by a two-dimensional
unstable shear layer that is maintained by an imposed body force, is examined
by decomposing it into dissipationless linear eigenmodes of the initial
profiles. The down-gradient momentum flux, as expected, originates from the
large-scale instability. However, continual up-gradient momentum transport by
large-scale linearly stable but nonlinearly excited eigenmodes is identified,
and found to nearly cancel the down-gradient transport by unstable modes. The
stable modes effectuate this by depleting the large-scale turbulent
fluctuations via energy transfer to the mean flow. This establishes a physical
mechanism underlying the long-known observation that coherent vortices formed
from nonlinear saturation of the instability reduce turbulent transport and
fluctuations, as such vortices are composed of both the stable and unstable
modes, which are nearly equal in their amplitudes. The impact of magnetic
fields on the nonlinearly excited stable modes is then quantified. Even when
imposing a strong magnetic field that almost completely suppresses the
instability, the up-gradient transport by the stable modes is at least two-
thirds of the down-gradient transport by the unstable modes, whereas for
weaker fields, this fraction reaches up to $98\%$. These effects are
persistent with variations in magnetic Prandtl number and forcing strength.
Finally, continuum modes are shown to be energetically less important, but
essential for capturing the magnetic fluctuations and Maxwell stress. A simple
analytical scaling law is derived for their saturated turbulent amplitudes. It
predicts the fall-off rate as the inverse of the Fourier wavenumber, a
property which is confirmed in numerical simulations.

##  I Introduction

Owing to their ubiquity in laboratory,[1] geophysical,[2, 3] and astrophysical
environments,[4, 5, 6, 7, 8, 9] shear layers have been extensively studied.
[10, 11] Observations and analyses from experiments and direct numerical
simulations have offered insights into the connection between large-scale
vortical structures formed from the instability of a shear layer and turbulent
transport across the layer.[12, 13, 14] Properties like shape and scale of the
nonlinearly saturated vortices, which dominate the transport, are generally
attributed to the linearly-unstable eigenmodes or closely related nonlinear
fluctuations. [15, 16] The nonlinear saturation of the instability, however,
can be more complex than just the finite-amplitude modifications of unstable
modes, as emerging understanding in fusion plasma instability
demonstrates.[17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]

Already in the late 1960s, using one of the early numerical simulations of
shear instability,[31] it was hinted that the nonlinear saturation of Kelvin-
Helmholtz instability involves, contrary to finite-amplitude modifications of
unstable modes, quasi-periodic oscillations in the fluctuations. Later, an
intuitive understanding of how such a phenomenon occurs in sheared fluids[32]
has been reported by invoking vortex nutation:[15] Fluctuation-amplitude
oscillations correlate with oscillations in the mean flow energy and lead to
vortex nutation. Fluctuations, however, are usually not decomposed into the
complete set of linear eigenmodes, and are commonly assumed[15] to be due to
unstable mode structures. But since unstable modes always drive a down-
gradient momentum transport, they cannot explain the increase of kinetic
energy in the mean flow.

Notably, occasional up-gradient momentum transport has been observed in
several experimental and numerical studies where an unstable shear layer
drives the turbulence.[33, 34, 35, 36, 37] Analyses of these transient events
[39, 38, 12, 40, 41] do not address the underlying conditions producing this
dynamics—whether the transient up-gradient transport is a part of an ongoing
subdominant process with occasional breakthroughs, or simply spontaneous
fluctuations. The laboratory and prior numerical experiments alone are not
sufficient to definitively answer this question. One way to expose the
underlying process is to examine the turbulent fluctuations using a complete
eigenmode decomposition, and assign roles and activities to each mode in the
transport phenomena. Indeed there can be modes other than the unstable modes
that are important in the turbulent phase, as an insightful study hints: the
dominant vortex in a turbulent background orients quasi-periodically against
(or towards) the mean flow and drives the down-gradient (or up-gradient)
momentum transport.[12] To understand such behaviors in detail, it is
instructive to also analyse how the instability saturates, a question that has
long been of interest[42] but for which understanding remains incomplete.

When turbulence is sustained via continuous energy injection from a large-
scale instability, there exist two primary candidates for instability
saturation. A common (but not necessarily justified) assumption is that energy
injected by the instability is transferred conservatively to increasingly
smaller scales in a forward, Kolmogorov-like cascade, where nonlinear
interactions move energy between linearly unstable or marginal modes until a
dissipation range is reached at small scales.[5] An alternative process
involves linearly stable eigenmodes at the large injection scales, which
absorb and remove significant energy from scales that launch the inertial
cascade. In several studies of microturbulence in fusion plasmas, linearly
stable modes have been found to be excited to significant levels via nonlinear
interactions and to drastically affect the saturated amplitudes and transport
characteristics of the system.[17, 18, 26, 27, 19, 20, 21, 29, 22, 28, 30]

Stable modes in shear flow turbulence, however, have been studied only
recently[43, 44, 45] and more remains to be understood, e.g., their role in
mixing and magnetic field evolution and how they might affect reduced models
of turbulence and transport. It was predicted in Ref. [43] that the Kelvin-
Helmholtz instability in its nonlinear evolution excites a linearly-stable
conjugate-root[46] of the inviscid instability, which affects the instability
saturation even when viscosity is finite. This was later verified in numerical
simulations of freely evolving shear layers. [44] However, the rapid
relaxation of the layer towards a stable profile on a time scale similar to
that of stable-mode excitation prevented general conclusions from being
reached, regarding how the turbulence and transport are affected by the stable
modes. The issue is aggravated by the addition of a flow-aligned magnetic
field, which causes the layer to relax even more rapidly. To circumvent this
challenge, one may drive the mean flow towards the unstable profile and thus
achieve a thorough statistical quantification of the stable modes. Note that
driven profiles are quite common in astrophysical shear flows, with forces
like gravity providing free energy for the drive. For these reasons, driven
shear flow is studied here.

The principal result of this study is that significant up-gradient momentum
transport is driven by nonlinearly excited (linearly-)stable modes, cancelling
a substantial portion of the down-gradient transport by unstable modes, and
notably this transport is present not just during turbulent momentum flux
reversals, but is continuously at work at a slightly lower level than that of
the unstable modes. This finding is robust even for variations of orders of
magnitude in background magnetic field strength, magnetic Prandtl number (or
resistivity), and forcing strength of the mean flow. Note that the stronger
background magnetic field tends to suppress the instability [46] and disrupt
the large-scale vortices, [47] while larger magnetic Prandtl number (weaker
resistivity for a fixed viscosity) extends the scale range of magnetic
fluctuations, compared to the flow fluctuations. [48] We also show, for
astrophysical applications, that a turbulent viscosity can be defined, with
the addition of stable modes, that can reliably capture the Reynolds stress:
Without stable modes, however, the stresses are greatly over-predicted by the
unstable modes.

This article is organized in the following manner. Section II entails the
magnetohydrodynamic (MHD) model of the shear flow and details the system set-
up. In Sec. III, the complete linear eigenspectrum is presented, along with a
discussion on the roles of different eigenmodes. Section IV shows the full
nonlinear evolution of MHD Kelvin-Helmholtz instability using direct numerical
simulations. A decomposition of the turbulent fluctuations onto linear
eigenmodes is performed in Sec. V, where a detailed study of imprints of
stable modes in turbulence and transport is presented. Section VI summarizes
the findings of this work.

##  II Model and simulation set-up

An incompressible magneto-fluid is considered in this study, and standard MHD
equations are adopted:

$\displaystyle\mathbf{\nabla}\cdot\mathbf{u}=0,$ (1a)
$\displaystyle\partial_{t}\mathbf{u}+\mathbf{u}\cdot\mathbf{\nabla}\mathbf{u}=-\frac{\nabla
P}{\rho}+\frac{\left(\nabla\times\mathbf{B}\right)\times\mathbf{B}}{4\pi\rho}+\nu\nabla^{2}\mathbf{u}+\mathbf{f},$
(1b) $\displaystyle\mathbf{\nabla}\cdot\mathbf{B}=0,$ (1c)
$\displaystyle\partial_{t}\mathbf{B}=\nabla\times\left(\mathbf{u}\times\mathbf{B}\right)+\eta\nabla^{2}\mathbf{B},$
(1d)

where $\mathbf{u}$, $\mathbf{B}$, $P$, $\rho$, $\nu$, $\eta$, and $\mathbf{f}$
respectively denote the fluid velocity, magnetic field, pressure, fluid
density, viscosity, ohmic diffusivity, and externally supplied acceleration to
the magneto-fluid.

###  II.1 Background flow, magnetic field, and forcing

A shear layer is examined on a two-dimensional $(x,z)$ plane with the initial
fluid velocity given by
$\mathbf{u}(x,z,t=0)=U_{0}\mathrm{tanh}(z/a)\hat{\mathbf{x}}$ and a flow-
aligned magnetic field, initially uniform, as
$\mathbf{B}(x,z,t=0)=B_{0}\hat{\mathbf{x}}$. The parameters
$a,U_{0},\mathrm{and\ }B_{0}$ signify the half-width of the flow-shear,
maximum initial fluid velocity, and initial magnetic field, respectively.
These parameters are exploited to non-dimensionalize all the variables
henceforth. Length, time, and energy (per unit mass) are hereafter measured in
units of $a$, $a/U_{0}$, and $U_{0}^{2}$, respectively. Thus the initial (or
reference) flow and magnetic field are represented by
$U_{\mathrm{ref}}(z)=\mathrm{tanh}(z)$ and $B_{\mathrm{ref}}(z)=1$ in the rest
of this article. The ratio of the maximum fluid speed $U_{0}$ to the Alfvén
speed can be written as the Alfvénic Mach number
$M_{\mathrm{A}}=U_{0}\sqrt{4\pi\rho}/B_{0}$. The viscosity and resistivity are
quantified via fluid Reynolds number $Re=U_{0}a/\nu$ and magnetic Reynolds
number $Rm=U_{0}a/\eta$, respectively.

In two dimensions, a more convenient formalism is available, using the
streamfunction $\phi$ and flux function $\psi$. Defining
$\mathbf{u}=\mathbf{\hat{{y}}}\times\nabla\phi$ and
$\mathbf{B}=\mathbf{\hat{{y}}}\times\nabla\psi$, the vorticity and the current
become $\nabla^{2}\phi\mathbf{\hat{{y}}}$ and
$\nabla^{2}\psi\mathbf{\hat{{y}}}$, respectively. Taking the curl of Eq. (1b),
and rewriting Eq. (1d) in terms of the stream- and flux-functions yields [49]

$\displaystyle\begin{split}&\partial_{t}\nabla^{2}\phi+\\{\nabla^{2}\phi,\phi\\}=M_{\mathrm{A}}^{-2}\\{\nabla^{2}\psi,\psi\\}+Re^{-1}\nabla^{4}\phi+\partial_{z}f(k_{x}\mathrm{=}0,z,t),\end{split}$
(2a) $\displaystyle\partial_{t}\psi=\\{\phi,\psi\\}+Rm^{-1}\nabla^{2}\psi,$
(2b)

where the Poisson bracket is
$\\{P,Q\\}=\partial_{x}P\cdot\partial_{z}Q-\partial_{z}P\cdot\partial_{x}Q$;
e.g., $\\{\phi,\psi\\}=-\mathbf{u}\cdot\mathbf{\nabla}\psi$. Here, $k_{x}$ is
the Fourier wavenumber along the $x$-axis. The parameters $Re=Rm=500$ are
chosen for all simulations unless mentioned otherwise (where $Rm$ is changed
to $50$ and $5\,000$ in different simulations). It should be emphasized that
these Reynolds numbers are defined using the initial scale $a$ of the sharpest
gradient in the flow as the characteristic length scale; however, as the
system evolves nonlinearly via vortex merging, despite the forced mean flow,
eddies of the size of the simulation box appear, which may be considered as
the characteristic length scale of motion.[11] When choosing this
normalization, non-dimensional numbers should be scaled accordingly, e.g.,
$Rm=5\,000$ becomes $Rm=5\,000\times L_{x}\approx 1.5\times 10^{5}$, where
$L_{x}$ represents the box-size along the mean flow direction. The external
body force, $\mathbf{f}=f(k_{x}\mathrm{=}0,z,t)\mathbf{\hat{x}}$, is applied
to the mean flow only, which is highlighted in Eq. 2a using the explicit
mention of $k_{x}\mathrm{=}0$. As in a recent study,[45] the forcing drives
the instantaneous mean flow towards the initial unstable profile
$U_{\mathrm{ref}}(z)$. A similar forcing mechanism exists for geo- and
astrophysical flows where gravitation[50] tends to build shear layers. We
assume here such a force, represented as a Krook operator, [51, 52, 53] as

$f=D_{\mathrm{Krook}}\left[U_{\mathrm{ref}}(z)-\langle
U(x,z,t)\rangle_{x}\right]+F_{0},$ (3)

where $D_{\mathrm{Krook}}$, sometimes also referred to as the profile
relaxation rate,[54] measures the forcing strength (in units of $U_{0}/a$);
and $\langle U(x,z,t)\rangle_{x}$ represents the instantaneous $x$-averaged
flow. If $D_{\mathrm{Krook}}=0$, the shear layer evolves freely and decaying
turbulence is realized as a result of the Kelvin-Helmholtz instability and the
turbulence it generates.

The time-independent force $F_{0}$ is implemented only to balance the viscous
diffusion of the initial shear layer:
$Re^{-1}\nabla^{2}U_{\mathrm{ref}}(z)+F_{0}=0$ ensures an initial equilibrium
state to which small-amplitude perturbations are added before the system is
evolved.

###  II.2 Initial and boundary conditions

As in the unforced study, [44] a simulation box of $L_{x}=10\pi$ is
considered, but double the size along the $z$-axis ($L_{z}=20\pi$), given that
the quasi-stationary turbulence simulated herein is run for much longer time,
which tends to create fully developed turbulent features that are larger in
size. Thus we adopt a larger domain to minimize their potential interactions
with the boundaries in the $z$-axis. Note the forcing applied to the mean flow
prevents profile relaxation and the turbulence remains mostly in the vicinity
of the shear layer. The numerical code Dedalus, [55] a pseudospectral solver,
is used in this study. Fourier modes along the $x$-axis and Chebyshev
polynomials along the $z$-axis are employed with $(N_{x},N_{z})=(2048,2048)$
spectral modes. We confirmed that the spectral energy density and dissipation
are converged at this resolution. Note also that these high resolutions
benefit the eigenmode projection of nonlinear data in the post-processing
analysis. Only for the simulation with magnetic Prandtl number of $10$, the
box size was changed to $(L_{x},L_{z})=(6\pi,8\pi)$ and the resolution was
increased to $(N_{x},N_{z})=(4096,4096)$; the same simulation was repeated
with $(N_{x},N_{z})=(4096,8192)$, but only for early times due to
computational cost, and found to reproduce, among others, the energy
evolution. All simulations use $3/2$ dealiasing rule, additionally. The
boundary conditions used in all simulations are periodic along the $x$-axis;
and perfectly conducting, no-slip, co-moving (with the initial flow) at the
top and bottom boundaries, $z=\pm L_{z}/2$.[44, 45]

The initial equilibrium state is seeded with small-amplitude perturbations
$(\widetilde{\phi},\widetilde{\psi})$ at all Fourier wavenumbers, as[44]

$\displaystyle\widetilde{\phi}(x,z,t=0)$
$\displaystyle=A_{\phi}\sum_{k_{x}\neq
0}k_{x}^{a}e^{ir_{\phi}(k_{x})}e^{-z^{2}/\sigma^{2}}e^{ik_{x}x},$ (4)

and

$\widetilde{\psi}(x,z,t=0)=A_{\psi}\sum_{k_{x}\neq
0}k_{x}^{a}e^{ir_{\psi}(k_{x})}e^{-z^{2}/\sigma^{2}}e^{ik_{x}x}.$ (5)

Here, $A_{\phi}$ and $A_{\psi}$ set the overall amplitudes of the
perturbations that have a Gaussian width controlled by $\sigma$ and the rate
at which they fall-off with the wavenumbers given by $a$. The random phases
$r_{\phi}(k_{x})$ and $r_{\psi}(k_{x})$, forming a uniform distribution in
$[0,2\pi)$, are issued for each different $k_{x}$ using a pseudo-random number
generator. Different choices of these initial conditions were investigated in
Ref. [44], motivating the choice here: $a=-1$, $\sigma=2$, and
$A_{\phi}=A_{\psi}=10^{-3}$. This set of parameters offers distinct linear and
nonlinear phases of evolution.

##  III Linear eigenmodes

Aiming to understand the nonlinear excitation of linear eigenmodes in the
turbulent phase, first the nonlinear initial-value problem is solved to
collect high-fidelity turbulent data. Afterward, a separate eigenvalue problem
is solved to obtain a complete linear eigenspectrum and eigenmodes, which are
used to expand the nonlinear data on this basis to track the amplitude of each
eigenmode. Such a basis is obtained by linearizing the governing equations
around the initial flow and magnetic field profiles, by dropping the
dissipative terms. The eigenmodes thus obtained are of a dissipationless
linear operator. Of course, the meaning and utility of this linear basis is a
priori unknown. Nevertheless, when a basis forms a complete set, one can
always expand an arbitrary fluctuation on that basis. As the non-dissipative
equations of motion preserve Parity-Time (PT-)reversal symmetry, such a system
is theoretically guaranteed to yield a complete basis as established recently
in PT-symmetric quantum mechanics.[56] Previous studies in gyrokinetic and MHD
plasmas have also revealed the usefulness of dissipationless linear eigenmodes
in interpreting dissipative nonlinear systems. [27, 44, 57]

###  III.1 Complete eigenspectrum

With the intent of obtaining dissipationless linear eigenmodes, the variables
$(\phi,\psi)$ in Eqs. (2a) and (2b) are decomposed into background and
perturbations,
$(\phi,\psi)=(\phi_{\mathrm{ref}},\psi_{\mathrm{ref}})+(\widetilde{\phi},\widetilde{\psi})$.
The linearized dissipationless equations for the evolution of perturbations
are

$\displaystyle{}\partial_{t}\nabla^{2}\widetilde{\phi}$
$\displaystyle=-\left[U_{\mathrm{ref}}\partial_{x}\nabla^{2}-(\partial_{z}^{2}U_{\mathrm{ref}})\cdot\partial_{x}\right]\widetilde{\phi}+\frac{1}{M_{\mathrm{A}}^{2}}\left[B_{\mathrm{ref}}\partial_{x}\nabla^{2}-(\partial_{z}^{2}B_{\mathrm{ref}})\cdot\partial_{x}\right]\widetilde{\psi},$
(6a) $\displaystyle\partial_{t}\widetilde{\psi}$
$\displaystyle=-U_{\mathrm{ref}}\partial_{x}\widetilde{\psi}+B_{\mathrm{ref}}\partial_{x}\widetilde{\phi}.$
(6b)

Fourier transforming along the $x$-axis and assuming time variation at each
Fourier wavenumber takes the form
$\hat{\phi}(k_{x},z,\omega)e^{i\omega(k_{x})t}$, Eqs. (6a)-(6b) become

$\displaystyle\omega\left(\partial_{z}^{2}-k_{x}^{2}\right)\hat{\phi}$
$\displaystyle=-k_{x}\left[U_{\mathrm{ref}}\left(\partial_{z}^{2}-k_{x}^{2}\right)-(\partial_{z}^{2}U_{\mathrm{ref}})\right]\hat{\phi}+\frac{1}{M_{\mathrm{A}}^{2}}k_{x}\left[B_{\mathrm{ref}}\left(\partial_{z}^{2}-k_{x}^{2}\right)-(\partial_{z}^{2}B_{\mathrm{ref}})\right]\hat{\psi},$
(7a) $\displaystyle\omega\hat{\psi}$
$\displaystyle=-U_{\mathrm{ref}}k_{x}\hat{\psi}+B_{\mathrm{ref}}k_{x}\hat{\phi}.$
(7b)

Figure 1: (a) Linear eigenspectrum of the MHD shear-flow system (at
$k_{x}=0.2$ with $M_{\mathrm{A}}=10$). The unstable and stable modes are shown
with thick (colored) crosses. Among the continuum modes that form a vertical
line, a zero-frequency continuum mode is displayed with a green-colored star.

(b)–(d) Eigenfunctions in $z$-space, with real (Re) and imaginary (Im) parts,
for unstable ($\phi_{1}$), stable ($\phi_{2}$), and one continuum ($\omega=0$)
mode. (e)–(g) Corresponding eigenmode structures in ($x,z$) space. Note that
the eigenmodes $\phi_{1}$ and $\phi_{2}$ are complex conjugate to each other.
Imaginary parts in their eigenfunctions induce relative tilt between them in
$(x,z)$ space, which will be consequential for momentum transport in Sec.
III.2. Each eigenmode is normalized to have unit total energy [following which
the maximum values of $\phi$ in (b)–(g) are chosen].

Solving Eqs. (7a)-(7b), the eigenvalues $\omega$ are found to be real except
when $0<|k_{x}|<1$, where two of the real eigenvalues coalesce to produce
imaginary eigenvalues, [58] as complex conjugate to each other. These are the
growth rates of the unstable eigenmode and its conjugate stable eigenmode,
which evolve in time as $\textrm{e}^{\gamma(k_{x})t}$ and
$\textrm{e}^{-\gamma(k_{x})t}$, respectively. This mode-pair is shown, for the
first Fourier wavenumber $k_{x}=2\pi/L_{x}=0.2$, in Fig. 1(a), along with all
the purely real eigenvalues. The latter constitute the eigenmode continuum
[59] and are theoretically infinite in number, although numerical
discretization always yields a finite but very large number of modes ($>3,000$
for each wavenumber in this study). These eigenvalues are given by the
relation $\omega/k_{x}+U_{\mathrm{ref}}(z)\pm v_{\mathrm{A,ref}}(z)=0$, where
$v_{\mathrm{A,ref}}(z)$ is the Alfvén speed along the reference magnetic field
at the vertical coordinate $z$.

The eigenfunctions, normalized to have unit total energy, are also shown in
Fig. 1: along the $z$-axis, see Figs. 1(b)–(d), and in $(x,z)$ space, see
Figs. 1(e)–(g). Note that complex conjugation transforms the unstable mode
$\phi_{1}(k_{x},z)$ into the stable mode $\phi_{2}(k_{x},z)$ and vice-versa.
This is a direct consequence of spontaneous PT-symmetry breaking in the ideal
shear-flow instability. [58] (The spontaneous symmetry breaking does not imply
that the equation of motion or the associated Hamiltonian breaks PT-symmetry;
it is rather some of the eigenfunctions of such a PT-symmetry-preserving
Hamiltonian that break PT-symmetry.)

A representative eigenfunction of a continuum mode, shown in Fig. 1(d),
exhibits sharp and narrow structure. To what physics each type of eigenmode
structure contributes will be explored in this article.

###  III.2 Competing roles of unstable and stable modes

Shown in Fig. 2 is a schematic diagram, illustrating how the relative tilts in
the eddies can transport momentum in opposite directions across the shear
layer.[14] It can be qualitatively observed from Figs. 1(e) and (f) that the
unstable and stable modes, drive down- and up-gradient momentum transport,
respectively. Precise quantative measures will be built and computed later in
Sec. V.4.

Figure 2: (a) The unstable mode of the flow transport momentum in the down-
gradient direction: $-x$-directed momentum at A is carried to
A$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$ and $+x$-directed
momentum at B is carried to
B$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$. Fluxes
A$\to$A$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$ and
B$\to$B$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$ act to relax
the mean flow gradient (shown with the long horizontal arrows). (b) Oppositely
tilted eddies, which correspond to a stable mode, transport momentum in the
up-gradient direction: $-x$-directed momentum at C is carried to
C$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$ and $+x$-directed
momentum at D is carried to
D$\mathbin{\rotatebox[origin={c}]{21.0}{${}^{\prime}$}}$. Both of these fluxes
replenish the mean flow. The direction of the streamlines (shown with grey
arrows on the elliptic eddies) does not alter these properties, but the tilt
does.

Since the unstable and stable modes compete with each other to transport
momentum in opposing directions, the excitation levels of these modes are
crucial. In the linear phase of instability evolution, the transport by the
unstable modes dominates over the transport by the stable modes. However, this
need not be the case in the nonlinear phase, as nonlinear processes can excite
the stable modes to appreciable levels. Whenever the stable modes surpass the
unstable modes in amplitudes, net momentum is transported in the up-gradient
direction. [60, 44] In extremely simplified models of transport, such as eddy
viscosity models, this contributes to negative eddy viscosity. Computing the
amplitude of each eigenmode in the nonlinear phase can thus be helpful to
build improved reduced transport models. A recent investigation also
demonstrated that this kind of competition between the two large-scale
eigenmodes alters the magnetic cascade substantially. [45]

##  IV Nonlinear evolution

Having provided a description of linear eigenmodes, we now turn to properties
of the nonlinear system, before discussing how expressions of linear modes may
be identified in turbulent fluctuations.

###  IV.1 Finte-amplitude Kelvin-Helmholtz instability

Small-amplitude perturbations in the flow and magnetic field evolve
exponentially fast in the linear regime of the instability, giving rise to a
chain of spiral vortices, as evident in Figs. 3(a) and (d). These structures
then interact nonlinearly with nearby vortices to yield even larger turbulent
structures as in Figs. 3(b) and (c). A contrast is to be made between forced
and unforced simulations. In the latter, the gradient of the mean flow
flattens out as the instability extracts energy. Decaying turbulence then
ensues. Forcing the mean flow, however, leads to a quasi-stationary
turbulence, as the energy in the gradient is replenished with the instability
drawing on its energy. In the saturated stage, energy input through the
unstable modes is balanced by energy removal via stable modes as well as
dissipative channels.

Figure 3: Time evolution of vorticity in (a)–(c) in a simulation with a forced
background flow, $D_{\mathrm{Krook}}=2$; and (d)–(f) in a simulation with a
freely evolving shear layer, $D_{\mathrm{Krook}}=0$; both for
$M_{\mathrm{A}}=10$. Panels (a)–(c) share the same colorbar, and (d)–(f) share
another colorbar. The instantaneous mean flow profile in each of the subplots
is shown with a black dotted curve, where the vertical axis represents the
$z$-coordinate and the horizontal direction corresponds to the $x$-velocity
$U_{0}(z,t)$, as exemplified in the inset of (a). Two arrows pointing in
opposite directions show the direction of the flow in the regions $z>0$ and
$z<0$. The initial flow profile $U_{0}(z,t=0)=\mathrm{tanh}(z)$ is shown with
a red dashed curve in (a) and (d). Rapid profile flattening is evident in (d).
While the instability dies out in the unforced case, quasi-stationary
turbulence is realized in the forced case in (c).

###  IV.2 Momentum transport

Figure 4: Time evolution of MHD stresses. (a) Reynolds stress $\tau_{u}(t,z)$.
(b) Reynolds stress $\tau_{u}(t,z=0)$ at the middle of the layer, $z=0$. (c)
Maxwell stress $\tau_{b}(z,t)$. (d) Maxwell stress $\tau_{b}(t,z=0)$ at $z=0$.
All data shown are for a single simulation with $M_{\mathrm{A}}=10$,
$D_{\mathrm{Krook}}=2$. The Reynolds stress in (b) reverses several times, in
contrast to the Maxwell stress in (d), which is almost always down-gradient.

It is now timely to discuss the turbulent transport of momentum in nonlinear
simulations. To derive the turbulent stresses, the evolution equation of the
mean flow can be written by $x$-averaging the momentum equation,

$\frac{\partial}{\partial t}\langle U\rangle_{x}=-\frac{\partial}{\partial
z}\left(\tau_{u}+\tau_{b}\right)+D_{\mathrm{Krook}}\left[U_{\mathrm{ref}}(z)-\langle
U\rangle_{x}\right]+\frac{1}{Re}\frac{\partial^{2}}{\partial
z^{2}}\left[\langle U\rangle_{x}-U_{\mathrm{ref}}(z)\right],$ (8)

where $U=U(x,z,t)$ represents the instantaneous flow,
$\langle\cdot\rangle_{x}$ signifies $x$-averaging operation, and $\tau_{u}$
and $\tau_{b}$ are the Reynolds and Maxwell stresses, arising from the
correlations of turbulent fluctuations of the flow and the magnetic fields,
respectively. Note that in Fraser et al.[44], a negative sign was
typographically missed in front of the first term on the right-hand side of
Eq. (8). With the sign displayed in Eq. (8) above, the turbulent stresses are
given by

$\displaystyle\tau_{u}$
$\displaystyle=\langle\widetilde{u}_{x}\widetilde{u}_{z}\rangle_{x}=-\langle\partial_{z}\widetilde{\phi}\cdot\partial_{x}\widetilde{\phi}\rangle_{x},$
(9a) $\displaystyle\tau_{b}$
$\displaystyle=-\frac{1}{M_{\mathrm{A}}^{2}}\langle\widetilde{b}_{x}\widetilde{b}_{z}\rangle_{x}=\frac{1}{M_{\mathrm{A}}^{2}}\langle\partial_{z}\widetilde{\psi}\cdot\partial_{x}\widetilde{\psi}\rangle_{x}.$
(9b)

These stresses are evaluated from nonlinear simulations and shown in Fig. 4.
Fluctuations of Reynolds stress are concentrated in the shear layer, near
$z\approx 0$. Time histories of the Reynolds and Maxwell stresses, at $z=0$,
where they are largest in magnitude, are compared in Figs. 4(b) and 4(d). Note
the recurring dominant up-gradient transport via the Reynolds stress. The
Maxwell stress, however, is almost always down-gradient. Figures 4(a) and 4(c)
also convey that the Maxwell stress is generally broader along the $z$-axis
than the Reynolds stress, which is more localized near the shear layer.

##  V Decomposition of nonlinear simulation onto linear modes

To probe the nonlinear simulation data, the turbulent fluctuations are
expanded on the linear eigenmode basis described in Sec. III. Consider an
arbitrary turbulent fluctuation
$\widetilde{\chi}_{\mathrm{turb}}=(\widetilde{\phi}_{\mathrm{turb}},\widetilde{\psi}_{\mathrm{turb}})$,
which is expanded as

$\widetilde{\chi}_{\mathrm{turb}}(x,z,t)=\sum_{k_{x}\neq
0}\textrm{e}^{ik_{x}x}\sum_{j}\beta_{j}(k_{x},t)\chi_{j}(k_{x},z)$ (10)

where the eigenmode basis $\chi_{j}(k_{x},z)$ is employed along the $z$-axis
at each wavenumber $k_{x}$ to decompose the fluctuations. The complex mode-
amplitude $\beta_{j}(k_{x},t)$, defined for each eigenmode $j$, can then be
computed using properties of the linear operator, described in the Appendix A,
even when the eigenmodes of the operator are non-orthogonal, as is the case
here.

Following earlier studies, [29, 43, 27, 44, 45] $j=1,2$ will be used to
represent unstable and stable modes, respectively. The computations herein
resolve as many as $3109$ eigenmodes at a particular $k_{x}$.

###  V.1 Nonlinear excitation of stable modes

The amplitudes of the unstable and stable modes are tracked in the nonlinear
simulations, and their time series are plotted in Fig. 5(a). As expected, the
unstable mode grows and the stable mode decays exponentially in the early
phase. However, as the fluctuations increase due to the growth of the unstable
modes, nonlinear interactions among them begin exciting the linearly stable
mode,[30] causing it to rise to almost the same level as the unstable mode at
that wavenumber, see Fig. 5(a). Later, in the fully nonlinear stage, all
eigenmodes can participate in the energy redistribution.

The energy in individual eigenmodes $|\beta_{j}|^{2}$, averaged over a
turbulent state ($t=150\textrm{--}1000$), is displayed in Fig. 5(b). It is
evident that the unstable and stable eigenmode pair contains a majority
($>70\%$) of the energy in the system. The remaining eigenmodes share a wide
spectrum of the remaining energy. This suggests that the turbulent system at
hand may be amenable to a substantial dimensionality reduction.[27, 19] For
the cases of the weaker magnetic fields, this finding is more prominent, as
evidenced in the Appendix B. In addition, the weaker fields support more
coherent amplitude-oscillations, unlike the large excursions in the amplitudes
observed with the stronger fields, e.g., $M_{\mathrm{A}}=10$ in Fig. 5(a). In
the latter case, the stronger Lorentz back-reaction acting on the large-scale
turbulent flow cause strong oscillations in the eigenmode amplitudes.

Figure 5: (a) Time traces of the eigenmode amplitudes are shown for
$k_{x}=0.2$ for a simulation with $M_{\mathrm{A}}=10$ and
$D_{\mathrm{Krook}}=2$. Note, in the inset, the nonlinear excitation of
linearly stable mode ($|\beta_{2}|$) in instability saturation ($|\beta_{1}|$
is the unstable mode amplitude). (b) All 3109 eigenmodes at $k_{x}=0.2$ are
plotted with their squared excitation levels in the nonlinearly saturated
phase, which represent the energy in each eigenmode. The diameter of each
circle shown corresponds to the energy in each eigenmode, and modes with lower
energy are plotted on top of more highly excited modes, such that all data
points are (partially) visible. Note that the total fluctuation energy is
composed of both the modal and non-modal energy because of the non-orthogonal
modes. Evaluating total energy at a wavenumber, $E=\int
dz\left[|\mathbf{u}|^{2}+|\mathbf{B}|^{2}/M_{\mathrm{A}}^{2}\right]/2=\int
dz\left[(\sum_{m}\beta_{m}\mathbf{u}_{m})\cdot(\sum_{n}\beta_{n}\mathbf{u}_{n})^{\ast}+(\sum_{m}\beta_{m}\mathbf{B}_{m})\cdot(\sum_{n}\beta_{n}\mathbf{B}_{n})^{\ast}/M_{\mathrm{A}}^{2}\right]/2=\sum_{m,n}E_{mn}$,
where $(\mathbf{u}_{m},\mathbf{B}_{m})$ represents the $m$-th eigenmode. When
$m$ and $n$ belong to discrete (d) modes, $E_{\mathrm{dd}}$ is, upon time-
averaging ($t=150\textrm{--}1000$), around $72\%$ of the total energy, whereas
when $m$ and $n$ belong to continuum (c) modes, $E_{\mathrm{cc}}$ is $\approx
22\%$; $E_{\mathrm{dc}}$ is $\approx 6\%$.

###  V.2 Reduced representation of the turbulent flow

To obtain a better understanding of turbulent dynamics, it is of interest to
compare different components of eigenmodes in the turbulent flow. An
approximate (reduced) representation of the turbulent flow
$\widetilde{\phi}_{\mathrm{approx}}$ can be constructed from a class of
eigenmodes at each wavenumber, e.g.,
$\widetilde{\phi}_{\mathrm{approx}}(x,z,t)$ can be written as a sum of an
unstable mode per wavenumber
$\sum_{k_{x}}\textrm{e}^{ik_{x}x}\beta_{1}(k_{x},t)\phi_{1}(k_{x},z)$, or as a
sum of an unstable and a stable mode per wavenumber
$\sum_{k_{x}}\textrm{e}^{ik_{x}x}\left[\beta_{1}(k_{x},t)\phi_{1}(k_{x},z)+\beta_{2}(k_{x},t)\phi_{2}(k_{x},z)\right]$.
Respective short-hand notations $\beta_{1}\phi_{1}$ and
$\beta_{1}\phi_{1}+\beta_{2}\phi_{2}$ will be used hereafter, i.e.,

$\displaystyle\beta_{1}\phi_{1}$
$\displaystyle\equiv\sum_{0<|k_{x}|<1}\textrm{e}^{ik_{x}x}\beta_{1}(k_{x},t)\phi_{1}(k_{x},z),$
(11a) $\displaystyle\beta_{1}\phi_{1}+\beta_{2}\phi_{2}$
$\displaystyle\equiv\sum_{0<|k_{x}|<1}\textrm{e}^{ik_{x}x}\left[\beta_{1}(k_{x},t)\phi_{1}(k_{x},z)+\beta_{2}(k_{x},t)\phi_{2}(k_{x},z)\right].$
(11b)

The nonlinear fluctuations of the flow are compared in Fig. 6, viewed at
different levels of truncation in the eigenmode expansion. The leftmost panel,
Fig. 6(a), shows the full turbulent fluctuations in the Kelvin-Helmholtz
(KH-)unstable wavenumbers $k_{x}=0.2,0.4,0.6,0.8$, which appear similar to the
full turbulent fluctuations that include all wavenumbers in the nonlinear
simulation (not shown); Fig. 6(b) displays the sum of unstable eigenmodes at
each of these KH-unstable wavenumbers; and Fig. 6(c) presents the sum of
unstable and stable eigenmodes at the same wavenumbers, while omitting all
continuum modes. Adding stable modes produces a substantial improvement in the
reconstruction. Note that such a reconstruction was found to deteriorate
quickly over time (i.e., a few instability $\mathrm{e}$-folding times where
one $\mathrm{e}$-folding time for the fastest growing mode $k_{x}=0.4$ is
$\gamma^{-1}\approx 5$) in the study of unforced shear layers, [44] as the
rapid relaxation of the layer towards a stable profile rendered the unstable
and stable eigenmodes of the system to be less representative of the decaying
turbulence. The turbulent fluctuation shown in Fig. 6 is at $t=702$, which
lies well within the nonlinear phase (the linear phase ends around $t\approx
30$). In this respect, the forced shear layer is markedly different from the
freely evolving layer.

Figure 6: (a) Full turbulent fluctuations in streamfunction in the Kelvin-
Helmholtz-unstable wavenumbers $0<|k_{x}|<1$ (thus called
$\widetilde{\phi}_{\mathrm{filtered}}$), observed in nonlinear simulations
with $M_{\mathrm{A}}=30$, $D_{\mathrm{Krook}}=2$. The shown plot of
fluctuations includes all types of eigenmodes—the unstable, stable and
continuum modes. (b) Reconstruction of the turbulent fluctuations by summing
only the unstable modes at the same wavenumber range. (c) Reconstruction by
adding stable and unstable modes, while omitting all continuum modes. The
reconstruction in (c) is clearly much alike the turbulent fluctuations in (a),
in contrast to the reconstruction in (b). Saturation theory of instability
that considers the unstable modes only, at best, can produce (b), but with
inclusion of the stable modes, substantial improvement can be achieved.

###  V.3 Performance of reduced representations

While the qualitative analysis of the turbulent-flow reconstruction in Sec.
V.2 is instructive, a quantitative measurement is desirable. To this end,
following Ref. [44], the reconstructive capability of reduced representations
is quantified, at each time step in the simulation, using the standard energy
norm that measures the fraction of kinetic energy lost when the eigenmode
basis is truncated, compared to the kinetic energy in the full turbulent flow
data—see the definition in Eq. (12). The energy norm is well-suited for
studying large-scale structures. Small-scale phenomena, however, may not be
amenable to such analysis, although one may be able to find ties between the
small- and large-scale pheonomena in some cases. This measure is also called a
“truncation error.” Note that this error arises not in the nonlinear
simulations but merely in the reduced representations of turbulent
fluctuations, when truncating the eigenmode basis in post-processing analyses.

Figure 7: Time traces of residuals, i.e., the fraction of energy missed in
truncated bases, normalized to the total energy in the turbulent flow at each
time step. The reconstruction uses truncated bases with unstable modes alone,
and unstable and stable modes together. The forcing strength is varied in
three different simulations with $M_{\mathrm{A}}=10$: (a)
$D_{\mathrm{Krook}}=1$, (b) $D_{\mathrm{Krook}}=0.1$, and (c)
$D_{\mathrm{Krook}}=0$. The unforced shear layer in (c) rapidly flattens out,
and thus instability no longer drives the turbulence. As long as the
turbulence is driven by the instability, the unstable and stable modes
together can reconstruct a large fraction of the turbulent flow features in
(a) and (b).

Using the energy norm, we define the relative truncation error, which may also
be called a normalized residual, in the following manner:

$\displaystyle\mathrm{Residual}=\frac{||\widetilde{\phi}_{\mathrm{exact}}-\widetilde{\phi}_{\mathrm{approx}}||_{\mathrm{KE}}^{2}}{||\widetilde{\phi}_{\mathrm{exact}}||_{\mathrm{KE}}^{2}}=\frac{||\widetilde{\phi}_{\mathrm{diff}}||_{\mathrm{KE}}^{2}}{||\widetilde{\phi}_{\mathrm{exact}}||_{\mathrm{KE}}^{2}}=\frac{\int
dxdz\left[\left(\partial_{x}\widetilde{\phi}_{\mathrm{diff}}\right)^{2}+\left(\partial_{z}\widetilde{\phi}_{\mathrm{diff}}\right)^{2}\right]}{\int
dxdz\left[\left(\partial_{x}\widetilde{\phi}_{\mathrm{exact}}\right)^{2}+\left(\partial_{z}\widetilde{\phi}_{\mathrm{exact}}\right)^{2}\right]},$
(12)

where $(\partial_{x}\phi)^{2}$ and $(\partial_{z}\phi)^{2}$ are the squared
$z$\- and $x$-components of velocities;
$\widetilde{\phi}_{\mathrm{diff}}=\widetilde{\phi}_{\mathrm{exact}}-\widetilde{\phi}_{\mathrm{approx}}$
with $\widetilde{\phi}_{\mathrm{exact}}$ and
$\widetilde{\phi}_{\mathrm{approx}}$ representing respectively the turbulent
streamfunction from nonlinear simulation and its reduced representation—either
a summation over the unstable modes alone or over the unstable and stable
modes together—both spanning fluctuations over a range of wavenumbers. Here,
this range, taken to be the same for both, is considered to be $0<|k_{x}|<1$,
which corresponds to the wavenumber range of the instability. If the residual
is less than unity, the truncation in the eigenmode expansion may be
considered as a representative of the full system and thus a candidate for
reduced-order model building. On the contrary, the residual being around unity
or more signfies the failure of the reduced representation in effectively
capturing the overall nonlinear fluctuations.

Figure 8: Shown are the time-averaged residuals for simulations different
$M_{\mathrm{A}}$. Note the residuals are the fractions of energy missed in the
truncated bases, compared to the total energy in the instantaneous full
turbulent flow. Time-averaging is performed over a quasi-stationary state of
turbulence $(t=600\text{--}1000)$. The reconstruction uses truncated bases
with unstable modes alone, and unstable and stable modes together, leaving all
the continuum modes. All simulations use $D_{\mathrm{Krook}}=2$. Note the
dramatic improvement with the inclusion of the stable modes. For
$M_{\mathrm{A}}=3$, the improvement is modest, as the stronger Lorentz force
back-reacts on the large-scale turbulent flow, producing more fluctuations in
the continuum modes.

The time evolution of the residuals is compared in Fig. 7 for varying forcing
strengths. As expected, the unstable modes entirely capture the fluctuations
in the linear phase (i.e., $t\lesssim 30$). In the nonlinear phase, however,
the unstable modes capture only a rather limited fraction of the turbulent
fluctuations. This is greatly improved when the stable modes are added. This
suggests that the success of quasilinear models in capturing key properties of
the turbulence can crucially depend on whether stable modes are considered
when constructing such models.

It is also interesting to note that the turbulence in the unforced shear
layer, see Fig. 7(c), is different from the forced cases. In the former, the
shear layer quickly flattens out and nearly shuts off the instability, leading
to a decaying turbulence. Regardless of whether the unstable and/or stable
modes are considered, the corresponding reconstructions fail to model the
turbulence with any degree of accuracy. By contrast, when the shear layer is
forced, a reduced representation of the turbulent flow with two modes (per
wavenumber) is found to perform well, recovering a substantial fraction of the
full nonlinear system.

A similar reconstruction is shown for various strengths of magnetic fields
$M_{\mathrm{A}}$ in Fig. 8, where residuals are time-averaged over a quasi-
stationary state of turbulence. With stronger magnetic fields (lower
$M_{\mathrm{A}}$), the vortices begin disrupting due to stronger Lorentz force
and consequently generate more fluctuations at scales beyond the Kelvin-
Helmholtz-instability (KHI) range. [47] This accounts for an increase of the
residual for low $M_{\mathrm{A}}$, although it remains below $0.2$ for
$M_{\mathrm{A}}=10$. For $M_{\mathrm{A}}=3$, the improvement with the
inclusion of the stable modes is modest. Momentum transport by large-scale
structures, formed from the unstable and stable modes, within the KHI range,
however, may still dominate over the transport contributed by much smaller
scales or the remaining continuum modes; hence, a quantitative analysis of
transport will be conducted next.

###  V.4 Competing up- and down-gradient momentum transport and their reduced
models

The Reynolds stress can be expressed in terms of the contribution from each
wavenumber, which can further be decomposed into the contribution from each
eigenmode. At a wavenumber $k_{x}$, the Reynolds stress from all the
fluctuations $\hat{\phi}_{k_{x}}$ read

$\tau_{u}(\mathrm{all\ modes})=2\mathrm{\
Im}[k_{x}\hat{\phi}_{k_{x}}\cdot\partial_{z}\hat{\phi}_{k_{x}}^{\ast}],$ (13)

whereas the contribution from an unstable mode alone, and from an unstable and
a stable mode alone, at that wavenumber are respectively given as

$\displaystyle\tau_{u}(\mathrm{unstable})$ $\displaystyle=2\mathrm{\
Im}[k_{x}(\beta_{1}\phi_{1,k_{x}})\cdot\partial_{z}(\beta_{1}{\phi}_{1,k_{x}})^{\ast}]$
(14) $\displaystyle=2|\beta_{1}|^{2}\mathrm{\
Im}[k_{x}\phi_{1,k_{x}}\cdot\partial_{z}{\phi}_{1,k_{x}}^{\ast}],$

and

$\displaystyle\tau_{u}(\mathrm{stable})$
$\displaystyle=2|\beta_{2}|^{2}\mathrm{\
Im}[k_{x}\phi_{2,k_{x}}\cdot\partial_{z}{\phi}_{2,k_{x}}^{\ast}]$ (15)
$\displaystyle=2|\beta_{2}|^{2}\mathrm{\
Im}[k_{x}\phi_{1,k_{x}}^{\ast}\cdot\partial_{z}{\phi}_{1,k_{x}}]$
$\displaystyle=-2|\beta_{2}|^{2}\mathrm{\
Im}[k_{x}\phi_{1,k_{x}}\cdot\partial_{z}{\phi}_{1,k_{x}}^{\ast}],$

where $\hat{\phi}_{k_{x}}$ is the Fourier transform of the streamfunction at
wavenumber $k_{x}$ and $\phi_{j,k_{x}}$ represents the $z$-dependent $j$-th
complex eigenmode: $j=1,2$ for unstable and stable modes, respectively. The
conjugate symmetry of unstable and stable modes, as shown in Figs. 1(b) and
1(c), is used in Eq. (15), i.e., $\phi_{2,k_{x}}(z)=\phi_{1,k_{x}}^{\ast}(z)$.
The negative sign of the last expression in Eq. (15) corresponds to the up-
gradient nature of momentum transport by stable modes, which was physically
analyzed in Sec. III.2 and visually demonstrated in Fig. 2.

The summed contributions of unstable and stable modes in transport, however,
can have cross-terms—quadratic correlations between unstable and stable
modes—that do not appear in Eqs. (14) and (15) where contribution from
individual modes are shown. But the cross-terms vanish when the unstable and
stable modes are exactly complex conjugates of each other, as is the case for
the ideal Kelvin-Helmholtz instability  (when this conjugate symmetry is
broken, e.g., in resistive tearing instability or in ion-temperature-gradient
instability,[60] the cross-terms can have non-zero contribution):

$\displaystyle\tau_{u}(\mathrm{unstable\ +\ stable})$
$\displaystyle=2\mathrm{\
Im}[k_{x}(\beta_{1}\phi_{1,k_{x}}+\beta_{2}\phi_{2,k_{x}})\cdot\partial_{z}(\beta_{1}{\phi}_{1,k_{x}}+\beta_{2}\phi_{2,k_{x}})^{\ast}]$
(16) $\displaystyle=2\left(|\beta_{1}|^{2}-|\beta_{2}|^{2}\right)\mathrm{\
Im}[k_{x}{\phi}_{1,k_{x}}\cdot\partial_{z}{\phi}_{1,k_{x}}^{\ast}]+\mathrm{cross\textnormal{-}terms},$

where

$\displaystyle\mathrm{cross\textnormal{-}terms}$ $\displaystyle=2\mathrm{\
Im}[k_{x}(\beta_{1}\phi_{1,k_{x}})\cdot\partial_{z}(\beta_{2}\phi_{2,k_{x}})^{\ast}]+2\mathrm{\
Im}[k_{x}(\beta_{2}\phi_{2,k_{x}})\cdot\partial_{z}(\beta_{1}\phi_{1,k_{x}})^{\ast}]$
(17) $\displaystyle=2\mathrm{\
Im}[\beta_{1}\beta_{2}^{\ast}k_{x}\phi_{1,k_{x}}\cdot\partial_{z}\phi_{2,k_{x}}^{\ast}]+2\mathrm{\
Im}[\beta_{2}\beta_{1}^{\ast}k_{x}\phi_{2,k_{x}}\cdot\partial_{z}\phi_{1,k_{x}}^{\ast}]$
$\displaystyle=2\mathrm{\
Im}[\beta_{1}\beta_{2}^{\ast}k_{x}\phi_{1,k_{x}}\cdot\partial_{z}\phi_{2,k_{x}}^{\ast}]+2\mathrm{\
Im}[\beta_{2}\beta_{1}^{\ast}k_{x}\phi_{1,k_{x}}^{\ast}\cdot\partial_{z}\phi_{2,k_{x}}]$
$\displaystyle=2\mathrm{\
Im}[\beta_{1}\beta_{2}^{\ast}k_{x}\phi_{1,k_{x}}\cdot\partial_{z}\phi_{2,k_{x}}^{\ast}]+2\mathrm{\
Im}[\left(\beta_{1}\beta_{2}^{\ast}k_{x}\phi_{1,k_{x}}\cdot\partial_{z}\phi_{2,k_{x}}^{\ast}\right)^{\ast}]$
$\displaystyle=0.$

Thus we obtain

$\displaystyle\tau_{u}(\mathrm{unstable\ +\ stable})$
$\displaystyle=2\left(|\beta_{1}|^{2}-|\beta_{2}|^{2}\right)\mathrm{\
Im}[k_{x}{\phi}_{1,k_{x}}\cdot\partial_{z}{\phi}_{1,k_{x}}^{\ast}]$ (18)
$\displaystyle=\tau_{u}(\mathrm{unstable})+\tau_{u}(\mathrm{stable}).$

These relations inform us about the $z$-profile of the Reynolds stress,
contributed by each wavenumber and each eigenmode. As largest momentum
transport happens in the region with the largest flow-gradient, it is
instructive to compute, in the forced shear layers, the turbulent stresses at
the middle of the layer at $z=0$, and compare the stress contributions from
different eigenmodes at various wavenumbers.

The total Reynolds stress from all modes and all wavenumbers in the
simulations is compared in Fig. 9 with the stress contributions from the
wavenumber range $0<|k_{x}|<1$, which is decomposed further into eigenmodes to
assess the contribution of the unstable modes, stable modes, and their sum.
The subplots demonstrate that the stable modes are highly efficient in
transporting momentum in the up-gradient direction, as compared to the down-
gradient transport by the unstable modes. Even for the strongest magnetic
field $M_{\mathrm{A}}=3$, close to the instability threshold, the stable modes
contribute significantly to a continuous reduction of the turbulent momentum
flux. In addition, the occasional breakthroughs in stable-mode activity cause
reversals of the transport direction. This reversal can be observed when the
total Reynolds stress in the system is computed, without decomposing the
stress into contributions by each eigenmode. However, when the stable modes
are not overtaking the unstable modes in transport, the resulting down-
gradient transport observed in simulations or experiments is difficult to
interpret, in regards to the contributions of stable modes in subdominantly
reducing the transport. An eigenmode decomposition of turbulent fluctuations,
however, uncovers a complete picture, as is shown here.

Figure 9: Time variations of Reynolds stress at the middle of the shear layer,
$z=0$. The stress contributions from unstable modes (blue), stable modes
(orange), their sum (red), and full nonlinear fluctuations, i.e., all modes
and all wavenumbers in the simulations (black), are compared, for varying
strengths of magnetic fields: (a) $M_{\mathrm{A}}=60$, (b)
$M_{\mathrm{A}}=10$, and (c) $M_{\mathrm{A}}=3$. Thin green lines represent
the zero level. All simulations use $D_{\mathrm{Krook}}=2$. Although, with
stronger magnetic fields, the up-gradient momentum transport by stable modes
are reduced, the up- and down-gradient transport nearly cancel each other
throughout all cases.

Similar variations of momentum transport across the middle of the shear layer
are compared in Fig. 10 for different forcing strengths. Note the unforced
case differs from the forced cases, as the nearly-flattened shear layer has
less momentum to be transported across the layer. As reported in Ref. [44],
despite the profile relaxation, the two eigenmodes per wavenumber describe
well the temporal variation of the Reynolds stress across the shear layer,
although the stress itself is very low (note its vertical scale). In all
cases, the stress captured via the sum of unstable and stable modes almost
completely follows the total stress from all modes.

Figure 10: Time variations of Reynolds stress at the middle of the shear
layer, $z=0$. The stress contributions from unstable modes (blue), stable
modes (orange), their sum (red), and full nonlinear fluctuations, i.e., all
modes and all wavenumbers in the simulations (black), are compared, for
varying forcing strengths: (a) $D_{\mathrm{Krook}}=25$, (b)
$D_{\mathrm{Krook}}=1$, and (c) $D_{\mathrm{Krook}}=0$. All simulations use
$M_{\mathrm{A}}=10$. Thin green lines represent the zero level. Qualitative
differences can be observed in unforced ($D_{\mathrm{Krook}}=0$) and forced
cases ($D_{\mathrm{Krook}}\neq 0$): as instability extracts energy from the
mean flow, the profile relaxation in the unforced layer leads to a decaying
turbulence, and the transport rates become very small [note the vertical axis
labels in (c)]. However, in all cases, the summed stable modes producing up-
gradient transport nearly cancel the down-gradient transport by unstable
modes. The addition of these two contributions produces a stress that is
almost identical to the stress from all modes.  Figure 11: Time variations of
Reynolds stress at the middle of the shear layer, $z=0$. The stress
contributions from unstable modes (blue), stable modes (orange), their sum
(red), and full nonlinear fluctuations, i.e., all modes and all wavenumbers in
the simulations (black), are compared, for varying magnetic Prandtl numbers
(resistivities): (a) $Pm=0.1$ and (b) $Pm=10$. All simulations use
$M_{\mathrm{A}}=10$, $D_{\mathrm{Krook}}=2$, and $Re=500$. Thin green lines
represent the zero level. It can be observed that the stable modes begin
driving up-gradient momentum transport at around $t\approx 30$ when the
nonlinear phase of evolution begins. By varying $Pm$ by two orders of
magnitude, around unity, the stable modes are found to substantially reduce
the down-gradient transport; note the case of $Pm=1$ is shown in Fig. 9(b).

In Fig. 11, the momentum transport by the unstable and stable modes is
presented as a function of magnetic Prandtl number $Pm=Rm/Re$. All simulations
until this point used $Rm=500$, which is now changed to $Rm=50$ and
$Rm=5\,000$. In both cases of $Pm=0.1$ and $Pm=10$, the stable modes still
substantially offset the turbulent momentum transport of the unstable modes.
The shorter time trace for $Pm=10$ is due to the higher simulation cost. It
should be noted that the quasi-stationary state in this simulation is still
undergoing changes, unlike in the case of $Pm=0.1$ in Fig. 11(a) or $Pm=1$ in
Fig. 9(b), all with the same $M_{\mathrm{A}}=10$, $D_{\mathrm{Krook}}=2$ and
$Re=500$.

The efficiency of time-averaged up-gradient momentum transport due to stable
modes is compared in Fig. 12 with the time-averaged down-gradient transport
due to unstable modes, via a measure, defined below:

$\mathrm{Transport\ reduction\
efficiency}=\frac{\langle\mathrm{Up\text{-}gradient\ transport\ by\ stable\
modes}\rangle_{t}}{\langle\mathrm{Down\text{-}gradient\ transport\ by\
unstable\ modes}\rangle_{t}},$ (19)

where $\langle A\rangle_{t}$ represents a time-averaging operation on $A$.

Figure 12: Parameter dependence of transport reduction efficiency, which is
defined as the ratio of time-averaged up-gradient Reynolds stress due to
stable modes and time-averaged down-gradient Reynolds stress due to unstable
modes. The stress is measured at the middle of the shear layer, $z=0$, where
the momentum transport is at its maximum. (a) Variations in
$M_{\mathrm{A}}=3,10,30,60,120$ with $D_{\mathrm{Krook}}=2$, $Pm=1$, and
linear $x$-scale. (b) Variations in $D_{\mathrm{Krook}}=0.1,1,2,6$ with
$M_{\mathrm{A}}=10$, $Pm=1$, and logarithmic $x$-scale. (c) Variations in
$\mathrm{Pm}=0.1,1,10$ (or, equivalent changes in resistivities) with
$M_{\mathrm{A}}=10$, $D_{\mathrm{Krook}}=2$, and logarithmic $x$-scale. All
plots have the same $y$-axis. The time-average for (a) and (b) is taken over a
long quasi-stationary state of turbulence $t=350\text{--}900$, while for (c),
it is $t=137\text{--}237$ where the quasi-stationary state is still undergoing
changes. In all cases, substantial reduction of transport by stable modes is
evident, which cancel, via their up-gradient transport, more than half of the
down-gradient transport by unstable modes, and this fraction reaches up to
$98\%$, see (a), for $M_{\mathrm{A}}=60$ and $M_{\mathrm{A}}=120$.

Variations in magnetic field strength, forcing strength, and magnetic Prandtl
number all demonstrate that the stable modes cancel an appreciable amount of
the turbulent momentum flux associated with the unstable modes. On average,
around $80\%$ of the down-gradient flux is offset in this manner.

A remark should be made now regarding the use of unstable and stable modes for
building a reliable reduced mode of transport for geo- and astro-physical
problems. One approach would be to relate the activity of these two modes with
a coefficient of diffusive flux (although the unstable and stable modes offer
spatial profiles of transport as well, with both diffusive and non-diffusive
fluxes, because they do not rely on an ad-hoc eddy-viscosity model, which is
an explicit diffusive-flux-based model). In the middle of the shear layer, the
diffusive flux, however, dominates because of the maximum in the flow-
gradient. The ad-hoc turbulent viscosity can thus be defined [27], more
importantly without a “free-parameter,” using Eq. (18) as

$\nu_{\mathrm{turb}}=\frac{\tau_{u}(\mathrm{unstable\ +\
stable})}{(dU_{0}/dz)\Big{|}_{z=0}}=\sum\limits_{0<k_{x}<1}2\left(|\beta_{1}|^{2}-|\beta_{2}|^{2}\right)\mathrm{\
Im}[k_{x}{\phi}_{1,k_{x}}\cdot\partial_{z}{\phi}_{1,k_{x}}^{\ast}]\Big{|}_{z=0}.$
(20)

Note that the denominator is unity for the shear-flow that has a linear
profile in the vicinity of $z=0$. To assess the importance of stable modes in
this construct, $|\beta_{1}|^{2}-|\beta_{2}|^{2}$ can be written as
$|\beta_{1}|^{2}\left(1-|\beta_{2}|^{2}/|\beta_{1}|^{2}\right)$. Since
$|\beta_{2}|^{2}$ has been found to on the same order of $|\beta_{1}|^{2}$,
e.g., see Fig. 12, where $|\beta_{2}|^{2}/|\beta_{1}|^{2}$ can range from
$\approx 0.8$ to $\approx 0.95$, yielding
$\left(1-|\beta_{2}|^{2}/|\beta_{1}|^{2}\right)\approx 0.05\mathrm{\ to\
}0.2$. Therefore, neglecting stable modes can overestimate the transport by a
factor of $5$ to $20$.

###  V.5 Reynolds vs. Maxwell stresses

With the above successful low-order representation of Reynolds stress above,
we now examine the fluctuations in the magnetic field that give rise of
Maxwell stress. The stress can be quantified as

$\tau_{b}(\mathrm{all\ modes})=-\frac{2}{M_{\mathrm{A}}^{2}}\mathrm{\
Im}[k_{x}\hat{\psi}_{k_{x}}\cdot\partial_{z}\hat{\psi}_{k_{x}}^{\ast}],$ (21)

$\tau_{b}(\mathrm{unstable\ +\
stable})=-\frac{2}{M_{\mathrm{A}}^{2}}\left(|\beta_{1}|^{2}-\beta_{2}|^{2}\right)\mathrm{\
Im}[k_{x}{\psi}_{1,k_{x}}\cdot\partial_{z}{\psi}_{1,k_{x}}^{\ast}],$ (22)

where $\hat{\psi}_{k_{x}}$ is the Fourier transform of the flux function at a
wavenumber $k_{x}$ and $\psi_{j,k_{x}}$ represents the $z$-dependent $j$-th
complex eigenmode ($j=1,2$ for unstable and stable modes, respectively).
Again, cross-terms arising from the correlation between the unstable and modes
can be shown to vanish, exactly as it was shown for the Reynolds stress in Eq.
(17).

Figure 13: Time-averaged turbulent stresses split into their contributions
from different wavenumbers. (a) Stresses on a linear scale. (b) Log-log
representation of the absolute value of the stresses. Note that only the
wavenumbers $k_{x}<1$ are Kelvin-Helmholtz-unstable. The simulation parameters
are $M_{\mathrm{A}}=60$ and $D_{\mathrm{Krook}}=2$; the time-average is taken
over a quasi-stationary state of turbulence, $t=350\text{--}1000$. The total
turbulent stress is dominated by the range $|k_{x}|<1$, which is captured by
the unstable and stable modes at those wavenumbers to a high precision. The
small amount of stresses that are contributed by smaller scales of
fluctuations span a broad range of wavenumbers, due to the smaller scales in
magnetic fields generated via straining by the flow.

As can be seen in Fig. 13, the Reynolds stress is dominated by large scales
while the Maxwell stress involves a large number of different scales. Figure
13(a), using axes with linear scales, shows the dominance of Reynolds stress
in the entire system, which the two-eigenmodes-per-wavenumber decomposition
(unstable and stable modes) captures, not only qualitatively, but also
quantitatively with great accuracy. In Fig. 13(b), a logarithmic scale is used
to expose the range of small scales that contribute significantly to the
magnetic fluctuations. Wavenumbers $k_{x}\lesssim 10$ have major
contributions, as opposed to $k_{x}<1$ for the fluctuations of the flow. The
fact that a large amount of flow energy resides at large scales suggests that
the shear-flow turbulence may be amenable to some form of quasilinear
modeling. Homogeneous isotropic turbulence, on the other hand, would not be
reliably captured with such models, as no scale separation exists therein.
Recent studies have highlighted that improved quasilinear models such as the
generalized quasilinear approximation are realizable in systems with length-
or time-scale separation. [61]

The magnetic fluctuations, on the other hand, span a broad range of scales.
This can be physically interpreted as a result of the straining of the
magnetic fields by the turbulent flow, which generates small scales in the
magnetic fields.[62, 63, 64, 48] The straining process by the large-scale
turbulent eddies converts the large-scale kinetic energy into the
intermediate-scale magnetic energy. [65] Magnetic fluctuations at such scales
can then, via Lorentz force, feed back on the flow, although mostly at smaller
scales. A comprehensive analysis of energy transfer for the present system
will be reported in a forthcoming publication where nonlinear mode-coupling
and energy transfer between fluctuations of discrete and continuum modes of
velocity and magnetic fields are also analyzed.

To model any aspect of magnetic fluctuations, one must thus rely on tools such
as statistical theories to obtain scaling laws that can offer insights into
these fluctuations. One such approach is detailed next.

###  V.6 Scaling law for continuum modes

Until this point, the discrete modes—unstable and stable modes—which describe
the turbulent flow well, have been our focus. The magnetic fluctuations, on
the other hand, result from the straining of field lines by the flow, exciting
the remaining continuum modes. Hence these modes are necessary for a
successful reconstruction of the magnetic fluctuations. Thus, we seek a simple
scaling law for the saturated turbulent amplitudes of the continuum modes.

####  V.6.1 Analytical prediction for continuum mode amplitudes

It is instructive to write the nonlinear MHD equations in the eigenmode basis,
arriving at what is also referred to as the mode-amplitude evolution
equation[30, 43, 45, 67, 66]

$\partial_{t}\beta_{j}(k_{x})=i\omega_{j}(k_{x})\beta_{j}(k_{x})+\sum_{\begin{subarray}{c}k_{x}^{\prime},k_{x}^{\prime\prime},m,n\\\
k_{x}^{\prime}+k_{x}^{\prime\prime}=k_{x}\end{subarray}}C_{jmn}(k_{x},k_{x}^{\prime})\beta_{m}(k_{x}^{\prime})\beta_{n}(k_{x}^{\prime\prime}),$
(23)

where $\beta_{j}(k_{x})$ represents the complex amplitude of an eigenmode $j$
at wavenumber $k_{x}$ with $\omega_{j}$ the associated mode-frequency; the
nonlinear mode coupling coefficient $C_{jmn}$ measures the three-wave overlap,
which dictates the strength of nonlinear beating between an eigenmode $m$ at
wavenumber $k_{x}^{\prime}$ and an eigenmode $n$ at wavenumber
$k_{x}^{\prime\prime}$, driving an eigenmode $j$ at wavenumber $k_{x}$ (with
the constraint $k_{x}^{\prime}+k_{x}^{\prime\prime}=k_{x}$).

For the continuum modes, as was mentioned in Sec. III.1, their frequencies
depend linearly on the wavenumber $k_{x}$ as [59]
$\omega/k_{x}+U_{\mathrm{ref}}(z)\pm v_{\mathrm{A,ref}}(z)=0$. This implies
$\omega_{j}=\omega\propto k_{x}$.

Heuristically, the scaling of the nonlinear mode coupling coefficient with
wavenumber can be obtained in the following manner: In Eqs. (2a) and (2b), the
separation of linear and nonlinear terms arises in Poisson brackets. Consider
a prototype equation,

$\displaystyle\partial_{t}\widetilde{P}$ $\displaystyle=\\{P,Q\\}+...$ (24)
$\displaystyle=\\{\widetilde{P},\widetilde{Q}\\}+\\{\widetilde{P},Q_{\mathrm{0}}\\}+\\{P_{\mathrm{0}},\widetilde{Q}\\}+...,$

where $P=P_{0}+\widetilde{P}$ and $Q=Q_{0}+\widetilde{Q}$ represent two fields
(e.g., $\nabla^{2}\phi$ or $\psi$ for the present problem), with $P_{0}$
representing the $x$-averaged mean component of $P$, and $\widetilde{P}$
standing for perturbations. The linear term, e.g.,
$ik_{x}\hat{P}\cdot\partial_{z}Q_{\mathrm{0}}$ which is in spectral space,
contains only one perturbed field, whereas the nonlinear term, e.g.,
$ik_{x}^{\prime}\hat{P}^{\prime}\cdot\partial_{z}\hat{Q}^{\prime\prime}$, has
two perturbed fields, with $\hat{P}^{\prime}$ and $\hat{Q}^{\prime\prime}$
representing the Fourier-transformed quantities at wavenumbers
$k_{x}^{\prime}$ and $k_{x}^{\prime\prime}$, respectively. It may be supposed
that the derivative $\partial_{z}$ on the perturbed quantities is roughly on
the scale of $|\partial_{z}|\sim k_{x}$. (This can be shown analytically for
all the eigenmodes, where the background flow is approximately uniform, see
Ref.[43].) Notice, however, that this argument applies only to the
perturbations: the operator $\partial_{z}$ acting on $Q_{0}$ clearly does not
produce a factor of $k_{x}$, which is zero for the mean component $Q_{0}$. We
now use this distinction to make a prediction for the amplitudes of
perturbations, in particular the continuum mode-amplitudes. The linear and
nonlinear terms thus assume the forms
$ik_{x}\hat{P}\cdot\partial_{z}Q_{\mathrm{0}}$ and
$ik_{x}^{\prime}\hat{P}^{\prime}\cdot
ik_{x}^{\prime\prime}\hat{Q}^{\prime\prime}$, respectively.

In Eq. 24, expanding the perturbations in the eigenmode basis, e.g.,
$\hat{P}=\sum_{l}\beta_{l}\hat{P}_{l}$ with $\hat{P}_{l}$ representing the
$l$-th eigenmode, and diagonalizing the linear terms (operator), one finds an
equation of the form given in Eq. 23. We can now attempt to understand the
behavior of the nonlinear mode coupling coefficients that drive the continuum
modes. Assuming nonlinear interactions between the continuum modes are local
in spectral space—interaction of three wavenumbers of similar scales—the
nonlinear term in Eq. 24 simplifies, e.g.,
$ik_{x}^{\prime}\hat{P}^{\prime}\cdot
ik_{x}^{\prime\prime}\hat{Q}^{\prime\prime}$ becomes
$-k_{x}^{2}\hat{P}\hat{Q}$; note the linear term has the form
$ik_{x}\hat{P}\cdot\partial_{z}Q_{\mathrm{0}}$.

In assuming local interaction between the continuum modes spectral space in
$k_{x}$, the involvement of unstable and stable modes in nonlinear
interactions is ignored, which otherwise could bring in non-local effects.
This may be a valid assumption for continuum modes at scales much above the
Kelvin-Helmholtz-unstable wavenumber range, i.e., $k_{x}>1$, as the wavenumber
convolution constraint of $k_{x}=k_{x}^{\prime}+k_{x}^{\prime\prime}$ does not
allow two (un-)stable modes to beat together to drive a continuum mode at
large $k_{x}$, e.g., $k_{x}>2$.

Continuing with the above assumption, the nonlinear term
$-k_{x}^{2}\hat{P}\hat{Q}$ has one extra $k_{x}$ compared to the linear term
$ik_{x}\hat{P}\cdot\partial_{z}Q_{\mathrm{0}}$. This implies that, for the
continuum modes, the nonlinear mode coupling coefficients $C$ are expected to
scale as

$C\propto k_{x}^{2},$ (25)

because the linear term for the continuum modes in Eq. 23 has the
eigenfrequency that depends linearly on $k_{x}$, i.e.,

$\omega\propto k_{x}.$ (26)

Such a property of nonlinear coupling coefficient is common in other
turbulence calculations, as well. [26]

In order to obtain a phenomenological scaling law, we now make no distinction
between different continuum modes, and thus balance the linear and nonlinear
terms of Eq. 23 in the quasi-stationary state as $\omega\beta\sim C\beta^{2}$.
Inserting their asymptotic dependences on $k_{x}$, the amplitudes of continuum
modes is found to follow

$\beta\sim k_{x}^{-1}.$ (27)

Note that the assumptions made in arriving at this simple scaling law are
crude. The next step will be to determine from nonlinear simulations whether
this scaling can be recovered or whether a number of assumptions made above
render the result inapplicable.

####  V.6.2 Numerical verification

Figure 14: (a) Dependence of time-averaged eigenmode amplitudes on $k_{x}$.
The indices $j$s of the eigenmodes are arranged in increasing order of their
real frequencies. Vertical lines signify the self-similar cascade of energy to
small scales in $k_{x}$. Mode amplitudes are averaged over a quasi-stationary
state of turbulence $t=400\text{--}800$ for a simulation with
$M_{\mathrm{A}}=120$ and $D_{\mathrm{Krook}}=2$. (b) $k_{x}$ spectra of mode
amplitudes in a nonlinear simulation (shown with empty circles) in a
$\log_{10}\mathrm{-}\log_{10}$ plot. Shown with a solid line is the analytical
prediction, made for the wavenumbers that lie beyond the Kelvin-Helmholtz-
unstable range, i.e., for $k_{x}>1$. The spectral index, predicted based upon
a number of simple assumptions, can be seen to fit the data reasonably well.

Time-averaged eigenmode amplitudes from nonlinear simulations, after
multiplying with $k_{x}$, are plotted in Fig. 14(a) as functions of $k_{x}$
and eigenmode index $j$, arranged in order of increasing real frequency of the
eigenmodes. The appearance of vertical near-equicontour lines signifies that
eigenmodes are excited in a similar pattern across a large range of scales.

The eigenmodes that lie within the yellow bands are localized in space
($z$-axis), but the band spans a range of heights, outside the shear layer
$|z|<1$. Empirically, we note that the center of the rightmost [leftmost] band
corresponds to $\omega/k_{x}=U_{00}$ [$\omega/k_{x}=-U_{00}$] where
$U_{00}=1$. These thick bands represent all eigenmodes that have phase speeds
$\omega/k_{x}=U_{00}+cv_{\mathrm{A},0}$
[$\omega/k_{x}=-U_{00}+cv_{\mathrm{A},0}$] where $-1<c<1$ and
$v_{\mathrm{A},0}=1/M_{\mathrm{A}}$; note that $|c|=1$ is not included in
these bands. All of these eigenmodes have peaks and oscillations in their
eigenfunctions outside of the shear layer. In the layer, the unstable and
stable modes maintain their dominance and thus these two discrete modes alone
almost completely regulate the momentum transport across the layer, as was
noted in Sec. V.4.

It is of interest now to compute from numerical simulation data how the
amplitude of each eigenmode $j$ falls off with $k_{x}$ and construct a
$j$-averaged spectral index. To this end, we note the amplitude
$\beta_{j}(k_{x}=0.2)$ for each mode $j$ at $k_{x}=0.2$ (the first wavenumber
in the simulation) and compute a scaled mode-amplitude
$\hat{\beta}_{j}(k_{x})$ as

$\hat{\beta}_{j}(k_{x})=\frac{\beta_{j}(k_{x})}{\beta_{j}{(k_{x}\mathrm{=}0.2)}},$
(28)

which is expected to fall-off with $k_{x}$ as $\sim k_{x}^{\alpha}$. In
principle, the spectral index $\alpha$ can depend on the eigenmode index $j$,
but a $j$-averaged spectral index is sought now, following the procedure

$\displaystyle\hat{\beta}_{j}(k_{x})$
$\displaystyle\propto(k_{x}/0.2)^{\alpha},$ (29)
$\displaystyle\mathrm{log\,}\hat{\beta}_{j}(k_{x})$
$\displaystyle\propto\alpha\,\mathrm{log\,}k_{x},$
$\displaystyle\langle\mathrm{log\,}\hat{\beta}_{j}(k_{x})\rangle_{j}$
$\displaystyle\propto\langle\alpha\rangle_{j}\,\mathrm{log\,}k_{x}.$

This $j$-averaged spectral distribution of the mode-amplitudes informs how, on
average, each eigenmode amplitude depend on $k_{x}$.

A plot of $\langle\mathrm{log\,}\hat{\beta}_{j}(k_{x})\rangle_{j}$ vs.
$\mathrm{log\,}k_{x}$ is shown in Fig. 14(b), along with the analytical
prediction of inverse-in-wavenumber fall-off of the mode amplitudes, at scales
above the Kelvin-Helmholtz-unstable wavenumber range. It should be highlighted
that the computation of all the eigenmode amplitudes at each wavenumber at
each simulation time is computationally demanding, as the process requires the
computation of modified left eigenmodes for each right eigenmode at each
wavenumber, apart from the mode projection calculation at each time step.
Therefore, only the first $24$ Fourier modes in $k_{x}$ are shown in Fig. 14.

A finding in Fig. 14 is the identification of self-similar cascade of energy
to smaller scales (larger $k_{x}$) in eigenmode space. This result also hints
that the interaction involving the continuum modes may be reasonably
simplified, and potentially be valuable in estimating the amplitudes of
unstable and stable modes, using closure theories (see Ref.[26] for a recent
example).

##  VI Conclusions

We have investigated MHD turbulence in two-dimensions, driven by a forced
unstable shear flow, using a complete eigenmode decomposition of fluctuations
in nonlinear simulations, which exposes the nonlinearly saturated excitation
level of each eigenmode and its role. Intrinsic to linear instability, the
unstable modes derive fluctuation energy from the mean flow gradient. The
linearly-decaying stable modes, however, contain almost the same amount of
energy as the unstable mode, which they receive via nonlinear excitation. This
truncated basis of two eigenmodes per wavenumber is found to reconstruct
essential large-scale features of turbulent flow and the associated momentum
transport via Reynolds stress. Quantifying transport due to unstable modes
alone shows an overestimation up to an order of magnitude higher relative to
the actual flux. The reduction in the flux is identified to be due to the
continuous up-gradient transport by the stable modes, which causes a near-
cancellation of down-gradient transport driven by unstable modes.

The continuum modes, on the other hand, describe small-scale fluctuations of
the flow and magnetic field, where the above large-scale unstable and stable
modes manifest themselves as a quasi-coherent vortex. To predict the mode
amplitudes of the continuum, a simple scaling law is derived from the
governing nonlinear equations and the predicted inverse-in-wavenumber fall-off
rate is found reasonably agree with the simulation data.

Although both the momentum transport and fluctuation energy are largely
described by the discrete modes, the former is more efficiently captured
[Figs. 9–13] as almost all the momentum transport occurs near $z=0$, which is
the region where the discrete modes dominate [Figs. 1(b) and 1(c)]. The
fluctuation energy, on the other hand, is related to fluctuations that are
scattered in and around the large-scale eddies; a portion of this energy is
claimed by the continuum modes, although a large fraction still belongs to the
discrete modes [Figs. 5(b), 6–8].

Transport reduction by stable modes can also be used to improve
phenomenological constructs like eddy viscosity, which are generally agnostic
as to the nonlinear excitation of stable modes. By predicting the turbulent
amplitudes of the unstable and stable modes, for astrophysically relevant
parameters, e.g., very large $Rm$, $Re$, $Pm$ compared to unity, the simple
relation between turbulent viscosity and eigenmode amplitudes [in Eq.(20)] can
be exploited to reliably model transport processes in astrophysical objects,
which otherwise cannot be solved using current state-of-the-art direct
numerical simulations. It should be noted that such a prediction for the mode
amplitudes $|\beta_{1}|$ and $|\beta_{2}|$ was recently made for ion-
temperature-gradient-driven fusion plasma turbulence [17] using statistical
closure theory. [68] Undertaking such a task for the present system is
interesting, but beyond the purview of this work and will thus be left for
future investigations.

The reduced representation of turbulent flow and transport presented here is
also useful for building sub-grid-scale models, which can allow performing
nonlinear simulations at extreme parameters with less-intensive computational
demands. Progress can thus be made in seeking models that reduce the number of
degrees of freedom while capturing essential features of the turbulent system.
Techniques like proper orthogonal decomposition, dynamic mode decomposition,
etc., also exist for such purposes, [69] but they operate on output from
nonlinear simulations, and it can be difficult to assign intuitive physical
meaning to the characteristic mode structures. Here, the truncated eigenmode
basis, composed of the unstable and stable modes, has been demonstrated to
reconstruct nonlinear fluctuations to an appreciable degree, thus suggesting
that they can be leveraged as a physically-motivated basis for extreme
parameter studies, without having to first perform a direct numerical
simulation. These modes may also be useful for generating, via their nonlinear
interactions with continuum modes, fluctuations associated with the continuum
modes. Such a test could be performed to analyze magnetic fluctuations. The
reduced basis, composed of unstable and stable mode alone, can also serve in
direct statistical simulations, [54] which have shown promises towards
simulating the slowly-evolving turbulent statistics, e.g., two-point two-time
correlations, three-point correlations between fluctuating fields, etc.,
rather than the fast-evolving field variables themselves, e.g., flow
velocities. Other improved forms of quasilinear models like the generalized
quasilinear approximation[61] may also benefit from using this truncated
basis. This possibility will be explored in a separate publication.

In the future, procedures similar to that employed here can be used to examine
the properties of other forms of instability-driven turbulence, such as
magneto-rotational-instability-driven [6] and stratified-shear-flow-driven
turbulence.[70] Building nonlinear energy transfer diagnostics in shear-flow
turbulence to study the physical processes and scales that impact the
difference in unstable- and stable-mode amplitudes is another possible avenue.
Such investigations constitute steps towards deployment in service to one-
dimensional stellar transport models. [5] Central to improved predictiveness
are the stable modes, whose properties will similarly require additional
studies.

###### Acknowledgements.

The authors appreciate K. Burns, E.H. Anders, and the Dedalus developers for
assistance in several aspects of leveraging the numerical code. Thanks also to
J. Fuller and other participants of the program “Transport in Stellar
Interiors, 2021” at the Kavli Institute of Theoretical Physics for useful
discussions. This material is based upon work funded by the Department of
Energy [DE-SC0022257] through the NSF/DOE Partnership in Basic Plasma Science
and Engineering. We also gratefully acknowledge support from NSF Grant Nos.
AST-1814327 and AST-1908338. The simulations reported herein were performed
using the XSEDE supercomputing resources via Allocation No. TG-PHY130027. The
data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Appendix A: Orthogonality of right eigenmode and modified left eigenmode

Due to the non-normality of the linear operator of the shear-flow instability,
the eigenmodes are not orthogonal. This presents a significant challenge in
the computation of mode amplitudes. An additional challenge is the generalized
eigenvalue nature of the problem at hand, when written in vorticity formalism,
as in Eq. (6a). This differs from the standard eigenvalue problem,
$L\xi_{j}=\lambda_{j}\xi_{j}$, where $L$ is a linear operator whose $j$-th
eigenmode is $\xi_{j}$ with eigenvalue $\lambda_{j}$. The generalized
eigenvalue problem that we encounter here is

$LX_{j}=\omega_{j}MX_{j},$ (A1)

where $L$ is a linear operator, $M$ is another linear operator containing the
Laplacian operation for our problem, $X_{j}$ is the $j$-th (right) eigenmode
with corresponding (right) eigenvalue $\omega_{j}$. [Often times, the
distinction between left and right eigenmodes of a linear operator is not made
as they happen to be the same; however this is not the case here for the non-
normal operator.] The right eigenmodes, although non-orthogonal to each other,
can be made orthogonal with an appropriate weight factor to the left
eigenmodes, which are solutions to another eigenvalue problem: [66, 67]

$Y_{j}^{\mathrm{T}}L=\sigma_{j}Y_{j}^{\mathrm{T}}M.$ (A2)

Here, $Y_{j}^{\mathrm{T}}$ is the transpose of the left eigenmode with its
left eigenvalue $\sigma_{j}$. A slight reformulation is possible to this
equation by taking Hermitian-transpose:

$L^{\dagger}Y_{j}^{\ast}=\sigma_{j}^{\ast}M^{\dagger}Y_{j}^{\ast}.$ (A3)

In the eigenvalue solver in Dedalus, the matrices $L^{\dagger}$ and
$M^{\dagger}$ for each wavenumber are passed, and their eigenmodes
$Y_{j}^{\ast}$ and eigenvalues $\sigma_{j}^{\ast}$ are found. It can be shown
that the eigenvalues $\sigma_{j}$ and $\omega_{j}$ are the same (i.e.,
$\sigma_{j}=\omega_{j}$), by analyzing Eqs. (A1) and (A2). A modified
orthogonality relation between the left and right eigenmodes can now be
derived: [66, 67]

$\displaystyle(Y_{j}^{\mathrm{T}}L)X_{i}=Y_{j}^{\mathrm{T}}(LX_{i})=Y_{j}^{\mathrm{T}}(\omega_{i}MX_{i})$
(A4) $\displaystyle\implies$
$\displaystyle(\sigma_{j}Y_{j}^{\mathrm{T}}M)X_{i}=Y_{j}^{\mathrm{T}}(\omega_{i}MX_{i})$
$\displaystyle\implies$ $\displaystyle
Y_{j}^{\mathrm{T}}MX_{i}(\sigma_{j}-\omega_{i})=0$ $\displaystyle\implies$
$\displaystyle Y_{j}^{\mathrm{T}}MX_{i}\propto\delta_{i,j},$

which means the left and right eigenmodes are orthogonal to each other with a
weight factor $M$, as long as their eigenvalues differ
($\sigma_{j}\neq\omega_{i}$). For numerical computation, it is convenient to
define
$Y_{j}^{\mathrm{T}}M=(M^{\mathrm{T}}Y_{j})^{\mathrm{T}}=Z_{j}^{\mathrm{T}}$
where $Z_{j}$ is the modified left eigenmode, which is—by construct—orthogonal
to the right eigenmode without any weight factor:
$Z_{j}^{\mathrm{T}}X_{i}\propto\delta_{i,j}$. Using this relation the
eigenmode coefficients $\beta_{j}({k_{x},t})$ in the eigenmode expansion of
turbulent fluctuations are computed at each wavenumber and at each time.

## Appendix B: Cyclic oscillations in mode-amplitudes for weak magnetic fields

Figure B1: Same as in Fig. 5, but for $M_{\mathrm{A}}=120$. The amplitude of
the nonlinearly excited stable mode is almost exactly the same as that of the
unstable mode. The nature of their oscillations are also similar, although the
oscillations in the stable-mode-amplitude lags behind that of the unstable
mode. The lag is likely an outcome of a time-delay in energy transfer from the
unstable to the stable mode at the same wavenumber, which thus requires a
series of nonlinear interactions with fluctuations at other wavenumbers.

## References

  * [1] E.C. Harding, J.F. Hansen, O.A. Hurricane, R.P. Drake, H.F. Robey, C.C. Kuranz, B.A. Remington, M.J. Bono, M.J. Grosskopf, and R.S. Gillespie, Observation of a Kelvin-Helmholtz Instability in a High-Energy-Density Plasma on the Omega Laser, Phys. Rev. Lett.  103, 045005 (2009).
  * [2] H. Hasegawa, M. Fujimoto, T.-D. Phan, H. Rème, A. Balogh, M.W. Dunlop, C. Hashimoto, and R. TanDokoro, Transport of solar wind into Earth’s magnetosphere through rolled-up Kelvin-Helmholtz vortices, Nature 430, 755 (2004).
  * [3] D.W. Waugh, A.H. Sobel, and L.M. Polvani, What is the polar vortex and how does it influence weather?, Bull. Am. Meteorol. Soc. 98, 37 (2017).
  * [4] P.L. Read, R.M.B. Young, and D. Kennedy, The turbulent dynamics of Jupiter’s and Saturn’s weather layers: order out of chaos?, Geosci. Lett. 7, 10 (2020).
  * [5] J. Fuller, A.L. Piro, and A.S. Jermyn, Slowing the spins of stellar cores, Mon. Not. R. Astron. Soc. 485, 3661 (2019).
  * [6] M.E. Pessah, C.-K. Chan, and D. Psaltis, The signature of the magnetorotational instability in the Reynolds and Maxwell stress tensors in accretion discs, Mon. Not. R. Astron. Soc. 372, 183 (2006).
  * [7] J. Goodman and G. Xu, Parasitic Instabilities in Magnetized, Differentially Rotating Disks, Astrophys. J.  432, 213 (1994).
  * [8] J. Alves, C. Zucker, A.A. Goodman, J.S. Speagle, S. Meingast, T. Robitaille, D.P. Finkbeiner, E.F. Schlafly, G.M. Green, A Galactic-scale gas wave in the solar neighbourhood, Nature 578, 237 (2020).
  * [9] R. Fleck, The ‘Radcliffe Wave’ as a Kelvin–Helmholtz instability, Nature 583, E24 (2020).
  * [10] A. Miura, Self-organization in the two-dimensional Kelvin-Helmholtz instability, Phys. Rev. Lett.  83, 1586 (1999).
  * [11] D. Lecoanet, M. McCourt, E. Quataert, K.J. Burns, G.M. Vasil, J.S. Oishi, B.P. Brown, J.M. Stone, and R.M. O’Leary, A validated non-linear Kelvin-Helmholtz benchmark for numerical hydrodynamics, Mon. Not. R. Astron. Soc. 455, 4274 (2016).
  * [12] C. Ho and P. Huerre, Perturbed free shear layers, Annu. Rev. Fluid Mech. 16, 365 (1984).
  * [13] F.K. Browand and C.M. Ho, The mixing layer, an example of quasi two-dimensional turbulence, J. de Mec. Theor. et Appl. 2, 99 (1983).
  * [14] V.P. Starr and N.E. Gaut, Negative viscosity, Sci. Am. 223, 72 (1970).
  * [15] A. Miura and T. Sato, Theory of vortex nutation and amplitude oscillation in an inviscid shear instability, J. Fluid Mech. 86, 33 (1978).
  * [16] W. Horton, T. Tajima, and T. Kamimura, Kelvin–Helmholtz instability and vortices in magnetized plasma, Phys. Fluids 30, 3485 (1987).
  * [17] P.W. Terry, P.-Y. Li, M.J. Pueschel, and G.G. Whelan, Threshold Heat-Flux Reduction by Near-Resonant Energy Transfer, Phys. Rev. Lett. 126, 025004 (2021).
  * [18] G.G. Whelan, M.J. Pueschel, and P.W. Terry, Nonlinear Electromagnetic Stabilization of Plasma Microturbulence, Phys. Rev. Lett.  120, 175002 (2018).
  * [19] M.J. Pueschel, B.J. Faber, J. Citrin, C.C. Hegna, P.W. Terry, and D.R. Hatch, Stellarator Turbulence: Subdominant Eigenmodes and Quasilinear Modeling, Phys. Rev. Lett.  116, 085001 (2016).
  * [20] K.D. Makwana, P.W. Terry, M.J. Pueschel, and D.R. Hatch, Subdominant Modes in Zonal-Flow-Regulated Turbulence, Phys. Rev. Lett.  112, 095002 (2014).
  * [21] D.R. Hatch, F. Jenko, A.B. Navarro, and V. Bratanov, Transition between saturation regimes of gyrokinetic turbulence, Phys. Rev. Lett.  111, 175001 (2013).
  * [22] D.R. Hatch, P.W. Terry, F. Jenko, F. Merz, and W.M. Nevins, Saturation of gyrokinetic turbulence through damped eigenmodes, Phys. Rev. Lett.  106, 115003 (2011).
  * [23] M.J. Pueschel, P.-Y. Li, and P.W. Terry, Predicting the critical gradient of ITG turbulence in fusion plasmas, Nucl. Fusion 61, 054003 (2021).
  * [24] P.-Y Li and P.W. Terry, Assessing physics of ion temperature gradient turbulence via hierarchical reduced-model representations, Phys. Plasmas 29, 042301 (2022).
  * [25] G.G. Whelan, M.J. Pueschel, P.W. Terry, J. Citrin, I.J. McKinney, W. Guttenfelder, and H. Doerk, Saturation and nonlinear electromagnetic stabilization of ITG turbulence, Phys. Plasmas 26, 082302 (2019).
  * [26] P.W. Terry, B.J. Faber, C.C. Hegna, V.V. Mirnov, M.J. Pueschel, and G.G. Whelan, Saturation scalings of toroidal ion temperature gradient turbulence, Phys. Plasmas 25, 012308 (2018).
  * [27] A.E. Fraser, M.J. Pueschel, P.W. Terry, and E.G. Zweibel, Role of stable modes in driven shear-flow turbulence, Phys. Plasmas 25, 122303 (2018).
  * [28] K.D. Makwana, P.W. Terry, J.-H. Kim, and D.R. Hatch, Damped eigenmode saturation in plasma fluid turbulence, Phys. Plasmas 18, 012302 (2011).
  * [29] K.D. Makwana, P.W. Terry, and J.-H. Kim, Role of stable modes in zonal flow regulated turbulence, Phys. Plasmas 19, 062310 (2012).
  * [30] P.W. Terry, D.A. Baver, and S. Gupta, Role of stable eigenmodes in saturated local plasma turbulence, Phys. Plasmas 13, 022307 (2006).
  * [31] R.H. Levy and R.W. Hockney, Computer Experiments on Low-Density Crossed-Field Electron Beams Phys. Fluids 11, 766 (1968).
  * [32] N.J. Zabusky and G.S. Deem, Dynamical Evolution of Two-Dimensional Unstable Shear Flows, J. Fluid Mech. 47, 353 (1971).
  * [33] L.-S. Huang and C.-M. Ho, Small-scale transition in a plane mixing layer, J. Fluid Mech. 210, 475 (1990).
  * [34] R.D. Moser and M.M. Rogers, The three-dimensional evolution of a plane mixing layer: Pairing and transition to turbulence, J. Fluid Mech. 247, 275 (1993).
  * [35] J.J. Riley and R.W. Metcalfe, Direct numerical simulation of a perturbed turbulent mixing layer, 18th Aerospace Sciences Meeting AIAA Paper 1980-0274, (1980).
  * [36] D. Oster and I. Wygnanski, The forced mixing layer between parallel streams, J. Fluid Mech. 123, 91 (1982).
  * [37] Y. Ito, K. Nagata, Y. Sakai, and O. Terashima, Momentum and mass transfer in developing liquid shear mixing layers, Exp. Therm. Fluid Sci. 51, 28 (2013).
  * [38] A. López Zazueta and L. Zavala Sansón, Self-oscillations of a two-dimensional shear flow with forcing and dissipation, Phys. Fluids 30, 044101 (2018).
  * [39] A. VanDine, H.T. Pham, and S. Sarkar, Turbulent shear layers in a uniformly stratified background: DNS at high Reynolds number, J. Fluid Mech. 916, A42 (2021).
  * [40] A.K.M.F. Hussain and K.B.M.Q. Zaman, An experimental study of organized motions in the turbulent plane mixing layer, J. Fluid Mech. 159, 85 (1985).
  * [41] A.K.M.F. Hussain, Coherent structures and turbulence, J. Fluid Mech. 173, 303 (1986).
  * [42] L. Landau, On the problem of turbulence, C.R. Acad. Sci. U.R.S.S. 44, 311 (1944).
  * [43] A.E. Fraser, P.W. Terry, E.G. Zweibel, and M.J. Pueschel, Coupling of damped and growing modes in unstable shear flow, Phys. Plasmas 24, 062304 (2017).
  * [44] A.E. Fraser, P.W. Terry, E.G. Zweibel, M.J. Pueschel, and J.M. Schroeder, The impact of magnetic fields on momentum transport and saturation of shear-flow instability by stable modes, Phys. Plasmas 28, 022309 (2021).
  * [45] B. Tripathi, A.E. Fraser, P.W. Terry, E.G. Zweibel, and M.J. Pueschel, Mechanism for Sequestering Magnetic Energy at Large Scales in Shear-Flow Turbulence, arXiv:2205.01298, submitted to Phys. Plasmas (2022).
  * [46] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Clarendon Press, Oxford 1961).
  * [47] J. Mak, S.D. Griffiths, and D.W. Hughes, Vortex disruption by magnetohydrodynamic feedback, Phys. Rev. Fluids 2, 113701 (2017).
  * [48] A.A. Schekochihin, J.L. Maron, S.C. Cowley, and J.C. McWilliams, The Small-Scale Structure of Magnetohydrodynamic Turbulence with Large Magnetic Prandtl Numbers, Astrophys. J.  576, 806 (2002).
  * [49] D. Biskamp, Magnetohydrodynamic Turbulence (Cambridge University Press, Cambridge 2003).
  * [50] F. Ebrahimi, S.C. Prager, and D.D. Schnack, Saturation of Magnetorotational Instability Through Magnetic Field Generation, Astrophys. J. 698, 233 (2009).
  * [51] M.J. Pueschel, D. Told, P.W. Terry, F. Jenko, E.G. Zweibel, V. Zhdankin, and H. Lesch, Magnetic Reconnection Turbulence in Strong Guide Fields: Basic Properties and Application to Coronal Heating, Astrophys. J., Suppl. Ser. 213, 30 (2014).
  * [52] J.B. Marston, E. Conover, and T. Schneider, Statistics of an Unstable Barotropic Jet from a Cumulant Expansion, J. Atmos. Sci. 65, 1955 (2008).
  * [53] K.M. Smith, C.P. Caulfield, and J.R. Taylor, Turbulence in forced stratified shear flows, J. Fluid Mech. 910, A42 (2021).
  * [54] A. Allawala, S.M. Tobias, and J.B. Marston, Dimensional reduction of direct statistical simulation, J. Fluid Mech. 898, A21 (2020).
  * [55] 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).
  * [56] C.M. Bender, PT symmetry: In quantum and classical physics (World Scientific Publishing, 2019).
  * [57] D.R. Hatch, F. Jenko, A.B. Navarro, V. Bratanov, P.W. Terry, and M.J. Pueschel, Linear signatures in nonlinear gyrokinetics: Interpreting turbulence with pseudospectra, New J. Phys. 18, 075018 (2016).
  * [58] Y. Fu and H. Qin, The physics of spontaneous parity-time symmetry breaking in the Kelvin-Helmholtz instability, New J. Phys. 22, 083040 (2020).
  * [59] K.M. Case, Stability of Inviscid Plane Couette Flow, Phys. Fluids 3, 143 (1960).
  * [60] P.W. Terry, D.A. Baver, and D.R. Hatch, Reduction of inward momentum flux by damped eigenmodes, Phys. Plasmas 16, 122305 (2009).
  * [61] J.B. Marston, G.P. Chini, and S.M. Tobias, Generalized quasilinear approximation: application to zonal jets, Phys. Rev. Lett.  116, 214501 (2016).
  * [62] G.K. Batchelor, On the spontaneous magnetic field in a conducting liquid in turbulent motion, Proc. Roy. Soc. London, Ser. A 201, 405 (1950).
  * [63] G.K. Batchelor and I. Proudman, The effects of rapid distortion of a fluid in turbulent motion, Q. J. Mech. Appl. Math 7, 83 (1954).
  * [64] A.A. Townsend, The Structure of Turbulent Shear Flow, $2^{\mathrm{nd}}$ Ed. (Cambridge University Press, Cambridge, 1976).
  * [65] A. Alexakis, P.D. Mininni, and A. Pouquet, Phys. Rev. E 72, 046301 (2005).
  * [66] K.J. Burns, Flexible spectral algorithms for simulating astrophysical and geophysical flows, Ph.D. thesis (Massachusetts Institute of Technology, 2018).
  * [67] A.E. Fraser, Role of stable eigenmodes in shear-flow instability saturation and turbulence, Ph.D. thesis (University of Wisconsin-Madison, 2020).
  * [68] S.A. Orszag, Analytical theories of turbulence, J. Fluid Mech. 41, 363 (1970).
  * [69] K. Taira, S.L. Brunton, S.T.M. Dawson, C.W. Rowley, T. Colonius, B.J. McKeon, O.T. Schmidt, S. Gordeyev, V. Theofilis, and L.S. Ukeiley, Modal Analysis of Fluid Flows: An Overview, AIAA J 55, 4013 (2017).
  * [70] P. Garaud, Double-Diffusive Convection at Low Prandtl Number, Annu. Rev. Fluid Mech. 50, 275 (2018).