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

Differentially rotating stars and planets transport angular momentum (AM) 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 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 AM 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.


INTRODUCTION
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 stelbtripathi@wisc.eduA.J. Barker@leeds.ac.uk lar 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-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 ux(x, z), uy(x, z), and uz(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., ux) grow the fastest in the linear phase (t=100), strong latitudinal jets uz are generated nonlinearly (t=6000).The color bar for t=100 is shared by ux, uy, and 3 uz; the color bar for t=6000 is shared by 3 ux, 3 uy, and uz.The turbulent transport of angular momentum, e.g., ⟨ ux uy⟩, is predicted in this paper using a jet-coupled turbulence closure.
Schubert-Fricke (GSF) instability1 (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. 2019Barker et al. , 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 selfconsistent 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.

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-Ω=Ωê 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 U 0 (x)= − Sxê y where x is the radial coordinate, y is the azimuthal coordinate, and S = −dΩ Shell (x)/d ln x is the local radial-shear rate, with Ω Shell (x) representing the "Shellular" rotation of the simplified star.A uniform gravity field with g = −gê x is directed radially inward (Fig. 1).A background radial temperature gradient ∇T 0 then stratifies the fluid density radially, with a thermal expansion coefficient α.In such a background state, any perturbations in velocity u and scaled temperature θ=αgT evolve (Barker et al. 2019) as where the variables p, ν, and κ are the fluid pressure (per unit density), the kinematic viscosity, and the thermal diffusivity, respectively.We also define the Prandtl number Pr=ν/κ.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 N , where N 2 êx =∇Θ 0 =αg∇T 0 .Henceforth, we non-dimensionalize all variables using the characteristic rotation time scale Ω −1 and length scale d, with d=(νκ/N 2 ) 1/4 , which is typically similar to the wavelengths of fastest-growing modes.Thus, N =N /Ω is the dimensionless buoyancy frequency and S=S/Ω the dimensionless shear rate (Rossby number).The GSF instability occurs in low-Pr fluids whenever r ∈ [0, 1], where r = Pr(1 + N 2 κ −2 ep )/(Pr − 1), with κ ep = 2(2 − S) representing the dimensionless epicyclic frequency (Barker et al. 2019).

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 L X=γ X, where X=[û x , ûy , ûz , θ] 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 three degrees of freedom at any given wavenumber-two components of velocity, and θ.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 GSFunstable [Re (γ) > 0, where 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 Im (γ) ̸ = 0, where Im denotes the imaginary part; Im (γ) corresponds to the frequency of inertial-gravity, or gravitoinertial, waves (IGWs), 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 3-component velocity and temperature fields.The dispersion relation then is a simple cubic polynomial in γ (Goldreich & Schubert 1967) as where γ ν = γ + νk 2 and γ κ = γ + κk 2 .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 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 L, we can decompose arbitrary incompressible fluctuation Xarb with k z ̸ =0 as Xarb = 3 j=1 β j Xj , where β j is the amplitude of the j th eigenvector Xj ; 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, Xarb =Eβ, where β is a (column-)vector of mode amplitudes and E is an eigenvector matrix, whose j th column is Xj .Thus, β=E −1 Xarb .
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 where the amplitudes of the eigenvectors are denoted by X , Y, Z, and Θ.We reserve the β-notation above for the amplitudes of eigenvectors with k z ̸ =0.

Initial value problem
We perform an ensemble of direct numerical simulations of Eqs.(1a) and (1b), by seeding a lowamplitude 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 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: , 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 xcomponent of the velocity ûx,m correspond to the m th eigenvector at k; and likewise for β n and ûy,n ; 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|β 1 | 2 Re ûx,1 û * y,1 .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 to-tal 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 =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.

Mode-amplitude evolution
To analyze nonlinear mode couplings, Eqs. ( 1a)-( 1b) are first spatially Fourier-transformed: , where X, X′ , and X′′ are state vectors at k, k ′ , and k ′′ , respectively, satisfying k = k ′ + k ′′ .Then, following Sec.2.1 we substitute X=Eβ, and likewise for X′ and X′′ .We multiply the obtained equation with E −1 and take the j th row of the resulting equation.This process yields an evolution equation for the j th eigenmode at k.Such an evolution equation for mode amplitude β j for k z ̸ =0 is different from that of the mode amplitude F ∈ {X , Y, Z, Θ} for k z =0, although both are coupled and nonlinear: where γ j is the complex-valued growth rate for the j th eigenmode with k z ̸ =0; the real-valued damping rate γ F is γ X when F is replaced with X in Eq. (3b); likewise for the replacement of F with Y, Z, and Θ; we note that measures the overlap of eigenmodes m with k ′ , n with k ′′ , and j with k.Such a mode coupling coefficient is found by applying E −1 to the (column) vector of N ( X′ m , X′′ n ), a process that incorporates all the nonlinearities of the system, thus making C (k,k ′ ) jmn ideal for a comprehensive instability-saturation analysis.) driven by an unstable mode at the linearly fastest-growing wavenumber k = (0, 0.63), and by an unstable mode at 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 ⟨ ux uy⟩ 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 kx=0, however, contributes negligibly (3%) to the transport (blue curve).Simulation parameters used are S=2.1,N 2 =10 and Pr=0.01.
On the right-hand side of Eq. (3a), the second term is the nonlinear coupling between eigenmodes with k ′ z ̸ =0 and k ′′ z ̸ =0 (hence the two βs), and the third term, with an F and a β, is the nonlinear coupling between eigenmodes with k ′ z =0 and k ′′ z ̸ =0 (see the last paragraph of Sec.2.1).Equations (3a) 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.

Mode-energy evolution
The energy evolution equation for each eigenmode can now be derived by multiplying Eq. (3a) by β * j and adding the complex conjugate of the resulting equation to arrive at where Q j =2Re γ j |β 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 j 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 ϵ ν in Fig. 3. From T jAA in Eq. ( 4), we separate out the nonlinear transfer T 1Z1 in a triad that involves a latitudinal flow (5) We now compare T 1Z1 with T 1AA in Fig. 4, for a GSFunstable 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 =2π/L x ; see Fig. 1, right-column snapshot at t=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 =0 are, however, linearly stable to the GSF instability, and, thus, must necessarily be excited nonlinearly by the interactions between the GSFunstable modes.This energy received is then viscously damped at k z =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 10 −4 0.00 0.02 0.04 0.06 0.08 γ Figure 3. Spectra of linear and nonlinear-saturation properties of the GSF instability.On the negative-kx domain, the colored square boxes (yellow-green-purple) display time-averaged energy extraction rates ⟨Q1⟩t by unstable modes from the mean gradients in a nonlinear simulation.On the positive-kx domain, colored square boxes (black-red-yellow) show the time-averaged viscous dissipation rates, which are pronounced at low kz.Over the entire (kx, kz)-plane, the non-square filled and line contours show the growth rates γ of unstable modes, with white dashed contour lines on the negative-kx domain and with bluish filled contours on the positive-kx domain.The fastest-growing mode resides at around k = (0, 0.63).The simulation parameters are S=2.1,N 2 =10 and Pr=0.01.
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.

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 (3a) 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. 2021Li et al. , 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 Z, with k z =0, dominates the nonlinear coupling, the complexity of the equation is significantly reduced (Terry et al. 2018).

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. (3a) has the structure: terms, e.g., βββ and ββX , have been dropped because the nonlinear energy transfer is almost entirely dominated by βZβ-the latitudinal-flow coupling (Fig. 4).Since βZβ also evolves, one can similarly obtain evolution equation for third-order correlator as ∂ t (βZβ) = ...βZβ + ...ββZZ.The closure solution then yields a relation βZβ = ...ββZZ.Although useful later, this relation does not predict the mode amplitude β, needed for the turbulent transport prediction.

Detailed Closure Model
To make quantitative predictions for transport, we take the EDQNM-closed evolution equation for the lat-itudinal flow energy [see Appendix B, Eq. (A20)], where, on the right-hand side, the second term contains a product of four amplitudes, but notably with |Z| 2 that also appears in the first term of the right-hand side.In Eq. ( 7), is the three-wave interaction time found from the EDQNM closure, and 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 . A more general expression for |β ′ 1 | 2 is found by using a standard Markovian assumption (Terry et al. 2021): |β ′ 1 | 2 is more weakly dependent on wavenumbers than the other factors in Eq. ( 7) arising from the coupling coefficients and τ 11Z .Such a consideration provides an expression for nonlinearly saturated squared-mode-amplitude with where we note that the coupling coefficients scale linearly with wavenumbers, and τ 11Z is the inverse of the sum of three growth rates of eigenmodes in a triad.The growth rates in τ 11Z 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, τ 11Z is approximated by using the linear growth rates (Terry et al. 2018;Terry et al. 2021).Using such, one identifies that the triplet interaction time τ 11Z is maximal when the triad involves a latitudinal flow (Z) and two GSF-unstable modes (j=1).Shorter triplet interaction times τ 12Z 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 τ 12Z via both the frequency and damping rate of the IGW.The largest interaction time τ 11Z dominates saturation.
The radial turbulent transport of angular momentum is measured by ⟨ũ is the unstable-mode amplitude at k ′′′ over which the summation is applied.Then, using |β ′′′ 1 | 2 from the above paragraph, The y-component û′′′ y,1 of velocity of the unstable eigenvector can be replaced with, e.g., its temperature perturbation θ′′′ 1 to predict the turbulent heat flux ⟨ũ x θ⟩.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.
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 =0, which Fig. 2 shows is inadequate.Hence, predictions of the quasilinear model that tend to reproduce the datavalidated closure model predictions are at best a fortuitous coincidence, occurring in a very limited parameter regime.

Impact of the new model in astrophysics
Since stellar interiors typically have extreme parameters such as Pr ≲ 10 −6 , current and anticipated nearfuture 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 Pr ≈ 10 −7 -1 and r ≈ 10 −5 -1 in Fig. 6.Noting that N 2 =2(S − 2) 1 + r(Pr −1 − 1) , these scans span N 2 ≲ 19 × 10 6 (in terms of Ω 2 ); this ratio is typically around 1 million for the Sun (Christensen-Dalsgaard et al. 1996).In Fig. 6, with S=3 (in terms of Ω), the Richardson number is as large as ≈ 2 × 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 ⟨ũ x ũy ⟩H −1 , where H≡U 0 (∂ x U 0 ) −1 is the scale 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 =10); the unfilled markers correspond to the cases where the squared Brunt-Väisälä frequency N 2 is varied (S=2.1).Both S-and N 2 -scan results collapse onto a single master curve, when ⟨ũx ũy⟩ 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, Pr).The GSF instability operates when r ∈ [0, 1) and Pr<1 (Pr=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.
Though Pr ∼ 10 −6 , the typical values of r in the solar tachocline and red giant stars are r ∼ 10 −3 -1 (varying with radius).
Then, using Fig. 6 where ⟨ũ x ũy ⟩ dimensionless is on average 0.5, we predict τ turb Ω ∼ 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 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 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.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 τ =d 2 /κ, and scaled temperature T in units of N 2 d [we label their x-coordinate with our z-coordinate and vice-versa]

Relation to the thermohaline instability
where the state vector [u B x , u B z , T B , µ B , p 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 /κ 2 ep may be recast as R 0 =1 + r(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 B x , u B z , T B , µ B ] are identical to Eqs. (1a) and (1b) for the GSF instability with the state vector [u x , u z , θ, u y ], when we note Hence, the transport rates with two kinds of nondimensionalizations-one using τ and d as the characteristic time scale and length scale, and another using Ω and d as the relevant scales-are related in the manner 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 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, 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.

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(Paxton et al. , 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.
script where the closure model is implemented.Other data used in this article will be shared on reasonable request to the corresponding authors.
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, stepby-step derivations.

C2. EIGENMODE-ENERGY EVOLUTION
To derive an evolution equation for energy in the j th eigenmode at k = (k x , k z ̸ =0), we multiply Eq. (A8) with β * j and add a complex conjugate of the resulting equation to arrive at Similar equations can be derived for fluctuation energy at GSF-stable wavenumbers k = (k x , k z =0) using Eqs.(A9a)-(A9c): Since numerical simulations inform us that the triplets with a latitudinal flow at (k x , 0), i.e., Z, dominates the nonlinear energy transfer, we may drop nonlinear terms on the right-hand side of Eq. (A10) that do not involve Z.In the resulting equation, because Y and Θ do not appear, Eqs. ( A11a) and (A11c) can also be removed, which allows us to write the following set of equations:

C3. TRIPLET CORRELATION EVOLUTION
To obtain evolution equations for the terms on the right-hand side of Eqs. ( A12a) and (A12b), we multiply Eq. (3a) with two amplitudes.For example, to determine the evolution of ⟨Z ′ β ′′ j β * j ⟩, first, we multiply Eq. (A8) for β j with Z ′ β ′′ j ; second, we multiply another equation, similar to Eq. ( A8), but for β ′′ j with Z ′ β * j ; and, finally, we multiply Eq. ( A9b), for Z ′ , with β ′′ j β * j , and add them all together.Using this, we provide below an example evolution equation for triplet correlations: . (A13) Since the numerical simulations show dominant coupling between the latitudinal flow 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 Y and 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 Y and Θ are unimportant, only the where the term inside the wavenumber-summation in Eq. (A19) has been symmetrized.

Figure 2 .
Figure 2. (a) Comparison of momentum transport (Reynolds Stress=2|β1| 2 Re ûx,1û * y,1 ) driven by an unstable mode at the linearly fastest-growing wavenumber k = (0, 0.63), and by an unstable mode at 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 ⟨ ux uy⟩ 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 kx=0, however, contributes negligibly (3%) to the transport (blue curve).Simulation parameters used are S=2.1,N 2 =10 and Pr=0.01.

∂Figure 4 .
Figure 4. Time evolution of the total nonlinear energy transfer T1AA to an unstable eigenmode at a wavenumber where the spectrum of ⟨ũx ũy⟩ peaks.T1Z1 is the energy transfer to the same unstable mode via interactions between the z-component of velocity (Z) with wavenumbers kz=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=2.1,N 2 =10 and Pr=0.01.

Figure 6 .
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 Pr ≪ 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 Pr < r and Pr > r found in (a).