# Three-dimensional shear-flow instability saturation via stable modes

B. Tripathi1<EMAIL_ADDRESS>P.W. Terry1 A.E. Fraser2,3,4 E.G. Zweibel1,5
M.J. Pueschel6,7 1Department of Physics, University of Wisconsin-Madison,
Madison, Wisconsin 53706, USA
2Department of Applied Mathematics, University of Colorado, Boulder, Colorado
80309, USA
3Department of Astrophysical and Planetary Sciences, University of Colorado,
Boulder, Colorado 80309, USA
4Laboratory for Atmospheric and Space Physics,, University of Colorado,
Boulder, Colorado 80309, USA
5Department of Astronomy, University of Wisconsin-Madison, Madison, Wisconsin
53706, USA
6Dutch Institute for Fundamental Energy Research, 5612 AJ Eindhoven, The
Netherlands
7Eindhoven University of Technology, 5600 MB Eindhoven, The Netherlands

###### Abstract

Turbulence in three dimensions ($3$D) supports vortex stretching that has long
been known to accomplish energy transfer to small scales. Moreover, net energy
transfer from large-scale, forced, unstable flow-gradients to smaller scales
is achieved by gradient-flattening instability. Despite such enforcement of
energy transfer to small scales, it is shown here not only that the shear-
flow-instability-supplied $3$D-fluctuation energy is largely inverse-
transferred from the fluctuation to the mean-flow gradient, but that such
inverse transfer is more efficient for turbulent fluctuations in $3$D than in
two dimensions ($2$D). The transfer is due to linearly stable eigenmodes that
are excited nonlinearly. The stable modes, thus, reduce both the nonlinear
energy cascade to small scales and the viscous dissipation rate. The vortex-
tube stretching is also suppressed. Up-gradient momentum transport by the
stable modes counters the instability-driven down-gradient transport, which
also is more effective in $3$D than in $2$D
($\mathrm{\approx}70\%\mathrm{\,\,vs.\,\,}\mathrm{\approx}50\%$). From
unstable modes, these stable modes nonlinearly receive energy via zero-
frequency fluctuations that vary only in the direction orthogonal to the plane
of $2$D shear flow. The more widely occurring $3$D turbulence is thus
inherently different from the commonly studied $2$D turbulence, despite both
saturating via stable modes.

##  I Introduction

The directionality of turbulent energy transfer has long been a topic of
interest.[1] Early on, energy transfer to small scales was posited in Navier-
Stokes inertial-range turbulence, as part of the scale-invariant property of
the nonlinearity.[2, 3] Later, inviscid dynamical invariants were shown to
govern energy transfer directions in two-dimensional ($2$D) turbulence.[4, 5]
In $2$D, because the vortex tubes cannot be stretched, two positively signed
invariants exist and demand oppositely directed cascades: enstrophy to small
and energy to large scales.[4, 5] Nonlinear energy transfer is also affected
by global symmetry breaking processes like rotation,[6, 7, 8]
stratification,[9] a mean shear flow,[10, 11, 12, 13, 14, 15] or an external
magnetic field via anisotropic linear physics.[16, 17, 18, 19] The important
effects observed are spectral condensation at large scales and anisotropic
turbulent structures that coincide with the anisotropic linear physics.[20,
21] The mechanism involves the correlation time of nonlinear interactions,
which depends on anisotropic frequencies of linear waves.[6, 7, 8, 18] Thus,
anisotropy is endowed to, otherwise isotropic, nonlinear energy transfer. It
is by this mechanism that the inverse energy cascades can be explained in $3$D
rotating turbulence,[6, 7, 8] and similar arguments apply in the creation of
large-scale zonal flows[18] in fusion plasmas where the mean magnetic field
breaks the symmetry. This understanding, however, must account for the physics
of the linear instability, when it drives the turbulence.

If the instability operates at large scales,[22] the tendency of instability
to relax its driving gradient implies energy transfer to smaller scales. This
traditional argument anticipates that instability-driven turbulence will
feature forward energy transfer. While true, the linear-instability-supplied
fluctuation energy need not necessarily, in its entirety, be available to a
nonlinear energy cascade to other scales. A discrepancy can arise because the
nonlinearity can excite stable-mode solutions of the operator of the linear
stability analysis.[22, 23, 24, 25, 26, 27, 28, 30, 29, 31, 32, 33, 34, 35,
36, 37] Stable-mode excitation has also been observed in laboratory
experiment.[38] When excited to large amplitudes, the stable modes directly
transfer turbulent energy from the instability-scale fluctuations to the
driving mean gradient,[27] a process that can be interpreted as an inverse
transfer of energy. Such direct fluctuation-to-mean transfer is distinct from
the nonlinear energy cascade. However, the former can drastically reduce the
latter,[27, 29] a phenomenon that has come to light in turbulence theory in
recent years.[18, 28]

All the effects described above can arise in shear flow-driven hydrodynamic
turbulence. Our chief concerns in this paper relate to energy transfer
directionality, spatial dimensionality, and instability saturation mechanism.
We investigate the turbulence driven by a forced, unstable shear flow in three
dimensions, specifically the Kelvin-Helmholtz (KH) instability. The KH
instability is an inviscid instability[39] whose eigenspectrum is composed of:
(i) discrete eigenfrequencies that are purely imaginary and that appear in
complex-conjugate pairs, and (ii) a theoretically infinite number of
eigenmodes that form a continuum in real frequency.[26] In this study, the
discrete modes are an unstable mode and a damped (conjugate-stable) mode. The
former extracts energy from the unstable flow-gradient and the latter reverses
it.[22, 26] Being a conjugate-stable mode, its mode structure is conceptually
similar to that of the unstable mode.[24] Because of the similarity in mode
structures and complex eigenfrequencies,[28] the conjugate-stable modes can
strongly nonlinearly couple to the unstable modes, and thus get excited to
large amplitudes.

In $2$D hydrodynamic[22] and magnetized shear flows,[24, 26, 25] the
conjugate-stable modes have been shown to be significantly excited in
saturation. The modes sequester energy at large scales, reduce the nonlinear
energy cascade, and lower the small-scale dissipation rate.[25] Driving up-
gradient momentum transport, the stable modes counter the instability-driven
down-gradient transport.[24, 26]

In $3$D fluid turbulence, the strong excitation of conjugate-stable modes
remains to be confirmed. Possible reasons for a lack of such a mechanism
relate to the added degree of freedom in $3$D, which may result in an
inefficient stable-unstable mode coupling, thus allowing the majority of
energy to cascade to small scales. However, with the additional dimension, the
number of new unstable and conjugate-stable modes (hereafter interchangeably
called stable modes) also increases.[39] Since the turbulence is inherently a
multi-scale nonlinear phenomenon, with complexity further enhanced by the
presence of large-scale energy sources and sinks in the form of unstable modes
and stable modes, reliable understanding can only be gained with quantitative
assessment. A detailed analysis of the energetics of $3$D turbulence is, thus,
crucially needed.

Instability saturation can also be influenced by the anisotropic dispersion
relation of the $3$D KH instability, where the anisotropy arises because the
perturbations that vary only in the spanwise direction (orthogonal to the
direction of the mean flow and the direction of shear) are distinct from the
perturbations that vary only in the direction of the mean flow. When the
spanwise perturbations are analyzed in terms of an eigenspectrum, one finds
neither the KH instability nor the waves forming a continuum in frequency: all
eigenfrequencies of the inviscid linear operator collapse to zero.[39] One
might, thus, assume that these spanwise fluctuations play no role in the
saturation of the KH instability, as the instability is largely caused by the
fastest-growing two-dimensional streamwise perturbations.[39] However,
nonlinear mode couplings can strongly drive the zero-frequency, purely
spanwise fluctuations, making the $3$D turbulence different from its $2$D
counterpart in essential regards. For example, the spanwise fluctuations, due
to their zero frequencies, can catalyze near-resonant[6, 7, 8, 18] energy
transfer between two waves of similarly high frequencies, as such a triad
maximizes the turbulent interaction time[30, 29]; in the context of the
inviscid-instability-driven dissipative turbulence, this translates to energy
transfer between eigenmodes of similar but opposite growth rates.[30, 29, 36,
37] Hence, this new channel involving spanwise fluctuation can participate in
the nonlinear saturation of the instability in $3$D. However, whether the
participation becomes dominant is a priori unknown.

Even in linear stability analysis, the $2$D and $3$D KH instabilities are
different, as shown already in 1933 by Squire, stating[40] that $2$D
perturbations are the least stable. Technically valid only for linear
analysis, the theorem is relied upon in insightful nonlinear studies of $2$D
turbulence, ignoring the possibility that the $3$D KH-unstable modes may
compete nonlinearly with the $2$D KH-unstable modes, and even dominate. A part
of the reason is the (computational) cost-effectiveness of the $2$D analysis.
Such studies commonly appeal to nonlinearly[7, 8] two-dimensionalizing effects
such as a strong rotation, although the rotation can introduce new waves.[41]
There have, however, been illuminating studies on $3$D KH-instability-driven
turbulence, where fully $3$D fluctuations are connected to the secondary and
tertiary instabilities of the fluctuations generated by the primary $2$D KH
instability (see, e.g, Refs. [43, 44, 45, 46, 47, 48, 42]). While such $3$D
higher-order instability analyses offer a unique perspective on the
development of the $3$D turbulence and characterize the “zoo” of secondary and
tertiary $3$D structures, it can be difficult to extract the importance of the
$3$D primary KH mode evolution because such slowly growing KH modes can take a
longer time to rise to appreciable levels. Before those are reached, the mean
shear flow can decay, causing a premature demise of the $3$D KH modes The
issue is, therefore, not just $3$D, but the time scale that the $3$D primary
KH-unstable modes need to be allowed to evolve self-consistently without
disadvantaging them by preventing the loss of the instability. A continuously
active free energy source for instability can arise in natural and laboratory
flows, including those in galaxies, stars, protoplanetary disks, atmospheres,
oceans, at the mouths of rivers and estuaries, and nuclear fusion devices,
where flow-shear is supplied externally; see, e.g., Refs. [48, 49, 50, 51, 52,
53, 54, 55] In such cases, how the instability saturates nonlinearly is not
understood.

Momentum transport by the KH instability is countered by the stable modes:[22,
23] the stable modes steepen the mean flow gradient that the KH instability
attempts to flatten. Understanding the efficiency of these competing processes
allows to construct reduced and truncated models of turbulent transport.[24,
26, 29] If the unstable and stable modes primarily contribute to momentum
transport in the $3$D KH-instability-driven turbulence, accurate transport
models may be achieved. Exploring such possibilities, therefore, benefits geo-
and astrophysical communities (see, e.g., Refs. [56, 57]) who wish to
understand the turbulence and predict the momentum and particle transport
rates due to the $3$D KH instability in nature.

One of the principal findings of this paper is that the $3$D KH instability
nonlinearly excites zero-frequency, spanwise fluctuations that near-resonantly
transfer energy to the conjugate-stable modes at KH-unstable wavenumbers.
These stable modes then return the turbulent energy to the mean flow. Such
stable modes inhibit the energy cascade to small scales, viscous dissipation
rate, momentum transport across the shear layer, and vortex stretching—a
purely $3$D phenomenon.

This paper is organized as follows. Section II presents the model of shear-
flow turbulence and simulation details. In Sec. III, a complete eigenmode
basis is derived using an inviscid linear operator. Stable-mode excitation is
demonstrated in Sec. IV. Section V shows significant counter-gradient momentum
transport by stable modes and presents a reduced transport model. Inverse
transfer of energy by the stable modes from the fluctuations to the mean flow
is investigated in Sec. VI. How such stable modes receive energy in the
nonlinear phase is explored in Sec. VII. The stable modes are shown in Sec.
VIII to retard $3$D vortex stretching and subdue small-scale dissipation.
Section IX presents broader issues and implications of the results, before
conclusions in Sec. X.

##  II Model and simulation setup

An incompressible hydrodynamic fluid evolves[39] according to

$\displaystyle\partial_{t}\mathbf{u}+\mathbf{u}\cdot\mathbf{\nabla}\mathbf{u}=-\nabla
p+\nu\nabla^{2}\mathbf{u}+\mathbf{f},$ (1a)
$\displaystyle\mathbf{\nabla}\cdot\mathbf{u}=0,$ (1b)

where $\mathbf{u}$ is the fluid velocity, $\nu$ the kinematic viscosity, and
$\mathbf{f}$ the forcing, or externally supplied acceleration. The pressure
(per unit density) $p$ imposes the incompressibility constraint Eq. (1b) on
the flow.

###  II.1 Background flow and forcing

We consider an initial flow profile
$\mathbf{u}(x,y,z,t\mathrm{=}0)=U_{0}\tanh(z/a)\hat{\mathbf{e}}_{x}$, which is
well-known to be unstable.[39] The flow amplitude $U_{0}$ and half-width $a$
of the flow-shear are the initial characteristic speed and length scale that
are used to non-dimensionalize all the variables henceforth. Thus, time is
measured in units of $a/U_{0}$ and energy (per unit mass) in terms of
$U_{0}^{2}$. The initial mean flow thus becomes
$\mathbf{u}(x,y,z,t\mathrm{=}0)=\tanh(z)\hat{\mathbf{e}}_{x}$. The viscosity
is non-dimensionally measured via the Reynolds number
$Re\mathrm{=}U_{0}a/\nu$.

In natural systems, the unstable mean shear flow can evolve in at least two
different qualitative ways: first, when the mean shear flow is allowed to
freely evolve following the response from the shear flow-driven turbulence,
the unstable shear flow relaxes to a stable configuration. Second, when the
mean shear flow is allowed to evolve in the presence of the external drive
that brings the shear flow in existence in the first place, the resulting
turbulence can be of quasi-stationary nature. Admittedly, the form of the
external drives can vary from system to system; avoiding such a system-
specific requirement, we use a general form of external forcing that
replenishes the initial, unstable equilibrium. Such type of restoring force
have previously been used to model external forces in natural
environments,[49, 48, 50, 58, 59] e.g., jets in stars, driven flows in geo-and
astrophysics, which can be related to various forms of external drives,[48,
49, 50, 51, 52, 53] for example, fingering convection that ceaselessly seeds
secondary instability with shear flow,[54, 55] and, in general, the
gravitational force that generates shear layers in accretion disks and stellar
interiors. We thus externally drive the mean flow using
$\mathbf{f}=f(k_{x}\mathrm{=}0,k_{y}\mathrm{=}0,z,t)\hat{\mathbf{e}}_{x}$,
where $k_{x}$ and $k_{y}$ are the wavenumbers along the $x$\- and $y$-axes,
respectively, with

$f=\frac{1}{t_{\mathrm{forcing}}}\left[U_{\mathrm{ref}}(z)-\left\langle
u_{x}(x,y,z,t)\right\rangle_{x,y}\right]+F_{0},$ (2)

where $U_{\mathrm{ref}}(z)\mathrm{=}\tanh(z)$ is the initial mean flow;
$\left\langle u_{x}(x,y,z,t)\right\rangle_{x,y}$ is the $(x,y)$-averaged
instantaneous flow, directed along the $x$-axis; and $1/t_{\mathrm{forcing}}$,
with units of $U_{0}/a$, is the rate at which the mean flow is forced toward
the initial equilibrium profile.[49, 48, 26] A time-independent forcing
$F_{0}$ helps to realize a true equilibrium of the initial mean flow:
$Re^{-1}\nabla^{2}U_{\mathrm{ref}}(z)+F_{0}=0$, allowing a formal linear
analysis.

On this equilibrium system, small-amplitude initial noise is seeded for the
present simulation studies.

###  II.2 Initial perturbations, boundary conditions, and simulation
parameters

To ensure incompressiblity $\mathbf{\nabla}\cdot\mathbf{u}\mathrm{=}0$ of the
flow in the initial random perturbations, careful steps are taken by first
writing $\mathbf{u}\mathrm{=}\nabla\times\mathbf{C}$, where
$\mathbf{C}\mathrm{=}\sum_{\mu\mathrm{=}\\{x,y,z\\}}C_{\mu}\hat{\mathbf{e}}_{\mu}$
is a $3$-dimensional vector field, with each component given as

$C_{\mu}(x,y,z)=\sum_{\begin{subarray}{c}0<|k_{x}|<k_{x}^{\mathrm{max}}\\\
0<|k_{y}|<k_{y}^{\mathrm{max}}\end{subarray}}\mathrm{e}^{ik_{x}x}\mathrm{e}^{ik_{y}y}\hat{C}_{\mu}(k_{x},k_{y},z).$
(3)

Initial perturbations are excited up to $k_{x}=k_{x}^{\mathrm{max}}$ and
$k_{y}=k_{y}^{\mathrm{max}}$. The Fourier components of such perturbations are

$\hat{C}_{\mu}(k_{x},k_{y},z)=\hat{\mathbf{e}}_{\mu}A\,\mathrm{e}^{-\frac{z^{2}}{\sigma^{2}}}\Bigg{(}|k_{x}|^{\alpha}\mathrm{e}^{i\theta_{\mu}(k_{x})}\delta_{k_{y},0}+|k_{y}|^{\alpha}\mathrm{e}^{i\theta_{\mu}(k_{y})}\delta_{k_{x},0}+|k_{x}k_{y}|^{\beta}\mathrm{e}^{i\theta_{\mu}(k_{x},k_{y})}\Bigg{)},$
(4)

where $A$ is an overall amplitude of the perturbation (measured non-
dimensionally with respect to the initial mean-flow amplitude $U_{0}$);
$\sigma$ is the width (measured non-dimensionally with respect to the half-
width of the flow-shear $a$) along the $z$-axis of the Gaussian envelope that
scales the initial perturbations; $\alpha$ and $\beta$ are the spectral slopes
of the perturbations in the wavenumber space $k_{x}$ and $k_{y}$,
respectively; each component of $\mathbf{C}$ is randomized in phase at every
wavenumber via random phase functions
$\theta_{\mu}(k_{x}),\theta_{\mu}(k_{y}),$ and $\theta_{\mu}(k_{x},k_{y})$. On
the right-hand side of Eq. (4), the first term inside the brackets excites the
$k_{y}\mathrm{=}0$ Fourier modes, which are the purely $2$D perturbations; the
second term corresponds to the fluctuations that vary only on the
$(y,z)$-plane, orthogonal to the streamwise direction; and the third term
excites perturbations that vary in all directions.

Motivated by earlier studies of $2$D shear-flow turbulence,[24, 25, 26, 27] we
use a small-amplitude perturbation, with $A\mathrm{=}10^{-3}$, to allow an
unambigious and decades-long linear evolution of the instability, before it
gets nonlinearly modified. The $2$D studies employed only the first term
inside the brackets on the right-hand side of Eq. (4). The initially excited
perturbations, up to $k_{x}^{\mathrm{max}}=128$ and
$k_{y}^{\mathrm{max}}=128$, have a near-flat spectrum of energy over the
($k_{x},k_{y}$)-plane, which we achieve by choosing $\alpha\mathrm{=}-1$ and
$\beta\mathrm{=}-1/2$ in Eq. (4). The perturbations peak at $z\mathrm{=}0$ and
have a Gaussian width of $\sigma\mathrm{=}2$.

To make the initial perturbations fully free of $\mathbf{u}$-divergences, the
pressure $p$ is obtained using
$p\mathrm{=}\nabla^{-2}\nabla\cdot\left[-(\mathbf{u}\cdot\nabla)\mathbf{u}\right]$
where $\nabla^{-2}$ stands for inversion of the Laplacian operator.

Boundary conditions along the $x$\- and $y$-axes are periodic, whereas, along
the $z$-axis, we impose no-slip, co-moving (with the initial flow) walls at
the top and bottom boundaries $z\mathrm{=}\pm L_{z}/2$.

While we have studied the current system with varying forms of initial
condition, box size, spectral resolution, and the Reynolds number, we present
in this paper results using a box size of
$L_{x}\mathrm{=}L_{y}\mathrm{=}L_{z}\mathrm{=}10\pi$, resolving a large number
of linearly KH-unstable wavenumbers, and using a spectral resolution of
$512^{3}$: Fourier modes on the $(x,y)$-plane and Chebyshev polynomials along
the $z$-axis. To avoid aliasing errors, we dealias the quadratic
nonlinearities of the governing equation via the $3/2$-rule in all directions.
Motivated by earlier studies,[25, 26, 27] $t_{\mathrm{forcing}}\mathrm{=}0.5$
is chosen throughout this study (the instability-growth time is $\sim 6$ and
the box-crossing time of the mean flow is $10\pi$). Such a choice adequately
replenishes the initial mean flow, thereby preventing the loss of the KH
instability. We use the pseudospectral solver Dedalus,[60] with third-order,
four-stage implicit-explicit Runge-Kutta time-integrator to advance the
system. At the cost of $\mathrm{\sim}10$ million CPU-hours, long time-
integrated $3$D simulations are studied for $Re=300,1500,$ and $5000$, in
addition to an ensemble of $10$ different three-dimensional numerical
simulation continuations (in Sec. VIII).

##  III Complete eigenmode basis

Before analyzing the nonlinear state of the shear-flow instability, we find
the inviscid eigenmodes by linearizing Eqs. (1a) and (1b), dropping the
viscous term, and Fourier-transforming the equations,

$\displaystyle\partial_{t}\hat{\mathbf{u}}+\mathbf{u}_{0}\cdot\mathbf{\nabla}\hat{\mathbf{u}}+\hat{\mathbf{u}}\cdot\mathbf{\nabla}\mathbf{u}_{0}=-\nabla\hat{p},$
(5a) $\displaystyle\mathbf{\nabla}\cdot\hat{\mathbf{u}}=0,$ (5b)

where $\mathbf{u}$ is decomposed into its mean
($k_{x}\mathrm{=}k_{y}\mathrm{=}0$) mode $\mathbf{u}_{0}$ and fluctuation
spectrum $\hat{\mathbf{u}}$ at $(k_{x},k_{y})$; likewise for the pressure
$\hat{p}$; the derivative operator $\nabla$ can be expressed as
$ik_{x}\hat{\mathbf{e}}_{x}+ik_{y}\hat{\mathbf{e}}_{y}+\hat{\mathbf{e}}_{z}\partial_{z}$.

With $\mathbf{u}_{0}\mathrm{=}U_{\mathrm{ref}}(z)\hat{\mathbf{e}}_{x}$, a
complete eigenspectrum can be found at any given wavenumber. When
$k_{y}\mathrm{=}0$, one obtains an eigenspectrum corresponding to $2$D
perturbations:[26] (i) unstable and conjugate-stable eigenmodes at
$|k_{x}|\lesssim 1$, frequencies of which are complex-conjugates, and (ii) a
large (theoretically infinite) number of eigenmodes that form a continuum in
real frequency and that exist at all $k_{x}\mathrm{\neq}0$. This set of
eigenmodes also exists for $3$D perturbations. To obtain the eigenmodes, we
discretize Eqs. (5a) and (5b) using $1024$ Chebyshev polynomials for each
variable ($\hat{u}_{x},\hat{u}_{y},\hat{u}_{z},\hat{p}$) along the $z$-axis,
and form a $4096\times 4096$-sized matrix eigenvalue problem at each
wavenumber, and then diagonalize the matrix to find the complete
eigenspectrum. The $3$D unstable and stable eigenmodes are visually very
similar as can be seen, for example, when examining their $z$-components of
velocity at the first Fourier mode numbers along the $x$\- and $y$-axes in
Fig. 1. The similarity appears because, in wavenumber space, they exhibit a
special complex-conjugate symmetry:
$(\hat{u}_{x},\hat{u}_{y},\hat{u}_{z},p)_{\mathrm{stable\
mode}}=(\hat{u}_{x}^{\ast},\hat{u}_{y}^{\ast},-\hat{u}_{z}^{\ast},p^{\ast})_{\mathrm{unstable\
mode}}$, where ∗ denotes the complex-conjugation operation.

Figure 1: For an unstable mean flow
$\mathbf{u}_{0}\mathrm{=}\tanh(z)\hat{\mathbf{e}}_{x}$, physical structures of
perturbed $\hat{u}_{z}$ of $3$D unstable and conjugate-stable modes are shown;
the colors represent the magnitude of the $\hat{u}_{z}$ (in arbitrary units).
The isocontours and the surface rendering highlight a symmetry that is
illustrated with a cartoon.

A growth rate spectrum of the unstable mode on the ($k_{x},k_{y}$)-plane is
shown in Fig. 2. At every wavenumber for which the unstable mode of the
inviscid linear operator exists, the conjugate-stable mode must necessarily,
as well exist—meaning that Fig. 2 also serves as a damping rate spectrum of
the conjugate-stable mode. For the special case of $k_{x}\mathrm{=}0$, the
entire eigenspectrum becomes infinitely degenerate, i.e., all eigenvalues
collapse to zero on the complex-frequency plane. In that case, the dropped
viscous term becomes singularly important, and is the sole active term that
damps the fluctuation spectrum.

Figure 2: Growth-rate spectrum of $3$D KH instability, showing an anisotropic
dispersion relation. The $2$D modes, $k_{y}\mathrm{=}0$, have the largest
growth rates. All wavenumbers with $k_{x}\mathrm{=}0$, shown with a vertical
blue dashed line, are linearly stable. The wavenumber $(0,0)$ corresponds to
the mean flow. As the Kelvin-Helmholtz instability is a large-scale
instability, the wavenumbers shown in white, whose boundary is approximated by
the black bounding curve, are linearly stable.

##  IV Nonlinear excitation of stable modes

Turning now to the nonlinear saturation of the instability, at any wavenumber
$\mathbf{k}\mathrm{=}(k_{x},k_{y})$, where discrete linear eigenfrequencies in
this system exist, an arbitrary Fourier-transformed state vector
$\hat{\chi}\mathrm{=}(\hat{u}_{x},\hat{u}_{y},\hat{u}_{z},\hat{p})$ can be
decomposed in the basis spanned by the complete set of eigenmodes in the
manner of $\hat{\chi}\mathrm{=}\sum_{j}\beta_{j}\hat{\chi}_{j}$, where
$\beta_{j}$ is the complex-valued amplitude of the $j$th eigenmode
$\hat{\chi}_{j}$. Since the eigenmodes at hand are non-orthogonal, the mode-
amplitudes $\beta_{j}$ in the turbulent state cannot be determined simply by
taking a Hermitian inner product between the eigenmodes. However, this is
merely a technical challenge; the $\beta_{j}$ can be uniquely determined by
employing the biorthogonal basis, found by solving an adjoint of Eqs. (5a) and
(5b), which produces the same eigenspectrum but a different set of eigenmodes,
usually known as left eigenmodes or adjoint solutions. To gain machine-
precision accuracy,[26, 61] we use adjoint solutions of discretized equations
rather than that of the analytical equations.

As the unstable modes grow exponentially in the linear phase, they attain
large amplitudes. Then the unstable modes, via nonlinear process, excite
eigenmodes that are linearly stable (i.e., decaying in the absence of
nonlinearity), as shown with the orange curves in Figs. 3(a) and 3(c). The
nonlinear interaction quickly involves multiple eigenmodes and multiple
wavenumbers, and, in the fully nonlinear phase, in Figs. 3(b) and 3(d), the
conjugate-stable mode attains an amplitude nearly identical to that of the
unstable mode. These figures show that the stable modes are excited even in
$3$D. Their excitation is universal, as long as the viscosity is not
enormously large to lead to a laminar flow. In support of this, we find, for
$Re=1500$ and other lower $Re$ cases (not shown), large excitation of stable
modes.

Figure 3: Amplitude evolution of unstable ($|\beta_{1}|$) and conjugate-stable
($|\beta_{2}|$) modes at $k_{x}\mathrm{=}k_{y}\mathrm{=}0.2$. In (a) and (b),
$Re\mathrm{=}5000$ is used in the initial value simulation; and in (c) and
(d), $Re\mathrm{=}300$ is used. In each case, the left panel shows the early
evolution using a logarithmic vertical scale, whereas the right panel displays
the full time evolution using linear scales on both axes. Stable modes are
nonlinearly excited much before the saturation of the unstable mode.

##  V Efficient $3$D counter-gradient momentum transport

We inquire if momentum is transported differently across the mean shear layer
due to the stable-mode-rich turbulence. Prior studies have shown counter-
gradient momentum transport due to stable modes, but only in $2$D systems,[24,
26, 62] leaving the question open whether such a result is generalizable to
$3$D systems. To address this issue, we measure the Reynolds stress by
decomposing the turbulent fluctuations at each KH-unstable wavenumber into
eigenmodes and by evaluating stress contribution from each eigenmode. Since
the momentum is transported maximally at the middle of the shear layer
$z\mathrm{=}0$, because of the strongest shear there, we compare the transport
rates at $z\mathrm{=}0$ from individual eigenmodes, along with the total
transport from mode-undecomposed fluctuations at $z\mathrm{=}0$.

The time-averaged spectrum of the Reynolds stress is shown in Fig. 4: in blue,
the stress from an unstable mode; in orange, the stress from a stable mode;
and in black, the total stress without performing the eigenmode decomposition.
The total stress $\tau_{\mathrm{all}}(z\mathrm{=}0)$ is large (the black bars
are taller) for low wavenumbers. Mode-decomposition shows that transport from
$3$D unstable modes is large, and correspondingly large is that from the $3$D
stable modes. This immediately challenges a prediction of transport based on
unstable modes only.[63, 64, 65, 54, 66] A correct transport level, however,
can be found by simply adding the contribution of stable modes.

To separate the $2$D and $3$D modes, we apply a wavenumber-summation to Fig. 4
and obtain a $2\mathrm{D-}3\mathrm{D}$ comparison in Fig. 5. It is clear that
the total stress, shown in black, is higher for the $3$D modes than for the
$2$D modes. This is easily explained with the aid of Fig. 4 where many $3$D
modes have comparable stress contribution to that of the $2$D modes, in
addition to the higher density of states of the $3$D modes, i.e., number of
wavenumbers $k_{y}$ per $k_{x}$.

Figure 4: Time-averaged spectrum of Reynolds stress due to: the unstable mode
$\tau_{1}$, stable mode $\tau_{2}$, and mode-undecomposed full fluctuations
$\tau_{\mathrm{all}}$. Where the peaks of the blue and black bars are at
similar heights, $\tau_{1}+\tau_{2}$ is a good approximation to
$\tau_{\mathrm{all}}$. The difference between $\tau_{1}+\tau_{2}$ and
$\tau_{\mathrm{all}}$ is small, and arises from the contribution of the
continuum modes to the total stress $\tau_{\mathrm{all}}$. All stresses are
averaged over $t\mathrm{=}200\textrm{--}708$ for $Re\mathrm{=}5000$. Figure 5:
Comparison of time-averaged Reynolds stresses at $z\mathrm{=}0$ for all KH-
unstable $3$D vs. $2$D modes. Further decomposition shows contributions from
the unstable mode $\tau_{1}$; stable mode $\tau_{2}$; the two added together,
$\tau_{1}+\tau_{2}$; and the mode-undecomposed full fluctuations
$\tau_{\mathrm{all}}$. The difference between $\tau_{\mathrm{all}}$ and
$\tau_{1}+\tau_{2}$ arises from the stress due to the continuum modes.
Stresses are averaged over $t\mathrm{=}200\textrm{--}708$ for
$Re\mathrm{=}5000$.

##  VI Inverse transfer of energy in $3$D

We now analyze the energy transfer rate $Q_{1}(\mathbf{k})$ from the mean flow
to the fluctuations at $\mathbf{k}$ by the unstable mode, and the rate of
inverse energy transfer $Q_{2}(\mathbf{k})$ from the fluctuation $\mathbf{k}$
to the mean flow by the stable mode. Both of these are linear processes, as
one of the involved wavenumbers is the mean flow—the reservoir of the free
energy for the instability and thus turbulence. Two additional linear
processes, however, exist: first, because the mean flow also evolves, the
instantaneous deviation of the mean from the initial profile can induce a
(typically small) linear coupling between the inviscid eigenmodes of the
initial profile, thus reducing $Q_{j}$ by an amount $R_{j}$. Second, the
viscous term introduces energy dissipation $D_{j}$. However, since the
unstable and stable eigenmodes that are analyzed here exist only at large
scales ($|k|\lesssim 1$), the dissipative effect at such scales is small,
compared to $Q_{j}$. This remains true as long as the viscosity is not
enormously large to stabilize completely the inviscid Kelvin-Helmholtz (KH)
instability; such a case, however, is irrelevant to the scenario we are
interested in, where the KH instability drives the turbulence.

To quantify the energy transfer rates described in the previous paragraph, we
Fourier-transform Eqs. (1a) and (1b), and project the equation onto an
eigenmode of interest, say the $j$th eigenmode, by multiplying the Fourier-
transformed equation with the $j$th modified left eigenmode $\hat{Z}_{j}$; see
Appendix A of Ref. [26] for a general overview of modified left eigenmodes. To
illustrate this here, Eqs. (1a) and (1b) can be structurally written, at
$\mathbf{k}\mathrm{\neq}(0,0)$, as

$\partial_{t}M\hat{\chi}=L\hat{\chi}+\sum_{\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime}:\mathbf{k}^{\prime}+\mathbf{k}^{\prime\prime}=\mathbf{k}}N(\hat{\chi}^{\prime},\hat{\chi}^{\prime\prime}),$
(6)

where the matrix $M$ of size $4\times 4$ has the form $[[\mathbb{I}_{3\times
3},0],[0,0]]$, with $\mathbb{I}_{3\times 3}$ as the identity matrix of size
$3\times 3$; furthermore, $L$ is a linear operator and $N$ a nonlinear
operator, with $\hat{\chi}^{\prime}$ and $\hat{\chi}^{\prime\prime}$
representing state vectors at $\mathbf{k}^{\prime}$ and
$\mathbf{k}^{\prime\prime}$, respectively.

Projecting Eq. (6) on the modified left eigenmode $\hat{Z}_{j}$, we obtain the
mode-amplitude evolution equation

$\partial_{t}\beta_{j}(\mathbf{k})=\gamma_{j}(\mathbf{k})\beta_{j}(\mathbf{k})+\sum_{l}B_{jl}(\mathbf{k})\beta_{l}(\mathbf{k})+\sum_{\begin{subarray}{c}\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime}:\\\
\mathbf{k}^{\prime}+\mathbf{k}^{\prime\prime}=\mathbf{k}\end{subarray}}\sum_{m,n}C_{jmn}(\mathbf{k},\mathbf{k}^{\prime})\beta_{m}^{\prime}\beta_{n}^{\prime\prime},$
(7)

where we have used the biorthogonality relation between the modified left
eigenmode $\hat{Z}_{j}$ and the right eigenmode $\hat{\chi}_{l}$ such that
$\langle\hat{Z}_{j},\hat{\chi}_{l}\rangle\propto\delta_{j,l}$. In Eq. (7),
$\gamma_{j}(\mathbf{k})$ is the complex-valued growth rate of the $j$th
eigenmode at a wavenumber $\mathbf{k}$; the coefficient $B_{jl}(\mathbf{k})$
measures the linear coupling (between the inviscid eigenmodes at $\mathbf{k}$)
induced by the viscous term and by the instantaneous deviation of the mean
flow from the initial profile; the coefficient
$C_{jmn}(\mathbf{k},\mathbf{k}^{\prime})$ measures the nonlinear coupling
between the eigenmode $m$ at $\mathbf{k}^{\prime}$ and the eigenmode $n$ at
$\mathbf{k}^{\prime\prime}$, driving the eigenmode $j$ at $\mathbf{k}$; and
$\beta_{m}^{\prime}$ and $\beta_{n}^{\prime\prime}$ represent the amplitudes
of modes $m$ at $\mathbf{k}^{\prime}$ and $n$ at $\mathbf{k}^{\prime\prime}$,
respectively. Now, we multiply Eq. (7) with $\beta_{j}^{\ast}$, and add the
complex-conjugate of the resulting equation to obtain the evolution equation
for the eigenmode energy, which reads

$\partial_{t}|\beta_{j}(\mathbf{k})|^{2}=Q_{j}+\left(R_{j}+D_{j}\right)+T_{j},$
(8)

where

$\displaystyle Q_{j}$
$\displaystyle=2\real\left[\gamma_{j}(\mathbf{k})\right]|\beta_{j}(\mathbf{k})|^{2},$
(9a) $\displaystyle R_{j}+D_{j}$
$\displaystyle=2\real\Big{[}\sum_{l}B_{jl}\beta_{l}\beta_{j}^{\ast}\Big{]},$
(9b) $\displaystyle T_{j}$
$\displaystyle=\sum_{\begin{subarray}{c}\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime}:\\\
\mathbf{k}^{\prime}+\mathbf{k}^{\prime\prime}=\mathbf{k}\end{subarray}}\sum_{m,n}T_{jmn}(\mathbf{k},\mathbf{k}^{\prime}),$
(9c) $\displaystyle T_{jmn}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle=2\real[C_{jmn}(\mathbf{k},\mathbf{k}^{\prime})\beta_{m}^{\prime}\beta_{n}^{\prime\prime}\beta_{j}^{\ast}].$
(9d)

In Eq. (8), $Q_{j}$ measures the energy transfer rate by the $j\mathrm{th}$
eigenmode to the fluctuation scale $\mathbf{k}$ from the initial mean flow;
$R_{j}$ is the energy transfer rate by the $j\mathrm{th}$ eigenmode to the
fluctuation scale $\mathbf{k}$, the energy source here being the instantaneous
deviations of the mean flow from the initial mean flow; $D_{j}$ is the viscous
dissipation rate by the $j\mathrm{th}$ eigenmode at the fluctuation scale
$\mathbf{k}$; $T_{j}$ is the nonlinear drive to the $j\mathrm{th}$ eigenmode.
Since the eigenmodes used for the basis expansion are the inviscid eigenmodes
of the initial mean flow, the terms $R_{j}$ and $D_{j}$ introduce a linear
coupling $B_{jl}$ between those eigenmodes, the effect of which we measure
next.

For both the $2$D and $3$D perturbations that the system generates, the energy
transfer rates $Q_{j}$ and $R_{j}+D_{j}$ are shown in Figs. 6(a) and 6(b) for
the energetically dominant wavenumbers. The energy transfer reduction term
$R_{j}+D_{j}$ is of opposite sign compared to $Q_{j}$. This physically means
that unstable modes are slightly less efficient in extracting energy from the
instantaneous mean flow when compared with their efficiency of energy
extraction from the initial mean flow; the same applies to the energy return
rate by the stable modes. This particular finding is similar to the results
reported for $2$D MHD shear-flow turbulence.[27]

Figure 6: Time traces of energy transfer to an unstable, $j\mathrm{=}1$, and a
stable, $j\mathrm{=}2$, mode at wavenumbers with (a) $2$D and (b) $3$D
perturbations, in a simulation with $Re\mathrm{=}5000$. The quantity $Q_{j}$
denotes the energy transfer rate from the initial mean flow to the
$j\mathrm{th}$ eigenmode at fluctuation scale $(k_{x},k_{y})$;
$R_{j}\mathrm{+}D_{j}$ collectively represents the effect of time-deviations
of the instantaneous mean flow from the initial profile and the effect of
viscous dissipation. In (c), $2$D and $3$D perturbations, decomposed by
unstable and stable modes, are compared in terms of the energy transfer rate
$Q_{j}$; the stable modes reverse the fluctuation energy to the mean flow at
all times, once the nonlinearity excites them.

At early times $t\lesssim 140$, the $2$D mode dominates over the $3$D mode;
Fig. 6(c). However, in the turbulent phase, the slowly growing $3$D wavenumber
competes with the $2$D wavenumber. In all cases, the stable modes are
significantly excited, nearly to the levels of the corresponding unstable
modes.

We note that, because of the added degree of freedom (the $y$-coordinate) in
the $3$D system, there exist more wavenumbers that are $3$D KH-unstable than
that are $2$D KH-unstable; for example, there is only one wavenumber that is
$2$D KH-unstable at $k_{x}\mathrm{=}0.2$, i.e., $(0.2,0)$, whereas there are
multiple wavenumbers that are $3$D KH-unstable at $k_{x}\mathrm{=}0.2$,
namely, $(0.2,\pm 0.2)$, $(0.2,\pm 0.4)$, etc. We avoid double-counting the
Hermitian conjugates ($-k_{x},-k_{y}$) of these wavenumbers ($k_{x},k_{y}$),
which are already included in the definition of $Q_{j}(\mathbf{k})$.

Time-averaged energy-transfer rates are compared in Fig. 7, where the rates
are decomposed by unstable and stable modes. Despite growing slowly, the $3$D
wavenumbers contribute significantly in the fully nonlinear phase in
extracting energy from the mean flow. However, the stable modes at the same
wavenumbers, whether $2$D or $3$D, always efficiently deplete fluctuation
energy and transfer energy in the reverse direction: from the fluctuation to
the mean flow.

The $3$D modes are more efficient than the $2$D modes in extracting energy
from the mean flow, as Fig. 8 shows. Additionally, the stable-to-unstable mode
fraction is higher for $3$D modes. The presence of the energy reversal by the
stable modes in $3$D confirms that the stable modes are not the manifestations
of the traditional inverse cascade of $2$D turbulence; the former is a direct
fluctuation-to-mean transfer of energy, whereas the latter is a fluctuation-
to-fluctuation transfer that takes multiple iterations and is a nonlinear
process. The quantities measured in Fig. 8 are the rates of energy transfer
directly between the mean and the fluctuations; we decompose such a linear
energy transfer rate into the contributions of the unstable and stable modes,
going beyond the simple traditional spectral transfer. The significant direct
energy reversal in $3$D by the stable modes is a somewhat surprising and
important finding of this paper.

Figure 7: Time-averaged rate of energy transfer from initial mean to
fluctuation $Q_{1}$ by unstable modes, and inverse transfer $Q_{2}$ by stable
modes, in a simulation with $Re\mathrm{=}5000$; the time-averaging interval is
$t\mathrm{=}200\textrm{--}708$. Because of the Hermitian conjugacy in the
perturbations over the wavenumbers, only a non-redundant half of the
$(k_{x},k_{y})$-plane is shown. The linearly fastest growing mode exists at
$\mathbf{k}\mathrm{=}(0.4,0)$. Figure 8: Energy transfer rates, classified by
$2$D and $3$D modes, with further decomposition into the unstable and the
stable modes. In the simulation with $Re\mathrm{=}5000$, time averaging is
taken over $t\mathrm{=}200\textrm{--}708$. The stable-to-unstable mode
fraction is larger for $3$D modes than for $2$D.

##  VII Nonlinear energy transfer to stable modes

Since the stable modes are found to contribute significantly to the energetics
and transport, what nonlinear process transfers energy to the stable modes?
Established[67, 68] diagnostics of energy transfer in wavenumber space fail
here, because they only inform about the energy transfer to a given
wavenumber, whereas, at any wavenumber, the total fluctuation is composed of
many different eigenmodes, making it desirable to put an eigenmode filter in
the transfer analysis.

The total nonlinear energy transfer $T_{2}$ from all the fluctuations to a
stable mode $j\mathrm{=}2$ at $\mathbf{k}\mathrm{=}(0.2,0)$ is shown in Fig.
9, using the expressions in Eqs. (8) and (9c). The overlaid curve is the
energy transfer rate $T_{2,\mathrm{Z}}$ based on triad interactions where one
of the triad members is a fluctuation with $k_{x}\mathrm{=}0$, channeling
energy to the stable mode. This is found by restricting the summation in Eq.
(9c) such that either $\mathbf{k}^{\prime}\mathrm{=}(0,k_{y}^{\prime})$ or
$\mathbf{k}^{\prime\prime}\mathrm{=}(0,k_{y}^{\prime\prime})$. We note that
$k_{x}\mathrm{=}0$ has zero growth rate and zero frequency. This wavenumber,
and, in particular, the first Fourier mode number along the $y$-axis, with
$k_{x}\mathrm{=}0$, is observed to have a large $x$-component in the turbulent
flow, and may be called a zonal mode, analog to secondary instabilities in
fusion plasmas.[69, 70, 30] A difference, however, exists: here, these zonal
modes contribute to the Reynolds stress, as opposed to those of fusion
plasmas, which do not transport momentum and, hence, the contribution of the
zonal modes to the stress is null.[29] The main similarity is that both of the
systems have qualitatively similar anisotropic dispersion relation and both
have zonal modes whose linear frequencies are zeros in the inviscid/ideal
limit (when one considers small-scale dissipation coefficient, as well, then
the complex frequencies feature damping rates).[30, 29]

We remark that the stable modes are excited in our system primarily via
nonlinear interaction with the $k_{x}\mathrm{=}0$ modes, once the initial
transient phase is complete by around $t\mathrm{=}100$ (Fig. 9). Since the
same triadic interaction continuously feeds energy to the stable mode in the
fully saturated phase, one may infer that the process of nonlinear formation
of the stable mode is rather coherent, and, to first order, does not involve
changing wavevectors; the triad driving the stable mode is largely the same.
Such an understanding of stable-mode formation is difficult to obtain in
spectral-only representation, as at a given KH-unstable wavenumber, there also
exists the stable mode, along with the continuum modes. This does not mean
that it is uninsightful to perform traditional diagnostic of energy transfer
in wavenumber space by integrating fluctuations along the direction of
inhomogeneity. In fact, some recent shear-flow studies have shown an
interesting finding of energy circulation over angles on the wavevector plane,
with the magnitude of the wavevector largely unchanged—a process that has been
dubbed as “transverse cascade.” [10, 11, 12, 13, 14, 15] One may investigate
the connection between such a cascade and the stable modes to learn if the
stable modes favor or impede the transverse cascade. Such a study, however, is
beyond the scope of the present paper.

The zero-frequency mode coupling leads to a near-resonant energy transfer to
the stable mode,[29, 30, 36, 37] see insets of Fig. 9, where both time-steady
and slow time-evolution of energy transfers are dominantly zonal-coupled. Such
a coupling allows a large amount of energy to be transferred coherently from
the unstable mode. This finding is of striking consequence, as it has a direct
potential to allow building a reduced, predictive model of saturation levels
of unstable and stable modes, as has recently been achieved in the context of
fusion microturbulence.[29, 30, 36, 71, 37, 72] The existence of zero-
frequency mode coupling and near-resonant energy transfer in $3$D hydrodynamic
turbulence, as found here, is surprising,[16] and may be considered a finding
of the foremost importance.

Figure 9: Nonlinear energy transfer $T_{2}$ to the stable mode, largely
channeled via spanwise-only (zero-frequency) fluctuation $T_{2,\mathrm{Z}}$,
which becomes dominant after $t\gtrsim 100$, shown in the yellow-shaded
region. The coupling to zero-frequency modes, labeled zonal modes, is further
probed using a power spectrum $P(\omega)$ of energy transfer in frequency
$\omega$, shown in the insets, with linear scales (left) and with logarithmic
scales (right), where
$T_{2,\mathrm{NZ}}\mathrm{=}T_{2}\mathrm{-}T_{2,\mathrm{Z}}$ represents a
purely non-zonal transfer. At $\omega\mathrm{=}0$, filled markers are used.
Zonal-coupling dominates at $|\omega|\mathrm{\approx}0$, implying that
nonlinear interactions feeding energy to stable modes are primarily coherent.

##  VIII Vortex stretching, cascade, and viscous dissipation

After quantifying the effect of stable modes in energetics and transport, we
now perform an ensemble of numerical experiments.

We randomly select five snapshots in the turbulent phase of the standard
simulation to construct an ensemble. In this ensemble, we project the state
vector at each of those times onto the stable eigenmodes at wavenumbers
$\mathbf{k}\mathrm{=}(0.2,0)$ and $\mathbf{k}\mathrm{=}(0.2,\pm 0.2)$, as
these wavenumbers dominate energetically (Fig. 7). Then, we delete these
stable modes from the chosen state vector, and use the resulting modified
state vector as an initial condition to start a new simulation. The restarted
simulations, therefore, produce responses to an impulsive zeroing of the
stable-mode amplitudes. For a second ensemble, we repeat this process,
projecting out this time the unstable modes from the same wavenumbers at the
same five snapshots. In both ensembles, the five responses of each impulsive
zeroing procedure are then averaged to produce two ensemble-averaged
responses, one from zeroing the stable modes and another from zeroing the
unstable modes.

We then evaluate the viscous dissipation rate, computed at wavenumbers other
than $(0.2,0)$ and $(0.2,\pm 0.2)$ so that the initial dissipation rates in
the mode-removed and standard simulations are identical. The time-evolution of
the dissipation rate $\Delta\epsilon_{\nu}$ is compared in Fig. 10 between the
two ensembles. Larger small-scale dissipation rates are observed when stable
modes are impulsively removed, similar to a recent finding in $2$D
turbulence.[25] Although removing the unstable modes initially causes a minor
positive change in dissipation rate, the real impact of the unstable-mode
removal is quickly seen as a reduced dissipation rate—lower than in the
standard simulation. This result is straightforward to explain: first, when an
eigenmode of significant amplitude is removed, the resulting large-scale flow
structure for a brief moment nonlinearly strains the other scales of the flow.
Second, nota bene, the linear interaction of the unstable modes (energy
source) or stable modes (energy sink) with the mean flow is impulsively
impaired, which immediately results in either a depletion of fluctuation
energy or a surplus. This second process dominates over the first one after
$\Delta t\approx 3$. Soon after $\Delta t\approx 5$, we observe different
behaviors in the dissipation rate—a set of curves lying in the positive
$\Delta\epsilon_{\nu}$ region, and the other set lying in the negative. The
negative is because the stable modes efficiently return the fluctuation energy
to the mean flow; the positive change is due to the enhanced forward cascade
of energy that is extracted from the mean flow by the unstable modes, which
are not countered by the stable modes here. Around $\Delta t\mathrm{=}5$, it
is observed that the stable-mode removed simulations, on average, feature the
highest dissipation rate, whereas the unstable-mode removed simulations, on
average, show a change in the dissipation rate that is nearly null.

Figure 10: Percentage changes in small-scale viscous dissipation rates. The
mean (solid curve) and one standard deviation (shaded) is found by simulating
an ensemble of $10$ different $3$D simulations, where either stable or
unstable modes are deleted instantaneously at a randomly selected time, and
the simulation resumed to measure the effect on the nonlinear cascade via the
small-scale dissipation rate. For all $\Delta t$ shown, stable-mode-removed
simulations (orange) show higher dissipation rates than that in a standard
simulation, corresponding to a positive $\Delta\epsilon_{\nu}$; lower
dissipation rates are observed after an initial impulsive transient in
unstable-mode-removed simulations (blue). All simulations were performed with
$Re\mathrm{=}300$.

Figure 11: Isocontours of squared vorticity in simulations at a time when the
viscous dissipation rate peaks. The same isocontour level is chosen in all
panels. In the standard simulation in (a), vortex stretching is suppressed
compared to a stable-mode-removed simulation in (b), where thin and elongated
filamentary vortices are prominent. In an unstable-mode-removed simulation in
(c), the long filamentary structures are not as pronounced as in (b). The
unstable modes, in the absence of stable modes, rapidly stretch and thin out
the vortex tubes; in contrast, the stable modes, in the absence of unstable
modes, deplete the fluctuation energy. Figure 12: Time evolution of enstrophy
$E_{\omega}=\langle|\hat{\mathbf{\omega}}|^{2}\rangle_{y,z}$ spectra in three
simulation-continuations: stable-mode removed, unstable-mode removed, and
standard. (a) All simulations are restarted with an identical initial
condition, except for the removed mode. The mentioned-modes are removed only
from the wavenumbers $k_{x}=0.2$ and $k_{y}=\\{0,\pm 0.2\\}$. The spectra are,
therefore, shown for $k_{x}>0.2$. (b)–(d) Impulsive responses to the mode-
removal are most pronounced when the unstable modes are removed, leading to
enhanced turbulence at all scales. The $k_{y}$-spectra of $(x,z)$-averaged
enstrophy are similar (not shown).

After quantitative analysis of dissipation rates, we show visualizations of
squared vorticity (whose volume integral is proportional to the dissipation
rate) and compare the physical structures of turbulence across the simulated
responses. At the time when the dissipation rate is maximal, the isocontours
of squared vorticity show, in a stable-mode impulsively-zeroed simulation, an
uninhibited stretching of vortex tubes by the unstable modes, compare Fig.
11(b) with Fig. 11(a). In the unstable-mode impulsively-zeroed simulation in
Fig. 11(c), a large-scale structure tends to develop instead of vortex
stretching that forms prominent small scales. We note that the vortex-tube
stretching is a purely $3$D phenomenon, and is associated with the forward
cascade of energy in $3$D hydrodynamic turbulence. A movie, found in the
online supplemental material of this paper, shows several snapshots of stable-
mode-removed simulations, where we observe prominent episodes of vortex
stretching.

To quantify the changes in vortex dynamics at smaller scales in response to
the large-scale-eigenmode removal, we measure enstrophy spectra in three
simulation-continuations, and show a comparison of their time evolution in
Fig. 12. Although started with an identical initial condition, the simulation
with stable-mode-removal rapidly evolves, as a nonlinear response to the mode-
deletion, to feature increased level of enstrophy at all scales. The unstable-
mode removed simulation also evolves; however, the increment in the small-
scale activity, early on, is relatively small, and, at later times, the
simulation, in fact, shows decreased enstrophy.

One may interpret the findings of the simulation experiments—increased
filamentary vortex structures due to the unstable modes, and reduced
filamentary vortex structures due to the stable modes—as the competing effects
of forward transfer of energy from the mean flow to the fluctuation by the
unstable mode and inverse transfer of energy by the stable mode in reversed
direction. The nonlinear process then responds to this linear physics.

##  IX Discussion

Some broader implications and interpretations of the results of this work
shall now be provided.

###  IX.1 Cascades, Couplings, and Competitions

We have shown with various quantitative measures that only a fraction of
unstable-mode energy is cascaded to small scales. A majority, near or greater
than $70\%$, of the instability-extracted energy from the mean flow to the
large-scale fluctuations is, in reality, nonlinearly transferred to the
conjugate-stable modes. This nonlinear mode-coupling at the instability-scale
controls the rate of the small-scale energy cascade right at its inception.

It remains true that larger Reynolds numbers push the Kolmogorov dissipation
length scale to smaller scales, thus allowing more scales to be dynamically
relevant in the turbulence. The small-scale turbulence is indeed affected by
the Reynolds number; however, the measures of interest, such as the energy
injection rate in the cascade channel where it is formed and the momentum
transport rate, are under the territorial jurisdiction of the fluctuations at
the instability scale. At such a large scale, the conjugate-stable modes are
excited, first via a parametric drive from the unstable modes through the
mode-coupling coefficient $C_{211}(\mathbf{k},\mathbf{k}^{\prime})$ in the
stable-mode evolution equation

$\partial_{t}\beta_{2}(\mathbf{k})=-\gamma(\mathbf{k})\beta_{2}(\mathbf{k})+\sum_{l}B_{2l}(\mathbf{k})\beta_{l}(\mathbf{k})+\sum_{\mathbf{k}^{\prime},\mathbf{k}^{\prime\prime}:\mathbf{k}^{\prime}+\mathbf{k}^{\prime\prime}=\mathbf{k}}C_{211}(\mathbf{k},\mathbf{k}^{\prime})\beta_{1}^{\prime}\beta_{1}^{\prime\prime}+\ldots,$
(10)

where $C_{211}(\mathbf{k},\mathbf{k}^{\prime})$ is independent of viscosity,
and, therefore, favors excitation of stable modes. Whether this excitation
becomes significant or not depends on the relative magnitudes of $\gamma$,
$B_{2l}$, and $C_{211}(\mathbf{k},\mathbf{k}^{\prime})$.[28]

At the large instability scale, with decreasing viscosity, $B_{2l}$ naturally
becomes smaller; this is different from the singular limit of small-scale
dissipation that does not vanish with decreasing viscosity. Therefore, the
large-scale-stable-mode excitation is only weakly impacted once the viscosity
is sufficiently low, as we have found in Fig. 3. In addition, the small-scale
dissipation rate, as well as the scales where intense viscous dissipation
occurs, are both controlled by the energy extracted from the mean flow-
gradient; the extraction rate, however, depends on the excitation levels of
the large-scale stable modes, as shown in Figs. 10 and 11. Thus, the Reynolds
number, although it spawns more small scales, only asymptotically and weakly
impacts the large-scale-stable-mode physics, as the stable modes are coupled
to the unstable modes via a nonlinear mode-coupling coefficient whose origin
lies in the overlap of the eigenmodes of an instability that is entirely
inviscid. Although a competition exists between small-scale cascades and
nonlinear coupling among the roots of the inviscid instability, the latter
dominates over the former as long as the Reynolds number is not so small that
the flow becomes near-laminar, strongly breaking the symmetry between the
unstable and stable modes—this latter case does not appear in high-Reynolds-
number turbulent flows in natural systems.

###  IX.2 Dimensionality, Density of states, Drives, and Dynamics

For a system with $n$ spatial dimensions (physical coordinate axes), we
conjecture that the number of KH-unstable wavenumbers grows following an
approximate power law $(k/k_{0})^{n-1}$, where $k$ is the magnitude of the
wavevector, and $k_{0}\mathrm{=}2\pi/L$ is the lowest wavenumber in the domain
of size $L$. The increment in the number of KH-unstable wavenumbers follows
from the density-of-states effect, associated with each dimension (Fig. 2);
one of the axes is not counted because that axis represents the direction of
inhomogeneity of the mean shear flow.

Such increased number of the KH-unstable wavenumbers put the mean shear-flow
under additional stress, all demanding energy and momentum redistribution.
Turbulence with increased spatial dimensions, hence, extracts more energy from
the mean flow and transports more momentum across the shear layer—i.e., the
turbulence features increased energy drive and dynamics. In decaying
turbulence, the slowly growing wavenumbers of higher spatial dimension
(greater than two) may not get an opportunity to attain high amplitudes in
comparison to the $2$D-mode amplitudes. However, for a sufficiently forced
mean-flow, such a discriminatory scenario is precluded; as evidenced in Figs.
5 and 8, the larger density of states of $3$D modes manifests in increased
levels of momentum transport and energy transfer.

The dynamics of shear-flow turbulence can be substantially different depending
on the linear process driving the turbulence. We shall remark here on two
important processes: the exponential instability growth and the transient non-
modal growth. When the Kelvin-Helmholtz instability is present, as in this
paper, the instability rapidly grows exponentially [Figs. 3(a) and (c), blue
curves]. Such a rapid exponential growth over longer times supersedes the
here-inconsequential transient growth due to non-modal effects.[73, 74] The
latter are crucial in modally stable non-normal flows,[75, 76, 77] e.g., shear
flows with no inflection point. With the fluctuations dominated by the
unstable modes, they nonlinearly excite conjugate-stable modes, with nearly
twice[28] the instability-growth rate [Figs. 3(a) and (c), orange curves].
Such stable modes, which are excited by the nonlinearity, are fundamentally
different from the linear, transient growth of fluctuations. We also note that
the conjugate-stable modes here continuously transfer energy from the
fluctuations to the mean flow, as seen in Fig. 6(c). This continuous process
differs from the transient non-modal growth.

##  X Conclusions

The ubiquitous Kelvin-Helmholtz instability is traditionally assumed to
saturate nonlinearly by cascading to small scales the unstable-mode energy, in
its entirety. Recently, a majority of the unstable-mode energy has been
reported to be returned to the mean flow from the fluctuation by the linearly
stable eigenmodes that are nonlinearly excited in two-dimensional
turbulence.[22, 25, 27] Whether such a consequential process exists in $3$D
fluid turbulence, as well, is the central question addressed in this work.
Here, the linearly stable $3$D modes, excited to significant amplitudes, are
shown to exhibit similarities with the reports of $2$D turbulence, as well as
new signatures that are fundamentally $3$D. The energy transfer path to stable
modes in $3$D is wholly different: in the turbulent phase, the transfer occurs
via modes, we label as zonal flows, that vary only in the direction orthogonal
to the $2$D shear flow. Such perturbations have zero linear frequencies due to
the anisotropic dispersion relation of the inviscid $3$D Kelvin-Helmholtz
instability. However, because of their zero frequencies, they allow near-
resonant energy transfer from unstable to conjugate-stable modes, saturating
the instability in $3$D.

The transfer of energy by the stable modes is directly from the instability-
scale fluctuation to the mean flow-gradient, which may be interpreted as an
inverse transfer of energy. This inversion is a linear process, different from
the traditional nonlinear energy cascade to small scales. We have shown here
not only the existence of such inverse transfer of energy in $3$D, but that
they are more efficient for the $3$D perturbations than for the $2$D
perturbations.

Vortex-tube stretching, an inherently $3$D process, is suppressed by the
stable modes. This has been demonstrated through a series of numerical
experiments, where we impulsively zeroed either stable or unstable modes in a
nonlinear simulation, at randomly selected times, and resumed the simulations
to obtain statistical responses in the small-scale activity. The viscous
dissipation rates and spectral amplitudes of enstrophy increase with stable-
mode removal, and decrease with unstable-mode removal. Further, $3$D
visualizations of vortex dynamics revealed that, in the absence of the stable
modes, the vortex tubes are significantly stretched and continuously thinned
out; the stable modes counter this mechanism and remove the thin filamentary
structures, attempting to form large-scale structures.

Transporting up-gradient momentum in $3$D, the stable modes counteract the
instability-driven down-gradient momentum transport more effectively than they
do in 2D ($\mathrm{\approx}70\%$ vs. $\mathrm{\approx}50\%$). This competition
occurs continuously in time and can occasionally reverse the net transport
direction, bearing important implications for shear-flow turbulence in geo-
and astrophysical environments.

For future work, this study raises important questions regarding the impact of
stable modes in other three-dimensional problems such as $3$D MHD shear flows,
where, analogous to the vortex stretching, there exists magnetic field
amplification, which, at small scales, can be inhibited by the stable modes,
favoring large-scale magnetic field generation[78]; $3$D stratified shear
flows, where the stable modes may serve as a zookeeper, managing a zoo of
instabilities[79, 44, 42]; the ideal, magneto-rotational instability
(MRI),[80] where it is possible that the transport measured in numerical
simulations of MRI-driven turbulence (see, e.g., Ref. [64]) has an unaccounted
contribution from the stable modes, and there is an opportunity to build
accurate transport models for astrophysical problems, magnetized[81, 82] or
not.

###### Acknowledgements.

We are grateful to Keaton Burns, Evan H. Anders, Daniel Lecoanet, and Michael
J. Gerard for useful discussions. This material is based upon work funded by
the Department of Energy (Grant No. DE-SC0022257) through the NSF/DOE
Partnership in Basic Plasma Science and Engineering. A.E.F. acknowledges
support from NASA HTMS Grant No. 80NSSC20K1280, and from the George Ellery
Hale Postdoctoral Fellowship in Solar, Stellar and Space Physics at the
University of Colorado, Boulder. The simulations used the ACCESS
supercomputing resources under Allocation No. TG-PHY130027. The data that
support the findings of this study are available from the corresponding author
upon reasonable request.

## References

  * [1] A. Alexakis and L. Biferale, Cascades and transitions in turbulent flows, Phys. Rep. 767, 1 (2018).
  * [2] A.N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers Dokl. Akad. Nauk S.S.S.R. 30, 299 (1941).
  * [3] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov (Cambridge University Press, Cambridge, England, 1995).
  * [4] R.H. Kraichnan, Inertial Ranges in Two-Dimensional Turbulence, Phys. Fluids 10, 1417 (1967).
  * [5] C.E. Leith, Diffusion Approximation for Two-Dimensional Turbulence, Phys. Fluids 11, 671 (1968).
  * [6] F. Waleffe, Inertial transfers in the helical decomposition, Phys. Fluids A 5, 677 (1993).
  * [7] L.M. Smith and F. Waleffe, Transfer of Energy to 2D Large Scales in Forced, Rotating 3D Turbulence, Phys. Fluids 11, 1608 (1999).
  * [8] L.M. Smith and F. Waleffe, Generation of Slow, Large Scales in Forced, Rotating, Stratified Turbulence, J. Fluid Mech. 451, 145 (2002).
  * [9] J.J. Riley and M.-P. Lelong, Fluid motions in the presence of strong stable stratification, Annu. Rev. Fluid Mech. 32, 613 (2000).
  * [10] W. Horton, J.-H. Kim, G.D. Chagelishvili, J.C. Bowman, and J.G. Lominadze, Angular redistribution of nonlinear perturbations: A universal feature of nonuniform flows, Phys. Rev. E 81, 066304 (2010).
  * [11] G.R. Mamatsashvili, D.Z. Gogichaishvili, G.D. Chagelishvili, and W. Horton, Nonlinear transverse cascade and two-dimensional magnetohydrodynamic subcritical turbulence in plane shear flows, Phys. Rev. E 89, 043101 (2014).
  * [12] G. Khujadze, G. Chagelishvili, S. Dong, J. Jiménez, and H. Foysi, Dynamics of homogeneous shear turbulence: A key role of the nonlinear transverse cascade in the bypass concept, Phys. Rev. E 94, 023111 (2016).
  * [13] D.Z. Gogichaishvili, G.R. Mamatsashvili, W. Horton, G.D. Chagelishvili, and G. Bodo, Nonlinear Transverse Cascade and Sustenance of MRI Turbulence in Keplerian Disks with an Azimuthal Magnetic Field, Astrophys. J.  845, 70 (2017).
  * [14] D.Z. Gogichaishvili, G.R. Mamatsashvili, W. Horton, and G.D. Chagelishvili, Active Modes and Dynamical Balances in MRI Turbulence of Keplerian Disks with a Net Vertical Magnetic Field, Astrophys. J.  866, 134 (2018).
  * [15] G.R. Mamatsashvili, G.D. Chagelishvili, M.E. Pessah, F. Stefani, and G. Bodo, Zero Net Flux MRI Turbulence in Disks: Sustenance Scheme and Magnetic Prandtl Number Dependence, Astrophys. J.  904, 47 (2020).
  * [16] D. Biskamp and A. Zeiler, Nonlinear Instability Mechanism in 3D Collisional Drift-Wave Turbulence, Phys. Rev. Lett. 74, 706 (1995).
  * [17] C. Ng and A. Bhattacharjee, Interaction of Shear-Alfven Wave Packets: Implication for Weak Magnetohydrodynamic Turbulence in Astrophysical Plasmas, Astrophys. J. 465, 845 (1996).
  * [18] P.W. Terry, Inverse energy transfer by near-resonant interactions with a damped-wave spectrum, Phys. Rev. Lett. 93, 235004 (2004).
  * [19] S. Du, H. Li, X. Fu, and Z. Gan, Anisotropic Energy Transfer and Conversion in Magnetized Compressible Turbulence, Astrophys. J. 948, 72 (2023).
  * [20] A.M. Rubio, K. Julien, E. Knobloch, and J.B. Weiss, Upscale energy transfer in three-dimensional rapidly rotating turbulent convection, Phys. Rev. Lett. 112, 144501 (2014).
  * [21] C. Guervilly, D. Hughes, and C.A. Jones, Large-scale vortices in rapidly rotating Rayleigh-Benard convection, J. Fluid Mech. 758, 407 (2014).
  * [22] 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).
  * [23] 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).
  * [24] 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).
  * [25] 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, Phys. Plasmas 29, 070701 (2022).
  * [26] B. Tripathi, A.E. Fraser, P.W. Terry, E.G. Zweibel, and M.J. Pueschel, Near-cancellation of up- and down-gradient momentum transports in magnetized shear flow turbulence due to stable modes, Phys. Plasmas 29, 092301 (2022).
  * [27] B. Tripathi, A.E. Fraser, P.W. Terry, E.G. Zweibel, M.J. Pueschel, and E.H. Anders, Nonlinear mode coupling and energetics of driven magnetized shear-flow turbulence, Phys. Plasmas 30, 072107 (2023).
  * [28] P.W. Terry, D.A. Baver, and S. Gupta, Role of stable eigenmodes in saturated local plasma turbulence, Phys. Plasmas 13, 022307 (2006).
  * [29] 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).
  * [30] 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).
  * [31] 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).
  * [32] G.G. Whelan, M.J. Pueschel, and P.W. Terry, Nonlinear Electromagnetic Stabilization of Plasma Microturbulence, Phys. Rev. Lett.  120, 175002 (2018).
  * [33] 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).
  * [34] 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).
  * [35] D.R. Hatch, P.W. Terry, F. Jenko, F. Merz, M.J. Pueschel, W.M. Nevins, and E. Wang, Role of subdominant stable modes in plasma microturbulence, Phys. Plasmas 18, 055706 (2011).
  * [36] P.-Y Li, P.W. Terry, G.G. Whelan, and M.J. Pueschel, Saturation physics of threshold heat-flux reduction, Phys. Plasmas 28, 102507 (2021).
  * [37] P.-Y Li and P.W. Terry, Assessing physics of ion temperature gradient turbulence via hierarchical reduced-model representations, Phys. Plasmas 29, 042301 (2022).
  * [38] T.M. Qian and M.E. Mauel, Observation of weakly damped modes using high resolution measurement of turbulence in a dipole confined plasma, Phys. Plasmas 27, 014501 (2020).
  * [39] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Clarendon Press, Oxford 1961).
  * [40] P. Drazin and W. Reid, Hydrodynamic Stability (Cambridge University Press, Cambridge, 2004).
  * [41] M. Buzzicotti, H. Aluie, L. Biferale, and M. Linkmann, Energy transfer in turbulence under rotation, Phys. Rev. Fluid 3, 034802 (2018).
  * [42] C. Liu, A. Kaminski, and W. Smyth, The butterfly effect and the transition to turbulence in a stratified shear layer, J. Fluid Mech. 953, A43 (2022).
  * [43] G.P. Klaassen and W.R. Peltier, The onset of turbulence in finite-amplitude Kelvin-Helmholtz billows, J. Fluid Mech. 155, 1 (1985).
  * [44] G.P. Klaassen and W.R. Peltier, The influence of stratification on secondary instability in free shear layers, J. Fluid Mech. 227, 71 (1991).
  * [45] A. Mashayek and W.R. Peltier, The ‘zoo’ of secondary instabilities precursory to stratified shear flow transition. Part 1 Shear aligned convection, pairing, and braid instabilities, J. Fluid Mech. 708, 5 (2012).
  * [46] A. Mashayek and W.R. Peltier, The ‘zoo’ of secondary instabilities precursory to stratified shear flow transition. Part 2 The influence of stratification, J. Fluid Mech. 708, 45 (2012).
  * [47] H. Salehipour, W. Peltier, and A. Mashayek, Turbulent diapycnal mixing in stratified shear flows: The influence of Prandtl number on mixing efficiency and transition at high Reynolds number, J. Fluid Mech. 773, 178 (2015).
  * [48] K.M. Smith, C.P. Caulfield, and J.R. Taylor, Turbulence in forced stratified shear flows, J. Fluid Mech. 910, A42 (2021).
  * [49] J.B. Marston, E. Conover, and T. Schneider, Statistics of an Unstable Barotropic Jet from a Cumulant Expansion, J. Atmos. Sci. 65, 1955 (2008).
  * [50] L. Cope, P. Garaud, and C. Caulfield, The dynamics of stratified horizontal shear flows at low Péclet number, J. Fluid Mech. 903, A1 (2020).
  * [51] P. Garaud and L. Kulenthirarajah, Turbulent Transport in a Strongly Stratified Forced Shear Layer with Thermal Diffusion, Astrophys. J.  821, 49 (2016).
  * [52] D. Lucas, C.P. Caulfield, and R.R. Kerswell, Layer formation in horizontally forced stratified turbulence: connecting exact coherent structures to linear instabilities, J. Fluid Mech 832, 409 (2017).
  * [53] C.P. Caulfield, Open questions in turbulent stratified mixing: Do we even know what we do not know?, Phys. Rev. Fluids 5, 110518 (2020).
  * [54] P. Garaud, Double-diffusive convection at low Prandtl number, Annu. Rev. Fluid Mech. 50, 275 (2018).
  * [55] A.E. Fraser, S.A. Reifenstein, and P. Garaud, Magnetized fingering convection in stars, arXiv (2023). arXiv:2302.11610.
  * [56] A. Masson and K. Nykyri, Kelvin-Helmholtz Instability: Lessons Learned and Ways Forward, Space Sci. Rev 214, 71 (2018).
  * [57] 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).
  * [58] M.J. Pueschel, F. Jenko, D. Told, and J. Büchner, Gyrokinetic simulations of magnetic reconnection, Phys. Plasmas 18, 112102 (2011).
  * [59] 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).
  * [60] 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).
  * [61] E. Paul, I. Abel, M. Landreman, and W. Dorland, An adjoint method for neoclassical stellarator optimization, J. Plasma Phys. 85, 795850501 (2019).
  * [62] P.W. Terry, D.A. Baver, and D.R. Hatch, Reduction of inward momentum flux by damped eigenmodes, Phys. Plasmas 16, 122305 (2009).
  * [63] J. Fuller, A.L. Piro, and A.S. Jermyn, Slowing the spins of stellar cores, Mon. Not. R. Astron. Soc. 485, 3661 (2019).
  * [64] 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).
  * [65] J. Goodman and G. Xu, Parasitic Instabilities in Magnetized, Differentially Rotating Disks, Astrophys. J.  432, 213 (1994).
  * [66] A.J. Barker, C.A. Jones, and S.M. Tobias, Angular momentum transport by the GSF instability: non-linear simulations at the equator, Mon. Not. R. Astron. Soc. 487, 1777 (2019).
  * [67] M.K. Verma, Energy Transfers in Fluid Flows: Multiscale and Spectral Perspectives (Cambridge University Press, Cambridge 2019).
  * [68] S. Gomé, L. Tuckerman, and D. Barkley, Patterns in transitional shear turbulence. Part 1. Energy transfer and mean-flow interaction, J. Fluid Mech. 964, A16 (2023).
  * [69] B.N. Rogers, W. Dorland, and M. Kotschenreuther, Generation and Stability of Zonal Flows in Ion-Temperature-Gradient Mode Turbulence, Phys. Rev. Lett.  85, 5336 (2000).
  * [70] M.J. Pueschel, T. G$\ddot{\mathrm{o}}$rler, F. Jenko, D.R. Hatch, and A.J. Cianciara, On secondary and tertiary instability in electromagnetic plasma microturbulence, Phys. Plasmas 20, 102308 (2013).
  * [71] 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).
  * [72] P.-Y. Li, M.J. Pueschel, P.W. Terry, and G.G Whelan, On the role of mode resonances in regulating zonal-flow-moderated plasma microturbulence, Nucl. Fusion 63, 026028 (2023).
  * [73] C. Arratia, C.P. Caulfield, and J.M. Chomaz, Transient perturbation growth in time-dependent mixing layers, J. Fluid Mech. 717, 90 (2013).
  * [74] A.K. Kaminski, C.P. Caulfield, and J.R. Taylor, Transient growth in strongly stratified shear layers, J. Fluid Mech. 758, R4 (2014).
  * [75] L.H. Gustavsson, Energy growth of three-dimensional disturbances in plane Poiseuille flow, J. Fluid Mech. 224, 241 (1991).
  * [76] G.D. Chagelishvili, A.G. Tevzadze, G. Bodo, and S.S. Moiseev. Linear mechanism of wave emergence from vortices in smooth shear flows, Phys. Rev. Lett.  79, 17 (1997).
  * [77] G.D. Chagelishvili, J-P. Zahn, A.G. Tevzadze, and J.G. Lominadze, On hydrodynamic shear turbulence in Keplerian disks: Via transient growth to bypass transition, Astron. Astrophys. 402, 401 (2003).
  * [78] S.M. Tobias and F. Cattaneo, Shear-driven dynamo waves at high magnetic Reynolds number, Nature 497, 463 (2013).
  * [79] C.P. Caulfield, Layering, Instabilities, and Mixing in Turbulent Stratified Flows, Annu. Rev. Fluid Mech. 53, 113 (2021).
  * [80] J.S. Oishi, G.M. Vasil, M. Baxter, A. Swan, K.J. Burns, D. Lecoanet, and B.P. Benjamin, The magnetorotational instability prefers three dimensions, Proc. R. Soc. A. 476, 20190622 (2020).
  * [81] S. Ji, J. Fuller, and D. Lecoanet, Magnetohydrodynamic simulations of the Tayler instability in rotating stellar interiors, Mon. Not. R. Astron. Soc. 521, 5372 (2022).
  * [82] R. Blandford and N. Globus, Ergomagnetosphere, ejection disc, magnetopause in M87 - I. Global flow of mass, angular momentum, energy, and current, Mon. Not. R. Astron. Soc. 514, 5141 (2022).