# Predicting the Slowing of Stellar Differential Rotation by Instability-
Driven Turbulence

B. Tripathi Department of Physics, University of Wisconsin–Madison, Madison,
Wisconsin 53706, USA<EMAIL_ADDRESS>A.J. Barker Department of Applied
Mathematics, School of Mathematics, University of Leeds, Leeds LS2 9JT, UK
<EMAIL_ADDRESS>A.E. Fraser Department of Applied Mathematics,
University of Colorado, Boulder, Colorado 80309, USA  Department of
Astrophysical and Planetary Sciences, University of Colorado, Boulder,
Colorado 80309, USA  Laboratory for Atmospheric and Space Physics, University
of Colorado, Boulder, Colorado 80309, USA  P.W. Terry Department of Physics,
University of Wisconsin–Madison, Madison, Wisconsin 53706, USA  E.G. Zweibel
Department of Physics, University of Wisconsin–Madison, Madison, Wisconsin
53706, USA  Department of Astronomy, University of Wisconsin–Madison, Madison,
Wisconsin 53706, USA

###### Abstract

Differentially rotating stars and planets transport angular momentum
internally due to turbulence at rates that have long been a challenge to
predict reliably. We develop a self-consistent saturation theory, using a
statistical closure approximation, for hydrodynamic turbulence driven by the
axisymmetric Goldreich–Schubert–Fricke (GSF) instability at the stellar
equator with radial differential rotation. This instability arises when fast
thermal diffusion eliminates the stabilizing effects of buoyancy forces in a
system where a stabilizing entropy gradient dominates over the destabilizing
angular momentum gradient. Our turbulence closure invokes a dominant three-
wave coupling between pairs of linearly unstable eigenmodes and a near-zero
frequency, viscously damped eigenmode that features latitudinal jets. We
derive turbulent transport rates of momentum and heat, and provide them in
analytic forms. Such formulae, free of tunable model parameters, are tested
against direct numerical simulations; the comparison shows good agreement.
They improve upon prior quasi-linear or ‘‘parasitic saturation” models
containing a free parameter. Given model correspondences, we also extend this
theory to heat and compositional transport for axisymmetric thermohaline
instability-driven turbulence in certain regimes.

Astrophysical fluid dynamics (101) — Solar differential rotation (1996) —
Stellar rotation (1629) — Stellar interiors (1606) — Hydrodynamics (1963) —
Extrasolar gaseous giant planets (509)

##  1 Introduction

Figure 1: (Left) A schematic diagram of a differentially rotating star with a
radial shear, gravity, and stable stratification. Such a system subject to the
GSF instability is studied using a local Cartesian model. (Right) Snapshots of
velocity components $\widetilde{u}_{x}(x,z)$, $\widetilde{u}_{y}(x,z)$, and
$\widetilde{u}_{z}(x,z)$ from the axisymmetric GSF instability-driven
turbulence; $x,y,$ and $z$ represent the local radial, azimuthal, and
latitudinal directions. Though finger-like horizontal structures (as shown by,
e.g., $\widetilde{u}_{x}$) grow the fastest in the linear phase
($t\mathrm{=}100$), strong latitudinal jets $\widetilde{u}_{z}$ are generated
nonlinearly ($t\mathrm{=}6000$). The color bar for $t\mathrm{=}100$ is shared
by $\widetilde{u}_{x}$, $\widetilde{u}_{y}$, and $3\widetilde{u}_{z}$; the
color bar for $t\mathrm{=}6000$ is shared by $3\widetilde{u}_{x}$,
$3\widetilde{u}_{y}$, and $\widetilde{u}_{z}$. The turbulent transport of
angular momentum, e.g., $\langle\widetilde{u}_{x}\widetilde{u}_{y}\rangle$, is
predicted in this paper using a jet-coupled turbulence closure.

Instability-driven turbulence is thought to play a major role in the transport
of angular momentum (AM), heat and composition in stellar and planetary
interiors (see, e.g., Garaud, 2018; Aerts et al., 2019; Spruit, 2002; Fuller
et al., 2019), as well as in astrophysical disks (e.g., Balbus & Hawley, 1998;
Lesur et al., 2023). Unfortunately, rates of turbulent transport are very
challenging to predict theoretically, and the lack of reliable theories has
hampered our understanding of the evolution of stellar and planetary internal
rotations and structures. For example, the AM redistribution in red giant
stars is currently poorly understood, and their core-envelope differential
rotations inferred from asteroseismology have not been adequately explained
(e.g., Beck et al., 2012; Eggenberger et al., 2012; Aerts et al., 2019). The
nearly solid-body rotation observed in the solar radiative interior also lacks
a robust explanation (e.g., Garaud & Garaud, 2008; Wood & McIntyre, 2011).

Differential rotation is known to drive a variety of hydrodynamic (and
hydromagnetic) instabilities. In this paper, we focus on modeling hydrodynamic
instabilities of differential rotation in stellar and planetary radiative
zones, and in particular on the Goldreich–Schubert–Fricke (GSF)
instability111It has also been referred to as the “Vertical Shear Instability”
(VSI) in accretion disks (e.g., Urpin & Brandenburg, 1998; Nelson et al.,
2013; Barker & Latter, 2015; Latter & Papaloizou, 2018) and as “inertial
instability” enabled by thermal diffusion in stellar interiors (Park et al.,
2020, 2021).(Goldreich & Schubert, 1967; Fricke, 1968). This is a double-
diffusive centrifugal instability in which rapid thermal diffusion (relative
to viscous momentum diffusion) enables instability by tempering the otherwise
stabilising effects of buoyancy forces. Prior work has studied the linear and
nonlinear properties of the instability, and the turbulence it drives
(Knobloch & Spruit, 1982; Knobloch, 1982; Korycansky, 1991; Rashid et al.,
2008; Barker et al., 2019, 2020; Dymott et al., 2023), but a reliable theory
for the resulting turbulent transport is lacking. This means that the effects
of the GSF instability on stellar rotational and chemical evolution have not
been modeled in a self-consistent manner. Instead, one typically invokes
unexplained ‘‘additional viscosities” or models that contain free parameters.
Such tunable parameters are intended to describe the effects of turbulence on
AM transport for which adequate knowledge is lacking.

A fully analytic model containing no free parameters is derived here for the
GSF instability-driven turbulence in $2.5$ dimensions ($2.5$-D), i.e., with
all three components of velocity but varying spatially only in two dimensions.
The predictions of our analytical model are in broad agreement with detailed
numerical simulations of turbulence, driven by the axisymmetric ($2.5$-D) GSF
instability at the equator of a star with radial differential rotation. Such
model of the instability is, for certain diffusivity ratios and in $2.5$-D,
formally and nonlinearly equivalent to the thermohaline, or salt-finger,
instability that transports heat and chemical elements (Knobloch, 1982; Barker
et al., 2019); thus, the turbulent transport arising from axisymmetric
fingering convection is also described by our theory.

The structure of this paper is as follows. In § 2, we present our model and
methods of analysis. Nonlinear mode coupling and saturation diagnostics of the
instability appear in § 3. Informed by such diagnostics, we present analytical
formulae, without any free parameters, to model the turbulence and its
transport properties in § 4. We discuss the astrophysical implications and
conclude in § 5. Details of the closure model are provided in the Appendices.

##  2 The GSF instability and inertial-gravity waves

To study a basic mechanism of angular momentum transport in a differentially
rotating star, we consider a local region inside the star near its equator,
where the rotation can be split into a uniform or mean
part—$\bm{\Omega}\mathrm{=}\Omega\hat{\bm{\mathrm{e}}}_{z}$, aligning with the
local latitudinal axis $z$—and a non-uniform part due to the radial
differential rotation. The latter is represented by a background linear shear
flow $\bm{U}_{0}(x)\mathrm{=}-\mathcal{S}x\hat{\bm{\mathrm{e}}}_{y}$ where $x$
is the radial coordinate, $y$ is the azimuthal coordinate, and
$\mathcal{S}=-d\Omega_{\mathrm{Shell}}(x)/d\ln x$ is the local radial-shear
rate, with $\Omega_{\mathrm{Shell}}(x)$ representing the ‘‘Shellular” rotation
of the simplified star. A uniform gravity field with
$\bm{g}=-g\hat{\bm{\mathrm{e}}}_{x}$ is directed radially inward (Fig. 1). A
background radial temperature gradient $\nabla T_{0}$ then stratifies the
fluid density radially, with a thermal expansion coefficient $\alpha$. In such
a background state, any perturbations in velocity $\bm{u}$ and scaled
temperature $\theta\mathrm{=}\alpha gT$ evolve (Barker et al., 2019) as

$\displaystyle
D\bm{u}+\bm{u}\cdot\nabla\bm{U}_{0}+2\bm{\Omega}\times\bm{u}=-\nabla
p+\theta\hat{\bm{\mathrm{e}}}_{x}+\nu\nabla^{2}\bm{u},$ (1a) $\displaystyle
D\theta+\bm{u}\cdot\nabla\Theta_{0}=\kappa\nabla^{2}\theta,$ (1b)
$\displaystyle\nabla\cdot\bm{u}=0,$ (1c) $\displaystyle
D\equiv\partial_{t}+(\bm{u}+\bm{U}_{0})\cdot\nabla,$ (1d)

where the variables $p,\nu,$ and $\kappa$ are the fluid pressure (per unit
density), the kinematic viscosity, and the thermal diffusivity, respectively.
We also define the Prandtl number Pr$\mathrm{=}\nu/\kappa$. Because the GSF
instability operates at length scales much smaller than a pressure scale
height in stars, the local approximation is valid; in such a case, when the
turbulence drives subsonic flows, the Boussinesq approximation is also
appropriate (Spiegel & Veronis, 1960). Assuming a uniform temperature
gradient, a fluid element perturbed radially oscillates with a constant Brunt-
Väisälä (buoyancy) frequency $\mathcal{N}$, where
$\mathcal{N}^{2}\hat{\bm{\mathrm{e}}}_{x}\mathrm{=}\nabla\Theta_{0}\mathrm{=}\alpha
g\nabla T_{0}$.

Henceforth, we non-dimensionalize all variables using the characteristic
rotation time scale $\Omega^{-1}$ and length scale $d$, with
$d\mathrm{=}(\nu\kappa/\mathcal{N}^{2})^{1/4}$, which is typically similar to
the wavelengths of fastest-growing modes. Thus,
$N\mathrm{=}\mathcal{N}/\Omega$ is the dimensionless buoyancy frequency and
$S\mathrm{=}\mathcal{S}/\Omega$ the dimensionless shear rate (Rossby number).
The GSF instability occurs in low-Pr fluids whenever $r\in[0,1]$, where
$r=\mathrm{Pr}(1+N^{2}\kappa_{\mathrm{ep}}^{-2})/(\mathrm{Pr}-1)$, with
$\kappa_{\mathrm{ep}}\mathrm{=}\sqrt{2(2-S)}$ representing the dimensionless
epicyclic frequency (Barker et al., 2019).

###  2.1 Eigenmode analysis

A linear analysis of Eqs. (1a)–(1b) for axisymmetric (uniform-in-$y$)
perturbations yields a simple matrix equation, which upon Fourier-transforming
becomes $\bm{L}\hat{\bm{X}}\mathrm{=}\gamma\hat{\bm{X}}$, where
$\hat{\bm{X}}\mathrm{=}[\hat{u}_{x},\hat{u}_{y},\hat{u}_{z},\hat{\theta}]^{\mathrm{T}}$,
with T as the transpose operation, is the state vector of spatially Fourier-
transformed components at wavevector $\mathbf{k}\mathrm{=}(k_{x},k_{z})$; the
matrix $\bm{L}$ is a linear operator, whose eigenvalues are the complex-valued
growth rates $\gamma$. The size of $\bm{L}$ demands four linearly independent
eigenvectors. Because of the additional constraint
$\nabla\mathrm{\cdot}\bm{u}\mathrm{=}0$, the system has only three degrees of
freedom at any given wavenumber—two components of velocity, and
$\hat{\theta}$. Hence, one eigenvector among the four eigenvectors does not
satisfy $\nabla\mathrm{\cdot}\bm{u}\mathrm{=}0$ and is rejected. We confirm
that this eigenvector is not excited within our incompressible Boussinesq
simulations. One among the remaining three eigenvectors at a given wavevector
becomes GSF-unstable [$\mathrm{Re}\,(\gamma)>0$, where $\mathrm{Re}\,$ denotes
the real part], whenever $r\in[0,1)$. The remaining two eigenvectors are
always stable, and their eigenvalues are complex conjugates of each other
whenever they satisfy $\mathrm{Im}\,(\gamma)\neq 0$, where $\mathrm{Im}\,$
denotes the imaginary part; $\mathrm{Im}\,(\gamma)$ corresponds to the
frequency of inertial-gravity, or gravito-inertial, waves (IGWs), modified by
the shear flow and damped by viscous and thermal diffusion.

The GSF instability grows dominantly via axisymmetric ($\partial_{y}\equiv 0$)
perturbations, therefore we focus upon the $(x,z)$-variations of the
$3$-component velocity and temperature fields. The dispersion relation then is
a simple cubic polynomial in $\gamma$ (Goldreich & Schubert, 1967) as

$\gamma_{\nu}^{2}\gamma_{\kappa}+\frac{\kappa_{\mathrm{ep}}^{2}k_{z}^{2}}{k^{2}}\gamma_{\kappa}+\frac{N^{2}k_{z}^{2}}{k^{2}}\gamma_{\nu}=0,$
(2)

where $\gamma_{\nu}=\gamma+\nu k^{2}$ and $\gamma_{\kappa}=\gamma+\kappa
k^{2}$. Equation (2) shows that, on the ($k_{x},k_{z}$)-plane, the growth rate
exhibits strong anisotropy: fluctuations with $k_{x}\mathrm{=}0$ (‘‘elevator
modes”) grow the fastest, whereas those with $k_{z}\mathrm{=}0$ are linearly
stable. This observation is critical for the nonlinear saturation of the GSF
instability because the anisotropy of the linear physics, in particular the
$k_{z}\mathrm{=}0$ fluctuation, can impose its anisotropy on the nonlinear
energy transfer, which is otherwise isotropic; such consequential effects have
been found in various systems such as $3$-D Kelvin-Helmholtz instability
(Tripathi et al., 2023b), rotating (Waleffe, 1993; Smith & Waleffe, 1999) and
stably stratified turbulence (Riley & Lelong, 2000), turbulence with an
external magnetic field in astrophysical (Ng & Bhattacharjee, 1996; Du et al.,
2023) and fusion plasmas (Biskamp & Zeiler, 1995; Terry, 2004).

Using the complete basis provided by the eigenvectors of the linear operator
$\bm{L}$, we can decompose arbitrary incompressible fluctuation
$\hat{\bm{X}}_{\mathrm{arb}}$ with $k_{z}\mathrm{\neq}0$ as
$\hat{\bm{X}}_{\mathrm{arb}}\mathrm{=}\sum_{j=1}^{3}\beta_{j}\hat{X}_{j}$,
where $\beta_{j}$ is the amplitude of the $j^{\mathrm{th}}$ eigenvector
$\hat{X}_{j}$; we reserve $j\mathrm{=}1$ for the GSF-unstable modes, and
$j\mathrm{=}2,3$ for the IGWs that are always linearly stable in this study.
More compactly, $\hat{\bm{X}}_{\mathrm{arb}}\mathrm{=}\bm{E}\bm{\beta}$, where
$\bm{\beta}$ is a (column-)vector of mode amplitudes and $\bm{E}$ is an
eigenvector matrix, whose $j^{\mathrm{th}}$ column is $\hat{X}_{j}$. Thus,
$\bm{\beta}\mathrm{=}\bm{E}^{-1}\hat{\bm{X}}_{\mathrm{arb}}$.

For the $k_{z}\mathrm{=}0$ modes, $\bm{E}$ turns out to be an identity matrix,
meaning that the three components of velocity, and the temperature,
individually form eigenvectors. In what follows, we therefore decompose an
arbitrary fluctuation with $k_{z}\mathrm{=}0$ into
$\hat{\bm{X}}_{\mathrm{arb}}^{\mathrm{T}}=\mathcal{X}[1,0,0,0]+\mathcal{Y}[0,1,0,0]+\mathcal{Z}[0,0,1,0]+\Theta[0,0,0,1]$,
where the amplitudes of the eigenvectors are denoted by
$\mathcal{X},\mathcal{Y},\mathcal{Z}$, and $\Theta$. We reserve the
$\bm{\beta}$-notation above for the amplitudes of eigenvectors with
$k_{z}\mathrm{\neq}0$.

###  2.2 Initial value problem

We perform an ensemble of direct numerical simulations of Eqs. (1a) and (1b),
by seeding a low-amplitude solenoidal random noise to $\bm{u}$, in a box of
size $(L_{x},L_{z})\mathrm{=}(100,100)$. To obtain numerically converged
results, a spatial resolution of up to $512^{2}$ grid points is used in the
pseudo-spectral solver SNOOPY (Lesur & Longaretti, 2005; Barker et al., 2019).

To determine the contribution of each eigenmode in, for example, the turbulent
momentum transport, we decompose the turbulent stress as:
$\langle\widetilde{u}_{x}\widetilde{u}_{y}\rangle\mathrm{=}\sum_{k_{x},k_{z}}\sum_{m,n}2\mathrm{Re}\left[\beta_{m}\hat{u}_{x,m}\beta_{n}^{\ast}\hat{u}_{y,n}^{\ast}\right]$,
where $\langle\cdot\rangle$ is an $(x,z)$-averaging operation; $m$ and $n$ are
summed from $1$ to $3$, corresponding to three excited eigenmodes at every
wavenumber $\bm{k}$; the amplitude $\beta_{m}$ and the $x$-component of the
velocity $\hat{u}_{x,m}$ correspond to the $m^{\mathrm{th}}$ eigenvector at
$\bm{k}$; and likewise for $\beta_{n}$ and $\hat{u}_{y,n}$; the operation ∗
denotes complex conjugation. Using such a decomposition, we obtain the
contribution of an unstable mode at $\bm{k}$ to the momentum transport rate,
which is
$2|\beta_{1}|^{2}\mathrm{Re}\left[\hat{u}_{x,1}\hat{u}_{y,1}^{\ast}\right]$.
This decomposition is performed for every wavenumber, hence allowing us to
trace evolution of transport contributions due to individual unstable modes
[see Fig. 2(a)].

The summed contributions of all eigenvectors from all wavenumbers reproduce,
to machine precision, the total transport rates found in the simulation before
performing mode decomposition, as we show in Fig. 2(b). The contributions of
the unstable modes are also compared across different wavenumber sums. Almost
identical results are found for heat transport (not shown). The unstable modes
from the linearly fastest-growing wavenumber branch $k_{x}\mathrm{=}0$
transport significantly less momentum than the other wavenumbers with
$k_{x}\mathrm{\neq}0$. This is our first surprising result, and it challenges
predictions of turbulent transport that rely on an unstable mode at the
fastest-growing wavenumber alone (e.g., Radko & Smith, 2012; Brown et al.,
2013; Barker et al., 2019). This finding also instructs us to investigate
nonlinear couplings between eigenmodes to understand the instability-
saturation mechanism.

Figure 2: (a) Comparison of momentum transport ($\mathrm{Reynolds\
Stress}\mathrm{=}2|\beta_{1}|^{2}\mathrm{Re}\left[\hat{u}_{x,1}\hat{u}_{y,1}^{\ast}\right]$)
driven by an unstable mode at the linearly fastest-growing wavenumber
$\mathbf{k}=(0,0.63)$, and by an unstable mode at $\mathbf{k}=(0.31,-0.13)$,
the wavenumber that has the largest contribution to the momentum transport in
the nonlinear phase. (b) Eigenmode decomposition of net Reynolds stress
$\langle\widetilde{u}_{x}\widetilde{u}_{y}\rangle$ in nonlinear simulation of
the GSF instability-driven turbulence, showing that the transport due to mode-
undecomposed fluctuations (red curve) and mode-decomposed all eigenmodes
(black curve) agree to machine precision. Transport is almost entirely
($88\%$) due to the unstable modes (green curve); the sum of fastest-growing
unstable modes at $k_{x}\mathrm{=}0$, however, contributes negligibly ($3\%$)
to the transport (blue curve). Simulation parameters used are
$S\mathrm{=}2.1,N^{2}\mathrm{=}10$ and Pr$\mathrm{=}0.01$. Figure 3: Spectra
of linear and nonlinear-saturation properties of the GSF instability. On the
negative-$k_{x}$ domain, the colored square boxes (yellow-green-purple)
display time-averaged energy extraction rates $\langle Q_{1}\rangle_{t}$ by
unstable modes from the mean gradients in a nonlinear simulation. On the
positive-$k_{x}$ domain, colored square boxes (black-red-yellow) show the
time-averaged viscous dissipation rates, which are pronounced at low $k_{z}$.
Over the entire $(k_{x},k_{z})$-plane, the non-square filled and line contours
show the growth rates $\gamma$ of unstable modes, with white dashed contour
lines on the negative-$k_{x}$ domain and with bluish filled contours on the
positive-$k_{x}$ domain. The fastest-growing mode resides at around
$\mathbf{k}=(0,0.63)$. The simulation parameters are
$S\mathrm{=}2.1,N^{2}\mathrm{=}10$ and Pr$\mathrm{=}0.01$.

##  3 Nonlinear saturation by coupling to latitudinal flow

###  3.1 Mode-amplitude evolution

To analyze nonlinear mode couplings, Eqs. (1a)–(1b) are first spatially
Fourier-transformed:
$\partial_{t}\hat{\bm{X}}=\bm{L}\hat{\bm{X}}+\sum_{\bm{k^{\prime}},\bm{k^{\prime\prime}}}\bm{N}(\hat{\bm{X}}^{\prime},\hat{\bm{X}}^{\prime\prime})$,
where $\hat{\bm{X}},\hat{\bm{X}}^{\prime}$, and $\hat{\bm{X}}^{\prime\prime}$
are state vectors at $\bm{k},\bm{k}^{\prime}$, and $\bm{k}^{\prime\prime}$,
respectively, satisfying $\bm{k}=\bm{k}^{\prime}+\bm{k}^{\prime\prime}$. Then,
following Sec. 2.1 we substitute $\hat{\bm{X}}\mathrm{=}\bm{E}\bm{\beta}$, and
likewise for $\hat{\bm{X}}^{\prime}$ and $\hat{\bm{X}}^{\prime\prime}$. We
multiply the obtained equation with $\bm{E}^{-1}$ and take the
$j^{\mathrm{th}}$ row of the resulting equation. This process yields an
evolution equation for the $j^{\mathrm{th}}$ eigenmode at $\bm{k}$. Such an
evolution equation for mode amplitude $\beta_{j}$ for $k_{z}\mathrm{\neq}0$ is
different from that of the mode amplitude
$F\in\\{\mathcal{X},\mathcal{Y},\mathcal{Z},\Theta\\}$ for $k_{z}\mathrm{=}0$,
although both are coupled and nonlinear:

$\displaystyle\partial_{t}\beta_{j}$
$\displaystyle=\gamma_{j}\beta_{j}+\sum_{\begin{subarray}{c}\bm{k}^{\prime},m,n\end{subarray}}C_{jmn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime}$
$\displaystyle\hskip
17.07182pt+\sum_{\begin{subarray}{c}\bm{k}^{\prime},F,n:k_{z}^{\prime}\mathrm{=}0\\\
F\in\\{\mathcal{X},\mathcal{Y},\mathcal{Z},\Theta\\}\end{subarray}}\left[C_{jFn}^{(\bm{k},\bm{k}^{\prime})}+C_{jnF}^{(\bm{k},\bm{k}^{\prime\prime})}\right]F^{\prime}\beta_{n}^{\prime\prime},$
(3a) $\displaystyle\partial_{t}F$
$\displaystyle=-\gamma_{F}F+\sum_{\bm{k}^{\prime},m,n:k_{z}\mathrm{=}0}C_{Fmn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime},$
(3b)

where $\gamma_{j}$ is the complex-valued growth rate for the $j^{\mathrm{th}}$
eigenmode with $k_{z}\mathrm{\neq}0$; the real-valued damping rate
$\gamma_{F}$ is $\gamma_{\mathcal{X}}$ when $F$ is replaced with $\mathcal{X}$
in Eq. (3b); likewise for the replacement of $F$ with
$\mathcal{Y},\mathcal{Z},$ and $\Theta$; we note that
$\gamma_{\mathcal{X}}=\gamma_{\mathcal{Y}}=\gamma_{\mathcal{Z}}=\nu
k_{x}^{2}$, and $\gamma_{\Theta}=\kappa k_{x}^{2}$. The nonlinear-coupling
coefficient, for example, $C_{jmn}^{(\bm{k},\bm{k}^{\prime})}$ measures the
overlap of eigenmodes $m$ with $\bm{k}^{\prime}$, $n$ with
$\bm{k}^{\prime\prime}$, and $j$ with $\bm{k}$. Such a mode coupling
coefficient is found by applying $\bm{E}^{-1}$ to the (column) vector of
$\bm{N}(\hat{\bm{X}}_{m}^{\prime},\hat{\bm{X}}_{n}^{\prime\prime})$, a process
that incorporates all the nonlinearities of the system, thus making
$C_{jmn}^{(\bm{k},\bm{k}^{\prime})}$ ideal for a comprehensive instability-
saturation analysis.

Figure 4: Time evolution of the total nonlinear energy transfer
$T_{1\mathrm{AA}}$ to an unstable eigenmode at a wavenumber where the spectrum
of $\langle\tilde{u}_{x}\tilde{u}_{y}\rangle$ peaks. $T_{1\mathcal{Z}1}$ is
the energy transfer to the same unstable mode via interactions between the
$z$-component of velocity ($\mathcal{Z}$) with wavenumbers $k_{z}\mathrm{=}0$,
and the other unstable modes. Comparison of two transfer functions reveals
that the dominant triad involves a latitudinal flow and two unstable modes.
The simulation parameters are $S\mathrm{=}2.1,N^{2}\mathrm{=}10$ and
Pr$\mathrm{=}0.01$.

On the right-hand side of Eq. (3), the second term is the nonlinear coupling
between eigenmodes with $k_{z}^{\prime}\mathrm{\neq}0$ and
$k_{z}^{\prime\prime}\mathrm{\neq}0$ (hence the two $\beta$s), and the third
term, with an $F$ and a $\beta$, is the nonlinear coupling between eigenmodes
with $k_{z}^{\prime}\mathrm{=}0$ and $k_{z}^{\prime\prime}\mathrm{\neq}0$ (see
the last paragraph of Sec. 2.1). Equations (3) and (3b) have the same number
of degrees of freedom as the original nonlinear equations in physical space,
Eqs. (1a)–(1b). These systems are completely equivalent, but one represents
dynamics in physical space and another in eigenmode space.

###  3.2 Mode-energy evolution

The energy evolution equation for each eigenmode can now be derived by
multiplying Eq. (3) by $\beta_{j}^{\ast}$ and adding the complex conjugate of
the resulting equation to arrive at

$\partial_{t}|\beta_{j}|^{2}=Q_{j}+T_{j\mathrm{A}\mathrm{A}},$ (4)

where $Q_{j}\mathrm{=}2\mathrm{Re}\,\gamma_{j}|\beta_{j}|^{2}$ is the linear
energy transfer rate to $\bm{k}$ from the mean gradients, and
$T_{j\mathrm{A}\mathrm{A}}$ is the total nonlinear energy transfer to the
$j^{\mathrm{th}}$ eigenmode from all possible nonlinear interactions; Eq.
(A10). We now show the spectrum of time-averaged $Q_{1}$, along with that of
the growth rate and time-averaged viscous dissipation rate $\epsilon_{\nu}$ in
Fig. 3.

From $T_{j\mathrm{A}\mathrm{A}}$ in Eq. (4), we separate out the nonlinear
transfer $T_{1\mathrm{\mathcal{Z}}1}$ in a triad that involves a latitudinal
flow $\mathcal{Z}$ at $k_{z}^{\prime}\mathrm{=}0$ and two GSF-unstable modes
($j\mathrm{=}1$) at $k_{z}^{\prime}\mathrm{\neq}0$:

$T_{1\mathrm{\mathcal{Z}}1}=\sum_{\bm{k}^{\prime}:k_{z}^{\prime}=0}2\mathrm{Re}\,\Bigg{\\{}\left[C_{1\mathcal{Z}1}^{(\bm{k},\bm{k}^{\prime})}+C_{11\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime\prime})}\right]\mathcal{Z}^{\prime}\beta_{1}^{\prime\prime}\beta_{1}^{\ast}\Bigg{\\}}.$
(5)

We now compare $T_{1\mathcal{Z}1}$ with $T_{1\mathrm{A}\mathrm{A}}$ in Fig. 4,
for a GSF-unstable mode with a wavenumber that contributes the largest to the
momentum transport. Repeating this transfer analysis at different wavenumbers
produces similar results. The two transfers are nearly identical, which
confirms the conjecture (Barker et al., 2019) that, in the fully nonlinear
phase, the GSF instability saturates via the formation of strong latitudinal
jets or flows. Such flows are $z$-directed, although with no $z$-variation,
and primarily have a wavenumber $k_{x}\mathrm{=}2\pi/L_{x}$; see Fig. 1,
right-column snapshot at $t\mathrm{=}6000$. These flows generalize to
meridional circulation in stars, and resemble zonal jets in planetary
atmospheres and fusion systems (Terry, 2019). Flows with $k_{z}\mathrm{=}0$
are, however, linearly stable to the GSF instability, and, thus, must
necessarily be excited nonlinearly by the interactions between the GSF-
unstable modes. This energy received is then viscously damped at
$k_{z}\mathrm{=}0$ and at low $k_{z}$, as seen in Fig. 3. To sum up, the mean
shear flow, destabilized by the thermal diffusion, lends energy to the
fluctuations via the GSF-unstable modes, which saturate by exciting
$k_{z}\mathrm{=}0$ latitudinal flows to a significant level. Such flows then
viscously dissipate the turbulent energy. This is the saturation mechanism of
the axisymmetric GSF instability found here.

##  4 Angular momentum transport model

The findings shown so far are sufficient to build a statistical closure model,
with no free parameters, and thus with predictive power.

Equation (3) has a quadratic nonlinearity, hence evolutionary equations for
mode energy contain triplet interactions [e.g., Eq. (5)]. To determine the
evolving triplet interaction terms, one can derive an equation with quadruplet
interactions [Eq. (A13)], and so on. To truncate this never-ending hierarchy
(the so-called ‘‘turbulence closure problem”), we invoke a standard turbulence
closure, the Eddy-Damped Quasi-Normal Markovian (EDQNM) approximation (see,
e.g., Orszag 1970; Terry et al. 2018; Hegna et al. 2018; Terry et al. 2021;
Pueschel et al. 2021; Li et al. 2021, 2023), that truncates the hierarchy at
fourth-order cumulants of the fluctuations, thereby assuming that the
statistics for the mode amplitudes are close to Gaussian. The resulting
equation, however, is still nonlinear and daunting. But when a latitudinal
flow $\mathcal{Z}$, with $k_{z}\mathrm{=}0$, dominates the nonlinear coupling,
the complexity of the equation is significantly reduced (Terry et al., 2018).

###  4.1 An outline of the Closure Model

We illustrate here the key steps involved to explain most simply our closure
model (by omitting details and treating all variables as real). First, we
observe that Eq. (3) has the structure:

$\partial_{t}\beta=...\beta+...\beta\beta+...\beta\mathcal{X}+...\beta\mathcal{Y}+...\beta\mathcal{Z}+...\beta\Theta\,,$
(6)

where the mode amplitudes are explicitly shown, and the dots ($...$) represent
terms such as the linear growth rate and the nonlinear coupling coefficients.
One can then obtain evolution equation for second-order correlator as
$\partial_{t}(\beta\beta)=...\beta\beta+...\beta\mathcal{Z}\beta$, where the
other nonlinear terms, e.g., $\beta\beta\beta$ and $\beta\beta\mathcal{X}$,
have been dropped because the nonlinear energy transfer is almost entirely
dominated by $\beta\mathcal{Z}\beta$—the latitudinal-flow coupling (Fig. 4).
Since $\beta\mathcal{Z}\beta$ also evolves, one can similarly obtain evolution
equation for third-order correlator as
$\partial_{t}(\beta\mathcal{Z}\beta)=...\beta\mathcal{Z}\beta+...\beta\beta\mathcal{Z}\mathcal{Z}$.
The closure solution then yields a relation
$\beta\mathcal{Z}\beta=...\beta\beta\mathcal{Z}\mathcal{Z}$. Although useful
later, this relation does not predict the mode amplitude $\beta$, needed for
the turbulent transport prediction.

To predict the mode amplitude $\beta$, we consider the latitudinal-flow
evolution equation, $\partial_{t}\mathcal{Z}=...\mathcal{Z}+...\beta\beta$,
and derive
$\partial_{t}(\mathcal{Z}\mathcal{Z})=...\mathcal{Z}\mathcal{Z}+...\beta\mathcal{Z}\beta$.
Then, $\beta\mathcal{Z}\beta$ can be replaced with a product of four
amplitudes using the closure solution in the previous paragraph. One thus
obtains
$\partial_{t}(\mathcal{Z}\mathcal{Z})=...\mathcal{Z}\mathcal{Z}+...\beta\beta\mathcal{Z}\mathcal{Z}$.
In quasi-stationary turbulence, $\partial_{t}\sim 0$, and thus
$\beta\beta\mathcal{Z}\mathcal{Z}=...\mathcal{Z}\mathcal{Z}$. The EDQNM
closure allows writing a fourth-order correlator
$\beta\beta\mathcal{Z}\mathcal{Z}$ as a sum of products of second-order
correlators, such as $|\beta|^{2}|\mathcal{Z}|^{2}$. Then, cancelling
$|\mathcal{Z}|^{2}$ from both sides of
$\beta\beta\mathcal{Z}\mathcal{Z}=...\mathcal{Z}\mathcal{Z}$, one predicts the
saturated mode amplitude (or energy): $\beta\beta=...\,$. Using this, one can
make predictions for turbulent transport rates as we shall show in the next
subsection.

###  4.2 Detailed Closure Model

To make quantitative predictions for transport, we take the EDQNM-closed
evolution equation for the latitudinal flow energy [see Appendix B, Eq.
(A20)],

$\displaystyle\partial_{t}|\mathcal{Z}|^{2}/2=$
$\displaystyle-\gamma_{\mathcal{Z}}|\mathcal{Z}|^{2}$ (7)
$\displaystyle+|\mathcal{Z}|^{2}\sum_{\bm{k}^{\prime}}|\beta_{1}^{\prime}|^{2}\,\mathrm{Re}\,\Bigg{.}\left[-\tau_{11\mathcal{Z}}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{\mathcal{Z}11}^{(\bm{k},\bm{k}^{\prime})}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{1\mathcal{Z}1}^{(\bm{k}^{\prime\prime},\bm{k})}\right],$

where, on the right-hand side, the second term contains a product of four
amplitudes, but notably with $|\mathcal{Z}|^{2}$ that also appears in the
first term of the right-hand side. In Eq. (7),

$\tau_{11\mathcal{Z}}=(\gamma_{1}^{\prime}+\gamma_{1}^{\prime\prime}-\gamma_{\mathcal{Z}}^{\ast})^{-1},$
(8)

is the three-wave interaction time found from the EDQNM closure, and
$\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{lmn}^{(\bm{p},\bm{q})}=C_{lmn}^{(\bm{p},\bm{q})}+C_{lnm}^{(\bm{p},\bm{p}-\bm{q})}$
is the symmetrized coupling coefficient. In quasi-stationary turbulence,
$\partial_{t}\sim 0$, and thus the linear and nonlinear terms must balance.
First, for simplicity, we consider a latitudinal flow at $(k_{x},0)$ that is
driven by two unstable modes at $(k_{x}^{\prime},k_{z}^{\prime})$ and
$(k_{x}\mathrm{-}k_{x}^{\prime},\mathrm{-}k_{z}^{\prime})$; then, using (7),
$|\beta_{1}^{\prime}|^{2}=\gamma_{\mathcal{Z}}\left(\mathrm{Re}\,[-\tau_{11\mathcal{Z}}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{\mathcal{Z}11}^{(\bm{k},\bm{k}^{\prime})}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{1\mathcal{Z}1}^{(\bm{k}^{\prime\prime},\bm{k})}]\right)^{-1}$.
A more general expression for $|\beta_{1}^{\prime}|^{2}$ is found by using a
standard Markovian assumption (Terry et al., 2021): $|\beta_{1}^{\prime}|^{2}$
is more weakly dependent on wavenumbers than the other factors in Eq. (7)
arising from the coupling coefficients and $\tau_{11\mathcal{Z}}$. Such a
consideration provides an expression for nonlinearly saturated squared-mode-
amplitude

$|\beta^{\prime}|^{2}=\gamma_{\mathcal{Z}}\times{\left(\frac{\gamma}{k^{2}}\right)}_{\mathrm{Closure}},$
(9)

with

${\left(\frac{\gamma}{k^{2}}\right)}_{\mathrm{Closure}}=\frac{1}{|\sum_{\bm{k^{\prime}}}\tau_{11\mathcal{Z}}\,\mathrm{Re}\,[-\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{\mathcal{Z}11}^{(\bm{k},\bm{k}^{\prime})}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{1\mathcal{Z}1}^{(\bm{k}^{\prime\prime},\bm{k})}]|},$ (10)

where we note that the coupling coefficients scale linearly with wavenumbers,
and $\tau_{11\mathcal{Z}}$ is the inverse of the sum of three growth rates of
eigenmodes in a triad.

The growth rates in $\tau_{11\mathcal{Z}}$ should, in principle, also have
amplitude-dependent eddy-damping rates as they become non-negligible, for
example, in homogeneous isotropic fluid turbulence; however, when waves or
instabilities exist, and when the turbulent transport spectrum is dominated by
low wavenumbers, as in this study, $\tau_{11\mathcal{Z}}$ is approximated by
using the linear growth rates (Terry et al., 2018, 2021). Using such, one
identifies that the triplet interaction time $\tau_{11\mathcal{Z}}$ is maximal
when the triad involves a latitudinal flow ($\mathcal{Z}$) and two GSF-
unstable modes ($j\mathrm{=}1$). Shorter triplet interaction times
$\tau_{12\mathcal{Z}}$ are expected for triads with, for example, the
latitudinal flow, an unstable mode, and a strongly damped IGW
($j\mathrm{=}2$), as such an interaction lowers $\tau_{12\mathcal{Z}}$ via
both the frequency and damping rate of the IGW. The largest interaction time
$\tau_{11Z}$ dominates saturation.

The radial turbulent transport of angular momentum is measured by
$\langle\tilde{u}_{x}\tilde{u}_{y}\rangle\mathrm{\approx}\sum_{\bm{k}^{\prime\prime\prime}}\hat{u}_{x,1}^{\prime\prime\prime}\hat{u}_{y,1}^{\prime\prime\prime\ast}|\beta_{1}^{\prime\prime\prime}|^{2}$,
where $\beta_{1}^{\prime\prime\prime}$ is the unstable-mode amplitude at
$\bm{k}^{\prime\prime\prime}$ over which the summation is applied. Then, using
$|\beta_{1}^{\prime\prime\prime}|^{2}$ from the above paragraph,

$\langle\tilde{u}_{x}\tilde{u}_{y}\rangle_{\mathrm{Closure}}={\left(\frac{\gamma}{k^{2}}\right)}_{\mathrm{Closure}}\times\gamma_{\mathcal{Z}}\sum_{\bm{k}^{\prime\prime\prime}}\hat{u}_{x,1}^{\prime\prime\prime}\hat{u}_{y,1}^{\prime\prime\prime\ast}.$
(11)

The $y$-component $\hat{u}_{y,1}^{\prime\prime\prime}$ of velocity of the
unstable eigenvector can be replaced with, e.g., its temperature perturbation
$\hat{\theta}_{1}^{\prime\prime\prime}$ to predict the turbulent heat flux
$\langle\tilde{u}_{x}\tilde{\theta}\rangle$.

In our simulations with radial differential rotation, the latitudinal momentum
flux is much lower than the radial flux, and, when time-averaged, it is nearly
null.

###  4.3 Tests of theoretical predictions

Figure 5: Tests of predictions of our closure model (red diamond) and a
quasilinear-type, parasitic-saturation model (blue inverted triangle) against
direct numerical simulations (DNS, black/gray circle). Variations of momentum
transport rates are shown in (a); the filled markers correspond to the cases
where the shear parameter $S$ is varied ($N^{2}\mathrm{=}10$); the unfilled
markers correspond to the cases where the squared Brunt-Väisälä frequency
$N^{2}$ is varied ($S\mathrm{=}2.1$). Both $S$\- and $N^{2}$-scan results
collapse onto a single master curve, when
$\langle\tilde{u}_{x}\tilde{u}_{y}\rangle$ is scaled by a factor shown on the
$y$-axis that transforms the governing equations of the GSF instability
studied here to depend on only two dimensionless parameters $(r,\mathrm{Pr})$.
The GSF instability operates when $r\in[0,1)$ and $\mathrm{Pr}\mathrm{<}1$
($\mathrm{Pr}\mathrm{=}0.01$ is chosen). The shown $y$-axis is precisely an
expression for the chemical transport rate for the thermohaline instability
[see Eq. (15a)]. Heat transport rates, shown in (b), display nearly identical
trends; [see Eq. (15b) for the scaling factor]. The closure prediction agrees
with full DNS better than the quasilinear prediction over the scanned range of
parameters.

A simple quasilinear model of the GSF-instability saturation was recently
proposed (Barker et al., 2019, 2020) by assuming that a secondary ‘‘parasitic”
instability feeds on the primary GSF-unstable mode. Such an assumption, also
called ‘‘parasitic saturation mechanism” (Goodman & Xu, 1994; Radko & Smith,
2012; Brown et al., 2013; Harrington & Garaud, 2019; Fraser et al., 2023), is
based on a single primary mode at the fastest-growing wavenumber
$\mathbf{k}^{\prime\prime\prime}$, which predicts the transport rate

$\langle\widetilde{u}_{x}\widetilde{u}_{y}\rangle_{\mathrm{QL}}=\frac{\gamma_{1}^{\prime\prime\prime
2}}{k_{z}^{\prime\prime\prime
2}}\hat{u}_{x,1}^{\prime\prime\prime}\hat{u}_{y,1}^{\prime\prime\prime\ast}f^{2}(\mathbf{k}^{\prime\prime\prime}),$
(12)

whose form is made manifestly similar to Eq. (11); the factor
$f(\mathbf{k}^{\prime\prime\prime})$, which is evaluated at
$\mathbf{k}^{\prime\prime\prime}\mathrm{=}(0,k_{z}^{\prime\prime\prime})$, is
the normalization factor of eigenmodes. Here,
$\gamma_{\nu}^{\prime\prime\prime}=\gamma_{1}^{\prime\prime\prime}+\nu
k^{\prime\prime\prime 2}$. To find
$\langle\widetilde{u}_{x}\widetilde{\theta}\rangle_{\mathrm{QL}}$, one can
replace $\hat{u}_{y,1}^{\prime\prime\prime}$ on the right-hand side of Eq.
(12) with $\theta_{1}^{\prime\prime\prime}$.

Predictions of Eqs. (11) and (12) are compared against transport rates from
direct numerical simulations in Figs. 5. Significant improvement in both the
momentum and heat transport predictions is observed with the statistical
closure model. The orders-of-magnitude variation in transport rates is
captured by the closure model.

In Fig. 5, for smaller values of $r$, though the transport rates of the
closure and the quasilinear models are similar, we emphasize that this
similarity is merely accidental: the physics the two models incorporate is
very different. The assumptions of the closure model are supported by detailed
numerical evidence (Figs. 2 and 3), including that of the dominant three-wave
coupling between two unstable modes and a latitudinal jet (Fig. 4). No jet
physics is considered in the quasilinear model. The quasilinear model predicts
transport rates based on only one fastest growing wavenumber $k_{z}$, with
$k_{x}\mathrm{=}0$, which Fig. 2 shows is inadequate. Hence, predictions of
the quasilinear model that tend to reproduce the data-validated closure model
predictions are at best a fortuitous coincidence, occurring in a very limited
parameter regime.

###  4.4 Impact of the new model in astrophysics

Since stellar interiors typically have extreme parameters such as
$\mathrm{Pr}\lesssim 10^{-6}$, current and anticipated near-future
computational resources are insufficient to permit direct numerical
simulations of realistic turbulence in them. In the face of such a challenge,
progress can be made by developing analytical theories, informed and tested by
numerical simulations at more accessible parameters. Thus, we now employ our
analytical theory to extrapolate and make predictions for realistic
astrophysical parameters. To achieve this, we derive fully analytic
expressions for all elements of the closure model, assuming that the coupling
of two GSF-unstable modes with the latitudinal jet remains dominant; see
Appendices A and B.

We then compare predictions of the closure model with those of the quasilinear
(QL) model, over a wide range of parameters $\mathrm{Pr}\approx
10^{-7}\textrm{--}1$ and $r\approx 10^{-5}\textrm{--}1$ in Fig. 6. Noting that
$N^{2}\mathrm{=}2(S-2)\left[1+r(\mathrm{Pr}^{-1}-1)\right]$, these scans span
$N^{2}\lesssim 19\times 10^{6}$ (in terms of $\Omega^{2}$); this ratio is
typically around 1 million for the Sun (Christensen-Dalsgaard et al., 1996).
In Fig. 6, with $S\mathrm{=}3$ (in terms of $\Omega$), the Richardson number
is as large as $\approx 2\times 10^{6}$. More extreme parameters can be easily
and quickly scanned with the analytic formula we have derived.

Now we predict transport efficiency of the GSF instability in stars. The
Reynolds stress is of order $\langle\tilde{u}_{x}\tilde{u}_{y}\rangle H^{-1}$,
where $H\mathrm{\equiv}U_{0}{(\partial_{x}U_{0})}^{-1}$ is the scale height of
the mean flow $U_{0}\mathrm{=}-\mathcal{S}x$. The time scale for modifying the
flow is $\tau_{\mathrm{turb}}\sim
U_{0}H/\langle\tilde{u}_{x}\tilde{u}_{y}\rangle\sim
Sx^{2}\Omega^{-1}d^{-2}/\langle\tilde{u}_{x}\tilde{u}_{y}\rangle_{\mathrm{dimensionless}}$
(see also Barker et al., 2019).

Though $\mathrm{Pr}\sim 10^{-6}$, the typical values of $r$ in the solar
tachocline and red giant stars are $r\sim 10^{-3}\textrm{--}1$ (varying with
radius). Then, using Fig. 6 where
$\langle\tilde{u}_{x}\tilde{u}_{y}\rangle_{\mathrm{dimensionless}}$ is on
average $0.5$, we predict $\tau_{\mathrm{turb}}\Omega\sim 2S(x/d)^{2}$. This
turbulent transport time scale is sufficiently short to be astrophysically
important, depending on the relative length scale $x/d$ of the mean flow and
shear strength $S$. For example, it can be as short as $\mathcal{O}(10)$ Myr
using values of $S$ and $x/d$ for the solar tachocline.

The turbulent transport rate depends sensitively also on the shear parameter
$S$ (and latitude and orientation of the shear, i.e., radial or mixed radial-
horizontal shear), and orders of magnitude faster turbulent transport is
possible. We highlight that, because the turbulent time scale for the GSF
instability can be shorter than $\mathcal{O}(10)$ Myr, incorporation of our
transport model for the GSF turbulence in stellar evolution codes is warranted
(particularly if extended to non-equatorial and 3-D GSF instabilities). Using
such, the long-term impact on the evolution of the rotation profile may be
assessed, informing us of the effects of the GSF instability in rapidly
rotating young stars. In this regard, the transport model built here for the
2.5D equatorial GSF instability (and the thermohaline instability) is
significant, as a reliable and reduced numerical treatment of the GSF
instability-driven turbulence is now available.

Figure 6: Predictions of (a) closure model and (b) quasilinear (QL) model for
turbulent momentum transport. The QL model, largely independent of Pr, fails
to reproduce the behavior of the closure model. For Pr closer to $1$, the
closure model predicts that the transport rate increases with decreasing
$r$—until an asymptotically large transport is attained. For $\mathrm{Pr}\ll
r$, the transport first increases and then decreases with $r$, in contrast to
the QL prediction. The white dashed-dotted line, with a unit slope, separates
distinct regimes of $\mathrm{Pr}<r$ and $\mathrm{Pr}>r$ found in (a).

###  4.5 Relation to the thermohaline instability

We emphasize that the equations describing 2.5-D thermohaline and GSF
instabilites at the equator are identical, when the compositional diffusivity
is equal to the kinematic viscosity—a case often realized in stars. We first
write Eqs. (7)–(9) of Brown et al. (2013), where they measure distance in
units of the characteristic length scales $d$ of the fingers, time in units of
the characteristic diffusion time scale $\tau\mathrm{=}d^{2}/\kappa$, and
scaled temperature $T$ in units of $N^{2}d$ [we label their $x$-coordinate
with our $z$-coordinate and vice-versa]

$\displaystyle\mathrm{Pr}^{-1}Du_{x}^{\mathrm{B}}$
$\displaystyle=-\partial_{x}p^{\mathrm{B}}+(T^{\mathrm{B}}-\mu^{\mathrm{B}})+\nabla^{2}u_{x}^{\mathrm{B}},$
(13a) $\displaystyle\mathrm{Pr}^{-1}Du_{z}^{\mathrm{B}}$
$\displaystyle=-\partial_{z}p^{\mathrm{B}}+\nabla^{2}u_{z}^{\mathrm{B}},$
(13b) $\displaystyle DT^{\mathrm{B}}$
$\displaystyle=-u_{x}^{\mathrm{B}}+\nabla^{2}T^{\mathrm{B}},$ (13c)
$\displaystyle D\mu^{\mathrm{B}}$
$\displaystyle=-\frac{u_{x}^{\mathrm{B}}}{R_{0}}+\Pr\nabla^{2}\mu^{\mathrm{B}},$
(13d)

where the state vector
$[u_{x}^{\mathrm{B}},u_{z}^{\mathrm{B}},T^{\mathrm{B}},\mu^{\mathrm{B}},p^{\mathrm{B}}]$
represents Brown-normalized $x$\- and $z$-velocities, scaled temperature,
chemical concentration, and fluid pressure, respectively. The so-called
density ratio $R_{0}\mathrm{=}-N^{2}/\kappa_{\mathrm{ep}}^{2}$ may be recast
as $R_{0}\mathrm{=}1+r(\mathrm{Pr}^{-1}-1)$. The solutions of Eqs. (13a)–(13d)
critically depend only on two parameters: Pr and $r$. Equations (13a)–(13d)
for the thermohaline instability with the state vector
$[u_{x}^{\mathrm{B}},u_{z}^{\mathrm{B}},T^{\mathrm{B}},\mu^{\mathrm{B}}]$ are
identical to Eqs. (1a) and (1b) for the GSF instability with the state vector
$[u_{x},u_{z},\theta,u_{y}]$, when we note

$\displaystyle u_{x}^{\mathrm{B}}$ $\displaystyle=\frac{u_{x}}{d/\tau},$ (14a)
$\displaystyle u_{z}^{\mathrm{B}}$ $\displaystyle=\frac{u_{z}}{d/\tau},$ (14b)
$\displaystyle T^{\mathrm{B}}$ $\displaystyle=\frac{\theta}{N^{2}d},$ (14c)
$\displaystyle\mu^{\mathrm{B}}$
$\displaystyle=-\frac{2\tau\Omega}{\mathrm{Pr}}\frac{u_{y}}{d/\tau}.$ (14d)

Hence, the transport rates with two kinds of non-dimensionalizations—one using
$\tau$ and $d$ as the characteristic time scale and length scale, and another
using $\Omega$ and $d$ as the relevant scales—are related in the manner

$\displaystyle-\langle\widetilde{u}_{x}^{\mathrm{B}}\widetilde{\mu}^{\mathrm{B}}\rangle=\langle\widetilde{u}_{x}\widetilde{u}_{y}\rangle\times
2\Omega\mathrm{Pr}^{1/2}(N^{2})^{{-3/2}}d^{-2},$ (15a)
$\displaystyle\langle\widetilde{u}_{x}^{\mathrm{B}}\widetilde{T}^{\mathrm{B}}\rangle=\langle\widetilde{u}_{x}\widetilde{\theta}\rangle\times\mathrm{Pr}^{1/2}(N^{2})^{{-3/2}}d^{-2},$
(15b)

where the variables with the superscripted ‘B’ are functions of two essential
parameters — $r$ and Pr only, whereas the variables without the superscript
are functions defined by the parameters $\mathrm{Pr},S,$ and $N^{2}$.

The expressions given in Eqs. (15a) and (15b) are plotted in Fig. 5. Thus our
Fig. 5 also represents a comparison of chemical and heat transport between
direct numerical simulations and analytical models for the thermohaline
instability-driven turbulence.

The reduction of the governing equations of the GSF instability to two
parameters ($r,\mathrm{Pr}$) is realized only in 2.5-D equatorial case, which
is where the analogy of the GSF instability with the thermohaline instability
becomes exact.

##  5 Discussion and conclusions

Turbulent transport in stellar interiors is a phenomenon too complex to
represent directly in stellar evolution models. It is often parametrized using
low-order models, such as mixing-length theories, or models that predict
transport rates based on the fastest-growing unstable mode (e.g., Denissenkov,
2010; Brown et al., 2013; Barker et al., 2019). The reliability of such models
can be compromised by several factors: first, the instability dispersion
relation is often anisotropic, a property that affects the nonlinear energy
transfer, thereby driving low-frequency fluctuations (Fig. 3). Second,
nonlinear mode coupling can strongly excite more weakly growing unstable modes
over a wide range of wavenumbers, presenting difficulties to single-mode
theory-based predictions. To circumvent such challenges, we build a nonlinear
mode-coupling theory, informed by detailed analyses of direct numerical
simulations, to arrive at a reliable and analytic transport model, that is
free of tunable parameters. We achieve this here for axisymmetric low-Pr
turbulence driven by centrifugally unstable differential rotation (the GSF
instability) at the equator in a stellar radiative zone.

Although 2.5-D turbulence driven by the GSF instability can differ from fully
3-D cases (Barker et al., 2019), strong secondary flows or jets have been
found in 3-D global systems, also. For example, recent simulations of fully
3-D spherical-shell non-rotating fingering convection have exhibited strong,
large-scale jets (Tassin et al., 2023). Such coherent jets are ubiquitious in
various settings such as in geo- and astrophysical observations, numerical
simulations of 3-D shear-flow instability-driven turbulence (Tripathi et al.,
2023b), and in laboratory fusion plasmas (Terry, 2019). The examples show the
large-scale flows can emerge even in global geometry. Hence, the success of
our theory offers a future possibility to extend the statistical closure
framework, presented here, to the more realistic 3-D simulations of the GSF
instability, ideally in a sphere, at a general latitude and for a range of Pr
(Barker et al., 2020; Dymott et al., 2023; Garaud & Brummell, 2015). It is
also possible that our closure-model framework can be adapted to magnetized
turbulence driven by unstable differential rotation.

Since our formulae are fully analytic, they are quick to implement in stellar
evolution codes such as MESA (Paxton et al., 2011, 2019) to reliably predict
axisymmetric GSF instability-driven turbulent transport rates that vary with
Pr and $r$ at different spatial grid points in an evolving star. It is
straightforward to compute the values of $r$ and Pr at a given spatial grid
point in a modeled star, and simply look up transport rates using our Fig.
6(a) to prescribe the rates of transport of angular momentum and heat; to
access the look up table of transport, see Data Availability. Since we have
also provided formulae and frameworks for the turbulent transport of the
axisymmetric GSF-analogous thermohaline instability, there now exists a
reliable chemical transport model, which employs DNS-confirmed key elements of
nonlinear saturation of the instability in stars.

## Acknowledgments

We thank the anonymous reviewer for constructive feedback. This work was
enabled by the support from the US Department of Energy Grant No. DE-
SC0022257, from the STFC Grant Nos. ST/S000275/1 and ST/W000873/1, and from
the NASA HTMS Grant No. 80NSSC20K1280. A.E.F. acknowledges support from the
George Ellery Hale Postdoctoral Fellowship in Solar, Stellar and Space Physics
at the University of Colorado. Computations were performed using ARC4, part of
the High Performance Computing facilities at the University of Leeds, and the
DiRAC Data Intensive service at Leicester, operated by the University of
Leicester IT Services, which forms part of the STFC DiRAC HPC Facility
(www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC
capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant
ST/R001014/1. DiRAC is part of the National e-Infrastructure. We also used
ACCESS computing resources via Allocation No. TG-PHY130027. The inception of
this work took place at the stimulating KITP program ‘‘Probes of Transport in
Stars” (supported in part by the NSF under Grant No. NSF PHY-1748958).

## Data Availability

The GitHub repository,
https://github.com/BindeshTripathi/GSF_transport, hosts the data presented in
Fig. 6(a), and a python script where the closure model is implemented. Other
data used in this article will be shared on reasonable request to the
corresponding authors.

## Appendix A: Nonlinear mode-coupling coefficients of the GSF instability

The $j^{\mathrm{th}}$ eigenvector of the linear operator $\mathbf{L}$ of the
GSF instability satisfies

$\hat{u}_{z,j}=-\frac{k_{x}}{k_{z}}\hat{u}_{x,j};\hskip
28.45274pt\hat{\theta}_{j}=-\frac{N^{2}}{\gamma_{\kappa,j}}\hat{u}_{x,j};\hskip
28.45274pt\hat{u}_{y,j}=-\frac{2\Omega-S}{\gamma_{\nu,j}}\hat{u}_{x,j},$ (A1)

using which the matrix $\mathbf{E}$ of eigenvectors can be created and
inverted as mentioned in the penultimate paragraph of Sec. 2.1. To expedite
analytic calculations, one may solve for the adjoint solutions
$\mathbf{Y}_{j}$ of $\mathbf{L}$, which are the eigenvectors of
$\mathbf{L}^{\dagger}$. Such adjoint solutions $\mathbf{Y}_{j}$ form a
biorthogonal basis with the eigenvectors $\mathbf{X}_{j}$ of $\mathbf{L}$
(Fraser et al., 2021; Tripathi et al., 2022a, 2023a, 2023b). That is,
$\langle\mathbf{Y}_{m},\mathbf{X}_{j}\rangle={\mathbf{Y}_{m}^{\mathrm{T}\ast}}\mathbf{X}_{j}\propto\delta_{j,m}$,
where $\mathrm{T}\ast$ is the transpose-conjugation operation; see Appendix A
of Tripathi et al. (2022b) for a general mathematical proof. It is known that
both $\mathbf{L}$ and $\mathbf{L}^{\dagger}$ have the same eigenvalues. The
expedient adjoint technique is equivalent to inverting the matrix $\mathbf{E}$
of size $4\times 4$, and they deliver completely identical results, which we
have verified. The $j^{\mathrm{th}}$ adjoint solution satisfies

$\hat{u}_{z,j}=0;\hskip
28.45274pt\hat{\theta}_{j}=\frac{k_{z}^{2}}{k^{2}\gamma_{\kappa,j}}\hat{u}_{x,j};\hskip
28.45274pt\hat{u}_{y,j}=\frac{2\Omega
k_{z}^{2}}{k^{2}\gamma_{\nu,j}}\hat{u}_{x,j}.$ (A2)

Using such, we write the nonlinear mode-coupling coefficient between the two
eigenmodes $m$ at wavenumber $\mathbf{k}^{\prime}$ and $n$ at
$\mathbf{k}^{\prime\prime}$, impacting the eigenmode $j$ at $\mathbf{k}$, as

$C_{jmn}^{(\mathbf{k},\mathbf{k}^{\prime})}=\frac{\langle\mathbf{Y}_{j}(\mathbf{k}),\mathbf{N}(\mathbf{X}_{m}^{\prime},\mathbf{X}_{n}^{\prime\prime})\rangle}{\langle\mathbf{Y}_{j}(\mathbf{k}),\mathbf{X}_{j}(\mathbf{k})},$
(A3)

where $\mathbf{N}(\mathbf{X}_{m}^{\prime},\mathbf{X}_{n}^{\prime\prime})$ is
the nonlinearity vector; for example, its $\theta$-component is
$-\mathbf{u}_{m}^{\prime}\cdot\nabla\theta_{n}^{\prime\prime}$.

Following this procedure, we have distilled the analytic coupling coefficients
for the GSF instability, and provide below their final expressions:

$\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{\mathcal{Z}11}^{(\mathbf{k},\mathbf{k}^{\prime})}=\frac{i\left(k_{x}^{\prime
2}-k_{x}^{\prime\prime 2}\right)}{k_{z}^{\prime}f^{\prime}f^{\prime\prime}},$
(A4)

where $f$ is the eigenmode normalization factor. With inverse dimensions of
$\hat{u}_{x,j}$, a suitable mode normalization can be
$f^{\prime}\mathrm{=}k^{\prime}/\gamma^{\prime}$ or
$f^{\prime}\mathrm{=}k^{\prime 2}d/(\gamma^{\prime}\gamma_{\nu}^{\prime}\tau)$
or their variants, with $\tau\mathrm{=}d^{2}/\kappa$ as the characteristic
diffusion time scale. Equation (A4) is simple as it represents the coupling
between two unstable modes that drive the latitudinal flow $\mathcal{Z}$. The
second coupling coefficient required in the closure model is

$\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{1\mathcal{Z}1}^{(\mathbf{k}^{\prime\prime},\mathbf{k})}=ik_{z}^{\prime}\frac{\left[\frac{k^{\prime\prime
2}}{N^{2}k_{z}^{\prime\prime
2}}\left(1-\frac{2k_{x}k_{x}^{\prime\prime}}{k^{\prime\prime
2}}\right)+\frac{1}{R_{0}\gamma_{\nu}^{\prime}\gamma_{\nu}^{\prime\prime}}-\frac{1}{\gamma_{\kappa}^{\prime}\gamma_{\kappa}^{\prime\prime}}\right]f^{\prime\prime}}{\left[\frac{k^{\prime\prime
2}}{N^{2}k_{z}^{\prime\prime 2}}+\frac{1}{R_{0}\gamma_{\nu}^{\prime\prime
2}}-\frac{1}{\gamma_{\kappa}^{\prime\prime 2}}\right]f^{\prime}}.$ (A5)

We note that the coupling coefficients depend only on wavenumbers and the
input parameters such as $N^{2}$, Pr,
$R_{0}\mathrm{=}-N^{2}/\kappa_{\mathrm{ep}}^{2}$ and the growth rates; the
growth rates in turn depend only on wavenumbers and the input parameters [Eq.
(2)]. Asymptotic approximation to the growth rate in the limit of, e.g., small
Pr is possible (Brown et al., 2013).

## Appendix B: Nonlinear mode-coupling coefficients of the Thermohaline
instability

Here we provide analytic expressions needed for the closure model applicable
to the thermohaline instability.

We find the $j^{\mathrm{th}}$ eigenvector of the linear operator $\mathbf{L}$
found from Eqs. (13a)–(13d) for the thermohaline instability satisfies

$\hat{u}^{\mathrm{B}}_{z,j}=-\frac{k_{x}}{k_{z}}\hat{u}^{\mathrm{B}}_{x,j};\hskip
28.45274pt\hat{T}^{\mathrm{B}}_{j}=-\frac{1}{\gamma_{\kappa,j}}\hat{u}^{\mathrm{B}}_{x,j};\hskip
28.45274pt\hat{\mu}^{\mathrm{B}}_{j}=-\frac{1}{R_{0}\gamma_{\nu,j}}\hat{u}^{\mathrm{B}}_{x,j}.$
(A6)

The $j^{\mathrm{th}}$ adjoint solution for the thermohaline instability is

$\hat{u}^{\mathrm{B}}_{z,j}=0;\hskip
28.45274pt\hat{T}^{\mathrm{B}}_{j}=\frac{k_{z}^{2}\mathrm{Pr}}{k^{2}\gamma_{\kappa,j}}\hat{u}^{\mathrm{B}}_{x,j};\hskip
28.45274pt\hat{\mu}^{\mathrm{B}}_{j}=-\frac{k_{z}^{2}\mathrm{Pr}}{k^{2}\gamma_{\nu,j}}\hat{u}^{\mathrm{B}}_{x,j}.$
(A7)

The coupling coefficients for the thermohaline instability are identically the
same as those for the GSF instability; Equations (A4) and (A5) require only a
minute modification: $N^{2}$ appearing twice in Eq. (A5) should be replaced
with Pr. We, conclude that our closure model and predictions are directly
applicable to the thermohaline instability, as well. This is significant
because the turbulent transport efficiencies of these two instabilites in
stars are not known, but are generally thought to be important. With the
closure model at hand, we can now make predictions reliably for both
instabilities.

## Appendix C: Details of Closure Model Calculations

To make the statistical closure model more accessible to a wide range of
readers, we provide below detailed, step-by-step derivations.

### C1. Amplitude evolution equation

The mode-amplitude $\beta_{j}$ evolution equation, given in Eq. (3), for a
wavenumber $\mathbf{k}=(k_{x},k_{z}\mathrm{\neq}0)$, with $j=1,2,$ or $3$, is

$\displaystyle\partial_{t}\beta_{j}=\gamma_{j}\beta_{j}+\sum_{\begin{subarray}{c}\bm{k}^{\prime},m,n\end{subarray}}C_{jmn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime}+\sum_{\bm{k}^{\prime},n}\Bigg{\\{}$
$\displaystyle\left[C_{j\mathcal{Y}n}^{(\bm{k},\bm{k}^{\prime})}+C_{jn\mathcal{Y}}^{(\bm{k},\bm{k}^{\prime\prime})}\right]\mathcal{Y}^{\prime}\beta_{n}^{\prime\prime}$
(A8)
$\displaystyle+\left[C_{j\mathcal{Z}n}^{(\bm{k},\bm{k}^{\prime})}+C_{jn\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime\prime})}\right]\mathcal{Z}^{\prime}\beta_{n}^{\prime\prime}$
$\displaystyle+\left[C_{j\Theta
n}^{(\bm{k},\bm{k}^{\prime})}+C_{jn\Theta}^{(\bm{k},\bm{k}^{\prime\prime})}\right]\Theta^{\prime}\beta_{n}^{\prime\prime}\Bigg{\\}},$

whereas, at $\mathbf{k}=(k_{x},k_{z}\mathrm{=}0)$, one finds that the
fluctuation amplitude evolves according to Eq. (3b), which is expanded below

$\displaystyle\partial_{t}\mathcal{Y}$
$\displaystyle=-\gamma_{\mathcal{Y}}\mathcal{Y}+\sum_{\bm{k}^{\prime},m,n:k_{z}\mathrm{=}0}C_{\mathcal{Y}mn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime},$
(A9a) $\displaystyle\partial_{t}\mathcal{Z}$
$\displaystyle=-\gamma_{\mathcal{Z}}\mathcal{Z}+\sum_{\bm{k}^{\prime},m,n:k_{z}\mathrm{=}0}C_{\mathcal{Z}mn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime},$
(A9b) $\displaystyle\partial_{t}\Theta$
$\displaystyle=-\gamma_{\Theta}\Theta+\sum_{\bm{k}^{\prime},m,n:k_{z}\mathrm{=}0}C_{\Theta
mn}^{(\bm{k},\bm{k}^{\prime})}\beta_{m}^{\prime}\beta_{n}^{\prime\prime},$
(A9c)

where $\gamma_{\mathcal{Y}}=\gamma_{\mathcal{Z}}=\nu k^{2}$ and
$\gamma_{\Theta}=\kappa k^{2}$ are the damping rates.

### C2. Eigenmode-energy evolution

To derive an evolution equation for energy in the $j^{\mathrm{th}}$ eigenmode
at $\mathbf{k}=(k_{x},k_{z}\mathrm{\neq}0)$, we multiply Eq. (A8) with
$\beta_{j}^{\ast}$ and add a complex conjugate of the resulting equation to
arrive at

$\displaystyle\partial_{t}|\beta_{j}|^{2}$
$\displaystyle=2\mathrm{Re}\,\gamma_{j}|\beta_{j}|^{2}+\sum_{\begin{subarray}{c}\bm{k}^{\prime},m,n\end{subarray}}2\mathrm{Re}\,\left[C_{jmn}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{m}^{\prime}\beta_{n}^{\prime\prime}\beta_{j}^{\ast}\rangle\right]$
(A10) $\displaystyle\hskip
71.13188pt+\sum_{\bm{k}^{\prime}}2\mathrm{Re}\,\Bigg{\\{}\left[C_{j\mathcal{Y}j}^{(\bm{k},\bm{k}^{\prime})}\langle\mathcal{Y}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle+C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime})}\langle\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle+C_{j\Theta
j}^{(\bm{k},\bm{k}^{\prime})}\langle\Theta^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle\right]\Big{|}_{k_{z}^{\prime}\mathrm{=}0}\Bigg{.}$
$\displaystyle\hskip
123.48485pt\Bigg{.}+\left[C_{jj\mathcal{Y}}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{j}^{\prime}\mathcal{Y}^{\prime\prime}\beta_{j}^{\ast}\rangle+C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{j}^{\prime}\mathcal{Z}^{\prime\prime}\beta_{j}^{\ast}\rangle+C_{jj\Theta}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{j}^{\prime}\Theta^{\prime\prime}\beta_{j}^{\ast}\rangle\right]\Big{|}_{k_{z}^{\prime}\mathrm{=}k_{z}}\Bigg{\\}}.$

Similar equations can be derived for fluctuation energy at GSF-stable
wavenumbers $\mathbf{k}=(k_{x},k_{z}\mathrm{=}0)$ using Eqs. (A9a)–(A9c):

$\displaystyle\partial_{t}|\mathcal{Y}|^{2}$
$\displaystyle=-2\mathrm{Re}\,\gamma_{\mathcal{Y}}|\mathcal{Y}|^{2}+\sum_{\bm{k}^{\prime},m}2\mathrm{Re}\,\left[C_{\mathcal{Y}mm}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{m}^{\prime}\beta_{m}^{\prime\prime}\mathcal{Y}^{\ast}\rangle\right]\Big{|}_{k_{z}\mathrm{=}0},$
(A11a) $\displaystyle\partial_{t}|\mathcal{Z}|^{2}$
$\displaystyle=-2\mathrm{Re}\,\gamma_{\mathcal{Z}}|\mathcal{Z}|^{2}+\sum_{\bm{k}^{\prime},m}2\mathrm{Re}\,\left[C_{\mathcal{Z}mm}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{m}^{\prime}\beta_{m}^{\prime\prime}\mathcal{Z}^{\ast}\rangle\right]\Big{|}_{k_{z}\mathrm{=}0},$
(A11b) $\displaystyle\partial_{t}|\Theta|^{2}$
$\displaystyle=-2\mathrm{Re}\,\gamma_{\Theta}|\Theta|^{2}+\sum_{\bm{k}^{\prime},m}2\mathrm{Re}\,\left[C_{\Theta
mm}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{m}^{\prime}\beta_{m}^{\prime\prime}\Theta^{\ast}\rangle\right]\Big{|}_{k_{z}\mathrm{=}0}.$
(A11c)

Since numerical simulations inform us that the triplets with a latitudinal
flow at $(k_{x},0)$, i.e., $\mathcal{Z}$, dominates the nonlinear energy
transfer, we may drop nonlinear terms on the right-hand side of Eq. (A10) that
do not involve $\mathcal{Z}$. In the resulting equation, because $\mathcal{Y}$
and $\Theta$ do not appear, Eqs. (A11a) and (A11c) can also be removed, which
allows us to write the following set of equations:

$\displaystyle\partial_{t}|\beta_{j}|^{2}$
$\displaystyle=2\mathrm{Re}\,\gamma_{j}|\beta_{j}|^{2}+\sum_{\bm{k}^{\prime}}2\mathrm{Re}\,\Bigg{\\{}\left[C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime})}\langle\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle\right]\Big{|}_{k_{z}^{\prime}\mathrm{=}0}+\left[C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{j}^{\prime}\mathcal{Z}^{\prime\prime}\beta_{j}^{\ast}\rangle\right]\Big{|}_{k_{z}^{\prime}\mathrm{=}k_{z}}\Bigg{\\}},$
(A12a) $\displaystyle\partial_{t}|\mathcal{Z}|^{2}$
$\displaystyle=-2\mathrm{Re}\,\gamma_{\mathcal{Z}}|\mathcal{Z}|^{2}+\sum_{\bm{k}^{\prime},m}2\mathrm{Re}\,\left[C_{\mathcal{Z}mm}^{(\bm{k},\bm{k}^{\prime})}\langle\beta_{m}^{\prime}\beta_{m}^{\prime\prime}\mathcal{Z}^{\ast}\rangle\right]\Big{|}_{k_{z}\mathrm{=}0}.$
(A12b)

### C3. Triplet correlation evolution

To obtain evolution equations for the terms on the right-hand side of Eqs.
(A12a) and (A12b), we multiply Eq. (3) with two amplitudes. For example, to
determine the evolution of
$\langle\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle$,
first, we multiply Eq. (A8) for $\beta_{j}$ with
$\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}$; second, we multiply another
equation, similar to Eq. (A8), but for $\beta_{j}^{\prime\prime}$ with
$\mathcal{Z}^{\prime}\beta_{j}^{\ast}$; and, finally, we multiply Eq. (A9b),
for $\mathcal{Z}^{\prime}$, with $\beta_{j}^{\prime\prime}\beta_{j}^{\ast}$,
and add them all together. Using this, we provide below an example evolution
equation for triplet correlations:

$\displaystyle\left[\partial_{t}-\left(-\gamma_{\mathcal{Z}}^{\prime}+\gamma_{j}^{\prime\prime}+\gamma_{j}^{\ast}\right)\right]\langle{\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime}\mathrm{=}0}=$
$\displaystyle\Bigg{\\{}\sum_{\bm{k}^{\prime\prime\prime},m}\left[C_{{\mathcal{Z}}mm}^{(\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\beta_{m}^{\prime\prime\prime}\beta_{m}(\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime})\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle\right]$
(A13)
$\displaystyle+\sum_{k_{x}^{\prime\prime\prime}}\left[C_{j\mathcal{Y}j}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\mathcal{Y}^{\prime\prime\prime}\beta_{j}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}\right.$
$\displaystyle\hskip
25.6073pt+C_{j\mathcal{Z}j}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\mathcal{Z}^{\prime\prime\prime}\beta_{j}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}$
$\displaystyle\hskip 25.6073pt+C_{j\Theta
j}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\Theta^{\prime\prime\prime}\beta_{j}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}$
$\displaystyle\hskip
25.6073pt+C_{jj\mathcal{Y}}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\beta_{j}^{\prime\prime\prime}\mathcal{Y}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}$
$\displaystyle\hskip
25.6073pt+C_{jj\mathcal{Z}}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\beta_{j}^{\prime\prime\prime}\mathcal{Z}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}$
$\displaystyle\hskip
25.6073pt\left.+C_{jj\Theta}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\beta_{j}^{\prime\prime\prime}\Theta(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}\right]$
$\displaystyle+\sum_{k_{x}^{\prime\prime\prime}}\left[C_{j\mathcal{Y}j}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\mathcal{Y}^{\prime\prime\prime\ast}\beta_{j}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}\right.$
$\displaystyle\hskip
25.6073pt+C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\mathcal{Z}^{\prime\prime\prime\ast}\beta_{j}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}$
$\displaystyle\hskip 25.6073pt+C_{j\Theta
j}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\Theta^{\prime\prime\prime\ast}\beta_{j}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}$
$\displaystyle\hskip
25.6073pt+C_{jj\mathcal{Y}}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\beta_{j}^{\prime\prime\prime\ast}\mathcal{Y}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}$
$\displaystyle\hskip
25.6073pt+C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\beta_{j}^{\prime\prime\prime\ast}\mathcal{Z}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}$
$\displaystyle\hskip
25.6073pt\Bigg{.}\left.\left.+C_{jj\Theta}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\beta_{j}^{\prime\prime\prime\ast}\Theta^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime}){\mathcal{Z}}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}\right]\Bigg{\\}}\right|_{k_{z}^{\prime}\mathrm{=}0}.$

Since the numerical simulations show dominant coupling between the latitudinal
flow $\mathcal{Z}$ and two unstable modes in a bath of turbulent interactions,
the interactions involving unstable modes only can be excluded. This
immediately implies that the first term on the right-hand side of Eq. (A13)
can be dropped in the face of remaining dominant terms. Among the remaining
terms, the fourth-order correlations that have different components of
velocity, e.g., the terms where $\mathcal{Y}$ and $\mathcal{Z}$ appear
together on the right-hand side of Eq. (A13), do not form terms that appear in
the definition of energy. Guided by numerical simulations where nonlinear
coupling to $\mathcal{Y}$ and $\Theta$ are unimportant, only the energy(-like)
terms with $\mathcal{Z}$ will be kept henceforth. The resulting equation then
reads

$\displaystyle\left[\partial_{t}-\left(-\gamma_{\mathcal{Z}}^{\prime}+\gamma_{j}^{\prime\prime}+\gamma_{j}^{\ast}\right)\right]\langle\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime}\mathrm{=}0}=$
$\displaystyle\Bigg{\\{}\sum_{k_{x}^{\prime\prime\prime}}\left[C_{j\mathcal{Z}j}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\mathcal{Z}^{\prime\prime\prime}\beta_{j}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime})\mathcal{Z}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}\right.$
(A14) $\displaystyle\hskip
25.6073pt+C_{jj\mathcal{Z}}^{(\bm{k}-\bm{k}^{\prime},\bm{k}^{\prime\prime\prime})}\langle\beta_{j}^{\prime\prime\prime}\mathcal{Z}(\bm{k}-\bm{k}^{\prime}-\bm{k}^{\prime\prime\prime})\mathcal{Z}^{\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}$
$\displaystyle\hskip
25.6073pt+C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\mathcal{Z}^{\prime\prime\prime\ast}\beta_{j}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime})\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}0}$
$\displaystyle\hskip
25.6073pt\left.\Bigg{.}\left.+C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime\prime\prime})\ast}\langle\beta_{j}^{\prime\prime\prime\ast}\mathcal{Z}^{\ast}(\bm{k}-\bm{k}^{\prime\prime\prime})\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\rangle|_{k_{z}^{\prime\prime\prime}\mathrm{=}k_{z}}\right]\Bigg{\\}}\right|_{k_{z}^{\prime}\mathrm{=}0}.$

The same procedure is then repeated to find evolution equations for the other
two triplet correlations that appear in Eqs. (A12a) and (A12b).

### C4: Quadruplet correlations and Statistical closure approximation

Equation (A14) can be solved using the technique of Green’s function inversion
and Markovianization, a standard step in EDQNM closure [although, here, we do
not modify the growth rate $\gamma$’s with the amplitude-dependent nonlinear
frequency, an approximation justifiable for the low-wavenumber regime; for
more details, see Terry et al. (2018)]. Such a solution yields

$\displaystyle\langle\mathcal{Z}^{\prime}\beta_{j}^{\prime\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime}\mathrm{=}0}=-\Bigg{\\{}\left(-\gamma_{\mathcal{Z}}^{\prime}+\gamma_{j}^{\prime\prime}+\gamma_{j}^{\ast}\right)^{-1}|\mathcal{Z}^{\prime}|^{2}$
$\displaystyle\Big{[}\left(C_{j\mathcal{Z}j}^{(\bm{k}^{\prime\prime},-\bm{k}^{\prime})}+C_{jj\mathcal{Z}}^{(\bm{k}^{\prime\prime},\bm{k})}\right)|\beta_{j}|^{2}+\left(C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime})\ast}+C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime\prime})\ast}\right)|\beta_{j}^{\prime\prime}|^{2}\Big{]}\Bigg{\\}}\Bigg{|}_{k_{z}^{\prime}\mathrm{=}0}.$
(A15)

Similarly, the other two triplet correlations that appear in Eqs. (A12a) and
(A12b) can also be solved to obtain

$\displaystyle\langle\beta_{j}^{\prime}\mathcal{Z}^{\prime\prime}\beta_{j}^{\ast}\rangle|_{k_{z}^{\prime}\mathrm{=}k_{z}}=-\Bigg{\\{}\left(\gamma_{j}^{\prime}-\gamma_{\mathcal{Z}}^{\prime\prime}+\gamma_{j}^{\ast}\right)^{-1}|\mathcal{Z}^{\prime\prime}|^{2}$
$\displaystyle\Big{[}\left(C_{j\mathcal{Z}j}^{(\bm{k}^{\prime},-\bm{k}^{\prime\prime})}+C_{jj\mathcal{Z}}^{(\bm{k}^{\prime},\bm{k})}\right)|\beta_{j}|^{2}+\left(C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime\prime})\ast}+C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime})\ast}\right)|\beta_{j}^{\prime}|^{2}\Big{]}\Bigg{\\}}\Bigg{|}_{k_{z}^{\prime}\mathrm{=}k_{z}},$
(A16)

and

$\displaystyle\langle\beta_{j}^{\prime}\beta_{j}^{\prime\prime}\mathcal{Z}^{\ast}\rangle|_{k_{z}\mathrm{=}0}=-\Bigg{\\{}\left(\gamma_{j}^{\prime}+\gamma_{j}^{\prime\prime}-\gamma_{\mathcal{Z}}^{\ast}\right)^{-1}|\mathcal{Z}|^{2}$
$\displaystyle\Big{[}\left(C_{j\mathcal{Z}j}^{(\bm{k}^{\prime},\bm{k})}+C_{jj\mathcal{Z}}^{(\bm{k}^{\prime},-\bm{k}^{\prime\prime})}\right)|\beta_{j}^{\prime\prime}|^{2}+\left(C_{j\mathcal{Z}j}^{(\bm{k}^{\prime\prime},\bm{k})}+C_{jj\mathcal{Z}}^{(\bm{k}^{\prime\prime},-\bm{k}^{\prime})}\right)|\beta_{j}^{\prime}|^{2}\Big{]}\Bigg{\\}}\Bigg{|}_{k_{z}\mathrm{=}0}.$
(A17)

### C5: A set of EDQNM-closed energy evolution equations

Solutions of the triplet correlations from Eqs. (A15)–(A17) are now
substituted into Eqs. (A12a) and (A12b), which results in the following set of
closed equations:

$\displaystyle\partial_{t}|\beta_{j}|^{2}=$ $\displaystyle
2\mathrm{Re}\,\gamma_{j}|\beta_{j}|^{2}-\sum_{\bm{k}^{\prime}}2\mathrm{Re}\,\Bigg{[}\frac{C_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime})}|\mathcal{Z}^{\prime}|^{2}}{\left(-\gamma_{\mathcal{Z}}^{\prime}+\gamma_{j}^{\prime\prime}+\gamma_{j}^{\ast}\right)}\left\\{\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k}^{\prime\prime},-\bm{k}^{\prime})}|\beta_{j}|^{2}+\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime})\ast}|\beta_{j}^{\prime\prime}|^{2}\right\\}\Bigg{]}\Bigg{|}_{k_{z}^{\prime}\mathrm{=}0}$
(A18) $\displaystyle\hskip
55.48277pt-\sum_{\bm{k}^{\prime}}2\mathrm{Re}\,\Bigg{[}\frac{C_{jj\mathcal{Z}}^{(\bm{k},\bm{k}^{\prime})}|\mathcal{Z}^{\prime\prime}|^{2}}{\left(\gamma_{j}^{\prime}-\gamma_{\mathcal{Z}}^{\prime\prime}+\gamma_{j}^{\ast}\right)}\left\\{\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k}^{\prime},-\bm{k}^{\prime\prime})}|\beta_{j}|^{2}+\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k},\bm{k}^{\prime\prime})\ast}|\beta_{j}^{\prime}|^{2}\right\\}\Bigg{]}\Bigg{|}_{k_{z}^{\prime}\mathrm{=}k_{z}},$

$\partial_{t}|\mathcal{Z}|^{2}=-2\mathrm{Re}\,\gamma_{\mathcal{Z}}|\mathcal{Z}|^{2}-|\mathcal{Z}|^{2}\sum_{\bm{k}^{\prime}}\sum_{j=1}^{3}2\mathrm{Re}\,\Bigg{.}\left[\frac{C_{\mathcal{Z}jj}^{(\bm{k},\bm{k}^{\prime})}}{\left(\gamma_{j}^{\prime}+\gamma_{j}^{\prime\prime}-\gamma_{\mathcal{Z}}^{\ast}\right)}\left\\{\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k}^{\prime},\bm{k})}|\beta_{j}^{\prime\prime}|^{2}+\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{j\mathcal{Z}j}^{(\bm{k}^{\prime\prime},\bm{k})}|\beta_{j}^{\prime}|^{2}\right\\}\right]\Bigg{|}_{k_{z}\mathrm{=}0},$
(A19)

where $\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{lmn}^{(\bm{p},\bm{q})}=C_{lmn}^{(\bm{p},\bm{q})}+C_{lnm}^{(\bm{p},\bm{p}-\bm{q})}$
is the symmetrized coupling coefficient.

Equation (A19) can be simplified by considering that it is the pairs of
unstable modes ($j=1$) that excite the latitudinal flow $\mathcal{Z}$:

$\partial_{t}|\mathcal{Z}|^{2}/2=-\gamma_{\mathcal{Z}}|\mathcal{Z}|^{2}+|\mathcal{Z}|^{2}\sum_{\bm{k}^{\prime}}|\beta_{1}^{\prime}|^{2}\,\mathrm{Re}\,\Bigg{.}\left[\frac{-\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{\mathcal{Z}11}^{(\bm{k},\bm{k}^{\prime})}\stackon[-4.2pt]{C}{\,\smash{\stackon[-1.95pt]{\mathchar
382\relax}{\rotatebox{180.0}{$\mathchar
382\relax$}}}}_{1\mathcal{Z}1}^{(\bm{k}^{\prime\prime},\bm{k})}}{\gamma_{1}^{\prime}+\gamma_{1}^{\prime\prime}-\gamma_{\mathcal{Z}}^{\ast}}\right],$
(A20)

where the term inside the wavenumber-summation in Eq. (A19) has been
symmetrized.

## References

  * Aerts et al. (2019) Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35, doi: 10.1146/annurev-astro-091918-104359
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1, doi: 10.1103/RevModPhys.70.1
  * Barker et al. (2019) Barker, A. J., Jones, C. A., & Tobias, S. M. 2019, MNRAS, 487, 1777, doi: 10.1093/mnras/stz1386
  * Barker et al. (2020) —. 2020, MNRAS, 495, 1468, doi: 10.1093/mnras/staa1327
  * Barker & Latter (2015) Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21, doi: 10.1093/mnras/stv640
  * Beck et al. (2012) Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55, doi: 10.1038/nature10612
  * Biskamp & Zeiler (1995) Biskamp, D., & Zeiler, A. 1995, Phys. Rev. Lett., 74, 706, doi: 10.1103/PhysRevLett.74.706
  * Brown et al. (2013) Brown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34, doi: 10.1088/0004-637X/768/1/34
  * Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286, doi: 10.1126/science.272.5266.1286
  * Denissenkov (2010) Denissenkov, P. A. 2010, ApJ, 723, 563, doi: 10.1088/0004-637X/723/1/563
  * Du et al. (2023) Du, S., Li, H., Fu, X., & Gan, Z. 2023, ApJ, 948, 72, doi: 10.3847/1538-4357/acc5e9
  * Dymott et al. (2023) Dymott, R. W., Barker, A. J., Jones, C. A., & Tobias, S. M. 2023, MNRAS, doi: 10.1093/mnras/stad1982
  * Eggenberger et al. (2012) Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4, doi: 10.1051/0004-6361/201219729
  * Fraser et al. (2023) Fraser, A. E., Reifenstein, S. A., & Garaud, P. 2023, arXiv e-prints, arXiv:2302.11610, doi: 10.48550/arXiv.2302.11610
  * Fraser et al. (2021) Fraser, A. E., Terry, P. W., Zweibel, E. G., Pueschel, M. J., & Schroeder, J. M. 2021, Physics of Plasmas, 28, 022309, doi: 10.1063/5.0034575
  * Fricke (1968) Fricke, K. 1968, ZAp, 68, 317
  * Fuller et al. (2019) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661, doi: 10.1093/mnras/stz514
  * Garaud (2018) Garaud, P. 2018, Annual Review of Fluid Mechanics, 50, 275, doi: 10.1146/annurev-fluid-122316-045234
  * Garaud & Brummell (2015) Garaud, P., & Brummell, N. 2015, ApJ, 815, 42, doi: 10.1088/0004-637X/815/1/42
  * Garaud & Garaud (2008) Garaud, P., & Garaud, J. D. 2008, MNRAS, 391, 1239, doi: 10.1111/j.1365-2966.2008.13930.x
  * Goldreich & Schubert (1967) Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571, doi: 10.1086/149360
  * Goodman & Xu (1994) Goodman, J., & Xu, G. 1994, ApJ, 432, 213, doi: 10.1086/174562
  * Harrington & Garaud (2019) Harrington, P. Z., & Garaud, P. 2019, ApJ, 870, L5, doi: 10.3847/2041-8213/aaf812
  * Hegna et al. (2018) Hegna, C. C., Terry, P. W., & Faber, B. J. 2018, Physics of Plasmas, 25, 022511, doi: 10.1063/1.5018198
  * Knobloch (1982) Knobloch, E. 1982, Geophys. Astrophys. Fluid Dyn., 22, 133, doi: 10.1080/03091928208221740
  * Knobloch & Spruit (1982) Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261
  * Korycansky (1991) Korycansky, D. G. 1991, ApJ, 381, 515, doi: 10.1086/170675
  * Latter & Papaloizou (2018) Latter, H. N., & Papaloizou, J. 2018, MNRAS, 474, 3110, doi: 10.1093/mnras/stx3031
  * Lesur & Longaretti (2005) Lesur, G., & Longaretti, P. Y. 2005, A&A, 444, 25, doi: 10.1051/0004-6361:20053683
  * Lesur et al. (2023) Lesur, G., Flock, M., Ercolano, B., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 465
  * Li et al. (2023) Li, P.-Y., Pueschel, M. J., Terry, P. W., & Whelan, G. G. 2023, Nuclear Fusion, 63, 026028, doi: 10.1088/1741-4326/acae37
  * Li et al. (2021) Li, P.-Y., Terry, P. W., Whelan, G. G., & Pueschel, M. J. 2021, Physics of Plasmas, 28, 102507, doi: 10.1063/5.0059335
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Ng & Bhattacharjee (1996) Ng, C. S., & Bhattacharjee, A. 1996, ApJ, 465, 845, doi: 10.1086/177468
  * Orszag (1970) Orszag, S. A. 1970, Journal of Fluid Mechanics, 41, 363, doi: 10.1017/S0022112070000642
  * Park et al. (2020) Park, J., Prat, V., & Mathis, S. 2020, A&A, 635, A133, doi: 10.1051/0004-6361/201936863
  * Park et al. (2021) Park, J., Prat, V., Mathis, S., & Bugnet, L. 2021, A&A, 646, A64, doi: 10.1051/0004-6361/202038654
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pueschel et al. (2021) Pueschel, M. J., Li, P. Y., & Terry, P. W. 2021, Nuclear Fusion, 61, 054003, doi: 10.1088/1741-4326/abf175
  * Radko & Smith (2012) Radko, T., & Smith, D. P. 2012, Journal of Fluid Mechanics, 692, 5, doi: 10.1017/jfm.2011.343
  * Rashid et al. (2008) Rashid, F. Q., Jones, C. A., & Tobias, S. M. 2008, A&A, 488, 819, doi: 10.1051/0004-6361:200810039
  * Riley & Lelong (2000) Riley, J. J., & Lelong, M.-P. 2000, Annual Review of Fluid Mechanics, 32, 613, doi: 10.1146/annurev.fluid.32.1.613
  * Smith & Waleffe (1999) Smith, L. M., & Waleffe, F. 1999, Physics of Fluids, 11, 1608, doi: 10.1063/1.870022
  * Spiegel & Veronis (1960) Spiegel, E. A., & Veronis, G. 1960, ApJ, 131, 442, doi: 10.1086/146849
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923, doi: 10.1051/0004-6361:20011465
  * Tassin et al. (2023) Tassin, T., Gastine, T., & Fournier, A. 2023, arXiv e-prints, arXiv:2309.01602, doi: 10.48550/arXiv.2309.01602
  * Terry (2004) Terry, P. W. 2004, Phys. Rev. Lett., 93, 235004, doi: 10.1103/PhysRevLett.93.235004
  * Terry (2019) Terry, P. W. 2019, Theory, ed. B. Galperin & P. L. Read (Cambridge University Press), 181–193, doi: 10.1017/9781107358225.010
  * Terry et al. (2018) Terry, P. W., Faber, B. J., Hegna, C. C., et al. 2018, Physics of Plasmas, 25, 012308, doi: 10.1063/1.5007062
  * Terry et al. (2021) Terry, P. W., Li, P. Y., Pueschel, M. J., & Whelan, G. G. 2021, Phys. Rev. Lett., 126, 025004, doi: 10.1103/PhysRevLett.126.025004
  * Tripathi et al. (2022a) Tripathi, B., Fraser, A. E., Terry, P. W., Zweibel, E. G., & Pueschel, M. J. 2022a, Physics of Plasmas, 29, 070701, doi: 10.1063/5.0096886
  * Tripathi et al. (2022b) —. 2022b, Physics of Plasmas, 29, 092301, doi: 10.1063/5.0101434
  * Tripathi et al. (2023a) Tripathi, B., Fraser, A. E., Terry, P. W., et al. 2023a, Physics of Plasmas, 30, 072107, doi: 10.1063/5.0156560
  * Tripathi et al. (2023b) Tripathi, B., Terry, P. W., Fraser, A. E., Zweibel, E. G., & Pueschel, M. J. 2023b, Physics of Fluids, 35, 105151, doi: 10.1063/5.0167092
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * Waleffe (1993) Waleffe, F. 1993, Physics of Fluids A, 5, 677, doi: 10.1063/1.858651
  * Wood & McIntyre (2011) Wood, T. S., & McIntyre, M. E. 2011, Journal of Fluid Mechanics, 677, 445, doi: 10.1017/jfm.2011.93