Transition from anti-solar to solar-like differential rotation: Dependence on Prandtl number

(abridged) Context: Late-type 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 thermal Prandtl number on the transition from anti-solar (slow equator, fast poles) to solar-like (fast equator, slow poles) differential rotation. Methods: Three-dimensional hydrodynamic and magnetohydrodynamic simulations in semi-global spherical wedge geometry are used to model convection zones of solar-like stars. Results: The overall convective velocity amplitude increases as the Prandtl number decreases in accordance with earlier studies. The transition from anti-solar to solar-like differential rotation is insensitive to the Prandtl number for Prandtl numbers below unity but for Prandtl numbers greater than unity, solar-like differential rotation becomes significantly harder to excite. Magnetic fields and more turbulent regimes with higher fluid and magnetic Reynolds numbers help in achieving solar-like differential rotation in near-transition cases. Solar-like differential rotation occurs only in cases with radially outward angular momentum transport at the equator. The dominant contribution to such 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 the current study. Magnetic fields have a greater effect on the differential rotation, although the inferred presence of a small-scale dynamo does not lead to drastically different results in the present study. The dominance of the thermal Rossby waves in the simulations is puzzling given the non-detection in the Sun. The current simulations are shown to be incompatible with the currently prevailing mean-field theory of differential rotation.


Introduction
The interplay of turbulent convection with the overall rotation of the Sun is the primary cause of differential rotation 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 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 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 tension between large-scale velocity amplitudes in simulations in comparison to the Sun (e.g. Hanasoge et al. 2012Hanasoge et al. , 2016Schumacher & 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, however, cannot be justified based on physical grounds 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: is is conceivable that sufficiently strong fields can suppress convection 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 magnetic field while Fan & Fang (2014) and Simitev et al. (2015) reported more positive outcomes. Nevertheless, these simulations most probably did not have high enough magnetic Reynolds numbers to excite a small-scale dynamo. This was addressed by 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 small-scale dynamo.
Another important parameter is the Prandtl number, Pr = ν/χ, where ν is the kinematic viscosity and χ is the thermal diffusivity. A notion that the solar convection zone is operating in a regime where the effective Prandtl number is large, has gained popularity recently (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 such set-ups, the prob-

The model
The simulation set-up is similar to those 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 are 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, 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 D ln p − c P D ln ρ, 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 = RρT , where R = 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 temperature-dependent according to Kramers opacity law (Weiss et al. 2004). The profile of K 1 is given by where K top = 2 3 F bot , with F bot = L/4πr 2 in 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 (Pencil Code Collaboration et al. 2021). In the present study the code employs third-order temporal and sixth-order spatial discretisation. Advective terms in Eqs. (1) to (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).

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, and the value of the modified Stefan-Boltzmann constant σ SB in the upper boundary condition σ SB T 4 surf = K∂T /∂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 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 = L/L ≈ 2.1 · 10 5 , where L is the dimensionless solar luminosity. The initial stratification is determined by the non-dimensional 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 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: 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 rotation is controlled by the Taylor number The fluid and magnetic Reynolds numbers and the Péclet number are given by respectively, where u rms = 3 2 (U 2 r + U 2 θ ) is the time-and volume averaged rms velocity where U 2 φ has been replaced by (U 2 r + U 2 θ )/2 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.1.
Mean quantities are denoted by overbars 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 are obtained by dividing the time series in three parts and computing averages over each one of them. The largest deviation of these sub-averages from the average over the whole time series is taken to represent the error.

Initial and boundary conditions
Initially the stratification is isentropic with 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 impenetrable and stress-free for the flow. On the bottom boundary, a fixed heat flux is prescribed while at the top a black body condition is applied. On 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. On 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 are initialized with random lowamplitude Gaussian noise fluctuations.

Results
Three sets of simulations were done 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 were remeshed to higher resolution and correspondingly higher Rayleigh, Péclet, and Reynolds numbers; see Table 1.
Only MHD variants of the P10 runs were run.

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 time-and azimuthal average: and the meridional flow is given by U mer = (U r , U θ , 0). In many simulations the latitudinal profiles of Ω are non-monotonic such that the rotation rate has a polar jet, a maximum at midlatitudes or sometimes several local minima and maxima as 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 Figures 1 to 3. Therefore the classification of AS and SL rotation profile is here based on the mean rotation profile at the equator where θ eq = π/2, and where the tildes refer to normalization by the rotation rate of the frame of reference, Ω 0 . If ∆Ω eq > 0 ( ∆Ω eq < 0) the run is classified as SL (AS) rotator. This  (29). PrM = 1 in all runs except in P01-4M and P01-6M where PrM = 1.67, and in P01-5M where PrM = 1.43 to obtain growing dynamos. The last column denotes the type of differential rotation with parenthesis indicating that the result is not statistically significant.
measure turns out to be a monotonic function of rotation and it is furthermore unaffected by equatorial asymmetries or latitudinal jets.
The current results indicate that the convective velocity increases when the Prandtl number is decreased. This is manifested by increasing fluid Reynolds number for decreasing SGS Prandtl number; see the eighth column of Table 1. Naively one could then expect that achieving SL differential rotation for low Pr SGS would be more difficult, that is, require faster rotation. Often the rotational influence on the flow is quantified by a simple definition of the Coriolis number where the convective length scale is assumed to be unchanged by rotation. Using this definition to characterize the results, the Coriolis number at which the rotation profile flips from AS to SL appears to decrease monotonically as Pr SGS decreases; see Fig. 4(a) where ∆Ω 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 Figures 1 and 2, whereas for Pr SGS = 10 the transition occurs at an even higher Coriolis number (Co ≈ 3); see Fig. 3.
However, the validity of this simplistic definition of the Coriolis number to characterize 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 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 when using this definition, essentially reducing the apparent dependence on Prandtl number significantly. Here this was tested by computing u according to where U p is the non-axisymmetric poloidal flow and the superscript refers to the corresponding 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 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 = π∆R/ u , and the Coriolis number based on this is given by where u rms is defined 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 per cent smaller 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), where ω = ∇ × u with u = U − 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 where 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, one 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 Figures 4(b) and (c), respectively. Ignoring the five runs at higher Reynolds and Péclet numbers 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 occurs still at a higher Co and Co ω than for the Pr SGS = 0.1 and 1 cases 2 . These re-2 Note that only the two most rapidly rotating runs with PrSGS = 10 have statistically significant SL differential rotation; see the 12th column of Table 1. sults are in accordance 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 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 horizontal error estimates 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 Figures 4(b) 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 , 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 where L = 4πr 2 in F bot is the luminosity, ρ is a reference density, and u is an estimate of an average convective velocity. Here 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 in determining 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. 4(d). 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 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 at Co = 0.46 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 are in accordance 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, if 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. (24), gives misleading results and should be avoided. Contrary to the SGS Prandtl number, the dependence on magnetic fields is clearer such that in MHD runs SL differential rotation is easier to excite.

Angular momentum transport
When contributions from molecular viscosity can be neglected, the angular momentum in the interior of the star is governed by the conservation equation where L = r 2 sin 2 θΩ is the specific angular momentum and is the total flux of angular momentum. 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 where b = B − B is the fluctuating magnetic field. Representative results of T are shown in Fig. 5 superimposed on top of the angular velocity in Runs P1-[2-4]M and P01-3Mh. The total angular momentum flux 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, T is directed predominantly equatorward irrespective 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 T at the equator: for positive (outward) flux the differential rotation is SL (Runs P1-4M and P01-3Mh in Fig. 5) whereas it is AS for negative (downward) flux (Run P1-2M). In the transitory case of Run P1-3M the flux at the equator converges at the local maxima of Ω similarly as in the SL cases of P1-4M and P01-3Mh. Therefore it appears sufficient to study the radial angular momentum flux to determine the difference between AS and SL differential rotation. Figure 6 shows all of the components of the radial angular momentum flux for runs P1-2M (P1-4M) with AS (SL) differential rotation. In both cases the Reynolds stress due to m = 0 contributions of the velocity is the dominant contributors while the Reynolds stress due to meridional circulation and Maxwell stresses are clearly subdominant. Moreover, in the AS case P1-2M the radial transport due to the Reynolds stresses is downward almost everywhere, while in the SL case P1-4M, both Q u rφ and Q U rφ are positive near the equator. The sign of the latter near the equator is determined by the sign of U φ because U r > 0 in all cases; see the meridional flow, for example, in Fig. 1.
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 U (36) Here the density fluctuations are assumed to be small and no Fourier filtering was applied to ρ. Representative results are shown in Fig. 7. In the AS case P1-2M, Fig. 7(a), the total Reynolds stress is negative everywhere except at the very base of the CZ. Contributions from the largest scale (m 1 , m 2 = 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 , m 2 = 3, 5) but the Reynolds stress for (m 1 , m 2 = 1, 5) is again very similar to the total stress. In the SL runs the largest scales (m 1 , m 2 = 1, 2) also contribute to a downward flux whereas the main contribution to the net outward flux comes from (m 1 , m 2 = 3, 5); see Fig. 7(b) for Run P1-4M. Similarly to the AS case, the contributions from m > 5 are very small. This indicates that practically all of the Reynolds stress at the equator is due to relatively large-scale structures which can be identified as the 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. 8 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 large-scale convective modes that appear essentially at a scale corresponding to forcing of turbulence.
The mechanism by which the differential rotation is generated in the current simulations is therefore different from that in Hotta et al. (2022) and Hotta et al. (2022) where the small-scale Maxwell stress that 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 modeled stars were rotating typically 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 no evidence currently 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 the large-scale convection modes such as the Busse columns do not appear. In the most commonly adopted approach of Kitchatinov & Rüdiger (2005), the radial angular momentum transport is downward for slow, and vanishing for rapid rotation at the equator. The SL differential rotation results in from a strong equatorward transport. Numerical simulations of isothermal homogeneous anisotropic turbulence also produce downward (slow rotation) or vanishing (rapid ro- tation) 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 meanfield theories avoid the problem of too prominent too strong thermal Rossby waves by simply neglecting them, whereas they fail to characterize both AS and SL cases in the current simulations. This tension 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.

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 also the cases Pr = 1 and Pr = 0.1 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 because 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 et al. 2016;Käpylä 2022) which also couple back to differential rotation. These aspects need to be revisited in future studies.
Many of the current simulations also included magnetic fields, albeit often in a parameter regime where 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 that likely also have small-scale dynamos. However, the effects of magnetic fields are likely to be more significant in more realistic higher-Re M systems as manifested by 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, which are still prominent in all current simulations, from the Sun still raises questions as to the generation mechanism of solar differential rotation. Finally, the tension between mean-field theories of differential rotation and 3D simulation results is pointed out as another aspect that requires further scrutiny in the future.