Issue 
A&A
Volume 669, January 2023



Article Number  A98  
Number of page(s)  12  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202244395  
Published online  17 January 2023 
Transition from antisolar to solarlike differential rotation: Dependence on Prandtl number
^{1}
GeorgAugustUniversität Göttingen, Institut für Astrophysik und Geophysik, FriedrichHundPlatz 1, 37077 Göttingen, Germany
email: pkaepyl@unigoettingen.de
^{2}
KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg 12, 10691 Stockholm, Sweden
Received:
1
July
2022
Accepted:
21
November
2022
Context. Latetype stars such as the Sun rotate differentially due to the interaction of turbulent convection and rotation.
Aims. The aim of the study is to investigate the effects of the effective thermal Prandtl number, which is the ratio of kinematic viscosity to thermal diffusivity, on the transition from antisolar (slow equator, fast poles) to solarlike (fast equator, slow poles) differential rotation.
Methods. Threedimensional hydrodynamic and magnetohydrodynamic simulations in semiglobal spherical wedge geometry were used to model the convection zones of solarlike stars.
Results. The overall convective velocity amplitude increases as the Prandtl number decreases, in accordance with earlier studies. The transition from antisolar to solarlike differential rotation is insensitive to the Prandtl number for Prandtl numbers below unity, but for Prandtl numbers greater than unity, solarlike differential rotation becomes significantly harder to excite. Magnetic fields and more turbulent regimes with higher fluid and magnetic Reynolds numbers help to achieve solarlike differential rotation in neartransition cases where antisolar rotation is found in more laminar simulations. Solarlike differential rotation occurs only in cases with radially outward turbulent angular momentum transport due to the Reynolds stress at the equator. The dominant contribution to this outward transport near the equator is due to prograde propagating thermal Rossby waves.
Conclusions. The differential rotation is sensitive to the Prandtl number only for large Prandtl numbers in the parameter regime explored in this study. Magnetic fields have a greater effect on the differential rotation, although the inferred presence of a smallscale dynamo did not lead to drastically different results. The dominance of the thermal Rossby waves in the simulations is puzzling because they are not detected in the Sun. The current simulations are shown to be incompatible with the currently prevailing meanfield theory of differential rotation.
Key words: turbulence / convection / dynamo / Sun: rotation
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
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). Threedimensional 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 antisolar (AS) differential rotation with equatorial deceleration (e.g., Fan & Fang 2014; Käpylä et al. 2014; Hotta et al. 2015), whereas solarlike (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 largescale 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 lowresolution 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 smallscale dynamo. This was addressed by the recent highresolution 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 smallscale 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 setups, 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 nonrotating 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 setup is similar to the setups 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, θ_{0} < θ < π − θ_{0} in colatitude, where θ_{0} = π/12, and 0 < ϕ < π/2 in longitude. Equations of fully compressible MHD were solved,
where A is the magnetic vector potential, U is the velocity, B = ∇ × A is the magnetic field, η is the magnetic diffusivity, μ_{0} is the permeability of vacuum, J = ∇ × B/μ_{0} 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. Ω_{0} = (cos θ, −sin θ, 0)Ω_{0} is the angular velocity vector, where Ω_{0} 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_{V}Dlnp − c_{P}Dlnρ, 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
where the semicolons refer to covariant derivatives (Mitra et al. 2009). The radiative flux is given by
where K is the heat conductivity. The latter consists of two parts, K = K_{1} + K_{2}, where K_{1} = K_{1}(r) is a fixed function of height and K_{2} = K_{2}(ρ, T) is density and temperaturedependent according to Kramers opacity law (Weiss et al. 2004). The profile of K_{1} is given by
where K_{top} = 2c_{P}F_{bot}/g(R), with where L is the luminosity of the star, and where d_{K} = 0.015R. The contribution K_{2} is given by
where ρ_{0} and T_{0} 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
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^{1} (Brandenburg 2021). The code in this study employs thirdorder temporal and sixthorder spatial discretisation. Advective terms in Eqs. (1)–(4) are written as fifthorder upwinding derivatives with a sixthorder hyperdiffusive correction, where the diffusion coefficient is flowdependent; 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_{0}, a, b, ρ_{0}, T_{0}, Ω_{0}, ν, η, χ_{SGS}, the profile of K_{1}, and the value of the modified StefanBoltzmann constant σ_{SB} in the upper boundary condition , 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 nondimensional luminosity is given by
where ρ_{0} 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^{5}, where ℒ_{⊙} is the dimensionless solar luminosity. The initial stratification is determined by the nondimensional pressure scale height at the surface,
where T_{1} = T(R, t = 0).
The relative strengths of viscosity, SGS diffusion, and magnetic diffusivity are given by the SGS and magnetic Prandtl numbers,
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
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,
where Δr = 0.3R is the depth of the layer, s_{hs} is the specific entropy in a onedimensional nonconvecting 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 Schwarzschildunstable 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,
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,
The fluid and magnetic Reynolds numbers and the Péclet number are given by
respectively, where is the time and volumeaveraged rms velocity, where has been replaced by to avoid contributions from differential rotation. The inverse of the wavenumber k_{1} = 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,
where t_{0} 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 subaverages from the average over the whole time series was taken to represent the error.
Summary of the runs.
2.2. Initial and boundary conditions
Initially, the stratification is isentropic with a polytropic index n = 1.5 and ξ_{0} = 0.02, resulting in an initial density contrast of 30. The value of K_{0} 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 stressfree for the flow. At the bottom boundary, a fixed heat flux is prescribed, while at the top, a blackbody 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
The azimuthal direction is periodic for all quantities. The velocity and magnetic fields were initialized with random lowamplitude 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 KelvinHelmholtz timescale, which is estimated as
where ℰ_{th} = ∫_{V}E_{th}dV is the total thermal energy of the initial nonconvecting state, with E_{th} = ρe, and where e = c_{V}T 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 KelvinHelmholtz timescale in the current simulations is , where is the freefall 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 timeaveraging, is typically of the order of τ_{KH}; see Fig. 1a for data from run P16H. Hydrodynamic simulations were considered thermally relaxed when the volumeaveraged (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} = ρϕ, , and . 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 P16H.
Fig. 1. Time series of total energy density and the rms velocity and magnetic fields. (a) Volumeaveraged total energy density ⟨E_{tot}⟩ normalized by its value at t = 0 from run P16H. The inset shows the volumeaveraged rms value of the azimuthal flow . The red (blue) part of the curve indicates the transient (statistically stationary) part of the time series, where indicates the starting time for the time averaging. τ_{KH} indicates the KelvinHelmholtz timescale. (b) Same as above, but for the corresponding MHD run P16M. The inset shows in addition to . The black (green) parts of the data for indicate transient (statistically stationary) parts. (c) Same as panel a, but for run P013H. 
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 twophase 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 volumeaveraged total energy density ⟨E_{tot}⟩ settles to a new statistically steady state. In some cases, such as in run P16M shown in Fig. 1b, the total energy density in the MHD state still shows temporal modulation due to quasicyclic largescale magnetism.
However, in a few of hydrodynamic cases (runs P012H, P013H, and P014H), a statistically stationary state was not reached, even though these models were run several KelvinHelmholtz 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 largescale flows that develop in rotating convective systems. The mean rotation profile is given by the time and azimuthal average,
and the meridional flow is given by . In many simulations, the latitudinal profiles of are nonmonotonic such that the rotation rate has a polar jet, a maximum at midlatitudes, 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
Fig. 2. Timeaveraged rotation profiles from MHD runs P01[25]M near the ASSL transition with Pr_{SGS} = 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 according to Eq. (24) are indicated in each panel. The white arrows indicate the meridional flow. 
where θ_{eq} = π/2, and where the tildes refer to normalization by the rotation rate of the frame of reference, Ω_{0}. If (), 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,
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 is shown for all runs as a function of Co. That is, in the lowresolution 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 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 Re_{M} (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 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 ASSL transition is reduced with this definition, essentially reducing the apparent dependence on Prandtl number significantly. Here this was tested by computing according to
where U_{p} is the nonaxisymmetric 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 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 , and the Coriolis number based on this is given by
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 lowresolution runs. Similarly, a Coriolis number based on fluctuating vorticity can be defined as (e.g., Brun et al. 2022)
where ω′ = ∇ × u with . 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^{2}. 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 higherRe runs because smaller convective scales are resolved and and ω′ pick these up. Furthermore, among the five higherRe runs, the single AS model (P011Mh) 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}, , 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
where is the luminosity, ρ_{⋆} is a reference density, and u_{⋆} is an estimate of an average convective velocity. We assume ρ_{⋆} = ρ_{0}^{3} and construct a stellar Coriolis number
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 Ω_{0}, to achieve SL differential rotation. However, the ASSL 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 higherresolution runs are SL or marginally SL for Co_{⋆} = 0.46, whereas the corresponding lowerresolution MHD and hydrodynamic runs are all AS. The current higherresolution runs with Re_{M} ≈ 94…113 also have smallscale 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 higherresolution 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
where is the specific angular momentum, with ϖ = r sin θ being the lever arm, and where
is the flux of angular momentum, where and are the fluctuating velocity and magnetic field, and where 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),
Additionally, the contributions to the stress due to molecular viscosity are given by
Although the viscous flux is negligibly small in the Sun, it is typically nonnegligible 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 P12M and P14M. 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 P12M and P14M. The contributions due to the turbulent Reynolds stress and due to the meridional flow 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 higherresolution runs, but 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 P12M and P14M. 
Fig. 7. Divergences of the different contributions to the angular momentum flux, Eqs. (33)–(37), from runs P12M (top row) and P14M (bottom). The rightmost panels show the divergence of the total stress, and the tildes refer to normalization by ρ_{0}GM. 
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 P12M and P14M 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 nonzero 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 P14M in Fig. 8), whereas it is AS for negative (downward) 𝒯^{u} (run P12M). Therefore, it appears sufficient to study the equatorial radial angular momentum flux to determine the difference between AS and SL differential rotation.
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 , where azimuthal orders ranging between and were retained^{4}. These flows are used to compute the radial Reynolds stress,
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 P12M, Fig. 9a, the total Reynolds stress is negative everywhere except at the base of the CZ. Contributions from the largest scale (, = 1, 2) nonaxisymmetric motions are statistically almost identical with the total stress. A weak positive contribution around the middle of the CZ is visible for (, = 3, 5), but the Reynolds stress for (, = 1, 5) is again very similar to the total stress. In the SL runs, the largest scales (, = 1, 2) also contribute to a downward flux, whereas the main contribution to the net outward flux comes from (, = 3, 5); see Fig. 9b for run P14M. 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 largescale 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 nonlinear progradepropagating thermal Rossby waves. It is therefore somewhat questionable to talk about turbulent Reynolds stress in this context since the Busse columns are essentially largescale convective modes.
Fig. 9. Fourierfiltered and total Reynolds stress component from runs P12M (a) and P14M (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 on the equatorial plane from run P016M (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 smallscale Maxwell stress is the dominant contribution to the radial angular momentum transport. While the highresolution runs in the present study have smallscale 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), largescale Busse columns can still be seen in the deep parts of their model; see for example their Fig. 6. If such largescale 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 solarlike in terms of the rotation profile, the conundrum with the too prominent largescale structures remains.
Furthermore, in hydrodynamic meanfield theories of differential rotation (e.g., Rüdiger 1989; Kitchatinov & Rüdiger 2005; Rogachevskii & Kleeorin 2015), the turbulence models are necessarily simplified, and largescale 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 meanfields 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 meanfield 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 ASSL 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 nonrotating 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 smallscale 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 higherresolution runs, which likely also have smallscale 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 nondetection 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 meanfield theories of differential rotation and 3D simulation results is pointed out as another aspect that requires further scrutiny in the future.
Only the two most rapidly rotating runs with Pr_{SGS} = 10 have statistically significant SL differential rotation; see the 12th column of Table 1.
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/41). This work was also partially supported by a grant from the Simons Foundation.
References
 Barekat, A., & Brandenburg, A. 2014, A&A, 571, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
 Brandenburg, A., Nordlund, A., & Stein, R. F. 2000, in Geophysical and Astrophysical Convection, Contributions from a workshop sponsored by the Geophysical Turbulence Program at the National Center for Atmospheric Research, October 1995, eds. P. A. Fox, & R. M. Kerr (The Netherlands: Gordon and Breach Science Publishers), 85 [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
 Brun, A. S., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1970a, ApJ, 159, 629 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1970b, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Camisassa, M. E., & Featherstone, N. A. 2022, ApJ, 938, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282 [CrossRef] [Google Scholar]
 Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [CrossRef] [Google Scholar]
 Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016, ApJ, 819, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Stejko, A. M., Kosovichev, A. G., Smolarkiewicz, P. K., & Strugarek, A. 2022, ApJ, 940, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
 Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [Google Scholar]
 Hotta, H., & Kusano, K. 2021, Nat. Astron., 5, 1100 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51 [Google Scholar]
 Hotta, H., Rempel, M., & Yokoyama, T. 2016, Science, 351, 1427 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., Kusano, K., & Shimada, R. 2022, ApJ, 933, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2019a, A&A, 622, A195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J. 2019b, A&A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J. 2021, A&A, 655, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J. 2022, ApJ, 931, L17 [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [CrossRef] [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., Olspert, N., Warnecke, J., & Brandenburg, A. 2017, A&A, 599, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Viviani, M., Käpylä, M. J., Brandenburg, A., & Spada, F. 2019, Geophys. Astrophys. Fluid Dyn., 113, 149 [CrossRef] [Google Scholar]
 Käpylä, P. J., Gent, F. A., Olspert, N., Käpylä, M. J., & Brandenburg, A. 2020, Geophys. Astrophys. Fluid Dyn., 114, 8 [CrossRef] [Google Scholar]
 Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&A, 576, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2004, Astron. Nachr., 325, 496 [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comp. Astrophys., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Matilsky, L. I., Hindman, B. W., & Toomre, J. 2020, ApJ, 898, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
 O’Mara, B., Miesch, M. S., Featherstone, N. A., & Augustson, K. C. 2016, AdSpR, 58, 1475 [Google Scholar]
 Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
 Brandenburg, A. 2021, J. Open Source Softw., 6, 2807 [CrossRef] [Google Scholar]
 Rogachevskii, I., & Kleeorin, N. 2015, J. Plasma Phys., 81, 395810504 [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solartype Stars (Berlin: Akademie Verlag) [Google Scholar]
 Rüdiger, G., & Hollerbach, R. 2004, The Magnetic Universe: Geophysical and Astrophysical Dynamo Theory (Weinheim: WileyVCH) [CrossRef] [Google Scholar]
 Rüdiger, G., Küker, M., Käpylä, P. J., & Strassmeier, K. G. 2019, A&A, 630, A109 [EDP Sciences] [Google Scholar]
 Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Schumacher, J., & Sreenivasan, K. R. 2020, Rev. Mod. Phys., 92, 041001 [NASA ADS] [CrossRef] [Google Scholar]
 Simitev, R. D., Kosovichev, A. G., & Busse, F. H. 2015, ApJ, 810, 80 [Google Scholar]
 Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiss, A., Hillebrandt, W., Thomas, H.C., & Ritter, H. 2004, Cox and Giuli’s Principles of Stellar Structure (Cambridge: Cambridge Scientific Publishers Ltd) [Google Scholar]
All Tables
All Figures
Fig. 1. Time series of total energy density and the rms velocity and magnetic fields. (a) Volumeaveraged total energy density ⟨E_{tot}⟩ normalized by its value at t = 0 from run P16H. The inset shows the volumeaveraged rms value of the azimuthal flow . The red (blue) part of the curve indicates the transient (statistically stationary) part of the time series, where indicates the starting time for the time averaging. τ_{KH} indicates the KelvinHelmholtz timescale. (b) Same as above, but for the corresponding MHD run P16M. The inset shows in addition to . The black (green) parts of the data for indicate transient (statistically stationary) parts. (c) Same as panel a, but for run P013H. 

In the text 
Fig. 2. Timeaveraged rotation profiles from MHD runs P01[25]M near the ASSL transition with Pr_{SGS} = 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 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[25]M in set P1. 

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

In the text 
Fig. 5. Normalized average angular velocity at the equator 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 Re_{M} (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 P12M and P14M. 

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

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

In the text 
Fig. 9. Fourierfiltered and total Reynolds stress component from runs P12M (a) and P14M (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 on the equatorial plane from run P016M (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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.