Open Access
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/0004-6361/202244395
Published online 17 January 2023

© The Authors 2023

Licence Creative CommonsOpen 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 Subscribe-to-Open 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). 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 rin < r < R in radius, where rin = 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,

(1)

(2)

(3)

(4)

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 = cVDlnp − cPDlnρ, where cV and cP are the specific heat capacities in constant volume and pressure, respectively. The gas is assumed to obey the ideal gas law, p = ℛρT, where ℛ = cP − cV is the gas constant. The rate of strain tensor is given by

(5)

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

(6)

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

(7)

where Ktop = 2cPFbot/|g(R)|, with where L is the luminosity of the star, and where dK = 0.015R. The contribution K2 is given by

(8)

where ρ0 and T0 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

(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 CODE1 (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, Fbot = −(KT/∂r)|r = rin, the values of K0, a, b, ρ0, T0, Ω0, ν, η, χSGS, the profile of K1, and the value of the modified Stefan-Boltzmann constant σSB in the upper boundary condition , where Tsurf 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

(10)

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

(11)

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

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

(12)

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

(13)

where χ = K/cPρ 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,

(14)

where Δr = 0.3R is the depth of the layer, shs is the specific entropy in a one-dimensional non-convecting hydrostatic model, evaluated near the top of the domain at rs = 0.95R, and where χs is the total thermal diffusivity K/cPρ from r = rs. 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 rs. 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,

(15)

where ⟨sθϕ is the time- and horizontal average of the specific entropy, and χtot = χSGS + ⟨χθϕ is the total thermal diffusivity. Rat is always significantly smaller than Ra because χSGS ≫ χ.

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

(16)

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

(17)

respectively, where is the time- and volume-averaged rms velocity, where has been replaced by to avoid contributions from differential rotation. The inverse of the wavenumber k1 = 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,

(18)

where t0 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.

Table 1.

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 K0 is chosen such that Frad = Ftot 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

(19)

(20)

(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 PrSGS = 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

(22)

where ℰth = ∫VEthdV is the total thermal energy of the initial non-convecting state, with Eth = ρe, and where e = cVT 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 , where 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 ⟨Etot⟩ reaches a statistically steady state, where ⟨Etot⟩=⟨Eth⟩+⟨Epot⟩+⟨Ekin⟩+⟨Emag⟩ with Epot = ρϕ, , 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 P1-6H.

thumbnail 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 . 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 Kelvin-Helmholtz timescale. (b) Same as above, but for the corresponding MHD run P1-6M. 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 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 ⟨Etot⟩ 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,

(23)

and the meridional flow is given by . In many simulations, the latitudinal profiles of 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. 24. Therefore the classification of the AS and SL rotation profile is here based on the mean rotation profile at the equator

(24)

thumbnail 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 according to Eq. (24) are indicated in each panel. The white arrows indicate the meridional flow.

thumbnail Fig. 3.

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

thumbnail 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, Ω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 PrSGS 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,

(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 PrSGS; see Fig. 5a where is shown for all runs as a function of Co. That is, in the low-resolution MHD runs with PrSGS = 0.1 (1), the transition occurs around Co ≈ 1.5 (Co ≈ 1.8); see Figs. 2 and 3, whereas for PrSGS = 10 the transition occurs at an even higher Coriolis number (Co ≈ 3); see Fig. 4.

thumbnail 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 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 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 according to

(26)

where Up 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 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

(27)

where urms 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)

(28)

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 PrSGS < 1. However, in both cases, the transition for PrSGS = 10 still occurs at a higher Co and Coω than for the PrSGS = 0.1 and 1 cases2. 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 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 urms, , 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

(29)

where is the luminosity, ρ is a reference density, and u is an estimate of an average convective velocity. We assume ρ = ρ03 and construct a stellar Coriolis number

(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 PrSGS = 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 PrSGS = 1 and 0.1. Notably, the MHD runs at PrSGS = 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 ReM ≈ 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 PrSGS = 10 than for PrSGS = 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

(31)

where is the specific angular momentum, with ϖ = r sin θ being the lever arm, and where

(32)

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),

(33)

(34)

(35)

(36)

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

(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 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 higher-resolution 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.

thumbnail Fig. 6.

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

thumbnail 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.

thumbnail 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 , where azimuthal orders ranging between and were retained4. These flows are used to compute the radial Reynolds stress,

(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 (,  = 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 (,  = 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 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.

thumbnail Fig. 9.

Fourier-filtered and total Reynolds stress component 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.

thumbnail Fig. 10.

Instantaneous normalized radial Reynolds stress component 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 (PrM = 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 (PrSGS). 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 PrSGS = 10 than for PrSGS = 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 PrSGS < 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 ReM 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.


2

Only the two most rapidly rotating runs with PrSGS = 10 have statistically significant SL differential rotation; see the 12th column of Table 1.

3

Using the average density of the star, , is another option, and in that case, the definition of Co is fully determined by the stellar parameters.

4

Due to the Δϕ = π/2 azimuthal extent of the simulation domain, m′ corresponds to 4m in a full sphere.

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

  1. Barekat, A., & Brandenburg, A. 2014, A&A, 571, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brandenburg, A. 2016, ApJ, 832, 6 [Google Scholar]
  4. 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]
  5. Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [Google Scholar]
  6. Brun, A. S., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21 [NASA ADS] [CrossRef] [Google Scholar]
  7. Busse, F. H. 1970a, ApJ, 159, 629 [NASA ADS] [CrossRef] [Google Scholar]
  8. Busse, F. H. 1970b, J. Fluid Mech., 44, 441 [NASA ADS] [CrossRef] [Google Scholar]
  9. Camisassa, M. E., & Featherstone, N. A. 2022, ApJ, 938, 65 [NASA ADS] [CrossRef] [Google Scholar]
  10. Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282 [CrossRef] [Google Scholar]
  11. Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 [NASA ADS] [CrossRef] [Google Scholar]
  12. Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
  13. Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76 [CrossRef] [Google Scholar]
  14. Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
  15. Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
  16. 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]
  17. Guerrero, G., Stejko, A. M., Kosovichev, A. G., Smolarkiewicz, P. K., & Strugarek, A. 2022, ApJ, 940, 151 [NASA ADS] [CrossRef] [Google Scholar]
  18. Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928 [Google Scholar]
  19. Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Ann. Rev. Fluid Mech., 48, 191 [Google Scholar]
  20. Hotta, H., & Kusano, K. 2021, Nat. Astron., 5, 1100 [NASA ADS] [CrossRef] [Google Scholar]
  21. Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51 [Google Scholar]
  22. Hotta, H., Rempel, M., & Yokoyama, T. 2016, Science, 351, 1427 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hotta, H., Kusano, K., & Shimada, R. 2022, ApJ, 933, 199 [NASA ADS] [CrossRef] [Google Scholar]
  24. Käpylä, P. J. 2019a, A&A, 622, A195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Käpylä, P. J. 2019b, A&A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Käpylä, P. J. 2021, A&A, 655, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Käpylä, P. J. 2022, ApJ, 931, L17 [CrossRef] [Google Scholar]
  28. Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [CrossRef] [Google Scholar]
  29. Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. 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]
  31. Käpylä, P. J., Viviani, M., Käpylä, M. J., Brandenburg, A., & Spada, F. 2019, Geophys. Astrophys. Fluid Dyn., 113, 149 [CrossRef] [Google Scholar]
  32. 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]
  33. 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]
  34. Karak, B. B., Miesch, M., & Bekki, Y. 2018, Phys. Fluids, 30, 046602 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kitchatinov, L. L., & Rüdiger, G. 2004, Astron. Nachr., 325, 496 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kitchatinov, L. L., & Rüdiger, G. 2005, Astron. Nachr., 326, 379 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comp. Astrophys., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
  38. Matilsky, L. I., Hindman, B. W., & Toomre, J. 2020, ApJ, 898, 111 [NASA ADS] [CrossRef] [Google Scholar]
  39. Miesch, M. S., & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
  40. Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923 [NASA ADS] [CrossRef] [Google Scholar]
  41. O’Mara, B., Miesch, M. S., Featherstone, N. A., & Augustson, K. C. 2016, AdSpR, 58, 1475 [Google Scholar]
  42. Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
  43. Brandenburg, A. 2021, J. Open Source Softw., 6, 2807 [CrossRef] [Google Scholar]
  44. Rogachevskii, I., & Kleeorin, N. 2015, J. Plasma Phys., 81, 395810504 [Google Scholar]
  45. Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solar-type Stars (Berlin: Akademie Verlag) [Google Scholar]
  46. Rüdiger, G., & Hollerbach, R. 2004, The Magnetic Universe: Geophysical and Astrophysical Dynamo Theory (Weinheim: Wiley-VCH) [CrossRef] [Google Scholar]
  47. Rüdiger, G., Küker, M., Käpylä, P. J., & Strassmeier, K. G. 2019, A&A, 630, A109 [EDP Sciences] [Google Scholar]
  48. Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
  49. Schumacher, J., & Sreenivasan, K. R. 2020, Rev. Mod. Phys., 92, 041001 [NASA ADS] [CrossRef] [Google Scholar]
  50. Simitev, R. D., Kosovichev, A. G., & Busse, F. H. 2015, ApJ, 810, 80 [Google Scholar]
  51. Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. 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

Table 1.

Summary of the runs.

All Figures

thumbnail 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 . 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 Kelvin-Helmholtz timescale. (b) Same as above, but for the corresponding MHD run P1-6M. 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 P01-3H.

In the text
thumbnail 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 according to Eq. (24) are indicated in each panel. The white arrows indicate the meridional flow.

In the text
thumbnail Fig. 3.

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

In the text
thumbnail Fig. 4.

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

In the text
thumbnail 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 ReM (Re) for MHD (hydrodynamic) runs.

In the text
thumbnail 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
thumbnail 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
thumbnail Fig. 8.

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

In the text
thumbnail Fig. 9.

Fourier-filtered and total Reynolds stress component 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
thumbnail Fig. 10.

Instantaneous normalized radial Reynolds stress component 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

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.