Predicting the Slowing of Stellar Differential Rotation by Instability-driven Turbulence

A linear analysis of Equations (1a)–(1b) for axisymmetric (uniform-in-y) perturbations yields a simple matrix equation, which upon Fourier-transforming becomes ${\boldsymbol{L}}\hat{{\boldsymbol{X}}}=\gamma \hat{{\boldsymbol{X}}}$, where $\hat{{\boldsymbol{X}}}={[{\hat{u}}_{x},{\hat{u}}_{y},{\hat{u}}_{z},\hat{\theta }]}^{{\rm{T}}}$, with ^T as the transpose operation, is the state vector of spatially Fourier-transformed components at wavevector k = (k_x , k_z ); the matrix L is a linear operator, whose eigenvalues are the complex-valued growth rates γ. The size of L demands four linearly independent eigenvectors. Because of the additional constraint ∇ · u = 0, the system has only 3 degrees of freedom at any given wavenumber—$\hat{\theta }$ and two components of velocity. Hence, one eigenvector among the four eigenvectors does not satisfy ∇ · u = 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 )\gt 0$, where $\mathrm{Re}\,$ denotes the real part), whenever r ∈ [0, 1). The remaining two eigenvectors are always stable, and their eigenvalues are complex conjugates of each other whenever they satisfy $\mathrm{Im}\,(\gamma )\ne 0$, where $\mathrm{Im}\,$ denotes the imaginary part; $\mathrm{Im}\,(\gamma )$ corresponds to the frequency of inertial gravity waves (IGWs; or gravito-inertial waves), modified by the shear flow and damped by viscous and thermal diffusion.

The GSF instability grows dominantly via axisymmetric (∂_y ≡ 0) perturbations, therefore we focus upon the (x, z)-variations of the three-component velocity and temperature fields. The dispersion relation is then a simple cubic polynomial in γ (Goldreich & Schubert 1967), as

$$\begin{eqnarray}&&{\gamma }_{\nu }^{2}{\gamma }_{\kappa }+\displaystyle \frac{{\kappa }_{\mathrm{ep}}^{2}{k}_{z}^{2}}{{k}^{2}}{\gamma }_{\kappa }+\displaystyle \frac{{N}^{2}{k}_{z}^{2}}{{k}^{2}}{\gamma }_{\nu }=0,\end{eqnarray} \tag{ 2 }$$

where γ_ν = γ + ν k² and γ_κ = γ + κ k². Equation (2) shows that on the (k_x , k_z )-plane, the growth rate exhibits strong anisotropy: fluctuations with k_x = 0 ("elevator modes") grow the fastest, whereas those with k_z = 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 = 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 3D Kelvin–Helmholtz instability (Tripathi et al. 2023b), rotating (Waleffe 1993; Smith & Waleffe 1999) and stably stratified turbulence (Riley & Lelong 2000), and 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 L , we can decompose the arbitrary incompressible fluctuation ${\hat{{\boldsymbol{X}}}}_{\mathrm{arb}}$ with k_z ≠ 0 as ${\hat{{\boldsymbol{X}}}}_{\mathrm{arb}}={\sum }_{j=1}^{3}{\beta }_{j}{\hat{X}}_{j}$, where β_j is the amplitude of the jth eigenvector ${\hat{X}}_{j}$; we reserve j = 1 for the GSF-unstable modes and j = 2, 3 for the IGWs that are always linearly stable in this study. More compactly, ${\hat{{\boldsymbol{X}}}}_{\mathrm{arb}}={\boldsymbol{E}}{\boldsymbol{\beta }}$, where β is a (column) vector of mode amplitudes and E is an eigenvector matrix, whose jth column is ${\hat{X}}_{j}$. Thus, ${\boldsymbol{\beta }}={{\boldsymbol{E}}}^{-1}{\hat{{\boldsymbol{X}}}}_{\mathrm{arb}}$.

For the k_z = 0 modes, 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 = 0 into ${\hat{{\boldsymbol{X}}}}_{\mathrm{arb}}^{{\rm{T}}}$ = ${ \mathcal X }[1,0,0,0]$ + ${ \mathcal Y }[0,1,0,0]$ + ${ \mathcal Z }[0,0,1,0]$ + Θ[0, 0, 0, 1], where the amplitudes of the eigenvectors are denoted by ${ \mathcal X },{ \mathcal Y },{ \mathcal Z }$, and Θ. We reserve the β -notation above for the amplitudes of eigenvectors with k_z ≠ 0.

We perform an ensemble of direct numerical simulations of Equations (1a) and (1b), by seeding a low-amplitude solenoidal random noise to u , in a box of size (L_x , L_z ) = (100, 100). To obtain numerically converged results, a spatial resolution of up to 512² 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 ={\sum }_{{k}_{x},{k}_{z}}{\sum }_{m,n}2\mathrm{Re}\left[{\beta }_{m}{\hat{u}}_{x,m}{\beta }_{n}^{* }{\hat{u}}_{y,n}^{* }\right]$, where 〈·〉 is an (x, z)-averaging operation; m and n are summed from 1 to 3, corresponding to three excited eigenmodes at every wavenumber k ; the amplitude β_m and the x-component of the velocity ${\hat{u}}_{x,m}$ correspond to the mth eigenvector at k , and likewise for β_n and ${\hat{u}}_{y,n}$; and the operation * denotes complex conjugation. Using such a decomposition, we obtain the contribution of an unstable mode at k to the momentum transport rate, which is $2| {\beta }_{1}{| }^{2}\mathrm{Re}\left[{\hat{u}}_{x,1}{\hat{u}}_{y,1}^{* }\right]$. This decomposition is performed for every wavenumber, hence allowing us to trace the evolution of transport contributions due to individual unstable modes (see Figure 2(a)).

Zoom In Zoom Out Reset image size

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 Figure 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 = 0 transport significantly less momentum than the other wavenumbers with k_x ≠ 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.

To analyze nonlinear mode couplings, Equations (1a)–(1b) are first spatially Fourier-transformed: ${\partial }_{t}\hat{{\boldsymbol{X}}}={\boldsymbol{L}}\hat{{\boldsymbol{X}}}\,+{\sum }_{{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}^{\prime\prime} }{\boldsymbol{N}}(\hat{{\boldsymbol{X}}}^{\prime} ,\hat{{\boldsymbol{X}}}^{\prime\prime} )$, where $\hat{{\boldsymbol{X}}},\hat{{\boldsymbol{X}}}^{\prime}$, and $\hat{{\boldsymbol{X}}}^{\prime\prime}$ are state vectors at ${\boldsymbol{k}},{\boldsymbol{k}}^{\prime}$, and k '', respectively, satisfying ${\boldsymbol{k}}={\boldsymbol{k}}^{\prime} +{\boldsymbol{k}}^{\prime\prime}$. Then, following Section 2.1, we substitute $\hat{{\boldsymbol{X}}}={\boldsymbol{E}}{\boldsymbol{\beta }}$, and likewise for $\hat{{\boldsymbol{X}}}^{\prime}$ and $\hat{{\boldsymbol{X}}}^{\prime\prime}$. We multiply the obtained equation with E ⁻¹ and take the jth row of the resulting equation. This process yields an evolution equation for the jth eigenmode at k . Such an evolution equation for mode amplitude β_j for k_z ≠ 0 is different from that of the mode amplitude $F\in \{{ \mathcal X },{ \mathcal Y },{ \mathcal Z },{\rm{\Theta }}\}$ for k_z = 0, although both are coupled and nonlinear:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\beta }_{j} & = & {\gamma }_{j}{\beta }_{j}+\displaystyle \sum _{{\boldsymbol{k}}^{\prime} ,m,n}{C}_{{jmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} }\\ & & +\,\displaystyle \sum _{\mathop{{\boldsymbol{k}}^{\prime} ,F,n:{k}_{z}^{{\prime} }=0}\limits_{F\in \{{ \mathcal X },{ \mathcal Y },{ \mathcal Z },{\rm{\Theta }}\}}}\left[{C}_{{jFn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}+{C}_{{jnF}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )}\right]F^{\prime} {\beta }_{n}^{{\prime\prime} },\end{array}\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&{\partial }_{t}F=-{\gamma }_{F}F+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m,n:{k}_{z}=0}{C}_{{Fmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} },\end{eqnarray} \tag{ 3b }$$

where γ_j is the complex-valued growth rate for the jth eigenmode with k_z ≠ 0; the real-valued damping rate γ_F is ${\gamma }_{{ \mathcal X }}$ when F is replaced with ${ \mathcal X }$ in Equation (3b); likewise for the replacement of F with ${ \mathcal Y },{ \mathcal Z },$ and Θ; and we note that ${\gamma }_{{ \mathcal X }}={\gamma }_{{ \mathcal Y }}={\gamma }_{{ \mathcal Z }}=\nu {k}_{x}^{2}$ and ${\gamma }_{{\rm{\Theta }}}=\kappa {k}_{x}^{2}$. The nonlinear coupling coefficient, for example, ${C}_{{jmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}$, measures the overlap of eigenmodes m with ${\boldsymbol{k}}^{\prime}$, n with k '', and j with k . Such a mode-coupling coefficient is found by applying E ⁻¹ to the (column) vector of ${\boldsymbol{N}}({\hat{{\boldsymbol{X}}}}_{m}^{{\prime} },{\hat{{\boldsymbol{X}}}}_{n}^{{\prime\prime} })$, a process that incorporates all the nonlinearities of the system, thus making ${C}_{{jmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}$ ideal for a comprehensive instability saturation analysis.

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

The energy evolution equation for each eigenmode can now be derived by multiplying Equation (3a) by ${\beta }_{j}^{* }$ and adding the complex conjugate of the resulting equation, to arrive at

$$\begin{eqnarray}&&{\partial }_{t}| {\beta }_{j}{| }^{2}={Q}_{j}+{T}_{j\mathrm{AA}},\end{eqnarray} \tag{ 4 }$$

where ${Q}_{j}=2\mathrm{Re}\,{\gamma }_{j}| {\beta }_{j}{| }^{2}$ is the linear energy transfer rate to k from the mean gradients, and T_jAA is the total nonlinear energy transfer to the jth eigenmode from all possible nonlinear interactions; Equation (C3). We now show the spectrum of time-averaged Q₁, along with that of the growth rate and time-averaged viscous dissipation rate _ν , in Figure 3.

Zoom In Zoom Out Reset image size

From T_jAA in Equation (4), we separate out the nonlinear transfer ${T}_{1{ \mathcal Z }1}$ in a triad that involves a latitudinal flow ${ \mathcal Z }$ at ${k}_{z}^{{\prime} }=0$ and two GSF-unstable modes (j = 1) at ${k}_{z}^{{\prime} }\ne 0$:

$$\begin{eqnarray}&&{T}_{1{ \mathcal Z }1}=\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }:{k}_{z}^{{\prime} }=0}2\mathrm{Re}\,\left\{\left[{C}_{1{ \mathcal Z }1}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}+{C}_{11{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )}\right]{ \mathcal Z }^{\prime} {\beta }_{1}^{{\prime\prime} }{\beta }_{1}^{* }\right\}.\end{eqnarray} \tag{ 5 }$$

We now compare ${T}_{1{ \mathcal Z }1}$ with T_1AA in Figure 4, for a GSF-unstable mode with a wavenumber that contributes the most 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 = 2π/L_x ; see the right-column snapshot at t = 6000 in Figure 1. These flows generalize to meridional circulation in stars and resemble zonal jets in planetary atmospheres and fusion systems (Terry 2019). Flows with k_z = 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 = 0 and at low k_z , as seen in Figure 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 = 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\partial }_{t}\beta =...\beta +...\beta \beta +...\beta { \mathcal X }+...\beta { \mathcal Y }+...\beta { \mathcal Z }+...\beta {\rm{\Theta }},\end{eqnarray} \tag{ 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 an evolution equation for the second-order correlator as ${\partial }_{t}(\beta \beta )=...\beta \beta +...\beta { \mathcal Z }\beta$, where the other nonlinear terms, e.g., β β β 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 (Figure 4). Since $\beta { \mathcal Z }\beta$ also evolves, one can similarly obtain an evolution equation for the 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 β, needed for the turbulent transport prediction.

To predict the mode amplitude β, 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, ∂_t ∼ 0, and thus $\beta \beta { \mathcal Z }{ \mathcal Z }=...{ \mathcal Z }{ \mathcal Z }$. The EDQNM closure allows us to write of a fourth-order correlator $\beta \beta { \mathcal Z }{ \mathcal Z }$ as a sum of the products of the second-order correlators, such as $| \beta {| }^{2}| { \mathcal Z }{| }^{2}$. Then, canceling $| { \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): β β = .... Using this, one can make predictions for turbulent transport rates, as we shall show in the next subsection.

To make quantitative predictions for transport, we take the EDQNM-closed evolution equation for the latitudinal flow energy (see Appendix B, Equation (C13)):

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| { \mathcal Z }{| }^{2}/2 & = & -{\gamma }_{{ \mathcal Z }}| { \mathcal Z }{| }^{2}\\ & & +\,| { \mathcal Z }{| }^{2}\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}| {\beta }_{1}^{{\prime} }{| }^{2}\,\mathrm{Re}\,\left[-{\tau }_{11{ \mathcal Z }}{\overleftrightarrow{C}}_{{ \mathcal Z }11}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\overleftrightarrow{C}}_{1{ \mathcal Z }1}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}\right],\end{array}\end{eqnarray} \tag{ 7 }$$

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 Equation (7),

$$\begin{eqnarray}&&{\tau }_{11{ \mathcal Z }}={({\gamma }_{1}^{{\prime} }+{\gamma }_{1}^{{\prime\prime} }-{\gamma }_{{ \mathcal Z }}^{* })}^{-1},\end{eqnarray} \tag{ 8 }$$

is the three-wave interaction time found from the EDQNM closure, and ${\overleftrightarrow{C}}_{{lmn}}^{({\boldsymbol{p}},{\boldsymbol{q}})}={C}_{{lmn}}^{({\boldsymbol{p}},{\boldsymbol{q}})}+{C}_{{lnm}}^{({\boldsymbol{p}},{\boldsymbol{p}}-{\boldsymbol{q}})}$ is the symmetrized coupling coefficient. In quasi-stationary turbulence, ∂_t ∼ 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}-{k}_{x}^{\prime} ,-{k}_{z}^{{\prime} })$; then, using Equation (7), $| {\beta }_{1}^{{\prime} }{| }^{2}={\gamma }_{{ \mathcal Z }}{\left(\mathrm{Re}\,[-{\tau }_{11{ \mathcal Z }}{\overleftrightarrow{C}}_{{ \mathcal Z }11}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\overleftrightarrow{C}}_{1{ \mathcal Z }1}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{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 Equation (7) arising from the coupling coefficients and ${\tau }_{11{ \mathcal Z }}$. Such a consideration provides an expression for the nonlinearly saturated squared mode amplitude:

$$\begin{eqnarray}&&| \beta ^{\prime} {| }^{2}={\gamma }_{{ \mathcal Z }}\times {\left(\displaystyle \frac{\gamma }{{k}^{2}}\right)}_{\mathrm{Closure}},\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\gamma }{{k}^{2}}\right)}_{\mathrm{Closure}}=\displaystyle \frac{1}{\left|{\sum }_{{\boldsymbol{k}}^{\prime} }{\tau }_{11{ \mathcal Z }}\,\mathrm{Re}\,[-{\overleftrightarrow{C}}_{{ \mathcal Z }11}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\overleftrightarrow{C}}_{1{ \mathcal Z }1}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}]\right|},\end{eqnarray} \tag{ 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 = 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 = 2), as such an interaction lowers ${\tau }_{12{ \mathcal Z }}$ via both the frequency and damping rate of the IGW. The largest interaction time τ_11Z dominates saturation.

The radial turbulent transport of AM is measured by $\langle {\tilde{u}}_{x}{\tilde{u}}_{y}\rangle \approx {\sum }_{{\boldsymbol{k}}\prime\prime\prime }{\hat{u}}_{x,1}^{\prime\prime\prime }{\hat{u}}_{y,1}^{\prime\prime\prime * }| {\beta }_{1}^{\prime\prime\prime }{| }^{2}$, where ${\beta }_{1}^{\prime\prime\prime }$ is the unstable mode amplitude at k ‴ over which the summation is applied. Then, using $| {\beta }_{1}^{\prime\prime\prime }{| }^{2}$ from the above paragraph,

$$\begin{eqnarray}&&\langle {\tilde{u}}_{x}{\tilde{u}}_{y}{\rangle }_{\mathrm{Closure}}={\left(\displaystyle \frac{\gamma }{{k}^{2}}\right)}_{\mathrm{Closure}}\times {\gamma }_{{ \mathcal Z }}\displaystyle \sum _{{\boldsymbol{k}}\prime\prime\prime }{\hat{u}}_{x,1}^{\prime\prime\prime }{\hat{u}}_{y,1}^{\prime\prime\prime * }.\end{eqnarray} \tag{ 11 }$$

The y-component ${\hat{u}}_{y,1}^{\prime\prime\prime }$ of the 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.

A simple quasi-linear (QL) 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 a "parasitic saturation mechanism" (Goodman & Xu 1994; Radko & Smith 2012; Brown et al. 2013; Harrington & Garaud 2019; Fraser et al. 2024), is based on a single primary mode at the fastest-growing wavenumber k ‴, which predicts the transport rate

$$\begin{eqnarray}&&\langle {\widetilde{u}}_{x}{\widetilde{u}}_{y}{\rangle }_{\mathrm{QL}}=\displaystyle \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 * }{f}^{2}({\boldsymbol{k}}\prime\prime\prime ),\end{eqnarray} \tag{ 12 }$$

whose form is manifestly similar to Equation (11); the factor f( k ‴), which is evaluated at ${\boldsymbol{k}}\prime\prime\prime =(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 Equation (12) with ${\theta }_{1}^{\prime\prime\prime }$.

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

Zoom In Zoom Out Reset image size

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

Since stellar interiors typically have extreme parameters, such as Pr ≲ 10⁻⁶, 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 the predictions of the closure model with those of the QL model, over a wide range of parameters Pr ≈ 10⁻⁷ − 1 and r ≈ 10⁻⁵ − 1 in Figure 6. Noting that ${N}^{2}=2(S-2)\left[1+r({\Pr }^{-1}-1)\right]$, these scans span N² ≲ 19 × 10⁶ (in terms of Ω²); this ratio is typically around 1 million for the Sun (Christensen-Dalsgaard et al. 1996). In Figure 6, with S = 3 (in terms of Ω), the Richardson number is as large as ≈2 × 10⁶. More extreme parameters can be easily and quickly scanned with the analytic formula we have derived.

Zoom In Zoom Out Reset image size

Now we predict the 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\equiv {U}_{0}{({\partial }_{x}{U}_{0})}^{-1}$ is the scale height of the mean flow ${U}_{0}=-{ \mathcal S }x$. The timescale for modifying the flow is ${\tau }_{\mathrm{turb}}\sim {U}_{0}H/\langle {\tilde{u}}_{x}{\tilde{u}}_{y}\rangle \sim {{Sx}}^{2}{{\rm{\Omega }}}^{-1}{d}^{-2}/\langle {\tilde{u}}_{x}{\tilde{u}}_{y}{\rangle }_{\mathrm{dimensionless}}$ (see also Barker et al. 2019).

Though Pr ∼ 10⁻⁶, the typical values of r in the solar tachocline and red giant stars are r ∼ 10⁻³ − 1 (varying with radius). Then, using Figure 6, where $\langle {\tilde{u}}_{x}{\tilde{u}}_{y}{\rangle }_{\mathrm{dimensionless}}$ is on average 0.5, we predict τ_turbΩ ∼ 2S(x/d)². This turbulent transport timescale 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 also depends sensitively 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 timescale for the GSF instability can be shorter than ${ \mathcal O }(10)$ Myr, the incorporation of our transport model for the GSF turbulence in stellar evolution codes is warranted (particularly if extended to nonequatorial and 3D 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.

We emphasize that the equations describing the 2.5D thermohaline and GSF instabilities at the equator are identical, when the compositional diffusivity is equal to the kinematic viscosity—a case often realized in stars. We first write Equations (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 timescale τ = d²/κ, and scaled temperature T in units of N² d (we label their x-coordinate with our z-coordinate and vice versa):

$$\begin{eqnarray}&&{\Pr }^{-1}{{Du}}_{x}^{{\rm{B}}}=-{\partial }_{x}{p}^{{\rm{B}}}+({T}^{{\rm{B}}}-{\mu }^{{\rm{B}}})+{{\rm{\nabla }}}^{2}{u}_{x}^{{\rm{B}}},\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&{\Pr }^{-1}{{Du}}_{z}^{{\rm{B}}}=-{\partial }_{z}{p}^{{\rm{B}}}+{{\rm{\nabla }}}^{2}{u}_{z}^{{\rm{B}}},\end{eqnarray} \tag{ 13b }$$

$$\begin{eqnarray}&&{{DT}}^{{\rm{B}}}=-{u}_{x}^{{\rm{B}}}+{{\rm{\nabla }}}^{2}{T}^{{\rm{B}}},\end{eqnarray} \tag{ 13c }$$

$$\begin{eqnarray}&&D{\mu }^{{\rm{B}}}=-\displaystyle \frac{{u}_{x}^{{\rm{B}}}}{{R}_{0}}+\Pr {{\rm{\nabla }}}^{2}{\mu }^{{\rm{B}}},\end{eqnarray} \tag{ 13d }$$

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

$$\begin{eqnarray}&&{u}_{x}^{{\rm{B}}}=\displaystyle \frac{{u}_{x}}{d/\tau },\end{eqnarray} \tag{ 14a }$$

$$\begin{eqnarray}&&{u}_{z}^{{\rm{B}}}=\displaystyle \frac{{u}_{z}}{d/\tau },\end{eqnarray} \tag{ 14b }$$

$$\begin{eqnarray}&&{T}^{{\rm{B}}}=\displaystyle \frac{\theta }{{N}^{2}d},\end{eqnarray} \tag{ 14c }$$

$$\begin{eqnarray}&&{\mu }^{{\rm{B}}}=-\displaystyle \frac{2\tau {\rm{\Omega }}}{\Pr }\displaystyle \frac{{u}_{y}}{d/\tau }.\end{eqnarray} \tag{ 14d }$$

Hence, the transport rates with two kinds of nondimensionalizations—one using τ and d as the characteristic timescale and length scale, and another using Ω and d as the relevant scales—are related in the manner

$$\begin{eqnarray}&&-\langle {\widetilde{u}}_{x}^{{\rm{B}}}{\widetilde{\mu }}^{{\rm{B}}}\rangle =\langle {\widetilde{u}}_{x}{\widetilde{u}}_{y}\rangle \times 2{\rm{\Omega }}{\Pr }^{1/2}{({N}^{2})}^{-3/2}{d}^{-2},\end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray}&&\langle {\widetilde{u}}_{x}^{{\rm{B}}}{\widetilde{T}}^{{\rm{B}}}\rangle =\langle {\widetilde{u}}_{x}\widetilde{\theta }\rangle \times {\Pr }^{1/2}{({N}^{2})}^{-3/2}{d}^{-2},\end{eqnarray} \tag{ 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 Pr, S, and N².

The expressions given in Equations (15a) and (15b) are plotted in Figure 5. Thus, our Figure 5 also represents a comparison of the 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, Pr) is realized only in the 2.5D equatorial case, which is where the analogy of the GSF instability with the thermohaline instability becomes exact.

The mode amplitude β_j evolution equation, given in Equation (3a), for a wavenumber k = (k_x , k_z ≠ 0), with j = 1, 2, or 3, is

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\beta }_{j} & = & {\gamma }_{j}{\beta }_{j}+\displaystyle \sum _{{\boldsymbol{k}}^{\prime} ,m,n}{C}_{{jmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} }\\ & & +\,\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },n}\left\{\left[{C}_{j{ \mathcal Y }n}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}+{C}_{{jn}{ \mathcal Y }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )}\right]{ \mathcal Y }^{\prime} {\beta }_{n}^{{\prime\prime} }\right.\\ & & +\,\left[{C}_{j{ \mathcal Z }n}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}+{C}_{{jn}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )}\right]{ \mathcal Z }^{\prime} {\beta }_{n}^{{\prime\prime} }\\ & & +\,\left.\left[{C}_{j{\rm{\Theta }}n}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}+{C}_{{jn}{\rm{\Theta }}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )}\right]{\rm{\Theta }}^{\prime} {\beta }_{n}^{{\prime\prime} }\right\},\end{array}\end{eqnarray} \tag{ C1 }$$

whereas, at k = (k_x , k_z = 0), one finds that the fluctuation amplitude evolves according to Equation (3b), which is expanded below:

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal Y }=-{\gamma }_{{ \mathcal Y }}{ \mathcal Y }+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m,n:{k}_{z}=0}{C}_{{ \mathcal Y }{mn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} },\end{eqnarray} \tag{ C2a }$$

$$\begin{eqnarray}&&{\partial }_{t}{ \mathcal Z }=-{\gamma }_{{ \mathcal Z }}{ \mathcal Z }+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m,n:{k}_{z}=0}{C}_{{ \mathcal Z }{mn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} },\end{eqnarray} \tag{ C2b }$$

$$\begin{eqnarray}&&{\partial }_{t}{\rm{\Theta }}=-{\gamma }_{{\rm{\Theta }}}{\rm{\Theta }}+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m,n:{k}_{z}=0}{C}_{{\rm{\Theta }}{mn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} },\end{eqnarray} \tag{ C2c }$$

where ${\gamma }_{{ \mathcal Y }}={\gamma }_{{ \mathcal Z }}=\nu {k}^{2}$ and γ_Θ = κ k² are the damping rates.

To derive an evolution equation for energy in the jth eigenmode at k = (k_x , k_z ≠ 0), we multiply Equation (C1) with ${\beta }_{j}^{* }$ and add a complex conjugate of the resulting equation to arrive at

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| {\beta }_{j}{| }^{2} & = & 2\mathrm{Re}\,{\gamma }_{j}| {\beta }_{j}{| }^{2}+\displaystyle \sum _{{\boldsymbol{k}}^{\prime} ,m,n}2\mathrm{Re}\,\left[{C}_{{jmn}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{m}^{{\prime} }{\beta }_{n}^{{\prime\prime} }{\beta }_{j}^{* }\rangle \right]\\ & & +\,\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}2\mathrm{Re}\,\left\{\left[{C}_{j{ \mathcal Y }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle { \mathcal Y }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle \right.\right.\\ & & +\,{\left.\left.{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle +{C}_{j{\rm{\Theta }}j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\rm{\Theta }}^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle \right]\right|}_{{k}_{z}^{{\prime} }=0}\\ & & +\,\left[{C}_{{jj}{ \mathcal Y }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{j}^{{\prime} }{ \mathcal Y }^{\prime\prime} {\beta }_{j}^{* }\rangle +{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{j}^{{\prime} }{ \mathcal Z }^{\prime\prime} {\beta }_{j}^{* }\rangle \right.\\ & & +\,\left.{\left.\left.{C}_{{jj}{\rm{\Theta }}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{j}^{{\prime} }{\rm{\Theta }}^{\prime\prime} {\beta }_{j}^{* }\rangle \right]\right|}_{{k}_{z}^{{\prime} }={k}_{z}}\right\}.\end{array}\end{eqnarray} \tag{ C3 }$$

Similar equations can be derived for the fluctuation energy at GSF-stable wavenumbers k = (k_x , k_z = 0) using Equations (C2a)–(C2c):

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| { \mathcal Y }{| }^{2} & = & -2\mathrm{Re}\,{\gamma }_{{ \mathcal Y }}| { \mathcal Y }{| }^{2}+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m}2\mathrm{Re}\,\\ & & \times \,{\left.\left[{C}_{{ \mathcal Y }{mm}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{m}^{{\prime} }{\beta }_{m}^{{\prime\prime} }{{ \mathcal Y }}^{* }\rangle \right]\right|}_{{k}_{z}=0},\end{array}\end{eqnarray} \tag{ C4a }$$

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| { \mathcal Z }{| }^{2} & = & -2\mathrm{Re}\,{\gamma }_{{ \mathcal Z }}| { \mathcal Z }{| }^{2}+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m}2\mathrm{Re}\,\\ & & \times \,{\left.\left[{C}_{{ \mathcal Z }{mm}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{m}^{{\prime} }{\beta }_{m}^{{\prime\prime} }{{ \mathcal Z }}^{* }\rangle \right]\right|}_{{k}_{z}=0},\end{array}\end{eqnarray} \tag{ C4b }$$

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| {\rm{\Theta }}{| }^{2} & = & -2\mathrm{Re}\,{\gamma }_{{\rm{\Theta }}}| {\rm{\Theta }}{| }^{2}+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m}2\mathrm{Re}\,\\ & & \times \,{\left.\left[{C}_{{\rm{\Theta }}{mm}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{m}^{{\prime} }{\beta }_{m}^{{\prime\prime} }{{\rm{\Theta }}}^{* }\rangle \right]\right|}_{{k}_{z}=0}.\end{array}\end{eqnarray} \tag{ C4c }$$

Since numerical simulations inform us that the triplets with a latitudinal flow at (k_x , 0), i.e., ${ \mathcal Z }$, dominate the nonlinear energy transfer, we may drop nonlinear terms on the right-hand side of Equation (C3) that do not involve ${ \mathcal Z }$. In the resulting equation, because ${ \mathcal Y }$ and Θ do not appear, Equations (C4a) and (C4c) can also be removed, which allows us to write the following set of equations:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| {\beta }_{j}{| }^{2} & = & 2\mathrm{Re}\,{\gamma }_{j}| {\beta }_{j}{| }^{2}+\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}2\mathrm{Re}\,\left\{{\left.\left[{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle \right]\right|}_{{k}_{z}^{{\prime} }=0}\right.\\ & & +\,\left.{\left.\left[{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{j}^{{\prime} }{ \mathcal Z }^{\prime\prime} {\beta }_{j}^{* }\rangle \right]\right|}_{{k}_{z}^{{\prime} }={k}_{z}}\right\},\end{array}\end{eqnarray} \tag{ C5a }$$

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| { \mathcal Z }{| }^{2} & = & -2\mathrm{Re}\,{\gamma }_{{ \mathcal Z }}| { \mathcal Z }{| }^{2}\\ & & +\,\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} },m}2\mathrm{Re}\,{\left.\left[{C}_{{ \mathcal Z }{mm}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}\langle {\beta }_{m}^{{\prime} }{\beta }_{m}^{{\prime\prime} }{{ \mathcal Z }}^{* }\rangle \right]\right|}_{{k}_{z}=0}.\end{array}\end{eqnarray} \tag{ C5b }$$

To obtain evolution equations for the terms on the right-hand side of Equations (C5a) and (C5b), we multiply Equation (3a) with two amplitudes. For example, to determine the evolution of $\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle$, we first multiply Equation (C1) for β_j with ${ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }$; second, we multiply another equation, similar to Equation (C1), but for ${\beta }_{j}^{{\prime\prime} }$ with ${ \mathcal Z }^{\prime} {\beta }_{j}^{* }$; and finally, we multiply Equation (C2b), for ${ \mathcal Z }^{\prime}$, with ${\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }$, and add them all together. Using this, we provide below an example evolution equation for triplet correlations:

$$\begin{eqnarray}&&\begin{array}{l}\left[{\partial }_{t}-\left(-{\gamma }_{{ \mathcal Z }}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }+{\gamma }_{j}^{* }\right)\right]\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{{\prime} }=0}\\ =\,\left\{\displaystyle \sum _{{\boldsymbol{k}}\prime\prime\prime ,m}\left[{C}_{{ \mathcal Z }{mm}}^{({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\beta }_{m}^{\prime\prime\prime }{\beta }_{m}({\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle \right]\right.\\ \ +\,\displaystyle \sum _{{k}_{x}^{\prime\prime\prime }}\left[{C}_{j{ \mathcal Y }j}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle { \mathcal Y }\prime\prime\prime {\beta }_{j}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\right.\\ \ +\,{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle { \mathcal Z }\prime\prime\prime {\beta }_{j}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\\ \ +\,{C}_{j{\rm{\Theta }}j}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\rm{\Theta }}\prime\prime\prime {\beta }_{j}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\\ \ +\,{C}_{{jj}{ \mathcal Y }}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\beta }_{j}^{\prime\prime\prime }{ \mathcal Y }({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\\ \ +\,{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\beta }_{j}^{\prime\prime\prime }{ \mathcal Z }({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\\ \ +\,\left.{C}_{{jj}{\rm{\Theta }}}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\beta }_{j}^{\prime\prime\prime }{\rm{\Theta }}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\right]\\ \ +\,\displaystyle \sum _{{k}_{x}^{\prime\prime\prime }}\left[{C}_{j{ \mathcal Y }j}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle { \mathcal Y }\prime\prime\prime * {\beta }_{j}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\right.\\ \ +\,{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle { \mathcal Z }\prime\prime\prime * {\beta }_{j}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\\ \ +\,{C}_{j{\rm{\Theta }}j}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle {\rm{\Theta }}\prime\prime\prime * {\beta }_{j}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\\ \ +\,{C}_{{jj}{ \mathcal Y }}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle {\beta }_{j}^{\prime\prime\prime * }{{ \mathcal Y }}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\\ \ +\,{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle {\beta }_{j}^{\prime\prime\prime * }{{ \mathcal Z }}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\\ \ +\,{\left.\left.\left.{C}_{{jj}{\rm{\Theta }}}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle {\beta }_{j}^{\prime\prime\prime * }{{\rm{\Theta }}}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\right]\right\}\right|}_{{k}_{z}^{{\prime} }=0}.\end{array}\end{eqnarray} \tag{ C6 }$$

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 Equation (C6) can be dropped in the face of the 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 Equation (C6), do not form terms that appear in the definition of energy. Guided by numerical simulations where nonlinear coupling to ${ \mathcal Y }$ and Θ are unimportant, only the energy(-like) terms with ${ \mathcal Z }$ will be kept henceforth. The resulting equation then reads

$$\begin{eqnarray}&&\begin{array}{l}\left[{\partial }_{t}-\left(-{\gamma }_{{ \mathcal Z }}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }+{\gamma }_{j}^{* }\right)\right]\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{{\prime} }=0}\\ =\,\left\{\displaystyle \sum _{{k}_{x}^{\prime\prime\prime }}\left[{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle { \mathcal Z }\prime\prime\prime {\beta }_{j}({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\right.\right.\\ \ +\,{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}}\prime\prime\prime )}\langle {\beta }_{j}^{\prime\prime\prime }{ \mathcal Z }({\boldsymbol{k}}-{\boldsymbol{k}}^{\prime} -{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\\ \ +\,{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle { \mathcal Z }\prime\prime\prime * {\beta }_{j}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }=0}\\ \ +\,{\left.\left.\left.{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}\prime\prime\prime )* }\langle {\beta }_{j}^{\prime\prime\prime * }{{ \mathcal Z }}^{* }({\boldsymbol{k}}-{\boldsymbol{k}}\prime\prime\prime ){ \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }\rangle {| }_{{k}_{z}^{\prime\prime\prime }={k}_{z}}\right]\right\}\right|}_{{k}_{z}^{{\prime} }=0}.\end{array}\end{eqnarray} \tag{ C7 }$$

The same procedure is then repeated to find evolution equations for the other two triplet correlations that appear in Equations (C5a) and (C5b).

Equation (C7) 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 γ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

$$\begin{eqnarray}&&\begin{array}{l}\langle { \mathcal Z }^{\prime} {\beta }_{j}^{{\prime\prime} }{\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{{\prime} }=0}\\ =\,-\left\{{\left(-{\gamma }_{{ \mathcal Z }}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }+{\gamma }_{j}^{* }\right)}^{-1}| { \mathcal Z }^{\prime} {| }^{2}\left[\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime\prime} ,-{\boldsymbol{k}}^{\prime} )}+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}\right)| {\beta }_{j}{| }^{2}\right.\right.\\ \ +{\left.\left.\left.\,\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )* }+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )* }\right)| {\beta }_{j}^{{\prime\prime} }{| }^{2}\right]\right\}\right|}_{{k}_{z}^{{\prime} }=0}.\end{array}\end{eqnarray} \tag{ C8 }$$

Similarly, the other two triplet correlations that appear in Equations (C5a) and (C5b) can also be solved to obtain

$$\begin{eqnarray}&&\begin{array}{l}\langle {\beta }_{j}^{{\prime} }{ \mathcal Z }^{\prime\prime} {\beta }_{j}^{* }\rangle {| }_{{k}_{z}^{{\prime} }={k}_{z}}\\ =\,-\left\{{\left({\gamma }_{j}^{{\prime} }-{\gamma }_{{ \mathcal Z }}^{\prime\prime} +{\gamma }_{j}^{* }\right)}^{-1}| { \mathcal Z }^{\prime\prime} {| }^{2}\left[\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime} ,-{\boldsymbol{k}}^{\prime\prime} )}+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}})}\right)| {\beta }_{j}{| }^{2}\right.\right.\\ \ +\ {\left.\left.\left.\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )* }+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )* }\right)| {\beta }_{j}^{{\prime} }{| }^{2}\right]\right\}\right|}_{{k}_{z}^{{\prime} }={k}_{z}},\end{array}\end{eqnarray} \tag{ C9 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}\langle {\beta }_{j}^{{\prime} }{\beta }_{j}^{{\prime\prime} }{{ \mathcal Z }}^{* }\rangle {| }_{{k}_{z}=0}\\ =\,-x\left\{{\left({\gamma }_{j}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }-{\gamma }_{{ \mathcal Z }}^{* }\right)}^{-1}| { \mathcal Z }{| }^{2}\left[\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}})}+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}^{\prime} ,-{\boldsymbol{k}}^{\prime\prime} )}\right)| {\beta }_{j}^{{\prime\prime} }{| }^{2}\right.\right.\\ \ +\ {\left.\left.\left.\left({C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}+{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}}^{\prime\prime} ,-{\boldsymbol{k}}^{\prime} )}\right)| {\beta }_{j}^{{\prime} }{| }^{2}\right]\right\}\right|}_{{k}_{z}=0}.\end{array}\end{eqnarray} \tag{ C10 }$$

The solutions of the triplet correlations from Equations(C8)–(C10) are now substituted into Equations (C5a) and (C5b), which results in the following set of closed equations:

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}| {\beta }_{j}{| }^{2}=2\mathrm{Re}\,{\gamma }_{j}| {\beta }_{j}{| }^{2}-\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}2\mathrm{Re}\,\\ \times \,{\left.\left[\displaystyle \frac{{C}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}| { \mathcal Z }^{\prime} {| }^{2}}{\left(-{\gamma }_{{ \mathcal Z }}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }+{\gamma }_{j}^{* }\right)}\left\{{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime\prime} ,-{\boldsymbol{k}}^{\prime} )}| {\beta }_{j}{| }^{2}+{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )* }| {\beta }_{j}^{{\prime\prime} }{| }^{2}\right\}\right]\right|}_{{k}_{z}^{{\prime} }=0}\\ -\,\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}2\mathrm{Re}\,\left[\displaystyle \frac{{C}_{{jj}{ \mathcal Z }}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}| { \mathcal Z }^{\prime\prime} {| }^{2}}{\left({\gamma }_{j}^{{\prime} }-{\gamma }_{{ \mathcal Z }}^{\prime\prime} +{\gamma }_{j}^{* }\right)}\right.\\ \times \,{\left.\left.\left\{{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime} ,-{\boldsymbol{k}}^{\prime\prime} )}| {\beta }_{j}{| }^{2}+{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime\prime} )* }| {\beta }_{j}^{{\prime} }{| }^{2}\right\}\right]\right|}_{{k}_{z}^{{\prime} }={k}_{z}},\end{array}\end{eqnarray} \tag{ C11 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}| { \mathcal Z }{| }^{2}=-2\mathrm{Re}\,{\gamma }_{{ \mathcal Z }}| { \mathcal Z }{| }^{2}-| { \mathcal Z }{| }^{2}\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}\displaystyle \sum _{j=1}^{3}2\mathrm{Re}\,\\ \ \times \,{\left.\left[\displaystyle \frac{{C}_{{ \mathcal Z }{jj}}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}}{\left({\gamma }_{j}^{{\prime} }+{\gamma }_{j}^{{\prime\prime} }-{\gamma }_{{ \mathcal Z }}^{* }\right)}\left\{{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}})}| {\beta }_{j}^{{\prime\prime} }{| }^{2}+{\overleftrightarrow{C}}_{j{ \mathcal Z }j}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}| {\beta }_{j}^{{\prime} }{| }^{2}\right\}\right]\right|}_{{k}_{z}=0},\end{array}\end{eqnarray} \tag{ C12 }$$

where ${\overleftrightarrow{C}}_{{lmn}}^{({\boldsymbol{p}},{\boldsymbol{q}})}={C}_{{lmn}}^{({\boldsymbol{p}},{\boldsymbol{q}})}+{C}_{{lnm}}^{({\boldsymbol{p}},{\boldsymbol{p}}-{\boldsymbol{q}})}$ is the symmetrized coupling coefficient.

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

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}| { \mathcal Z }{| }^{2}/2 & = & -{\gamma }_{{ \mathcal Z }}| { \mathcal Z }{| }^{2}+| { \mathcal Z }{| }^{2}\displaystyle \sum _{{{\boldsymbol{k}}}^{{\prime} }}| {\beta }_{1}^{{\prime} }{| }^{2}\,\mathrm{Re}\,\\ & & \times \,\left[\displaystyle \frac{-{\overleftrightarrow{C}}_{{ \mathcal Z }11}^{({\boldsymbol{k}},{\boldsymbol{k}}^{\prime} )}{\overleftrightarrow{C}}_{1{ \mathcal Z }1}^{({\boldsymbol{k}}^{\prime\prime} ,{\boldsymbol{k}})}}{{\gamma }_{1}^{{\prime} }+{\gamma }_{1}^{{\prime\prime} }-{\gamma }_{{ \mathcal Z }}^{* }}\right],\end{array}\end{eqnarray} \tag{ C13 }$$

where the term inside the wavenumber summation in Equation (C12) has been symmetrized.