1. Introduction
Turbulence in the strongly stratified regions of planetary oceans, atmospheres and the interiors of stars and gas giants provides an important source of vertical transport of chemicals and momentum, thereby playing a critical role in their long-term evolution (Zahn Reference Zahn1974, Reference Zahn1992; Fernando Reference Fernando1991; Pinsonneault Reference Pinsonneault1997; Maeder & Meynet Reference Maeder and Meynet2000; Ivey, Winters & Koseff Reference Ivey, Winters and Koseff2008; Ferrari & Wunsch Reference Ferrari and Wunsch2009; Aerts, Mathis & Rogers Reference Aerts, Mathis and Rogers2019; Garaud Reference Garaud2021). However, current understanding of the extremely low thermal-Prandtl-number regime of astrophysical turbulence remains disjoint from the order-unity Prandtl-number regime of geophysical turbulence. This is despite identical asymptotic limits for the Reynolds number (
$Re\gg 1$), Froude number (
$Fr\ll 1$) and buoyancy Reynolds number (
$Rb=Re\,Fr^2\gg 1$). The Prandtl number 
$Pr=\nu /\kappa$ measures the ratio of the microscopic viscosity 
$\nu$ to the thermal diffusivity 
$\kappa$, and is extremely small in stellar plasma, in particular. Photons rapidly diffuse heat compared to the much slower momentum diffusion by ion–ion collisions, leading to 
$\kappa \gg \nu$. For reference, stellar radiative zones have Prandtl numbers that can range from 
$Pr=O(10^{-9})$ to at most 
$Pr=O(10^{-5})$, in contrast to Earth's fluids, which range from 
$Pr\simeq 0.7$ in the atmosphere to 
$Pr\simeq 10$ in the ocean (Garaud Reference Garaud2021). A Prandtl number as high as 
$Pr\simeq 700$ can be reached in parts of the ocean dominated by salt stratification (Thorpe Reference Thorpe2005; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018; Gregg Reference Gregg2021) – where the salt diffusivity replaces thermal diffusivity. Compared to the 
$Pr=1$ case, a small 
$Pr$ can have significant effects on large-scale hydrodynamic instabilities and the resulting turbulence (Garaud Reference Garaud2021), while a large 
$Pr$ influences scales comparable to and smaller than the viscous scale and can have significant effects on the buoyancy flux, mixing efficiency and properties of shear-induced turbulence (Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015; Legaspi & Waite Reference Legaspi and Waite2020; Okino & Hanazaki Reference Okino and Hanazaki2020).
Many important features of the unity Prandtl-number regime are becoming better understood from a combination of numerical experiments (Waite & Bartello Reference Waite and Bartello2004; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Riley & Lindborg Reference Riley and Lindborg2010; Maffioli & Davidson Reference Maffioli and Davidson2016; Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2017; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019), observational data (Lindborg & Brethouwer Reference Lindborg and Brethouwer2007; Riley & Lindborg Reference Riley and Lindborg2008; Falder, White & Caulfield Reference Falder, White and Caulfield2016; Lefauve & Linden Reference Lefauve and Linden2022) and theoretical developments (Billant & Chomaz Reference Billant and Chomaz2001; Lindborg Reference Lindborg2006; Chini et al. Reference Chini, Michel, Julien, Rocha and Colm-cille2022). Strongly stratified turbulence with 
$Pr=O(1)$ forced at some horizontal length scale leads to emergent vertical scales set by the stratification and exhibits two distinct inertial ranges that transfer energy from the large forcing scales to the small viscous and thermal scales where it is dissipated. The transition length scale between the two ranges is known as the Ozmidov scale (Ozmidov Reference Ozmidov1992). The turbulence is highly anisotropic from the outer forcing scale down to the Ozmidov scale, with the down-scale energy transfer likely containing contributions from both a local energy cascade and non-local energy transfer mechanisms, e.g. shear instabilities, wave–wave interactions (Waite Reference Waite2011; Augier, Billant & Chomaz Reference Augier, Billant and Chomaz2015; Maffioli & Davidson Reference Maffioli and Davidson2016; Khani Reference Khani2018). Below the Ozmidov scale, the turbulence is isotropic down to the diffusive scales. The properties of the energy cascade, relevant dimensionless parameters and scaling relations for the emergent vertical scales are now fairly well understood, although many further questions remain, such as on the efficiency of mixing (Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018; Monismith, Koseff & White Reference Monismith, Koseff and White2018; Legg Reference Legg2021) and the origin of self-organized criticality (Smyth & Moum Reference Smyth and Moum2013; Salehipour, Peltier & Caulfield Reference Salehipour, Peltier and Caulfield2018; Smyth, Nash & Moum Reference Smyth, Nash and Moum2019; Chini et al. Reference Chini, Michel, Julien, Rocha and Colm-cille2022; Lefauve & Linden Reference Lefauve and Linden2022).
The main aim of this study is to extend the theoretical arguments used in the 
$Pr=O(1)$ regime to make analogous predictions for the emergent vertical scales and energy cascades in the asymptotically low 
$Pr$ regime. If only the thermal diffusivity is increased while keeping all other parameters constant (thereby decreasing 
$Pr$), one should expect a smooth transition to occur between two asymptotic regimes when thermal diffusion shifts from being important only on the smallest viscous scales to playing an important role on the mesoscales, i.e. on scales comparable to or larger than the Ozmidov scale.
To understand this transition, we use the critical balance framework proposed by Nazarenko & Schekochihin (Reference Nazarenko and Schekochihin2011) for anisotropic wave systems. Critical balance argues that linear wave 
$\omega ^{-1}$ and nonlinear interaction 
$\tau _{{NL}}$ time scales are comparable, 
$\omega ^{-1} \sim \tau _{{NL}}$, on a scale-by-scale basis throughout a local energy cascade, giving a prediction for the anisotropy of the turbulence as a function of scale. Originally applied in mean-field magnetohydrodynamic (MHD) turbulence (Goldreich & Sridhar Reference Goldreich and Sridhar1995) (see also Schekochihin (Reference Schekochihin2022) and references therein), critical balance successfully predicts power laws of the anisotropic energy spectra and associated transition scales in rotating turbulence and unity Prandtl-number strongly stratified turbulence (Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011).
We propose that critical balance should naturally extend to low thermal-Prandtl-number strongly stratified turbulence with a modification in its physical interpretation. As 
$Pr$ is decreased (by increasing 
$\kappa$ while keeping 
$\nu$ small and fixed), the full internal gravity wave (IGW) dispersion 
$\omega$ smoothly transitions from the asymptotic dispersion for adiabatic, inviscid, propagating IGWs when 
$Pr=O(1)$ (i.e. 
$\omega \sim \omega _{{N}}$) to an asymptotic dispersion for overdamped, inviscid, IGWs modified by the interaction of buoyancy and fast thermal diffusion when 
$Pr\ll 1$ (i.e. 
$\omega \sim \omega _{{lPe}}$, where ‘lPe’ refers to the ‘low turbulent Péclet-number’ limit, as discussed in more detail in §§ 2 and 4). As a result, scaling laws for the emergent vertical scales and two cascades predicted by critical balance will likewise smoothly change as 
$Pr$ is decreased. We note that, because critical balance assumes a local energy cascade, the relative role of non-local energy transfer mechanisms in the low 
$Pr$ limit remains to be understood.
The paper layout is as follows. The Boussinesq equations used to model stably stratified turbulence are defined in § 2, followed by an argument for when a transition of turbulence regimes should occur. Critical balance arguments in the unity 
$Pr$ regime are reviewed in § 3, which can then be compared with the new critical balance arguments in the low 
$Pr$ regime given in § 4. Astrophysical applications of the new scaling laws are discussed in § 5.
2. The Boussinesq approximation and the 
$Pr\,Rb=1$ transition
2.1. The Boussinesq approximation
Stably stratified regions of stars and gas giants below their surfaces typically sustain turbulent motions with vertical length scales that are small compared with the local scale height and velocity fluctuations that are small compared to the local sound speed. (The local scale height is the local e-folding scale of a background thermodynamic variable. For example, the pressure scale height in a stellar interior is 
$H_P\equiv (\partial _r \ln P)^{-1}$.) In this limit, the Spiegel–Veronis–Boussinesq equations are a rigorous approximation for the fluctuations of the velocity field and thermodynamic variables of a compressible fluid on top of a stably stratified background (Spiegel & Veronis Reference Spiegel and Veronis1960). The approximation effectively filters out the high frequencies of sound waves compared to the lower frequencies of IGWs and fluid motions of interest.
The governing equations for the velocity field 
$\boldsymbol {u}'$ and the buoyancy variable 
$\theta '=\alpha g T'=-\rho 'g/\rho _m$ are
\begin{gather} \partial_t \boldsymbol{u}'+\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}' ={-}\frac{1}{\rho_m}\boldsymbol{\nabla} p'+\theta'\hat{z}+\nu \nabla^2 \boldsymbol{u}', \end{gather}
\begin{gather}\partial_t \theta'+\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla} \theta' ={-}N^2u_z'+\kappa \nabla^2\theta', \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'=0 . \end{gather}
Here 
$\alpha$ is the coefficient of thermal expansion, 
$g$ is the local gravitational constant, 
$T'$ is the temperature perturbation, 
$\rho '$ is the density perturbation, 
$\rho _m$ is the mean density of the region, 
$p'$ is the pressure perturbation and 
$N$ is the local Brunt–Väisälä frequency. Primed variables here denote dimensional quantities and unprimed variables will later denote dimensionless quantities. Note that the Brunt–Väisälä frequency captures all the relevant local thermodynamics of the medium in the Boussinesq approximation of any fluid (Bois Reference Bois1991). As a result, the above equations are formally equivalent to those used in geophysical fluid studies of the Earth's oceans, but with a different definition of 
$N$. For example, in the case of an ideal gas one has 
$N^2=\alpha g (\partial _z \bar {T}-\partial _z T_{ad})$, while in the case of a liquid, one has 
$N^2=-g\partial _z \bar {\rho }/\rho _m$, where 
$\partial _z\bar {\rho }$, 
$\partial _z\bar {T}$ and 
$\partial _z T_{ad}$ are the background density, background temperature and adiabatic temperature gradients, respectively.
2.2. Geophysical regime
A stably stratified fluid forced on a horizontal length 
$L$ and velocity scale 
$U$ leads to turbulence with an emergent outer vertical length 
$l_z$ and velocity scale 
$u_z$ set by the physical parameters of the fluid 
$\{N, \nu, \kappa \}$. From dimensional analysis, only three dimensionless parameters characterize the fluid: the Reynolds number 
$Re=UL/\nu$, the Froude number 
$Fr=U/NL$ and the Prandtl number 
$Pr=\nu /\kappa$. Understanding the scaling of the emergent outer vertical scales as well as the structure of the subsequent (anisotropic) energy cascade are important theoretical and experimental goals.
In the geophysical fluid regime where 
$Pr=O(1)$, evidence from theoretical arguments, simulations and experimental data strongly suggests that the outer vertical length and vertical scales are directly set by the Froude number when the viscosity is sufficiently small: 
$l_z\sim Fr\,L$ and 
$u_z\sim Fr\,U$, where 
$l_z$ is often called the buoyancy length scale. Strong stratification (
$Fr\ll 1$) thus leads to highly anisotropic structures of the large-scale eddies characterized by long horizontal scales and short vertical scales as well as significantly more energy in the horizontal compared to the vertical velocity components. The injected energy undergoes an anisotropic forward energy cascade at large length scales until the Ozmidov scale, an intermediate scale given by 
$l_O=Fr^{3/2}L$ (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). For scales smaller than 
$l_O$, the effects of buoyancy are negligible on the fast turnover times of the eddies, and an isotropic Kolmogorov cascade operates down to the dissipation scales, which are the viscous (
$l_\nu =Re^{-3/4}L$) and thermal (
$l_\kappa \simeq l_\nu$) scales. The range of the isotropic cascade is set by the buoyancy Reynolds number 
$Rb=Re\,Fr^2$ (i.e. 
$l_O/l_\nu =Rb^{3/4}$) and needs to be sufficiently large in order to support the two cascades characteristic of strongly stratified turbulence (Bartello & Tobias Reference Bartello and Tobias2013).
The horizontal and vertical scales from above in the limits 
$Fr\ll 1$ and 
$Rb\gg 1$ suggest a consistent rescaling of the Boussinesq equations using the following dimensionalization (identical to that of Billant & Chomaz (Reference Billant and Chomaz2001)):
\begin{gather} \boldsymbol{u}_h'=U \boldsymbol{u}_h, \quad u_z'=Fr\,Uu_z, \quad \theta'=\frac{1}{Fr}\frac{U^2}{L}\theta, \quad p'=\rho_m U^2p, \end{gather}
\begin{gather}x'=Lx, \quad y'=L y, \quad z'=Fr\,L z, \quad t'=\frac{L}{U} t, \end{gather}
where the scaling of 
$\theta '$ is determined by a balance of 
$\boldsymbol {u}'\boldsymbol {\cdot }\boldsymbol {\nabla }\theta '\sim N^2u_z'$, since thermal diffusivity is considered to be sufficiently small. Substituting the above, the Boussinesq equations become
\begin{gather} \partial_t \boldsymbol{u}_h+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_h ={-}\boldsymbol{\nabla}_h p+ \left[\frac{1}{Re}{\nabla}_h^2 +\frac{1}{Rb}{\nabla}_z^2\right]\boldsymbol{u}_h, \end{gather}
\begin{gather}Fr^2[\partial_t u_z+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}u_z] ={-}\boldsymbol{\nabla}_zp+\theta+Fr^2\left[\frac{1}{Re}{\nabla}_h^2 +\frac{1}{Rb}{\nabla}_z^2\right]u_z, \end{gather}
\begin{gather}\partial_{t}\theta+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta ={-}u_z+\left[\frac{1}{Pr\,Re}{\nabla}_h^2 +\frac{1}{Pr\,Rb}{\nabla}_z^2\right]\theta, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
where 
$\boldsymbol {u}_h$ and 
$\boldsymbol {\nabla }_h$ denote the horizontal components of the velocity and gradient. Note that 
$Rb$ and 
$Pr\,Rb$ act as effective Reynolds numbers in the vertical part of the momentum and thermal diffusion terms, respectively. Examination of the dominant balance to lowest order is helpful:
\begin{gather} \partial_t \boldsymbol{u}_h+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_h ={-}\boldsymbol{\nabla}_h p, \end{gather}
\begin{gather}0={-}\boldsymbol{\nabla}_z p+\theta, \end{gather}
\begin{gather}\partial_t \theta+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta ={-}u_z, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0. \end{gather}Advection in the horizontal momentum equation (2.4a) is balanced by horizontal pressure gradients, while the dominant balance in the vertical momentum equation (2.4b) is instead between the vertical pressure gradient and buoyancy fluctuations. These balances will change if the viscosity is increased and the vertical gradients of the momentum diffusion term become important. Further, the buoyancy equation (2.4c) is a balance between temperature advection and displacement of the fluid against the background stratification gradient. It is this latter balance that will change if thermal diffusion is increased and the vertical gradients of the thermal diffusion term become important, as we show in the next section.
2.3. Transitions from the geophysical regime
We now aim to understand when transitions occur from the regime of geophysical fluid turbulence where 
$Fr\ll 1$, 
$Rb\gg 1$ and 
$Pr=O(1)$. First, let us consider the better understood transition to the viscosity-affected stratified flow regime as 
$Re$ is decreased with fixed 
$Fr$ and 
$Pr=O(1)$ (Godoy-Diana, Chomaz & Billant Reference Godoy-Diana, Chomaz and Billant2004). The case of decreasing 
$Pr$ turns out to behave in an analogous manner. From the perspective of the turbulent cascades, as viscosity is increased, the viscous scale will grow until it is comparable to the Ozmidov scale, 
$l_\nu \sim l_O$, at which point 
$Rb\sim 1$ and the isotropic cascade at small scales disappears (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). This appears as a shift of the dominant balance in the horizontal momentum equation (2.3a) between advection and the vertical gradient of the momentum diffusivity,
\begin{equation} \frac{\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_h'}{\nu{\nabla}_z^2\boldsymbol{u}_h'}\sim \frac{u_zl_z}{\nu}\sim Re\,Fr^2=Rb, \end{equation}
where vertical and horizontal advection are comparable due to the incompressibility constraint (i.e. 
$u_z/l_z\sim U/L$) and 
$l_z/L\sim u_z/U\sim Fr$ is used to estimate the vertical scales near 
$Rb\sim 1$. Heuristically, 
$Rb$ is the ratio of the eddy turnover rate to the viscous diffusion rate at the outer vertical scales, i.e. 
$Rb\sim (u_z/l_z)/(\nu /l_z^2)$. For 
$Rb<1$, a change of dominant balance (
$\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}_h'\sim \nu {\nabla }_z^2\boldsymbol {u}_h'$) leads to an alternative scaling 
$l_z/L=Fr/Rb^{1/2}$ (usually written as 
$l_z/L\sim Re^{-1/2}$) where viscosity dominates the coupling of adjacent vertical layers. The same shift occurs in the buoyancy equation (
$N^2 u_z' \sim \kappa {\nabla }_z^2\theta '$), and the new vertical velocity scale becomes 
$u_z/U\sim Fr\, Rb^{1/2}$ (derived with 
$\theta '\sim \boldsymbol {\nabla }_z p'/\rho _m$ from the unchanged dominant balance in the vertical momentum equation). This regime is often reached by simulations because computational constraints limit how small the viscosity can be set. Note that the vertical scales smoothly transition from 
$l_z/L\sim u_z/L\sim Fr$ at 
$Rb=1$ as 
$Rb$ is decreased.
Returning to the physically interesting limit 
$Fr\ll 1$ and 
$Rb\gg 1$ of strongly stratified turbulence, we now consider the effect of decreasing 
$Pr$ at fixed 
$Re$ and 
$Fr$, equivalent to increasing the thermal diffusivity while keeping the viscosity fixed. From the perspective of the turbulent cascades, as thermal diffusivity is increased, the thermal scale 
$l_\kappa \sim (Pr\,Re)^{-3/4}L$ will grow until it is comparable to the Ozmidov scale 
$l_\kappa \sim l_O$, at which point 
$Pr\,Rb\sim 1$ and thermal diffusion becomes important on mesoscales that are influenced by buoyancy forces (Lignières Reference Lignières2019). This appears as a shift of the dominant balance only in the buoyancy equation (2.3c) between advection and the vertical gradient of thermal diffusion:
\begin{equation} \frac{\boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} \theta'}{\kappa{\nabla}_z^2\theta'}\sim \frac{u_zl_z}{\kappa}\equiv Pe_t. \end{equation}
Here 
$Pe_t$ is the turbulent Péclet number, which can be interpreted heuristically in a similar way to 
$Rb$ as the ratio of the eddy turnover rate to the thermal diffusion rate at the outer vertical scales, i.e. 
$Pe_t=(u_z/l_z)/(\kappa /l_z^2)$. If we use the scalings 
$l_z/L\sim u_z/U\sim Fr$ from the 
$Pr=O(1)$ regime to estimate 
$Pe_t$ in the geophysical regime, we see that 
$Pe_t\sim Pr\,Rb$ and so the transition from 
$Pe_t>1$ to 
$Pe_t<1$ occurs around 
$Pr\,Rb\sim 1$, exactly like the 
$l_\kappa \sim l_O$ transition discussed above. Thus, thermal diffusion will cause a transition in turbulent regimes if 
$Pr\,Rb<1$. The scalings for 
$u_z$ and 
$l_z$ from the 
$Pr\,Rb>1$ regime will then change (the scaling for 
$Pe_t$ will correspondingly change as well).
The emergent turbulent Péclet number is a more important parameter than the standard Péclet number 
$Pe=Pr\,Re$ (Zahn Reference Zahn1992; Lignières Reference Lignières2019; Cope, Garaud & Caulfield Reference Cope, Garaud and Caulfield2020), which measures the ratio of advection to the horizontal gradient of thermal diffusion: 
$(\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } \theta ')/(\kappa {\nabla }_h^2\theta ')\sim UL/\kappa$. This is because 
$Pe_t\ll Pe$ – thermal diffusion will always be more important in the vertical than the horizontal direction. In astrophysical systems, the extremely large Reynolds numbers often keep 
$Pe\gg 1$ despite small Prandtl numbers. As a canonical example, turbulence from horizontal shear instabilities in the solar tachocline approximately sustain 
$Re=O(10^{14})$ and 
$Pe=O(10^{8})$ using 
$Pr=O(10^{-6})$ (Garaud Reference Garaud2020). Thus horizontal thermal diffusion is likely to be less important on outer scales in other stars as well. On the other hand, we find 
$Pr\,Rb=O(10^{1})$ (and hence 
$Pe_t=O(10^1)$) by using 
$Fr\sim 3\times 10^{-4}$ and 
$Rb\sim O(10^7)$. The near unity value of 
$Pe_t$ shows that vertical thermal diffusion can easily become relevant in the astrophysical case (Garaud Reference Garaud2021), in particular, in stars with much lower 
$Pr$ (down to 
$Pr=O(10^{-9})$ in some stars), stronger background stratification (higher 
$N$) or weaker driving (lower 
$U$ or larger 
$L$). Interestingly, these values of 
$Re$, 
$Fr$ and 
$Rb$ are similar to those found in Earth's atmosphere (Lilly Reference Lilly1983; Waite, von Larcher & Williams Reference Waite, von Larcher and Williams2014).
An important asymptotic model known as the low Péclet-number approximation is often used to study the 
$Pe_t<1$ regime (Lignières Reference Lignières1999). It can be derived from the shift in dominant balance in the buoyancy equation. If 
$Pe_t\ll 1$, then the 
$-N^2u_z'$ term balances 
$\kappa \nabla ^2\theta '$ instead of 
$\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } \theta '$. In other words, to lowest order, buoyancy fluctuations are generated by vertical advection of the mean background profile while advection of the buoyancy fluctuations is unimportant because of rapid thermal diffusion. As a result, the buoyancy fluctuations can be solved for directly, 
$\theta '=({N^2}/{\kappa }) \nabla ^{-2}u_z'$. This is substituted back into (2.1a) to get a closed momentum equation:
\begin{equation} \partial_t\boldsymbol{u}'+\boldsymbol{u}'\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u}' ={-}\frac{1}{\rho_m}\boldsymbol{\nabla} p' +\frac{N^2}{\kappa}\nabla^{{-}2}u_z'\hat{z}+\nu \nabla^2\boldsymbol{u}', \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}'=0. \end{equation} There are now only two dimensional parameters, 
$N^2/\kappa$ and 
$\nu$, which, alongside the imposed 
$U$ and 
$L$, imply that the scaling in this limit can only be a function of the dimensionless parameters 
$Pe\,Fr^{-2}=Pr\,Rb/Fr^{4}$ and 
$Re$ (Lignières Reference Lignières2019; Cope et al. Reference Cope, Garaud and Caulfield2020).
The discussion so far has argued that the geophysical regime with 
$Fr\ll 1$, 
$Rb\gg 1$ and 
$Pr\,Rb>1$ can transition either to a flow dominated by viscous effects when 
$Rb<1$ or to stratified turbulence modified by thermal diffusion when 
$Pr\,Rb<1$. These transitions correspond to the effective Reynolds numbers of the vertical part of the momentum and thermal diffusion terms becoming smaller than unity, respectively. Our aim now is to derive a formal scaling for the astrophysically motivated 
$Pr\,Rb<1$ regime with the help of the critical balance hypothesis, but first we review critical balance arguments for the geophysical 
$Pr\,Rb>1$ regime.
3. Critical balance and scaling for $Pr\,Rb>1$
We derive here the standard scaling relations for the limit of 
$Fr\ll 1$, 
$Rb\gg 1$ and 
$Pr\,Rb>1$ using critical balance arguments. Consider a stably stratified fluid where energy is injected with power 
$P$ at a low wavenumber 
$k_f=2{\rm \pi} /L$ and sustains turbulence with a kinetic energy dissipation rate 
$\epsilon \sim P\sim U^3/L$, where 
$U$ and 
$L$ are the outer horizontal velocity and length scales. Viewing anisotropic structures in the turbulence as a function of horizontal scale 
$k_\perp ^{-1}$, the goal is to find the characteristic vertical scale 
$k_\parallel ^{-1}$ associated with each 
$k_\perp$. Here 
$k_\parallel$ and 
$k_\perp$ are wavenumbers parallel and perpendicular to the direction of gravity, respectively. Such a structure will have an associated linear time scale 
$\omega ^{-1}$ related to the wave dispersion and a nonlinear time scale 
$\tau _{NL}$ related to the self-straining time scale. We discuss estimation of both time scales for the 
$Pr\,Rb>1$ case before applying critical balance arguments that connect 
$k_\perp$ and 
$k_\parallel$.
The Boussinesq system supports linear motions with a dispersion relation given by 
$\omega (k_\parallel,k_\perp )$. In the limit of arbitrarily small viscosity and thermal diffusivity, wave damping is negligible (
$\nu k^2,\kappa k^2\ll \omega$) and the linear motions are propagating waves closely approximated by adiabatic, inviscid IGWs with dispersion 
$\omega \approx \omega _{{N}}= N k_\perp /k$. Because large-scale vertical motion is strongly restricted, the vertical scales are much finer than the horizontal scales, with an anisotropy at large scales quantified by 
$k_\perp /k_\parallel \ll 1$. The linear wave frequency is then approximately
\begin{equation} \omega_{{N}}\approx \frac{N k_\perp }{k_\parallel}. \end{equation} On the other hand, nonlinear interactions break up eddies and transfer energy from larger to smaller scales, setting up a cascade from the forcing to the diffusive scales where the energy is dissipated. Nonlinear interactions occur via the advection term 
$\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}$. Incompressibility, 
$\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}=0$, requires 
$\boldsymbol {\nabla }_\parallel u_\parallel \sim \boldsymbol {\nabla }_\perp u_\perp$ and allows the estimate 
$\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}\simeq \boldsymbol {u}_\perp \boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}_\perp \simeq \boldsymbol {u}_\parallel \boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {u}_\perp$, where 
$u_\parallel$ and 
$u_\perp$ are the scale-dependent vertical and horizontal velocities, respectively. The nonlinear interaction time scale using the perpendicular nonlinearity is given by
\begin{equation} \tau_{{NL}}^{{-}1}\sim k_\perp u_\perp(k_\perp). \end{equation} The scaling of 
$\tau _{NL}$ with 
$k_\perp$ can be found if a separate relation can connect 
$u_\perp$ and 
$\tau _{NL}$. This comes from assuming that a local (in scale) cascade brings energy down from larger to smaller scales with a ‘cascade time’ 
$\tau _{{cas}}(k_\perp )$ that determines the energy spectrum 
$E(k_\perp )$:
\begin{equation} k_\perp E(k_\perp)\sim u_\perp^2(k_\perp)\sim \epsilon\tau_{{cas}}(k_\perp). \end{equation}
A relation between 
$\tau _{cas}$ and 
$\tau _{NL}$ would provide the desired 
$\tau _{NL}(k_\perp )$. Note that the effects of non-local energy transfer mechanisms as well as energy irreversibly lost to the buoyancy flux are not included in this cascade, but would be needed for a more detailed theory. We expect that the presented arguments should be reasonable as long as the local kinetic energy cascade comprises an 
$O(1)$ fraction of the total energy cascade.
In homogeneous isotropic turbulence, the only dimensional option is to set 
$\tau _{{cas}}\sim \tau _{{NL}}$, which results in the standard Kolmogorov scalings 
$u_\perp (k_\perp )\sim (\epsilon /k_\perp )^{1/3}$, 
$\tau _{NL}\sim (L/U)(k_\perp L)^{-2/3}$ and 
$E(k_\perp )\sim \epsilon ^{2/3}k_\perp ^{-5/3}$. In anisotropic turbulence, the two additional dimensionless parameters 
$k_\parallel /k_\perp$ and 
$\omega \tau _{{NL}}$ no longer constrain the system, and the energy spectrum and cascade time can be an arbitrary functions of these dimensionless groups. A further physically motivated constraint is needed.
Nazarenko & Schekochihin (Reference Nazarenko and Schekochihin2011) propose critical balance as a universal scaling conjecture for strong turbulence in anisotropic wave systems. Critical balance states that the linear propagation 
$\omega ^{-1}$ and the nonlinear interaction time scales 
$\tau _{{NL}}$ are approximately equal, 
$\omega \tau _{{NL}}\sim 1$, on a scale-by-scale basis at all scales where the source of anisotropy is important. Physically, in essence, this is a causality argument in a system where the perpendicular nonlinearity dominates (e.g. rotating turbulence, MHD with a mean field): a fluctuation with some 
$k_\perp$ cannot maintain an extent longer than 
$k_\parallel ^{-1}$ set by requiring the linear propagation time in the parallel direction to be comparable to the nonlinear breakup time in the horizontal direction. However, the causality argument is more complicated in stably stratified turbulence (Nazarenko & Schekochihin Reference Nazarenko and Schekochihin2011) because the nonlinearity has equal strength in the parallel and perpendicular directions, 
$\tau _{NL}\sim k_\perp u_\perp \sim k_\parallel u_\parallel$, so one could equally argue a balance between linear propagation time in the perpendicular direction and nonlinear breakup time in the vertical direction. In either case, since the group velocity sets the propagation speed of information, the linear time scale across either the parallel or perpendicular extent of a fluctuation is 
$(k_\perp v_{g,\perp })^{-1}\sim (k_\parallel v_{g,\parallel })^{-1}\sim \omega _N^{-1}$. Several physical mechanisms are known to be consistent with critical balance, including zigzag (Billant & Chomaz Reference Billant and Chomaz2000a,Reference Billant and Chomazb) and shear instabilities (see § 4.1); however, a complete physical picture is still an area of investigation (see further discussions in Lindborg (Reference Lindborg2006)). Going forward, we assume the critical balance hypothesis and that IGWs effectively set the linear propagation time scale in the 
$Pr\,Rb>1$ limit.
Critical balance removes the ambiguity in determining the cascade time scale (i.e. 
$\tau _{{cas}}\sim \tau _{{NL}}$), which results in a Kolmogorov spectrum for the horizontal spectrum at all scales. The vertical spectrum can then be easily determined from 
$E(k_\parallel )k_\parallel \sim E(k_\perp )k_\perp$ once 
$k_\perp$ and 
$k_\parallel$ are related. Applying critical balance 
$\omega _N \tau _{NL}\sim 1$ and rearranging gives the following relation between 
$k_\perp$ and 
$k_\parallel$:
\begin{equation} k_\perp\sim \left(\frac{\epsilon}{N^3}\right)k^3_\parallel=l_O^2k_\parallel^3, \end{equation}
where 
$l_O=(\epsilon /N^3)^{1/2}=k_O^{-1}$ is the Ozmidov scale. The anisotropy 
$k_\perp /k_\parallel \sim (l_Ok_\parallel )^2$ decreases at smaller scales until the Ozmidov scale 
$k^{-1}\sim l_O$, where the turbulence returns to isotropy, 
$k_\perp \sim k_\parallel$ and 
$\tau _{{NL}}^{-1}\sim N$. The Ozmidov scale is thus the largest horizontal scale that can overturn before restoration by buoyancy forces becomes significant. The horizontal and vertical spectra for 
$k< k_O$ are then
\begin{equation} \frac{E(k_\perp)}{U^2L}\sim (k_\perp L)^{{-}5/3},\quad \frac{E(k_\parallel)}{U^2L}\sim Fr(k_\parallel l_z)^{{-}3}. \end{equation}
These energy spectra agree with the theoretical scaling predictions in the geophysical literature (Dewan Reference Dewan1997; Billant & Chomaz Reference Billant and Chomaz2001; Lindborg Reference Lindborg2006), where the parallel energy spectrum is often written in a dimensional form as 
$E(k_\parallel )\sim N^2k_\parallel ^{-3}$.
The turbulence no longer feels the large-scale stratification gradients below the Ozmidov scale since 
$\omega _{{N}} \tau _{{NL}}\ll 1$ for 
$k\gg k_O$. Consequently, an isotropic Kolmogorov cascade results at small scales because 
$\tau _{{NL}}$ becomes the only dimensionally available time scale (i.e. the slow buoyancy restoration time scale is irrelevant). The isotropic cascade extends from 
$k_O$ to the viscous wavenumber 
$k_\nu$, where 
$\tau _{{NL}}^{-1}\sim \nu k^2\rightarrow l_\nu \sim (\nu /\epsilon ^{1/3})^{3/4}$. As a result, the vertical spectrum has a break at 
$k\sim k_O$ from 
$k_\parallel ^{-3}$ to 
$k_\parallel ^{-5/3}$, while the horizontal spectrum remains 
$k_\perp ^{-5/3}$ throughout. The temperature thus plays an important role providing buoyancy for scales 
$k< k_O$ but is simply advected as a passive scalar for scales 
$k>k_O$ until it is dissipated at thermal diffusion scales.
For clarity, it is useful to summarize the energy cascade in terms of the dimensionless parameters 
$Fr$ and 
$Re$. Energy injected at large horizontal scales 
$L$ undergoes an anisotropic cascade until the Ozmidov scale, 
$l_O=Fr^{3/2}L$. Following the anisotropic cascade is an isotropic cascade down to the viscous scale, which can be written as 
$l_\nu =Re^{-3/4}L$. The size of the isotropic cascade is 
$l_O/l_\nu =Rb^{3/4}$. Here 
$Rb$ essentially plays the role of the effective Reynolds number for the isotropic cascade with outer length scale 
$l_O$ and velocity scale 
$u_{\parallel }(k_O)=(\epsilon l_O)^{1/3}=Fr^{1/2}U$ (i.e. 
$Rb=u_{\parallel }(k_O)l_O/\nu$). A schematic of the energy cascade in terms of dimensionless parameters is shown in figure 1.
Figure 1. Energy cascade in strongly stratified turbulence for 
$Pr\,Rb\gg 1$ relevant to the 
$Pr=O(1)$ regime of Earth's atmosphere and ocean. Here 
$E(k_\perp )$ is the horizontal energy spectrum of the velocity field. Energy is injected at low wavenumber 
$k_f$ (large scales) with power 
$P\sim \epsilon$ and undergoes a forward cascade down to dissipation scales. An anisotropic cascade results across horizontal wavenumbers in the range 
$k_f\lesssim k_\perp \lesssim k_O$, with associated vertical wavenumbers in the range 
$2{\rm \pi} /l_z\lesssim k_\parallel \lesssim k_O$. An isotropic cascade follows for wavenumbers larger than 
$k_O$ up to the dissipation wavenumber 
$k_\nu$. The strength of the anisotropy is quantified by 
$k_\perp /k_\parallel$.
Critical balance has provided a prediction for the anisotropy in the energy cascade for 
$Rb\gg 1$ and can therefore also predict the scaling of the outer vertical length and velocity scales of the system by considering the largest scales of the cascade. Substituting 
$k_\perp \sim 1/L$ and 
$k_\parallel \sim 1/l_z$ into (3.4) predicts that 
$l_z\sim Fr\,L$ for the outer vertical length scales. Enforcing incompressibility, 
$u_\parallel \sim (k_\perp /k_\parallel ) u_\perp$, subsequently predicts that 
$u_z=u_\parallel (2{\rm \pi} /L)\sim Fr\,U$ for the outer vertical velocity scale. These are exactly the scaling relations discussed in § 2.2 – critical balance successfully reproduces the known scaling relations in the 
$Pr\,Rb>1$ regime.
4. Critical balance and scaling for $Pr\,Rb<1$
We turn to deriving scaling relations in the 
$Pr\,Rb<1$ regime by extending the critical balance arguments presented in the previous section. The idea is to replace the adiabatic, inviscid IGW frequency 
$\omega _{{N}}$ with the corresponding frequency in the low (turbulent) Péclet-number limit 
$\omega _{{lPe}}$ for the estimate of the linear time scale. A linear analysis of the Boussinesq equations with 
$\nu =0$ and 
$\kappa \neq 0$ modelling the 
$Pr\ll 1$ limit gives the dispersion relation
\begin{equation} \omega^2-{\rm i}\omega\gamma_\kappa-\omega_{{N}}^2=0,\quad \gamma_\kappa=\kappa k^2, \quad \omega^2_N=N^2\frac{k^2_\perp}{k^2}. \end{equation} In the limit 
$Pr\,Rb\ll 1$ (i.e. 
$l_\kappa \gg l_O$), thermal diffusion is faster than the restoring buoyancy time scale, so 
$\gamma _\kappa \gg \omega _{{N}}$ and the two roots of (4.1a–c) become 
$\omega \sim {\rm i}\omega ^2_N/\gamma _\kappa$ and 
$\omega \sim {\rm i}\gamma _\kappa$. The latter is an uninteresting rapid thermal diffusion rate, while the former is an effective damping rate, 
$\omega _{{lPe}}={\rm i}\omega ^2_N/\gamma _\kappa ={\rm i}N^2k_\perp ^2/\kappa k^4$ (Lignières Reference Lignières1999, Reference Lignières2019). Thus, the linear response frequency is no longer a real frequency of a restoring oscillation, but instead a damping rate 
$\gamma _{{lPe}}$ (i.e. 
$\omega _{{lPe}}={\rm i}\gamma _{{lPe}}$) corresponding to the effective rate at which restoring buoyancy and strong thermal diffusion operate. A direct linear analysis of the low turbulent Péclet-number equations (2.7) also gives the eponymous damping rate 
$\omega ={\rm i}\gamma _{{lPe}}$ – the two approaches nicely agree. Using the expectation of strong anisotropy at large scales 
$k_\perp /k_\parallel \ll 1$, the linear damping time scale can be estimated as
\begin{equation} \gamma_{{lPe}}\sim \frac{N^2 k_\perp^2}{\kappa k_\parallel^4}. \end{equation} Before the critical balance hypothesis can be applied, a new justification is needed because the causality argument in the 
$Pr\,Rb>1$ regime no longer applies: waves are overdamped rather than propagate. The instantaneous propagation of information is a peculiarity of the low Péclet-number equations (2.7), since the temperature and vertical velocity fields are coupled by an elliptic equation to lowest order. Critical balance instead becomes an argument for selective decay. The dependence between the damping rate of a fluctuation and its vertical extent, 
$\gamma _{lPe}\sim l_\parallel ^{4}$ according to (4.2), means that longer vertical structures overdamp faster. As a result, any fluctuation at some 
$k_\perp$ with parallel extent longer than the 
$l_\parallel \sim k_\parallel ^{-1}$ set by 
$\gamma _{lPe}\tau _{NL}\sim 1$ will rapidly decay away before nonlinear effects can become significant. Critical balance in the thermally diffusive regime thus determines the longest parallel structure for a given 
$k_\perp$ that can sustain before nonlinear breakup. A sketch of the physical argument is shown in figure 2, alongside a comparison with the causality argument in the 
$Pr\,Rb>1$ regime.
Figure 2. A sketch of the physical arguments used for critical balance for eddies with perpendicular extent 
$l_\perp \sim k_\perp ^{-1}$. The notation 
$l_{\parallel,CB}$ refers to the 
$l_\parallel \sim k_\parallel ^{-1}$ that satisfies critical balance in each regime. Panel (a) shows the causality argument for the geophysical regime (
$Pr\,Rb>1$). Since the group velocity, 
$\boldsymbol {v}_{g}$, of the associated IGW sets the speed at which information can propagate, any eddy with 
$l_\parallel >l_{\parallel, CB}$ would be causally disconnected before it nonlinearly breaks up within time 
$t\sim \tau _{NL}$. Panel (b) shows the selective decay argument for the thermally diffusive regime (
$Pr\,Rb<1$). Since the damping rate increases with 
$\gamma _{lPe}\sim l_\parallel ^{4}$, any eddy with 
$l_\parallel > l_{\parallel, CB}$ would rapidly decay away before it could evolve.
With the modified physical interpretation, application of critical balance 
$\gamma _{{lPe}}\tau _{{NL}}\sim 1$ and rearrangement gives the new relationship between 
$k_\perp$ and 
$k_\parallel$ as
\begin{equation} \frac{k_\perp}{k_\parallel}\sim \left(\frac{\kappa\epsilon^{1/3}}{N^2}\right)^{3/4}k^2_\parallel=l_{OM}^2k_\parallel^2, \end{equation}
where the modified Ozmidov scale is defined as 
$l_{OM}=(\kappa \epsilon ^{1/3}/N^2)^{3/8}$. One can check that the critical balance prediction self-consistently maintains 
$\omega _{{N}}/\gamma _\kappa \ll 1$ all the way to the largest scales since 
$\omega _{{N}}/\gamma _\kappa \sim N k_\perp /\kappa k_\parallel ^3\sim (Pr\,Rb)^{1/4}\ll 1$. The turbulence now returns to isotropy at the modified Ozmidov scale where the overdamping rate for a fluctuation with 
$k_\perp \sim k_\parallel$ is comparable to its eddy turnover time, 
$N^2l_{OM}^2/\kappa \sim \tau _{NL}^{-1}$. In analogy to the Ozmidov scale, the modified Ozmidov scale is the largest horizontal scale that can overturn before overdamping becomes significant. Note that the modified Ozmidov is larger than the Ozmidov scale, as would be expected:
\begin{equation} l_{OM}=\frac{Fr^{3/2}}{(Pr\,Rb)^{3/8}}L,\quad \frac{l_{OM}}{l_O}=(Pr\,Rb)^{{-}3/8}>1. \end{equation} The outer vertical scale (
$k_\parallel \sim l_z^{-1}$) at the largest horizontal scale where 
$k_\perp \sim L^{-1}$ can again be found. Subsequent enforcement of incompressibility 
$(\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}'=0)$ then gives 
$u_z\sim l_z U/L$. The result is shown below:
\begin{equation} l_z\sim \frac{Fr}{(Pr\,Rb)^{1/4}}L, \quad u_z\sim \frac{Fr}{(Pr\,Rb)^{1/4}}U. \end{equation} These scalings self-consistently predict a small turbulent Péclet number 
$Pe_t\sim (Pr\,Rb)^{1/2}\ll 1$. At this point it becomes suggestive to define a modified Froude number as 
$Fr_M\equiv Fr/(Pr\,Rb)^{1/4}$ so that
\begin{equation} l_z=Fr_ML,\quad l_{OM}=Fr_M^{3/2}L. \end{equation}
These are the same exponents as in the geophysical regime. Using the new definition, the corresponding horizontal and vertical spectra for 
$k< k_{OM}$ are
\begin{equation} \frac{E(k_\perp)}{U^2L}\sim (k_\perp L)^{{-}5/3},\quad \frac{E(k_\parallel)}{U^2L}\sim Fr_M(k_\parallel l_z)^{{-}3}. \end{equation} The dimensional form of the parallel energy spectrum, 
$E(k_\parallel )\sim N(\epsilon /\kappa )^{1/2}k_\parallel ^{-3}$, now depends on the dissipation and thermal diffusivity, unlike in the geophysical regime, where 
$E(k_\parallel )\sim N^2k_\parallel ^{-3}$. A schematic of the energy cascade is shown in figure 3, and a comparison of the cascade path in the 
$k_\perp$–
$k_\parallel$ plane with the 
$Pr\,Rb>1$ regime is shown in figure 4.
Figure 3. Energy cascade for strongly stratified turbulence for 
$Pr\,Rb\ll 1$ relevant to the 
$Pr\ll 1$ regime of stellar radiative zones. Similar to figure 1, an anisotropic cascade results in the range 
$k_f\lesssim k_\perp \lesssim k_{OM}$ followed by an isotropic cascade from 
$k_{OM}$ to 
$k_\nu$. Note that the beginning of the isotropic subrange has moved to larger scales 
$k_{OM}< k_O$.
Figure 4. Comparison of the cascade path in the 
$k_\perp$–
$k_\parallel$ plane between the geophysical (teal) and thermally diffusive (orange) regimes. Energy is injected at low perpendicular wavenumber 
$k_\perp \sim L^{-1}$, follows the critical balance path until the respective Ozmidov scale (
$k_{O}$ or 
$k_{OM}$), enters the isotropic cascade and then is dissipated at the viscous scales.
It is now clear that the transition from the 
$Pr\,Rb>1$ to the 
$Pr\,Rb<1$ regime simply corresponds to a replacement 
$Fr\rightarrow Fr_M$. How can this modified Froude number be physically understood? If the Froude number is reinterpreted to define the ratio of the emergent vertical length scale to the imposed horizontal scale (i.e. 
$l_z/L\equiv Fr$), then critical balance at the largest scales simply sets the Froude number. By substituting 
$k_\perp \sim L^{-1}$ and 
$k_\parallel \sim l_z^{-1}$ into the linear wave frequencies 
$\omega _{{N}}\sim N l_z/L$ and 
$\gamma _{{lPe}}\sim N^2l_z^4/\kappa L^2$ and then comparing both with the corresponding nonlinear frequency scale 
$\tau _{{NL}}^{-1}\sim U/L$, the two Froude numbers emerge:
\begin{gather} Pr\,Rb>1{:} \quad \frac{l_z}{L}\sim \frac{U}{LN}\equiv Fr, \end{gather}
\begin{gather}Pr\,Rb<1{:} \quad \frac{l_z}{L}\sim \left(\frac{\kappa U}{N^2L^3}\right)^{1/4}\equiv Fr_M. \end{gather} As should be expected, the outer vertical scales smoothly transition from 
$l_z/L=Fr$ to 
$l_z/L=Fr/(Pr\,Rb)^{1/4}$ at 
$Pr\,Rb=1$ when 
$Pr$ is decreased. This is analogous to the smooth transition from 
$l_z/L=Fr$ to 
$l_z/L=Fr/Rb^{1/2}$ at 
$Rb=1$ when 
$Re$ is decreased, as discussed in § 2.3. Additionally, we note that 
$Fr_M$ can be rewritten in terms of 
$Pe\,Fr^{-2}$ and 
$Re$ (i.e. 
$Fr_M=(Pe\,Fr^{-2})^{-1/4}$), as was required by the low turbulent Péclet-number approximation.
With the acquired scaling relations from critical balance above, we suggest the following dimensionalization for rescaling the Boussinesq equations in the 
$Pr\,Rb\ll 1$ limit:
\begin{gather} \boldsymbol{u}_h'=U\boldsymbol{u}_h, \quad u_z'=Fr_MUu_z, \quad \theta'=\frac{1}{Fr_M}\frac{U^2}{L}\theta, \quad p'= \rho_m U^2p, \end{gather}
\begin{gather}x'=Lx, \quad y'=L y, \quad z'=Fr_ML z, \quad t'=\frac{L}{U} t. \end{gather}
Aside from 
$z'$ and 
$u_z'$, the only other variable whose scaling changed is 
$\theta '$, which is now determined by the new dominant balance between 
$N^2u_z'\sim \kappa {\nabla }_z^2\theta '$ as discussed in § 2.3. The Boussinesq equations become
\begin{gather} \partial_t \boldsymbol{u}_h+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_h={-}\boldsymbol{\nabla}_h p+ \left[\frac{1}{Re}{\nabla}_h^2+\frac{1}{Rb_M}{\nabla}_z^2\right]\boldsymbol{u}_h, \end{gather}
\begin{gather}Fr_M^2[\partial_t u_z+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}u_z]={-}\boldsymbol{\nabla}_z p +\theta+Fr_{M}^2\left[\frac{1}{Re}{\nabla}_h^2+\frac{1}{Rb_M}{\nabla}_z^2\right]u_z, \end{gather}
\begin{gather}(Pr\,Rb)^{1/2}[\partial_{t}\theta+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\theta]={-}u_z+ \left[\frac{(Pr\,Rb)^{1/2}}{Pr\,Re}{\nabla}_h^2+{\nabla}_z^2\right]\theta, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather} We see that 
$Rb_M\equiv Re\,Fr_M^2=Rb/(Pr\,Rb)^{1/2}$ acts as the new effective Reynolds number and nicely matches with the scale separation between the modified Ozmidov and the viscous scale, 
$l_{OM}/l_\nu =(Rb_M)^{3/4}$. To lowest order,
\begin{gather} \partial_t \boldsymbol{u}_h+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_h={-}\boldsymbol{\nabla}_h p, \end{gather}
\begin{gather}0={-}\boldsymbol{\nabla}_z p+\theta, \end{gather}
\begin{gather}0={-}u_z+{\nabla}_z^2\theta, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather} Comparing with the 
$Pr\,Rb>1$ equation set (2.4), the only difference appears in the dominant balance of the buoyancy equation, where 
$u_z\sim {\nabla }_z^2\theta$ instead of 
$u_z\sim \boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla }\theta$, as expected. The vertical momentum equation remains a balance of buoyancy fluctuations with the vertical pressure gradient. Advection of the vertical momentum is now suppressed by a factor of 
$Fr_M^2$ instead of 
$Fr^2$. In principle, this term can become of order unity at sufficiently low 
$Pr$ if 
$Pr< Fr^2/Re$; or, equivalently, when thermal diffusion is so efficient that 
$k_{OM}< k_f$ (or 
$\kappa >N^2L^3/U$), at which point the turbulence is simply isotropic.
Looking back at all the arguments so far in both the 
$Pr\,Rb<1$ and 
$Pr\,Rb>1$ regimes, we note the parallelism between the critical balance arguments for 
$\omega \tau _{{NL}}$ and the dominant balance arguments for the Boussinesq equations. The time scale 
$\tau _{{NL}}$ represented the nonlinear advection terms (i.e. 
$\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}_h$), while 
$\omega$ captured the combined effects of all the remaining linear terms. As dominant balance shifted between linear terms in various equations, 
$\omega$ shifted correspondingly to the dominant root of the full IGW dispersion relation (and vice versa). For example, 
$\omega _{{N}}\rightarrow {\rm i}\gamma _\nu \equiv {\rm i}\nu k_z^2$ as 
$Re$ was decreased in § 2, and 
$\omega _{{N}}\rightarrow \omega _{{lPe}}$ as 
$Pr$ was decreased in § 4. Critical balance and the corresponding dominant balance arguments are therefore one and the same.
4.1. Role of vertical shear instabilities
Examining the role of vertical shear instabilities in the turbulence offers an alternative physical interpretation of the scaling results. Many experimental and numerical studies in the geophysical regime find evidence suggesting that turbulent structures are organized to be marginally unstable to local vertical shear instabilities, a behaviour thought to be closely related to self-organized criticality (Smyth & Moum Reference Smyth and Moum2013; Salehipour et al. Reference Salehipour, Peltier and Caulfield2018; Smyth et al. Reference Smyth, Nash and Moum2019; Chini et al. Reference Chini, Michel, Julien, Rocha and Colm-cille2022; Lefauve & Linden Reference Lefauve and Linden2022). A heuristic argument for such behaviour is that the vertical shear of a turbulent eddy will generally be driven towards the marginal value of the shear for instability. An eddy with a strong vertical shear will go unstable within a dynamical time and reduce its shear to the marginal value, while one with a weak shear can become amplified to the marginal value before going unstable. As a result, vertical gradients of the horizontal velocity of eddies may be maintained at the marginal condition for vertical shear instability. As a concrete example, vertically adjacent quasi-horizontal eddy motions (e.g. ‘pancake’ eddies or modes) have been suggested and observed to exhibit such behaviour in early studies, such as Lin & Pao (Reference Lin and Pao1979), Lilly (Reference Lilly1983), Spedding, Browand & Fincham (Reference Spedding, Browand and Fincham1996) and Riley & de BruynKops (Reference Riley and de BruynKops2003) (and references therein). We now examine how the critical balance scalings in both the geophysical and thermally diffusive regimes marginally satisfy the corresponding vertical shear instability criteria.
In the 
$Pr\,Rb>1$ (i.e. 
$Pe_t\sim Pr\,Rb>1$) regime, the nonlinear criterion to sustain turbulence in a vertical shear instability requires 
$J=N^2/S^2\lesssim 1$ (Richardson Reference Richardson1920; Howard & Maslowe Reference Howard and Maslowe1973) (i.e. Richardson's criterion), which is essentially based on an energy argument requiring more kinetic energy release than potential energy cost upon instability of the flow. To satisfy marginal stability, eddies with horizontal velocity scale 
$u_\perp (k_\perp )$ would need to maintain a vertical separation 
$l_\parallel$ such that 
$J(k_\perp )=N^2/S^2\sim 1$ holds with 
$S\sim u_\perp /l_\parallel$. Indeed, substituting the critical balance results (i.e. 
$k_\perp E(k_\perp )\sim u_\perp ^2(k_\perp )$ alongside (3.4)) into 
$J(k_\perp )\sim N^2l_\parallel ^2/u_\perp ^2(k_\perp )$ shows that 
$J(k_\perp )\sim 1$ is satisfied on a scale-by-scale basis from 
$k_f\leq k_\perp \leq k_O$. Interpretation of strongly stratified turbulence in terms of marginal instability of local vertical shear instabilities is therefore consistent with critical balance scalings in the geophysical regime.
In the 
$Pr\,Rb<1$ (i.e. 
$Pe_t\sim (Pr\,Rb)^{1/2}<1$) regime, the nonlinear instability criterion is relaxed due to the weakening of the potential energy cost by the significant thermal diffusion rate at large vertical scales (
$u_z/l_z<\kappa /l_z^2$ since 
$Pe_t<1$). The criterion proposed on phenomenological grounds by Zahn (Reference Zahn1992) requires 
$J Pe_t\lesssim 1$ for instability, which smoothly connects across the 
$Pe_t=1$ transition and has been approximately verified (in accessible parameter regimes) in recent numerical studies (Prat & Lignières Reference Prat and Lignières2013, Reference Prat and Lignières2014; Garaud, Gallet & Bischoff Reference Garaud, Gallet and Bischoff2015; Prat et al. Reference Prat, Guilet, Viallet and Müller2016; Garaud, Gagnier & Verhoeven Reference Garaud, Gagnier and Verhoeven2017). To satisfy marginal stability, eddies with horizontal velocity scale 
$u_\perp (k_\perp )$ would need to maintain a vertical separation 
$l_\parallel$ such that 
$J(k_\perp )Pe_t(k_\perp )=N^2u_\parallel l_\parallel / S^2\kappa \sim 1$. Similar to before, substituting the critical balance results (i.e. 
$k_\perp E(k_\perp )\sim u_\perp ^2(k_\perp )$ alongside (4.3)) into 
$J(k_\perp )Pe_t(k_\perp )$ shows that 
$J(k_\perp )Pe_t(k_\perp )\sim 1$ is satisfied on a scale-by-scale basis from 
$k_f\leq k_\perp \leq k_{OM}$. Thus the consistent interpretation between the marginal instability of local vertical shear instabilities and critical balance scalings is maintained across the 
$Pe_t=1$ transition into the thermally diffusive regime.
We note that some studies add the additional requirement 
$Re_t=u_zl_z/\nu >Re_{crit}\sim 10^3$ to ensure that viscous effects play no role in suppressing the nonlinear criterion for instability. Combining this assumption with 
$JPe_t\lesssim 1$ gives a weaker condition for instability, 
$JPr\lesssim (JPr)_{crit}$ (necessary, but not sufficient, since this is equivalent to 
$JPe_t\lesssim Re_t/Re_{crit}$). Reporting simulations in terms of 
$JPr$ can be useful since it may be a proxy of 
$Re_t\sim (JPr)^{-1}$ (Prat et al. Reference Prat, Guilet, Viallet and Müller2016). However, requiring 
$Re_t\gg 1$ as an additional assumption is unnecessary for strongly stratified turbulence because 
$Re_t\sim Rb_M\gg 1$ is already built in to the asymptotics. This is physically reasonable because the dynamics of scales larger than the modified Ozmidov scale cannot be affected by viscous effects (i.e. since 
$l_\nu \ll l_{OM}$).
4.2. Comparison with previous work
This work follows in the light of a series of simulations and an evolving discussion in Cope et al. (Reference Cope, Garaud and Caulfield2020) and Garaud (Reference Garaud2020) attempting to understand the complicated parameter space of 
$Pr<1$ stratified turbulence driven by horizontal shear instabilities. Cope et al. (Reference Cope, Garaud and Caulfield2020) examine the low Péclet-number limit where 
$Pe<1$ and 
$Pe_t\ll 1$, which falls squarely into our 
$Pr\,Rb<1$ regime and allows for comparison. Predictions of the vertical outer scales from critical balance (
$u_z\sim Fr_MU$ and 
$l_z\sim Fr_ML$) are in conflict with the theoretical predictions of Cope et al. (Reference Cope, Garaud and Caulfield2020) shown below:
\begin{equation} u_z\sim Fr_M^{2/3}U,\quad l_z\sim Fr_M^{4/3}L, \end{equation}
where their notation has been translated using 
$Fr_M=(BPe)^{-1/4}$.
A major source of the discrepancy arises due to different assumptions in Cope et al. (Reference Cope, Garaud and Caulfield2020) for the dominant balance of the vertical momentum equation. Examining the proposed balances directly on the low turbulent Péclet-number equations gives a clear way to compare. The horizontal component of the momentum equation (2.7) requires the pressure to be of order unity, 
$p\sim U^2/L$. For the vertical component, our arguments in this work essentially balance 
$\boldsymbol {\nabla }_z p'\sim (N^2/\kappa ) \nabla ^{-2}u_z'$, which directly implies 
$l_z\sim Fr_M L$ by using 
$\nabla ^{-2}\sim l_z^2$ and enforcing incompressibility, 
$u_z/l_z\sim U/L$, at the outer scales. The vertical advection terms are then smaller by 
$O(Fr_M^2)$ (see (4.10b)) and result in a consistent balance similar to the 
$Pr\,Rb>1$ regime (see (2.3b)). Cope et al. (Reference Cope, Garaud and Caulfield2020) alternatively assume a balance of vertical advection and the thermally modified buoyancy term: 
$\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } u_z'\sim (N^2/\kappa ) \nabla ^{-2}u_z'$. This assumption leads to an inconsistency because it results in a vertical pressure gradient that is too large (
$\boldsymbol {\nabla }_z p'\gg \boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla } u_z'$), as discussed in § 5.2.2 of Cope et al. (Reference Cope, Garaud and Caulfield2020). It is unclear how (4.12a,b) would smoothly transition from the 
$Pr\,Rb>1$ regime, perhaps requiring an additional intermediate turbulence regime. We also note that (4.12a,b) do not satisfy the incompressibility constraint, which stems from their assumption that the horizontal and vertical length scales of the spatial derivatives are similar. This assumption may be justified in their case, where turbulence is driven in a wide band of horizontal modes by the horizontal shear instability, while our arguments suppose an idealized forcing at a single horizontal scale.
The large set of 
$Pe<1$ simulations in Cope et al. (Reference Cope, Garaud and Caulfield2020) (e.g. their figure 8) should allow differentiation between the two theoretical predictions in principle; however, the comparison is difficult because of the small differences between the fractional exponents. A detailed follow-up analysis of the simulations in light of the new critical balance predictions hopefully should resolve the question of the correct scaling relations, and will be explored in forthcoming work. It is promising that Garaud (Reference Garaud2020) finds a sharp transition at 
$Pe_t\sim 1$ (as shown in her figure 3), which we interpret as the 
$Pr\,Rb\sim 1$ transition.
We note that Garaud (Reference Garaud2020) predicts an additional turbulent regime for 
$Pr<1$, 
$Pe\gg 1$ and 
$Pe_t>1$. According to the critical balance framework, this falls into the 
$Pr\,Rb>1$ regime, and scaling from the geophysical literature (
$l_z\sim FrL$) should simply apply. Garaud (Reference Garaud2020) instead theoretically proposes, and finds numerical evidence for, scaling relations given by 
$u_z/U\sim l_z/L\sim Fr^{2/3}$. Similar to our 
$Pr\,Rb>1$ regime, Garaud (Reference Garaud2020) argues that 
$\boldsymbol {\nabla }_zp'\sim \theta '$ and 
$\boldsymbol {u}'\boldsymbol {\cdot } \boldsymbol {\nabla }\theta '\sim -N^2u_z'$, but, instead of assuming incompressibility to get 
$l_z/L\sim Fr$, Garaud argues for a constant scaling of the time-dependent mixing efficiency to arrive at 
$l_z/L\sim Fr^{2/3}$, with motivation and support from her numerical results (see her figure 4). As discussed in § 4.4.4 of Garaud (Reference Garaud2020), several factors may be at play. First, 
$l_z$ may be sensitive to the method used to extract the vertical scale from simulations. For example, Lindborg (Reference Lindborg2006) use a weighted average of 
$k_z$ with the energy spectrum to find a 
$l_z/L\sim Fr$ scaling in 
$Pr=1$ simulations, while Garaud (Reference Garaud2020) uses the vertical correlation length in her 
$Pr=O(10^{-1})$ simulations. The different methods correspond to different weighting functions when integrating with the energy spectrum. It would be useful to compare both diagnostics at the same 
$Pr$ to resolve this possible issue. Another issue is that reaching asymptotic values of 
$Rb\gg 1$ while keeping 
$Fr\ll 1$ is computationally expensive due to a requirement for two large-scale separations between 
$L$ and 
$l_{OM}$ as well as 
$l_{OM}$ and 
$l_{\nu }$. Measuring scaling exponents in simulations with 
$Fr=O(10^{-1})$ and 
$Rb=O(10^1)$ may lead to transitional scaling relations (Bartello & Tobias Reference Bartello and Tobias2013). Indeed, the arguments presented in this paper and the theoretical literature, strictly speaking, only rigorously hold for the asymptotic limits of 
$Pr\,Rb\gg 1$ or 
$Pr\,Rb\ll 1$.
5. Astrophysical applications
5.1. Diffusion coefficients
Mixing of chemical elements in stellar interiors by stratified turbulence can have important consequences for stellar evolution, and quantifying its efficiency is essential for comparison with stellar observations (see Maeder & Meynet Reference Maeder and Meynet2000; Salaris & Cassisi Reference Salaris and Cassisi2017; and references therein). One-dimensional stellar evolution models require an effective vertical diffusion coefficient, 
$D_v$, that is typically estimated by the product of some characteristic vertical velocity and length scales of the turbulence, i.e. 
$D_v\sim u_z l_z$. Which scales to choose remains an important problem and may depend on the instability that drives the turbulence. Typical choices include using the outer vertical scales or the largest isotropic scales. The scale-dependent anisotropy obtained from critical balance allows an estimate of the contribution to 
$D_v$ from eddies of different scales in stratified turbulence with horizontal length scale 
$L$ and vertical scale 
$U$ set by some instability (e.g. shear instabilities of differential rotation). We define the scale-dependent diffusion coefficient 
$\tilde {D}_v(k_\parallel )=u_\parallel (k_\parallel )k_\parallel ^{-1}$ from eddies of vertical length scale 
$k_\parallel ^{-1}$ and corresponding vertical velocity scale 
$u_\parallel (k_\parallel )$. Using the incompressibility relation and the anisotropy relations for 
$k_\perp /k_\parallel$ gives
\begin{gather} k_\parallel\leq k_{O(M)}{:} \quad \frac{\tilde{D}_v(k_\parallel)}{UL}=Fr_{(M)}^2, \end{gather}
\begin{gather}k\geq k_{O(M)}{:} \quad \frac{\tilde{D}_v(k)}{UL}=\frac{Fr_{(M)}^2}{(kl_{O(M)})^{4/3}}\quad (k_\parallel\sim k), \end{gather}
where 
$k_{O(M)}$ is defined to be 
$k_O$ for 
$Pr\,Rb>1$ or 
$k_{OM}$ for 
$Pr\,Rb<1$ (and similarly for 
$l_{O(M)}$ and 
$Fr_{(M)}$). The turbulent diffusion coefficient in principle has contributions from turbulence at all scales, but is often argued to be dominated by either the outer scales, 
$D_v\sim \tilde {D}_v(l_z^{-1})$, or the Ozmidov scales, 
$D_v\sim \tilde {D}_v(l_{O(M)}^{-1})$. The constant value of 
$\tilde {D}_v$ in the large-scale anisotropic subrange for 
$k_\parallel < k_{O(M)}$ means that one can equally use either choice, i.e. 
$D_v\sim u_\parallel (l_z^{-1})l_z\sim u_\parallel (l_{O(M)}^{-1})l_{O(M)}$. Thus scaling relations based on critical balance predict that 
$D_v\sim Fr^2UL$ in the 
$Pr\,Rb>1$ regime and 
$D_v\sim Fr_M^2UL$ in the 
$Pr\,Rb<1$ regime. The estimated turbulent diffusion coefficient in the 
$Pr\,Rb<1$ limit is larger than in the 
$Pr\,Rb>1$ limit by a factor of
\begin{equation} \frac{D_v^{[Pr\,Rb\,<\,1]}}{D_v^{[Pr\,Rb\,>\,1]}}\sim \left(\frac{Fr_M}{Fr}\right)^2=\frac{1}{(Pr\,Rb)^{1/2}}, \end{equation}
which is significant only if 
$Pr\,Rb\ll 1$.
The foundational work of Zahn (Reference Zahn1992) (see also Lignières Reference Lignières2019) uses the modified Ozmidov scales, 
$D_v\sim u_\parallel (l_{OM}^{-1})l_{OM}\sim (\epsilon \kappa /N^2)^{1/2}$, in agreement with the prediction above. Our estimates for 
$D_v$, therefore, do not affect the results of previous astrophysical studies based on the estimate for 
$D_v$ by Zahn (Reference Zahn1992). Lastly, we note that the diffusion coefficient in the 
$Pr\,Rb<1$ limit agrees with the prediction from Cope et al. (Reference Cope, Garaud and Caulfield2020) (using the outer vertical scales for estimation of 
$D_v\sim u_\parallel (l_z^{-1})l_z$), although the agreement may be coincidental given their different predictions for 
$u_z$ and 
$l_z$ as discussed in § 4.2.
5.2. Small-scale dynamo instability criterion
Turbulence in the fully ionized and highly conductive plasma of a stellar radiative zone may be able to generate and sustain magnetic fields on scales smaller than the forcing scale through the small-scale dynamo (SSD). Anisotropy in the velocity field caused by stable stratification makes dynamo action less efficient compared to isotropic turbulence at the same magnetic Reynolds number, 
$Rm=UL/\eta$, where 
$\eta$ is the resistivity. Indeed, the limit of infinitely strong stratification leads to a planar velocity field (i.e. 
$\boldsymbol {u}=(u_x(x,y,z,t),u_y(x,y,z,t),0)$), which is a known anti-dynamo flow (Zeldovich & Ruzmaikin Reference Zeldovich and Ruzmaikin1980). Using a large set of direct numerical simulations, Skoutnev, Squire & Bhattacharjee (Reference Skoutnev, Squire and Bhattacharjee2021) found that the SSD is unstable in the 
$Pr=O(1)$ regime if the scale separation of the Ozmidov scale and the magnetic resistive scale is sufficiently large. Quantitatively, this translates to requiring the magnetic buoyancy Reynolds number 
$Rb_m=Pm\,Rb$ to be larger than a critical value (i.e. 
$Rb_m>Rb_m^c$), where 
$Pm=\nu /\eta$ is the magnetic Prandtl number. The direct analogy to isotropic turbulence is the requirement that 
$Rm$ be larger than a critical value 
$Rm>Rm^c$.
We propose that, in the 
$Pr\,Rb<1$ regime, the SSD criterion switches to requiring a sufficient scale separation between the modified Ozmidov scale and the resistive scale. The criterion in the two regimes is then given by
\begin{gather} Pr\,Rb>1{:} \quad Rb_m=Pm\,Rb>Rb_m^c \end{gather}
\begin{gather}Pr\,Rb<1{:} \quad Rb_{m,M}=Pm\,Rb_M=\frac{Pm\,Rb}{(Pr\,Rb)^{1/2}}>Rb_{m,M}^c. \end{gather}
This extension naturally carries through the conjecture in Skoutnev et al. (Reference Skoutnev, Squire and Bhattacharjee2021) that the SSD only operates within the isotropic portion of the turbulent cascade. The new criterion predicts that the SSD will become more unstable as 
$Pr$ is decreased at fixed 
$Rb$ into the 
$Pr\,Rb<1$ regime since 
$Rb_M>Rb$. Qualitatively, dynamo action becomes more efficient as the anisotropy caused by stratification is decreased by increased thermal diffusion.
We expect that both 
$Rb_m^c=Rb_m^c(Pm)$ and 
$Rb_{m,M}^c=Rb_{m,M}^c(Pm)$ should have a weak 
$Pm$ dependence that is strongest around 
$Pm=O(1)$, similar to the dependence of 
$Rm^c(Pm)$ in the case of isotropic turbulence (Iskakov et al. Reference Iskakov, Schekochihin, Cowley, McWilliams and Proctor2007; Schekochihin et al. Reference Schekochihin, Iskakov, Cowley, McWilliams, Proctor and Yousef2007). It would be reasonable that they are of comparable magnitudes, 
$Rb_{m,M}^c\sim Rb_{m}^c$, since both act as effective outer magnetic Reynolds numbers for their respective isotropic subranges.
6. Conclusion
Critical balance is a theory for strong turbulence in anisotropic wave systems that argues for a balance of linear wave and nonlinear interaction time scales to predict the scale-by-scale structure of the turbulent cascade. We have proposed that critical balance can unify the unity 
$Pr$ geophysical regime and the extremely low 
$Pr$ astrophysical regime of strongly stratified turbulence because both support anisotropic linear wave motions in different asymptotic limits of the IGW dispersion relation. The dispersion reduces to adiabatic inviscid IGWs in the 
$Pr=O(1)$ limit and to IGWs overdamped on buoyancy-modified time scales in the 
$Pr\ll 1$ limit. We find that a smooth transition between the two regimes occurs at 
$Pr\,Rb=O(1)$ as 
$Pr$ is decreased, or equivalently when the turbulent Péclet number 
$Pe_t=u_zl_z/\kappa$ drops below order unity and shifts dominant balance in the buoyancy equation. Application of critical balance in the 
$Pr\,Rb<1$ regime predicts an anisotropic cascade and scaling relations for the outer vertical scales that are identical to the geophysical regime if the Froude number, 
$Fr$, is simply replaced by a modified Froude number, 
$Fr_M\equiv Fr/(Pr\,Rb)^{1/4}$ (e.g. the outer vertical length scale is 
$l_z\sim Fr_ML$). Indeed, a smooth transition occurs at 
$Pr\,Rb=1$ between the two regimes.
The scaling relations from critical balance offer a theoretical framework for understanding the properties of strongly stratified turbulence in stellar radiative zones. Application to estimating the vertical turbulent diffusion coefficient in the 
$Pr\,Rb<1$ regime gives the same scaling result as the original estimates of Zahn (Reference Zahn1992), thereby validating their dimensional arguments. However, a complete understanding of strongly stratified turbulence in stellar radiative zones would need also to take into account the effects of non-local energy transfer, rotation and magnetic fields (large and small scale). In analogy with the geophysical regime, energy at large scales in the thermally diffusive regime is likely also transferred through non-local transfer mechanisms, not only through a local energy cascade as assumed in our critical balance arguments. Rotation may have a significant effect when the driving horizontal scales are sufficiently large to be influenced by Coriolis forces, which is thought to be the case in the upper solar radiative zone (Garaud Reference Garaud2020), for example. Similarly, large-scale magnetic fields will modify the wave dispersion relation and support magneto-internal waves, whose properties will change the nature of the turbulence presented here depending on the magnitude and direction of the large-scale field. Lastly, small-scale magnetic fields can grow in the isotropic portion of the turbulence through the small-scale dynamo if the stratification is not too strong in the 
$Pr\,Rb>1$ regime (Skoutnev et al. Reference Skoutnev, Squire and Bhattacharjee2021). We proposed a modified instability criterion for the SSD in the 
$Pr\,Rb<1$ limit based on the new scaling laws. The modified instability criterion predicts a more easily excited dynamo, since thermal diffusion ameliorates the effect of stratification. Stellar stratified turbulence is, therefore, at the very least, in a saturated state of the SSD, containing small-scale magnetic fields in near equipartition with the isotropic scales of the velocity field. The effects of a saturated SSD on mixing and transport are unknown.
Acknowledgements
The author would like to thank M. Kunz for an early suggestion to explore more generally the ideas of critical balance, A. Bhattacharjee for providing intuition of their application to MHD turbulence, and P. Garaud for many insightful discussions on stratified turbulence that eventually led to this work. We also thank the referees for improvements to the physical interpretation of the results. We thank K. Whitman for illustration of the graphical abstract.
Funding
V.A.S. was supported by the Max-Planck/Princeton Center for Plasma Physics (NSF grant PHY-1804048).
Declaration of interests
The author reports no conflict of interest.
