© The Authors 2023

1. Introduction

The interplay of turbulent convection with the overall rotation of the Sun is the primary cause of differential rotation that is observed at the solar surface and in the interior (e.g., Rüdiger 1989; Miesch & Toomre 2009). Three-dimensional numerical simulations solving the equations of magnetohydrodynamics (MHD) capture the essence of this process and routinely produce solutions that are qualitatively similar to those for the Sun with equatorial acceleration (e.g., Gilman 1983; Brun et al. 2004; Guerrero et al. 2013; Käpylä et al. 2014). However, it has become increasingly clear recently that even the most sophisticated current simulations are missing something essential. The most striking manifestation of this is that simulations using a nominal solar luminosity and rotation rate often produce anti-solar (AS) differential rotation with equatorial deceleration (e.g., Fan & Fang 2014; Käpylä et al. 2014; Hotta et al. 2015), whereas solar-like (SL) differential rotation is achieved only with significantly more rapid rotation (e.g., Viviani et al. 2018; Matilsky et al. 2020).

This is related to the convective conundrum (O’Mara et al. 2016), which is essentially the difference between large-scale velocity amplitudes in simulations in comparison to the Sun (e.g., Hanasoge et al. 2012, 2016; Schumacher & Sreenivasan 2020). Until recently, the most common way to ensure SL differential rotation in simulations with solar luminosity and rotation rate has been to lower the convective velocities by artificially enhancing the radiative diffusivity (e.g., Fan & Fang 2014; Käpylä et al. 2014; Hotta et al. 2016). This cannot be justified based on physical grounds, however, since convection is thought to carry practically all of the energy flux through the solar convection zone (CZ), with the exception of very deep layers. Another more plausible effect is due to magnetic fields: it is conceivable that sufficiently strong fields can suppress convection on large scales to a degree where the differential rotation flips from AS to SL. Early results with relatively low-resolution simulations were mixed: Karak et al. (2015) found essentially no dependence on the magnetic field, while Fan & Fang (2014) and Simitev et al. (2015) reported more positive outcomes. Nevertheless, the magnetic Reynolds numbers of these simulations were most probably not high enough to excite a small-scale dynamo. This was addressed by the recent high-resolution simulations of Hotta & Kusano (2021) and Hotta et al. (2022), which suggest that SL differential rotation can indeed be achieved with the help of an efficient small-scale dynamo.

Another important parameter is the Prandtl number, Pr = ν/χ, where ν is the kinematic viscosity and χ is the thermal diffusivity. The notion that the solar convection zone is operating in a regime in which the effective Prandtl number is large has recently gained popularity (e.g., O’Mara et al. 2016; Bekki et al. 2017; Karak et al. 2018). While these studies indicate that the overall velocity amplitudes are decreased in these set-ups, the problem with the differential rotation is even exacerbated (Karak et al. 2018). This is because it is not only the velocity amplitude that is sensitive to Pr, but turbulent transport of angular momentum and heat are also affected (e.g., Cattaneo et al. 1991; Käpylä 2021). Furthermore, theoretical arguments suggest that Pr ≪ 1 in the solar CZ (e.g., Ossendrijver 2003; Schumacher & Sreenivasan 2020).

Prandtl numbers that strongly deviate from unity are challenging numerically, and therefore most simulations are conducted in the Pr ≈ 1 regime. It is commonly acknowledged that reaching realistic parameter regimes in terms of, for example, Prandtl, Reynolds, and Rayleigh numbers with current or foreseeable simulations of stellar convection is not feasible (e.g., Kupka & Muthsam 2017). The main aim of the present study is to vary the Prandtl number within the range that can be achieved with numerical simulations with values above and below unity. The current study is also inspired by recent results from hydrodynamic non-rotating convection in Cartesian geometry (Käpylä 2021), where the convective energy transport and velocity statistics were found to be sensitive to the effective Prandtl number.

2. The model

The simulation set-up is similar to the set-ups used in Käpylä et al. (2019) and Käpylä et al. (2020). The simulation domain is a spherical wedge that spans r_in < r < R in radius, where r_in = 0.7R and R is the radius of the star, θ₀ < θ < π − θ₀ in colatitude, where θ₀ = π/12, and 0 < ϕ < π/2 in longitude. Equations of fully compressible MHD were solved,

$$\begin{array}{cc}& \frac{\partial \mathit{A}}{\partial t}=\mathit{U}\times \mathit{B}-\eta {\mu }_0\mathit{J},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\mathrm{D}ln\rho }{\mathrm{D}t}=-\mathbf{\nabla }\mathbf{·}\mathit{U},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\mathrm{D}\mathit{U}}{\mathrm{D}t}=\mathit{g}-2{\mathbf{\Omega }}_0\times \mathit{U}-\frac{1}{\rho }(\mathbf{\nabla }p-\mathit{J}\times \mathit{B}-\mathbf{\nabla }\mathbf{·}2\nu \rho \mathsf{S}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& T\frac{\mathrm{D}s}{\mathrm{D}t}=\frac{1}{\rho }[\eta {\mu }_0{\mathit{J}}^2-\mathbf{\nabla }\mathbf{·}({\mathit{F}}^{\mathrm{rad}}+{\mathit{F}}^{\mathrm{SGS}})]+2\nu {\mathsf{S}}^2,\hfill \end{array}$$(4)

where A is the magnetic vector potential, U is the velocity, B = ∇ × A is the magnetic field, η is the magnetic diffusivity, μ₀ is the permeability of vacuum, J = ∇ × B/μ₀ is the current density, D/Dt = ∂/∂t + U⋅∇ is the advective time derivative, ρ is the density, and g = −∇ϕ is the acceleration due to gravity, where ϕ = −GM/r is a fixed spherically symmetric gravitational potential, with G and M being the universal gravitational constant and the stellar mass, respectively. Ω₀ = (cos θ, −sin θ, 0)Ω₀ is the angular velocity vector, where Ω₀ is the rotation rate of the frame of reference, p is the pressure, ν is the kinematic viscosity, T is the temperature, and s is the specific entropy with Ds = c_VDlnp − c_PDlnρ, where c_V and c_P are the specific heat capacities in constant volume and pressure, respectively. The gas is assumed to obey the ideal gas law, p = ℛρT, where ℛ = c_P − c_V is the gas constant. The rate of strain tensor is given by

$$\begin{array}{c}\hfill {\mathsf{S}}_{\mathit{ij}}=\frac{1}{2}(U_{i;j}+U_{j;i})-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }\mathbf{·}\mathit{U},\end{array}$$(5)

where the semicolons refer to covariant derivatives (Mitra et al. 2009). The radiative flux is given by

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{rad}}=-K\mathbf{\nabla }T,\end{array}$$(6)

where K is the heat conductivity. The latter consists of two parts, K = K₁ + K₂, where K₁ = K₁(r) is a fixed function of height and K₂ = K₂(ρ, T) is density- and temperature-dependent according to Kramers opacity law (Weiss et al. 2004). The profile of K₁ is given by

$$\begin{array}{c}\hfill K_1=K_{\mathrm{top}}[1+tanh\left(\frac{r-R}{d_K}\right)],\end{array}$$(7)

where K_top = 2c_PF_bot/|g(R)|, with $F_{\text{bot}}=L/4\pi r_{\text{in}}^2$ where L is the luminosity of the star, and where d_K = 0.015R. The contribution K₂ is given by

$$\begin{array}{c}\hfill K_2(\rho ,T)=K_0{(\rho /{\rho }_0)}^{-(a+1)}{(T/T_0)}^{3-b},\end{array}$$(8)

where ρ₀ and T₀ are reference values of density and temperature, and the values a = 1 and b = −7/2 correspond to the Kramers opacity law. This formulation was first used in convection simulations by Brandenburg et al. (2000).

The subgrid scale (SGS) flux is given by

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{SGS}}=-{\chi }_{\mathrm{SGS}}\rho \mathbf{\nabla }{s}^{\prime },\end{array}$$(9)

where χ_SGS is the (constant) SGS diffusion coefficient for the entropy fluctuation s′(r, θ, ϕ) = s − ⟨s⟩_θϕ, where ⟨s⟩_θϕ is the spherically symmetric part of the specific entropy. The SGS flux does not contribute to the net radial energy transport because it is decoupled from the mean stratification, and therefore changing χ_SGS does not lead to drastic changes in the boundary layer thickness near the surface.

The simulations were made using the PENCIL CODE¹ (Brandenburg 2021). The code in this study employs third-order temporal and sixth-order spatial discretisation. Advective terms in Eqs. (1)–(4) are written as fifth-order upwinding derivatives with a sixth-order hyperdiffusive correction, where the diffusion coefficient is flow-dependent; see Appendix B of Dobler et al. (2006).

2.1. System parameters and diagnostics quantities

The simulations are defined by the energy flux imposed at the bottom boundary, F_bot = −(K∂T/∂r)|_r = rin, the values of K₀, a, b, ρ₀, T₀, Ω₀, ν, η, χ_SGS, the profile of K₁, and the value of the modified Stefan-Boltzmann constant σ_SB in the upper boundary condition ${\sigma }_{\text{SB}}{T}_{\text{surf}}^4=-K\partial T/\partial r$, where T_surf is the (unconstrained) surface temperature. The current models use a significantly enhanced luminosity in comparison to real stars to bring the thermal and dynamical timescales close enough to be resolved in the simulations. This leads to correspondingly higher convective velocities, and therefore, the rotation rate is increased accordingly to capture a similar rotational influence on the flow in the simulations in comparison to real stars; see Appendix A of Käpylä et al. (2020).

The non-dimensional luminosity is given by

$$\begin{array}{c}\hfill \mathcal{L}=\frac{L_0}{{\rho }_0{(GM)}^{3/2}{R}^{1/2}},\end{array}$$(10)

where ρ₀ is the initial density at the base of the convection zone. The degree of luminosity enhancement is given by the ratio L_ratio = ℒ/ℒ_⊙ ≈ 2.1 × 10⁵, where ℒ_⊙ is the dimensionless solar luminosity. The initial stratification is determined by the non-dimensional pressure scale height at the surface,

$$\begin{array}{c}\hfill {\xi }_0=\frac{\mathcal{R}{T}_1}{GM/R},\end{array}$$(11)

where T₁ = T(R, t = 0).

The relative strengths of viscosity, SGS diffusion, and magnetic diffusivity are given by the SGS and magnetic Prandtl numbers,

$$\begin{array}{c}\hfill {\mathrm{Pr}}_{\mathrm{SGS}}=\frac{\nu }{{\chi }_{\mathrm{SGS}}},\mathrm{Pm}=\frac{\nu }{\eta }.\end{array}$$(12)

We use Pm = 1 in most of the runs and vary Pr_SGS between 0.1 and 10. The thermal Prandtl number related to the radiative conductivity is given by

$$\begin{array}{c}\hfill \mathrm{Pr}=\frac{\nu }{\chi },\end{array}$$(13)

where χ = K/c_Pρ is the radiative diffusivity, which in general varies as a function of radius, latitude, and time. In the current simulations, χ_SGS ≫ χ almost everywhere. The efficiency of convection is quantified by the Rayleigh number,

$$\begin{array}{c}\hfill \mathrm{Ra}=\frac{GM{(\mathrm{\Delta }r)}^4}{\nu {\chi }_{\mathrm{s}}{R}^2}{(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}{s}_{\mathrm{hs}}}{\mathrm{d}r})}_{r_{\mathrm{s}}},\end{array}$$(14)

where Δr = 0.3R is the depth of the layer, s_hs is the specific entropy in a one-dimensional non-convecting hydrostatic model, evaluated near the top of the domain at r_s = 0.95R, and where χ_s is the total thermal diffusivity K/c_Pρ from r = r_s. The hydrostatic solution is Schwarzschild-unstable only in a thin layer near the surface (see e.g., Barekat & Brandenburg 2014; Brandenburg 2016) which is why the Rayleigh number is evaluated at r_s. Moreover, χ_SGS does not contribute to Ra because it only acts on deviations from the spherically symmetric specific entropy. Additionally, a turbulent Rayleigh number is quoted,

$$\begin{array}{c}\hfill {\mathrm{Ra}}_{\mathrm{t}}=\frac{GM{(\mathrm{\Delta }r)}^4}{\nu {\chi }_{\mathrm{tot}}{R}^2}{(-\frac{1}{c_{\mathrm{P}}}\frac{\mathrm{d}{⟨s⟩}_{\theta \varphi }}{\mathrm{d}r})}_{r_{\mathrm{s}}},\end{array}$$(15)

where ⟨s⟩_θϕ is the time- and horizontal average of the specific entropy, and χ_tot = χ_SGS + ⟨χ⟩_θϕ is the total thermal diffusivity. Ra_t is always significantly smaller than Ra because χ_SGS ≫ χ.

The magnitude of the rotation is controlled by the Taylor number,

$$\begin{array}{c}\hfill \mathrm{Ta}=\frac{4{\mathrm{\Omega }}_0^2{(\mathrm{\Delta }r)}^4}{{\nu }^2}.\end{array}$$(16)

The fluid and magnetic Reynolds numbers and the Péclet number are given by

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_1},{\mathrm{Re}}_{\mathrm{M}}=\frac{u_{\mathrm{rms}}}{\eta k_1},\mathrm{Pe}=\frac{u_{\mathrm{rms}}}{{\chi }_{\mathrm{SGS}}{k}_1},\end{array}$$(17)

respectively, where $u_{\mathrm{rms}}=\sqrt{\frac{3}{2}(U_r^2+U_{\theta }^2)}$ is the time- and volume-averaged rms velocity, where $U_{\varphi }^2$ has been replaced by $(U_r^2+U_{\theta }^2)/2$ to avoid contributions from differential rotation. The inverse of the wavenumber k₁ = 2π/Δr ≈ 21/R is used to characterize the radial extent of the convection zone. Several definitions of the Coriolis number that describes the rotational influence on the flow are discussed in Sect. 3.2.

Mean quantities are denoted by overbars and are defined by the time- and azimuthal average,

$$\begin{array}{c}\hfill \overline{f}(r,\theta )=\frac{1}{\mathrm{\Delta }\varphi \mathrm{\Delta }t}{\displaystyle {\int }_{t_0}^{t_0+\mathrm{\Delta }t}{\int }_0^{\mathrm{\Delta }\varphi }f(r,\theta ,\varphi ,t)\mathrm{d}\varphi \mathrm{d}t,}\end{array}$$(18)

where t₀ and Δt are the beginning and the length of the statistically steady part of the simulation, and where Δϕ = π/2 is the azimuthal extent of the simulation domain. Error estimates were obtained by dividing the time series into three parts and by computing averages over each one of them. The largest deviation of these sub-averages from the average over the whole time series was taken to represent the error.

2.2. Initial and boundary conditions

Initially, the stratification is isentropic with a polytropic index n = 1.5 and ξ₀ = 0.02, resulting in an initial density contrast of 30. The value of K₀ is chosen such that F_rad = F_tot at the bottom of the domain.

The radial and latitudinal boundaries are assumed to be impenetrable and stress-free for the flow. At the bottom boundary, a fixed heat flux is prescribed, while at the top, a black-body condition is applied. At the latitudinal boundaries, the gradients of thermodynamic quantities are set to zero; see Käpylä et al. (2013). For the magnetic field, we apply a radial field condition at the upper and a perfect conductor condition at the lower boundary. At the latitudinal boundaries, the field is assumed to be tangential to the boundary. These conditions are given in terms of the magnetic vector potential by

$$\begin{array}{cc}& A_r=0,\frac{\partial A_{\theta }}{\partial r}=-\frac{A_{\theta }}{r},\frac{\partial A_{\varphi }}{\partial r}=-\frac{A_{\varphi }}{r}(r=R),\hfill \end{array}$$(19)

$$\begin{array}{cc}& \frac{\partial A_r}{\partial r}=A_{\theta }=A_{\varphi }=0(r=r_{\mathrm{in}}),\hfill \end{array}$$(20)

$$\begin{array}{cc}& A_r=\frac{\partial A_{\theta }}{\partial \theta }=A_{\varphi }=0(\theta ={\theta }_0,\pi -{\theta }_0).\hfill \end{array}$$(21)

The azimuthal direction is periodic for all quantities. The velocity and magnetic fields were initialized with random low-amplitude Gaussian noise fluctuations.

3. Results

Three sets of simulations were run, where Pr_SGS = 0.1 (set P01), 1 (P1), and 10 (P10), respectively. The first two sets contain hydrodynamic and MHD runs, and a subset of the MHD runs was remeshed to higher resolution and to correspondingly higher Rayleigh, Péclet, and Reynolds numbers; see Table 1. Only MHD variants of the P10 runs were run.

3.1. Thermal relaxation and convergence

Before discussing the main results, we briefly describe the measures we undertook to ensure that the simulations are thermally relaxed and that the results converged. The relevant timescale for thermal relaxation is the Kelvin-Helmholtz timescale, which is estimated as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{{\mathcal{E}}_{\mathrm{th}}}{L},\end{array}$$(22)

where ℰ_th = ∫_VE_thdV is the total thermal energy of the initial non-convecting state, with E_th = ρe, and where e = c_VT is the specific internal energy. In general, the runtime of the simulations, Δt, needs to be of the order of τ_KH to ensure thermal relaxation. A shorter runtime suffices if the initial condition is close to the final thermally relaxed convecting state. However, this can typically be concluded only a posteriori.

The normalized Kelvin-Helmholtz timescale in the current simulations is ${\stackrel{\sim }{\tau }}_{\mathrm{KH}}={\tau }_{\mathrm{KH}}/{\tau }_{\mathrm{ff}}\approx 4.1\times {10}^3$, where ${\tau }_{\mathrm{ff}}=\sqrt{GM/R^3}$ is the free-fall time. By virtue of the enhanced luminosity method used here (see e.g., Käpylä et al. 2020), τ_KH can be resolved in the simulations; see the 13th column of Table 1. In the hydrodynamical simulations, the initial transient, data from which is not used in the time-averaging, is typically of the order of τ_KH; see Fig. 1a for data from run P1-6H. Hydrodynamic simulations were considered thermally relaxed when the volume-averaged (denoted by angle brackets) total energy density ⟨E_tot⟩ reaches a statistically steady state, where ⟨E_tot⟩=⟨E_th⟩+⟨E_pot⟩+⟨E_kin⟩+⟨E_mag⟩ with E_pot = ρϕ, $E_{\mathrm{kin}}=\frac{1}{2}\rho {\mathit{U}}^2$, and $E_{\mathrm{mag}}=\frac{1}{2}{\mu }_0^{-1}{\mathit{B}}^2$. The flows often reach a statistically steady state already significantly earlier; see the inset of Fig. 1a for representative results for the time evolution of the rms value of U_ϕ from run P1-6H.

Fig. 1. Time series of total energy density and the rms velocity and magnetic fields. (a) Volume-averaged total energy density ⟨Etot⟩ normalized by its value at t = 0 from run P1-6H. The inset shows the volume-averaged rms value of the azimuthal flow ${\stackrel{\sim }{U}}_{\varphi }^{\mathrm{rms}}=U_{\varphi }^{\mathrm{rms}}/\sqrt{GM/R}$. The red (blue) part of the curve indicates the transient (statistically stationary) part of the time series, where ${\stackrel{\sim }{t}}_0=t_0/{\tau }_{\mathrm{ff}}$ indicates the starting time for the time averaging. τKH indicates the Kelvin-Helmholtz timescale. (b) Same as above, but for the corresponding MHD run P1-6M. The inset shows ${\stackrel{\sim }{B}}_{\varphi }^{\mathrm{rms}}=B_{\varphi }^{\mathrm{rms}}/({\rho }_0{\mu }_0\sqrt{GM/R})$ in addition to ${\stackrel{\sim }{U}}_{\varphi }^{\mathrm{rms}}$. The black (green) parts of the data for ${\stackrel{\sim }{B}}_{\varphi }^{\mathrm{rms}}$ indicate transient (statistically stationary) parts. (c) Same as panel a, but for run P01-3H.

The MHD simulations were typically started from thermally relaxed hydrodynamic cases. In these cases, the transient due to thermal relaxation is significantly shorter than in the hydrodynamic progenitor run. However, the growth of the magnetic field is slow at the modest magnetic Reynolds numbers presented here. Therefore, another two-phase transient follows, in which the magnetic field first grows from a weak seed field to become dynamically important, which is followed by a period in which the thermal state relaxes toward a new statistically steady MHD state; see Fig. 1b. Again, the criterion for convergence is taken to be that the sum of the volume-averaged total energy density ⟨E_tot⟩ settles to a new statistically steady state. In some cases, such as in run P1-6M shown in Fig. 1b, the total energy density in the MHD state still shows temporal modulation due to quasi-cyclic large-scale magnetism.

However, in a few of hydrodynamic cases (runs P01-2H, P01-3H, and P01-4H), a statistically stationary state was not reached, even though these models were run several Kelvin-Helmholtz times; see Fig. 1c. These runs exhibit excursions between different flow states where the mean azimuthal flow also shows sign changes, leading to a large uncertainty in the estimates of differential rotation. The most probable explanation for this behaviour is that these models lie in a bistable parameter regime in which different differential rotation profiles, including AS and SL states, are possible, depending on the history of the run (Gastine et al. 2014; Käpylä et al. 2014). Whether these models finally settle to one or the other stable solution is unclear and would possibly require significantly longer integration times. None of the corresponding MHD runs exhibits similar excursions, thus being consistent with the absence of hysteresis in the MHD regime (Karak et al. 2015).

3.2. Differential rotation and meridional circulation

The main focus of the current study is to explore the effects of the Prandtl number for the large-scale flows that develop in rotating convective systems. The mean rotation profile is given by the time and azimuthal average,

$$\begin{array}{c}\hfill \overline{\mathrm{\Omega }}(r,\theta )={\mathrm{\Omega }}_0+\frac{{\overline{U}}_{\varphi }(r,\theta )}{rsin\theta },\end{array}$$(23)

and the meridional flow is given by ${\overline{\mathit{U}}}_{\mathrm{mer}}=({\overline{U}}_r,{\overline{U}}_t,0)$. In many simulations, the latitudinal profiles of $\overline{\mathrm{\Omega }}$ are non-monotonic such that the rotation rate has a polar jet, a maximum at mid-latitudes, or sometimes several local minima and maxima as a function of latitude. Furthermore, equatorial asymmetries can occur, rendering the amplitude of the latitudinal shear an unreliable diagnostic of the overall sense of differential rotation; see representative examples in Figs. 2–4. Therefore the classification of the AS and SL rotation profile is here based on the mean rotation profile at the equator

$$\begin{array}{c}\hfill ⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩=\frac{{\int }_{r_{\mathrm{in}}}^R{r}^2[\stackrel{\sim }{\overline{\mathrm{\Omega }}}(r,{\theta }_{\mathrm{eq}})-1]\mathrm{d}r}{{\int }_{r_{\mathrm{in}}}^R{r}^2\mathrm{d}r},\end{array}$$(24)

Fig. 2. Time-averaged rotation profiles from MHD runs P01-[2-5]M near the AS-SL transition with PrSGS = 0.1 with the rotation rate increasing from left to right. Coriolis numbers according to Eqs. (25) and (28) and the mean differential rotation at the equator $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩$ according to Eq. (24) are indicated in each panel. The white arrows indicate the meridional flow.

Fig. 3. Same as Fig. 2 but for runs P1-[2-5]M in set P1.

Fig. 4. Same as Fig. 2 but for runs P10-[6-9]M in set P10.

where θ_eq = π/2, and where the tildes refer to normalization by the rotation rate of the frame of reference, Ω₀. If $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩>0$ ($⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩<0$), the run is classified as SL (AS) rotator. This measure turns out to be a monotonic function of rotation and is furthermore much less affected by equatorial asymmetries or polar jets. See Camisassa & Featherstone (2022) for a more granular classification making use of more diagnostics to describe the transition.

The current results indicate that the convective velocity increases when the Prandtl number is decreased. This is manifest when the fluid Reynolds number increases for a decreasing SGS Prandtl number; see the eighth column of Table 1. Naively, it might then be expected that achieving SL differential rotation for low Pr_SGS would be more difficult, that is, require faster rotation. The rotational influence on the flow is often quantified by a simple definition of the Coriolis number,

$$\begin{array}{c}\hfill \mathrm{Co}=\frac{2{\mathrm{\Omega }}_0}{u_{\mathrm{rms}}{k}_1},\end{array}$$(25)

where the convective length scale is assumed to be unchanged by rotation. When this definition is used to characterize the results, the Coriolis number at which the rotation profile flips from AS to SL decreases monotonically with Pr_SGS; see Fig. 5a where $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩$ is shown for all runs as a function of Co. That is, in the low-resolution MHD runs with Pr_SGS = 0.1 (1), the transition occurs around Co ≈ 1.5 (Co ≈ 1.8); see Figs. 2 and 3, whereas for Pr_SGS = 10 the transition occurs at an even higher Coriolis number (Co ≈ 3); see Fig. 4.

Fig. 5. Normalized average angular velocity at the equator $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩$ for all runs as a function of Co (a), Coℓ(b), Coω(c), and Co⋆(d). Circles (crosses) denote MHD (hydrodynamic) runs, and the colour of the symbols indicates the SGS Prandtl number. The sizes of the symbols are proportional to ReM (Re) for MHD (hydrodynamic) runs.

However, the validity of this simplistic definition of the Coriolis number in characterizing the simulations can be questioned based on its very crude estimate of the convective length scale and velocity amplitude. For example, Gastine et al. (2014) argued that a local Rossby (inverse Coriolis) number based on the length scale from the mean spherical harmonic degree ${\overline{\ell }}_u$ of the m ≠ 0 poloidal flows gives a more accurate estimate (see also Schrinner et al. 2012). Furthermore, they showed that the scatter near the AS-SL transition is reduced with this definition, essentially reducing the apparent dependence on Prandtl number significantly. Here this was tested by computing ${\overline{\ell }}_u$ according to

$$\begin{array}{c}\hfill {\overline{\ell }}_u={\displaystyle \sum _{\ell }\ell \frac{⟨{\mathit{U}}_{\mathrm{p}}^{\ell }\mathbf{·}{\mathit{U}}_{\mathrm{p}}^{\ell }⟩}{⟨{\mathit{U}}_{\mathrm{p}}\mathbf{·}{\mathit{U}}_{\mathrm{p}}⟩},}\end{array}$$(26)

where U_p is the non-axisymmetric poloidal flow, and where ℓ is the spherical harmonic degree. Data from a varying number of horizontal slices from near the base, at the middle, and near the top of the CZ were analyzed for each run, and the resulting ${\overline{\ell }}_u$ is an average over the depths and time. The number of time slices per run varies between 7 and roughly 60. The corresponding length scale is $d_u=\pi \mathrm{\Delta }R/{\overline{\ell }}_u$, and the Coriolis number based on this is given by

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\ell }=\frac{2{\mathrm{\Omega }}_0{d}_u}{u_{\mathrm{rms}}},\end{array}$$(27)

where u_rms is defined in the same way as in Eq. (17). The numbers in Table 1, fifth column, indicate that the value of Co_ℓ is sensitive to the Reynolds number such that in the runs with the highest Re, the values of Co_ℓ are roughly 30% lower than those of the low-resolution runs. Similarly, a Coriolis number based on fluctuating vorticity can be defined as (e.g., Brun et al. 2022)

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\omega }=2{\mathrm{\Omega }}_0/{\omega }_{\mathrm{rms}}^{\prime },\end{array}$$(28)

where ω′ = ∇ × u with $\mathit{u}=\mathit{U}-\overline{\mathit{U}}$. This quantity shows a similar sensitivity to the fluid Reynolds number as Co_ℓ; see the sixth column of Table 1. Both of these definitions pick up smaller length scales at the highest Reynolds number runs, resulting in lower Coriolis numbers. This is likely an indication that the simulations are still far away from an asymptotic regime in which the results would be independent of the diffusivities. Nevertheless, Co_ℓ and Co_ω characterize the rotational influence on the flow more accurately than Co by being sensitive to the actual dominant length scale. However, due to the Reynolds number sensitivity, we should only compare results from runs with comparable Re when these definitions are used to characterize the results.

The average radial differential rotation at the equator is shown as functions of Co_ℓ and Co_ω in Figs. 5b and c, respectively. When the five runs at higher Reynolds and Péclet numbers are ignored for the time being, these results suggest that the dependence of the differential rotation transition as a function of the Prandtl number all but vanishes for Pr_SGS < 1. However, in both cases, the transition for Pr_SGS = 10 still occurs at a higher Co_ℓ and Co_ω than for the Pr_SGS = 0.1 and 1 cases². These results agree with those reported by Karak et al. (2018), who also found that a Prandtl number above unity promotes AS differential rotation due to enhanced downward angular momentum transport. Both definitions give much lower Coriolis numbers for the higher-Re runs because smaller convective scales are resolved and ${\overline{\ell }}_u$ and ω′ pick these up. Furthermore, among the five higher-Re runs, the single AS model (P01-1Mh) appears to have a marginally larger Coriolis number than the SL counterparts in both cases, although the error estimates for Co_ℓ and Co_ω are large in both cases. Whether this is a real effect or due to insufficient statistics remains open at this point. As a side note, the resemblance of Figs. 5b and c, or alternatively, the correlation between Co_ω and Co_ℓ, suggests that Co_ω captures the rotation dependence of the convective length scale almost as well as Co_ℓ without having to perform numerically expensive spherical harmonic decomposition.

Each of the definitions of the Coriolis number discussed so far rely on diagnostic quantities such as u_rms, ${\overline{\ell }}_u$, and ω′, that are sensitive to other details, such as the fluid Reynolds number, of the system. Yet another alternative is to define a Coriolis number that depends only on stellar input parameters such as the luminosity, mass, and rotation rate of the modeled star. This can be constructed by assuming that

$$\begin{array}{c}\hfill L={\rho }_{\star }{u}_{\star }^3{R}^2,\end{array}$$(29)

where $L=4\pi r_{\text{in}}^2{F}_{\text{bot}}$ is the luminosity, ρ_⋆ is a reference density, and u_⋆ is an estimate of an average convective velocity. We assume ρ_⋆ = ρ₀³ and construct a stellar Coriolis number

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\star }=\frac{2{\mathrm{\Omega }}_{0}R}{u_{\star }}=\frac{2{\mathrm{\Omega }}_0{R}^{5/3}{\rho }_0^{1/3}}{L^{1/3}}.\end{array}$$(30)

While this definition is imperfect in the sense that the actual convective flow speed or scale do not enter, it is useful to determine whether a model with a given rotation rate and luminosity is an AS or SL rotator. This makes particularly sense in the homogeneous set of simulations considered here, where the stellar mass, luminosity, and radius are all fixed. The results are shown in Fig. 5d. It is apparent that the runs with Pr_SGS = 10 require a significantly higher Co_⋆, corresponding here to higher Ω₀, to achieve SL differential rotation. However, the AS-SL transition occurs at the same Co_⋆ for Pr_SGS = 1 and 0.1. Notably, the MHD runs at Pr_SGS = 0.1 and 1 for Co_⋆ = 0.50 are at least marginally SL in comparison to the corresponding AS hydrodynamic cases. Similarly, the higher-resolution runs are SL or marginally SL for Co_⋆ = 0.46, whereas the corresponding lower-resolution MHD and hydrodynamic runs are all AS. The current higher-resolution runs with Re_M ≈ 94…113 also have small-scale dynamos, which was tested with separate runs where the axisymmetric (m = 0) magnetic fields were suppressed. Therefore, it is plausible that the main contribution to the earlier appearance of SL differential rotation in these cases is due to the growing importance of the magnetic fields. However, no corresponding higher-resolution hydrodynamic runs were made to confirm this. The current results agree with the results of Hotta & Kusano (2021) and Hotta et al. (2022), who argue in favour of the magnetic field being the decisive factor in the transition. Finally, when the simulations are scaled to physical units as in Käpylä et al. (2020), the lowest value in the present study, Co_⋆ = 0.38, corresponds to the case of solar rotation rate at solar luminosity.

The conclusion of comparing the results characterized in terms of the four variants of the Coriolis number is that SL differential rotation is substantially more difficult to obtain for Pr_SGS = 10 than for Pr_SGS = 1 and 0.1, whereas in the latter two cases, there is no clear difference. Furthermore, the use of the simplistic Coriolis number, Eq. (25), gives misleading results and should be avoided. In contrast to the SGS Prandtl number, the dependence on magnetic fields is clearer, such that in MHD runs, SL differential rotation is easier to excite.

3.3. Angular momentum transport

The angular momentum in the interior of the star is governed by the conservation equation

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(\overline{\rho }\mathcal{L})+\mathbf{\nabla }\mathbf{·}\mathcal{T}=0,\end{array}$$(31)

where $\mathcal{L}={\varpi }^2\overline{\mathrm{\Omega }}$ is the specific angular momentum, with ϖ = r sin θ being the lever arm, and where

$$\begin{array}{cc}\hfill \mathcal{T}& =\varpi [\overline{\rho \mathit{u}{u}_{\varphi }}+\overline{\rho }\overline{\mathit{U}}({\overline{U}}_{\varphi }+\varpi {\mathrm{\Omega }}_0)\hfill \\ \hfill & -(\overline{\mathit{b}{b}_{\varphi }}+\overline{\mathit{B}}{\overline{B}}_{\varphi })/{\mu }_0-2\nu \rho \overline{\mathsf{S}}\mathbf{·}\widehat{\mathit{\varphi }}],\hfill \end{array}$$(32)

is the flux of angular momentum, where $\mathit{u}=\mathit{U}-\overline{\mathit{U}}$ and $\mathit{b}=\mathit{B}-\overline{\mathit{B}}$ are the fluctuating velocity and magnetic field, and where $\widehat{\mathit{\varphi }}$ is the unit vector in the azimuthal direction. The internal angular velocity profile is thus determined by the spatial distribution of the angular momentum fluxes. The main transporters are due to the Reynolds and Maxwell stress due to fluctuating and mean flows and fields (see Rüdiger 1989; Rüdiger & Hollerbach 2004),

$$\begin{array}{cc}& {\mathcal{T}}^{\mathit{u}}=\varpi {\overline{\mathit{Q}}}^{\mathit{u}}=\varpi \overline{\rho \mathit{u}{u}_{\varphi }},\hfill \end{array}$$(33)

$$\begin{array}{cc}& {\mathcal{T}}^{\mathit{b}}=\varpi {\overline{\mathit{M}}}^{\mathit{b}}=-\varpi \overline{\mathit{b}{b}_{\varphi }}/{\mu }_0,\hfill \end{array}$$(34)

$$\begin{array}{cc}& {\mathcal{T}}^{\overline{\mathit{U}}}=\varpi {\overline{\mathit{Q}}}^{\overline{\mathit{U}}}=\varpi \overline{\rho }\overline{\mathit{U}}({\overline{U}}_{\varphi }+\varpi {\mathrm{\Omega }}_0),\hfill \end{array}$$(35)

$$\begin{array}{cc}& {\mathcal{T}}^{\overline{\mathit{B}}}=\varpi {\overline{\mathit{M}}}^{\overline{\mathit{B}}}=-\varpi \overline{\mathit{B}}{\overline{B}}_{\varphi }/{\mu }_0.\hfill \end{array}$$(36)

Additionally, the contributions to the stress due to molecular viscosity are given by

$$\begin{array}{c}\hfill {\mathcal{T}}^{\nu }=\varpi {\overline{\mathit{Q}}}^{\nu }=-2\varpi \nu \rho \overline{\mathsf{S}}\mathbf{·}\widehat{\mathit{\varphi }}.\end{array}$$(37)

Although the viscous flux is negligibly small in the Sun, it is typically non-negligible in simulations (e.g., Guerrero et al. 2022).

Representative results of 𝒯 are shown in Fig. 6 superimposed on top of the angular velocity in runs P1-2M and P1-4M. The total angular momentum flux closely follows the pattern of the meridional circulation; see also Fig. 3. To see which terms dominate the angular momentum transport, we proceed to compute the divergences of the various terms that are encapsulated in 𝒯. Figure 7 show the divergences of the various contributions to the angular momentum flux given by Eqs. (33)–(37) for runs P1-2M and P1-4M. The contributions due to the turbulent Reynolds stress ${\overline{\mathit{Q}}}^{\mathit{u}}$ and due to the meridional flow ${\overline{\mathit{Q}}}^{\overline{\mathit{U}}}$ are by far the largest and approximately balance out in all of the cases considered here. The viscous and Maxwell stress contributions are clearly subdominant. The former (latter) decreases (increases) in the higher-resolution runs, but ${\overline{\mathit{M}}}^{\mathit{b}}$ remains subdominant in these cases as well. The rightmost panels of Fig. 7 show the residual divergence ∇⋅𝒯, which has relatively high peak values, but primarily on small scales.

Fig. 6. Total angular momentum flux 𝒯 (arrows) superimposed on the mean angular velocity profile (colour contours) from runs P1-2M and P1-4M.

Fig. 7. Divergences of the different contributions to the angular momentum flux, Eqs. (33)–(37), from runs P1-2M (top row) and P1-4M (bottom). The rightmost panels show the divergence of the total stress, and the tildes refer to normalization by ρ0GM.

Because the divergence of the stress is very similar in the AS and SL cases, the difference in the angular momentum transport between these cases is rather subtle. Furthermore, the stress from meridional circulation is radially outward near the equator in both AS and SL cases; see for example Fig. 3, and the viscous and magnetic contributions are weak in all of the current runs. Therefore the most plausible explanation is a change in the turbulent Reynolds stress 𝒯^u, which is shown in Fig. 8 for runs P1-2M and P1-4M superimposed on the mean angular velocity. The turbulent Reynolds stress is nearly radially downward at high latitudes in the bulk of the CZ. The region of downward radial flux is roughly confined inside the tangent cylinder in all cases. Outside the tangent cylinder and near the surface at all latitudes, 𝒯^u is directed predominantly equatorward regardless of whether the rotation profile is AS or SL. Outside the tangent cylinder, the angular momentum transport is increasingly axial, although the radial flux at the equator remains non-zero in all cases. The sense of the differential rotation is related to the sign of the radial component of 𝒯^u at the equator: for positive (outward) 𝒯^u, the differential rotation is SL (run P1-4M in Fig. 8), whereas it is AS for negative (downward) 𝒯^u (run P1-2M). Therefore, it appears sufficient to study the equatorial radial angular momentum flux to determine the difference between AS and SL differential rotation.

Fig. 8. Same as Fig. 6 but showing the turbulent Reynolds stress 𝒯u.

To study the equatorial Reynolds stress in more detail, a number of slices of the m ≠ 0 velocity field U(r, θ_eq, ϕ) were analyzed. Azimuthal Fourier filtering was applied to produce filtered velocity fields ${\mathit{U}}^{\text{f}}(m_1^{\prime },m_2^{\prime })$, where azimuthal orders ranging between $m_1^{\prime }$ and $m_2^{\prime }$ were retained⁴. These flows are used to compute the radial Reynolds stress,

$$\begin{array}{c}\hfill {\overline{Q}}_{r\varphi }^{\mathit{u}}(m_1^{\prime },m_2^{\prime })=\frac{1}{\mathrm{\Delta }\varphi \mathrm{\Delta }t}{\displaystyle {\int }_{t_0}^{t_0+\mathrm{\Delta }t}{\int }_0^{\mathrm{\Delta }\varphi }\rho U_r^{\mathrm{f}}(m_1^{\prime },m_2^{\prime })U_{\varphi }^{\mathrm{f}}(m_1^{\prime },m_2^{\prime })\mathrm{d}\varphi \mathrm{d}t.}\end{array}$$(38)

Here the density fluctuations are assumed to be small, and no filtering was applied to ρ. Representative results are shown in Fig. 9. In the AS case P1-2M, Fig. 9a, the total Reynolds stress is negative everywhere except at the base of the CZ. Contributions from the largest scale ($m_1^{\prime }$, $m_2^{\prime }$ = 1, 2) non-axisymmetric motions are statistically almost identical with the total stress. A weak positive contribution around the middle of the CZ is visible for ($m_1^{\prime }$, $m_2^{\prime }$ = 3, 5), but the Reynolds stress for ($m_1^{\prime }$, $m_2^{\prime }$ = 1, 5) is again very similar to the total stress. In the SL runs, the largest scales ($m_1^{\prime }$, $m_2^{\prime }$ = 1, 2) also contribute to a downward flux, whereas the main contribution to the net outward flux comes from ($m_1^{\prime }$, $m_2^{\prime }$ = 3, 5); see Fig. 9b for run P1-4M. Similarly to the AS case, the contributions from m′> 5 are small. This indicates that practically all of the Reynolds stress at the equator is due to relatively large-scale structures that can be identified as Busse columns (Busse 1970a,b), which are also often referred to as banana cells. Such features are often prominently visible in snapshots of the velocity field; see Fig. 10 for a representative example. The Busse cells are manifestations of non-linear prograde-propagating thermal Rossby waves. It is therefore somewhat questionable to talk about turbulent Reynolds stress in this context since the Busse columns are essentially large-scale convective modes.

Fig. 9. Fourier-filtered and total Reynolds stress component ${\stackrel{\sim }{\overline{Q}}}_{r\varphi }^{\mathit{u}}$ from runs P1-2M (a) and P1-4M (b) as indicated in the legend. The shaded areas indicate error estimates according to the definition in Sect. 2.1.

Fig. 10. Instantaneous normalized radial Reynolds stress component ${\stackrel{\sim }{\overline{Q}}}_{r\varphi }^{\mathit{u}}$ on the equatorial plane from run P01-6M (colour contours). The m ≠ 0 flows are indicated by the arrows, the width of which is proportional to the local flow amplitude.

The mechanism by which the differential rotation is generated in the current simulations is therefore different from that in Hotta et al. (2022), where the small-scale Maxwell stress is the dominant contribution to the radial angular momentum transport. While the high-resolution runs in the present study have small-scale dynamos and show an increased tendency for SL differential rotation, the Reynolds stress due to the thermal Rossby waves is still the dominant contribution to the angular momentum flux in all of the runs considered here. In an earlier study (Käpylä et al. 2017), the Maxwell stresses were found to be comparable to the Reynolds stress at the highest magnetic Reynolds numbers, but in that study, the magnetic Prandtl number was higher (Pr_M = 5), and the modelled stars typically rotated three to four times faster than in the current study. Although the Maxwell stress dominates the angular momentum transport in the simulations of Hotta et al. (2022), large-scale Busse columns can still be seen in the deep parts of their model; see for example their Fig. 6. If such large-scale convective patters were as prominent in the Sun, they should have been detected by helioseismology, but there is currently no evidence to this effect. Therefore, it seems that although highly magnetized simulations are more solar-like in terms of the rotation profile, the conundrum with the too prominent large-scale structures remains.

Furthermore, in hydrodynamic mean-field theories of differential rotation (e.g., Rüdiger 1989; Kitchatinov & Rüdiger 2005; Rogachevskii & Kleeorin 2015), the turbulence models are necessarily simplified, and large-scale convection modes such as Busse columns do not appear. In the most commonly adopted approach of Kitchatinov & Rüdiger (2005), the radial angular momentum transport at the equator is downward for slow rotation and vanishes for rapid rotation. The SL differential rotation results from a strong equatorward transport. Numerical simulations of isothermal homogeneous anisotropic turbulence also produce downward (slow rotation) or vanishing (rapid rotation) radial angular momentum flux at the equator (Käpylä 2019a) in qualitative accordance with Kitchatinov & Rüdiger (2005). Hydrodynamic mean-fields models based on these concepts do not typically produce AS solutions unless strong magnetic fields are present (Kitchatinov & Rüdiger 2004), although more recently, a hydrodynamic mechanism has also been discussed (Rüdiger et al. 2019). This latter process relies on poleward horizontal angular momentum flux at slow rotation, which was also found from local simulations, but which appears to be absent in global models such as those presented here. The mean-field theories avoid the problem of too prominent thermal Rossby waves by simply neglecting them, whereas they fail to characterize both AS and SL cases in the current simulations. This difference is yet another facet of the convective conundrum, the resolution of which is likely to require further critical assessment of both theoretical and simulation approaches.

4. Conclusions

The transition from AS to SL differential rotation was studied as a function of the SGS Prandtl number (Pr_SGS). Four definitions of the Coriolis number were used to quantify the exact point of transition from simulations where the rotation of the star was varied. While this transition occurs at a higher Coriolis number for Pr_SGS = 10 than for Pr_SGS = 1 and 0.1, no statistically relevant difference was found between the last two cases. This suggests that the Prandtl number dependence of the AS-SL transition is weak for Pr_SGS < 1, whereas a high Prandtl number makes it significantly more difficult to achieve SL differential rotation.

These results are puzzling because earlier non-rotating local simulations (Käpylä 2021) suggested that the cases Pr = 1 and Pr = 0.1 also differ significantly in many respects. A notable difference to the study of Käpylä (2021) is that the current simulations do not include a radiative layer below the CZ. This can explain why no subadiabatic layers develop at the base of the CZs in the current simulations: the effects of overshooting are absent. The latter was found to be particularly sensitive to the Prandtl number in the local simulations Käpylä (2019b) and Käpylä (2021). The inclusion of a radiative layer has also consequences for the dynamo solutions (e.g., Guerrero & Smolarkiewicz 2016; Käpylä 2022), which also couple back to differential rotation. Furthermore, rotation was not taken into account in the earlier Cartesian studies. These aspects need to be revisited in future studies.

Many of the current simulations also included magnetic fields, but often in a parameter regime in which the small-scale dynamo is not excited. Thus the influence of magnetic fields is relatively weak in most of the current runs. Nevertheless, the magnetic fields make it easier to excite SL differential rotation especially in the current higher-resolution runs, which likely also have small-scale dynamos. However, the effects of magnetic fields are likely to be more significant in higher Re_M systems, as shown in the recent results of Hotta et al. (2022). Therefore, magnetism appears to be the most promising candidate to explain the discrepancy between solar observations and global simulations. Nevertheless, the non-detection of thermal Rossby waves from the Sun, which remain prominent in all current simulations, still raises questions about the generation mechanism of solar differential rotation. Finally, the difference between mean-field theories of differential rotation and 3D simulation results is pointed out as another aspect that requires further scrutiny in the future.

Acknowledgments

The anonymous referee is acknowledged for their constructive comments. I acknowledge the hospitality of Nordita during the program “The Shifting Paradigm of Stellar Convection: From Mixing Length Concepts to Realistic Turbulence. Modelling”. I also thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme DYT2 where part of the work on this paper was undertaken. This work was supported by EPSRC grant EP/R014604/1. The simulations were made within the Gauss Center for Supercomputing project “Cracking the Convective Conundrum” in the Leibniz Supercomputing Centre’s SuperMUC–NG supercomputer in Garching, Germany. This work was supported by the Deutsche Forschungsgemeinschaft Heisenberg programme (grant No. KA 4825/4-1). This work was also partially supported by a grant from the Simons Foundation.

References

All Tables

All Figures

Fig. 1. Time series of total energy density and the rms velocity and magnetic fields. (a) Volume-averaged total energy density ⟨Etot⟩ normalized by its value at t = 0 from run P1-6H. The inset shows the volume-averaged rms value of the azimuthal flow ${\stackrel{\sim }{U}}_{\varphi }^{\mathrm{rms}}=U_{\varphi }^{\mathrm{rms}}/\sqrt{GM/R}$. The red (blue) part of the curve indicates the transient (statistically stationary) part of the time series, where ${\stackrel{\sim }{t}}_0=t_0/{\tau }_{\mathrm{ff}}$ indicates the starting time for the time averaging. τKH indicates the Kelvin-Helmholtz timescale. (b) Same as above, but for the corresponding MHD run P1-6M. The inset shows ${\stackrel{\sim }{B}}_{\varphi }^{\mathrm{rms}}=B_{\varphi }^{\mathrm{rms}}/({\rho }_0{\mu }_0\sqrt{GM/R})$ in addition to ${\stackrel{\sim }{U}}_{\varphi }^{\mathrm{rms}}$. The black (green) parts of the data for ${\stackrel{\sim }{B}}_{\varphi }^{\mathrm{rms}}$ indicate transient (statistically stationary) parts. (c) Same as panel a, but for run P01-3H.

In the text

Fig. 2. Time-averaged rotation profiles from MHD runs P01-[2-5]M near the AS-SL transition with PrSGS = 0.1 with the rotation rate increasing from left to right. Coriolis numbers according to Eqs. (25) and (28) and the mean differential rotation at the equator $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩$ according to Eq. (24) are indicated in each panel. The white arrows indicate the meridional flow.

In the text

	Fig. 3. Same as Fig. 2 but for runs P1-[2-5]M in set P1.
In the text

	Fig. 4. Same as Fig. 2 but for runs P10-[6-9]M in set P10.
In the text

	Fig. 5. Normalized average angular velocity at the equator $⟨\mathrm{\Delta }{\stackrel{\sim }{\mathrm{\Omega }}}_{\mathrm{eq}}⟩$ for all runs as a function of Co (a), Coℓ(b), Coω(c), and Co⋆(d). Circles (crosses) denote MHD (hydrodynamic) runs, and the colour of the symbols indicates the SGS Prandtl number. The sizes of the symbols are proportional to ReM (Re) for MHD (hydrodynamic) runs.
In the text

	Fig. 6. Total angular momentum flux 𝒯 (arrows) superimposed on the mean angular velocity profile (colour contours) from runs P1-2M and P1-4M.
In the text

	Fig. 7. Divergences of the different contributions to the angular momentum flux, Eqs. (33)–(37), from runs P1-2M (top row) and P1-4M (bottom). The rightmost panels show the divergence of the total stress, and the tildes refer to normalization by ρ0GM.
In the text

	Fig. 8. Same as Fig. 6 but showing the turbulent Reynolds stress 𝒯u.
In the text

	Fig. 9. Fourier-filtered and total Reynolds stress component ${\stackrel{\sim }{\overline{Q}}}_{r\varphi }^{\mathit{u}}$ from runs P1-2M (a) and P1-4M (b) as indicated in the legend. The shaded areas indicate error estimates according to the definition in Sect. 2.1.
In the text

	Fig. 10. Instantaneous normalized radial Reynolds stress component ${\stackrel{\sim }{\overline{Q}}}_{r\varphi }^{\mathit{u}}$ on the equatorial plane from run P01-6M (colour contours). The m ≠ 0 flows are indicated by the arrows, the width of which is proportional to the local flow amplitude.
In the text