© M. Viviani and M. J. Käpylä 2021

1. Introduction

The solar surface differential rotation has been known for a long time (Scheiner 1630; Carrington 1863): the equator completes a turn in around 25 days, while the poles take roughly 30 days. Helioseismic inferences also allowed the subsurface rotation to be uncovered (Schou et al. 1998), and revealed that the lines of constant angular velocity are radial. This was somewhat unexpected, as in a uniform, incompressible flow, the Taylor–Proudman theorem (Chandrasekhar 1961) states that the horizontal components of the velocity field cannot vary in the direction of the rotation axis, and the flow is forced to move in vertical columns, in which case constant angular-velocity contours on cylinders would be observed. Therefore, the Sun is able to break the Taylor–Proudman balance by some means. Another surprising observational result came from time-distance helioseismology (Hanasoge et al. 2012), which revealed a lack of power in the kinetic energy spectrum at large scales, where the peak for giant cells should be located. Such a peak would be expected from mixing-length theory (MLT; Vitense 1953; Böhm-Vitense 1958): in its original formulation, MLT predicts convection at all possible scales, which would also correspond to cells of the diameter of the entire convective layer. Also, more recent measurements (e.g., Rincon et al. 2017) suggest that supergranulation may indeed be the largest scale excited in the Sun.

Theoretical explanations as to how the Sun breaks the Taylor–Proudman balance include a “thermal wind”, generating a clockwise meridional circulation pattern. This circulation results from a latitudinal temperature gradient, which is such that the pole is warmer than the equator by a small amount, namely of the order of a few Kelvin (Rüdiger 1989). This temperature difference is comparable to the error of current instruments, although Rast et al. (2008) reported an enhancement of ∼2.5 K at the Sun’s poles. One possible theoretical explanation of such a temperature gradient evokes the importance of turbulent effects, such as latitudinally anisotropic heat flux, which has been shown to be able to lead to a temperature difference of ∼4 K (Kitchatinov & Rüdiger 1995). Also, the presence of a weakly subadiabatic layer at the base of the convection zone has been shown to generate a thermal wind and sustain the necessary temperature gradient in a mean-field hydrodynamic model (Rempel 2005). Furthermore, small-scale dynamo action (see, e.g., Kazantsev 1968) could be the reason for the Sun not to be in Taylor–Proudman balance, as has been shown recently by Hotta (2018).

Modeling efforts of stellar convection in spherical or semi-spherical shells still struggle to produce solutions in which the Taylor–Proudman balance is self-consistently broken, and thus still tend to show cylindrical isocontours for the differential rotation (e.g., Guerrero et al. 2013; Gastine et al. 2014; Käpylä et al. 2014; Augustson et al. 2015). This is commonly interpreted as suggesting too strong a rotational influence. Moreover, such models are unable to reproduce an accelerated equator when using the solar rotation rate (Gastine et al. 2014; Käpylä et al. 2014; Karak et al. 2015). Most present-day numerical setups use fixed heat conduction profiles and depths of convection zones (e.g., Brun et al. 2011; Käpylä et al. 2013) motivated by MLT. Although MLT-designed setups have been successful in reproducing the pattern of granulation and supergranulation in surface convection (e.g., Brandenburg et al. 2005; Nordlund et al. 2009), models simulating deeper parts of the convection zone (CZ) produce far more power in the velocity spectrum at large scale than measured in observations (Gizon & Birch 2012). All of the above-mentioned discrepancies between observations and numerical models are collectively known as the “convective conundrum” (see, e.g., O’Mara et al. 2016) and solving it is one of the major challenges of contemporary solar physics.

One proposed way to crack the convective conundrum is to hypothesize that the actual convectively unstable layer in the Sun, according to the Schwarzschild criterion (Chandrasekhar 1961), is shallower than expected. Spruit (1997) described convection as being driven by cool threads descending from the surface into deeper layers, overwhelming the convection driven by heating from below. Such a phenomenon is now denoted as entropy rain (Brandenburg 2016) and describes surface-driven convection that would excite only small to medium length scales. By extending MLT to include entropy fluctuations, Brandenburg (2016) identified the presence of a Schwarzschild stable, sub-adiabatic layer in which the convective flux is still positive. Such a layer was first identified in the Earth’s atmosphere (Deardorff 1961, 1966), and was therefore termed the “Deardorff layer”.

Although recent studies (e.g., Hotta et al. 2019) have not found evidence to support their existence, the formation of such sub-adiabatic layers was reported in the hydrodynamic studies of for example Hotta (2017), Korre et al. (2017), and Käpylä et al. (2017). The study by Käpylä et al. (2017) is especially relevant here, as these authors demonstrated the emergence of a substantial sub-adiabatic layer and the existence of nonlocal surface-driven convection using Kramers opacity law in a Cartesian model. Consequently, they redefined the convection zone as the sum of the convection zone in the traditional sense –now called the buoyancy zone– plus the sub-adiabatic part –now denoted the Deardorff layer. The depth of the layers was not determined a priori, but was rather an outcome of the simulations.

However, these studies did not investigate large-scale dynamo action. The effect of sub-adiabatic layers in global MHD simulations was investigated by Käpylä et al. (2019). The formation of a stably stratified layer at the bottom of the domain allowed for the storage of magnetic field beneath it, also found in an earlier study by Browning et al. (2006), but these strong fields were also seen to be capable of suppressing the oscillating magnetic field at the surface. Käpylä et al. (2019) also considered the effect of sub-adiabatic layers on the convective velocity spectra, but found that the decrease in power at large scales was not enough to solve this part of the conundrum.

Another mechanism to reduce the overly high rotational influence on convection in simulations, studied first in a Cartesian model by Hotta et al. (2015) and then in fully spherical models by Karak et al. (2018), could be provided by the Lorentz force feedback from the magnetic to the velocity field. Such feedback could result from strong magnetic fluctuations, originating for example from the action of a small-scale dynamo instability operating in the CZ. Thus-generated magnetic fluctuations could suppress the turbulent velocity field through the Lorentz force, hence acting as an enhanced viscosity, and increasing the magnetic Prandtl number, that is, the ratio of the viscosity and resistivity of the fluid. Karak et al. (2018) investigated such a situation numerically; their simulations developed an overshoot zone at the base of the domain, and also showed a decrease in the convective power at large scales due to downward-directed plumes. These results, although arising for a different reason, are consistent with the results of Käpylä et al. (2017, 2019). Another finding of Karak et al. (2018) was that the plumes, carrying their angular momentum inward, caused the rotation profile to switch to anti-solar.

Observations of rapidly rotating stars, younger and more active than the Sun, indicate concentrations of magnetic activity at high latitudes persisting for a long time (e.g., Berdyugina & Tuominen 1998). A common configuration is two activity patches on two “active” longitudes, separated by roughly 180° in longitude (e.g., Jetsu 1996). Active longitudes usually migrate in the rotational frame of the star, forming azimuthal dynamo waves (ADW; Berdyugina & Tuominen 1998). The direction of migration of these structures can follow the plasma rotation, and in this case we discuss prograde ADWs; these can also drift in the opposite direction (retrograde ADWs), or can appear stationary to the observer (standing ADWs). These ADWs can persist for time-spans extending to ten years (e.g., Lindborg et al. 2011), or their appearance can be more erratic (e.g., Olspert et al. 2015), with a short-lived ADW reappearing after some time. Lehtinen et al. (2016) and See et al. (2016) reported on a threshold in activity above which stars show active longitudes. In the study of Lehtinen et al. (2016), the active longitudes found were mostly migrating in the prograde direction. The appearance of active longitudes has been attributed to nonaxisymmetric dynamo modes operating in these stars (Tuominen et al. 2002), in contrast to the axisymmetric dynamo operating in less active stars. The transition from nonaxisymmetric to axisymmetric dynamos has also been studied numerically (Cole et al. 2014; Viviani et al. 2018), but these latter studies reported a majority of retrograde ADWs, and a transition from axi- to nonaxisymmetric solutions at overly low rotation rates in comparison to observations. Both studies used prescribed and MLT-motivated profiles for heat conduction, resulting in a priori fixed depth of the CZ.

The aim of this paper is to extend the study of Viviani et al. (2018) to include a dynamically adaptable heat conduction. In order to do this, we use a Kramers-like opacity law, as was done in Käpylä et al. (2019) for semi-spherical wedge simulations. We use computational domains extending over the full longitudinal extent in order to be able to study both axi- and nonaxisymmetric dynamo solutions.

2. Setup and model

We apply a similar setup as in Käpylä et al. (2013, 2019), representing the outer envelope of a solar-like star, 0.7R ≤ r ≤ R (with R the radius of the star), in a semi-spherical domain, 0 ≤ ϕ ≤ 2π and θ₀ ≤ θ ≤ π − θ₀ (θ₀ = 15°). We solve the system of MHD equations numerically:

$$\begin{array}{cc}& \frac{\mathrm{D}ln\rho }{\mathrm{D}t}+\mathrm{\nabla }·\mathit{U}=0,\hfill \\ \hfill & \frac{\mathrm{D}\mathit{U}}{\mathrm{D}t}=\mathit{g}-2{\mathbf{\Omega }}_0\times \mathit{U}+\frac{1}{\rho }(\mathit{J}\times \mathit{B}-\mathrm{\nabla }p+\mathrm{\nabla }·2\nu \rho \mathsf{S}),\hfill \\ \hfill & T\frac{\mathrm{D}s}{\mathrm{D}t}=\frac{1}{\rho }[-\mathrm{\nabla }·({\mathit{F}}^{\mathrm{rad}}+{\mathit{F}}^{\mathrm{SGS}})+{\mu }_0\eta {\mathit{J}}^2-{\mathrm{\Gamma }}_{\mathrm{cool}}]+2\nu {\mathsf{S}}^2,\hfill \\ \hfill & \frac{\partial \mathit{A}}{\partial t}=\mathit{U}\times \mathit{B}-{\mu }_0\eta \mathit{J},\hfill \end{array}$$(1)

where ρ and U are the density and the velocity field, g = −GM/r³ is the gravitational acceleration, with G being the gravitational constant and M the mass of the star, Ω₀ = Ω₀(cos θ,− sin θ,0) is the bulk rotation, J, B, and A are the electric current, the magnetic field, and the vector potential, respectively, p, ν, and μ₀ are the pressure, the viscosity, and the magnetic permeability in vacuum, while η is the magnetic diffusivity. Here, S is the rate-of-strain tensor, and F^rad and F^SGS are the radiative and sub-grid scale (SGS) fluxes, expressed by:

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{rad}}=-K\mathrm{\nabla }T\mathrm{and}{\mathit{F}}^{\mathrm{SGS}}=-{\chi }_{\mathrm{SGS}}\rho T\mathrm{\nabla }s\prime ,\end{array}$$(2)

where K is the radiative heat conductivity, χ_SGS is the SGS heat diffusivity, assumed to be constant, and s′ is the fluctuating entropy, $s\prime =s-{\langle \overline{s}\rangle }_{\theta }$, where the overbar denotes longitudinal average and the brackets express averaging over the variable in the subscript. Finally, Γ_cool is a term acting near the surface and cooling towards a reference temperature, and its flux is expressed as:

$$\begin{array}{c}\hfill {\overline{F}}_{\mathrm{cool}}={\displaystyle {\int }_{0.7R}^R{\mathrm{\Gamma }}_{\mathrm{cool}}\mathrm{d}r.}\end{array}$$(3)

In our previous studies with fixed K-profile in spherical geometry (e.g., Käpylä et al. 2013) and in the Cartesian study adopting Kramers opacity law, Käpylä et al. (2017), we employed a SGS scheme that acted on the mean entropy gradient rather than on its fluctuations in order to take the role of the cooling layer applied here, namely to transport the flux near the surface. In addition to that, Käpylä et al. (2017) introduced an additional SGS scheme that acts on the entropy fluctuations to ensure numerical stability in the bulk of the convection zone, where the thermal diffusivity χ = K/(ρc_P) varies over several orders of magnitude when Kramers opacity law is used. This same scheme was used in Käpylä et al. (2019) and is also adopted here.

The initial velocity and magnetic fields are Gaussian seeds. The initial stratification is isentropic. The radiative heat conductivity, K, follows from Kramers opacity law for free–free and bound–free transitions (used also in Barekat & Brandenburg 2014; Käpylä et al. 2017, 2019, 2020; Käpylä 2019):

$$\begin{array}{c}\hfill K=K_0{\left(\frac{\rho }{{\rho }_0}\right)}^{-2}{\left(\frac{T}{T_0}\right)}^{13/2},\end{array}$$(4)

where ρ₀ and T₀ are reference values for density and temperature taken at the bottom of the convection zone. The constant K₀ is defined via:

$$\begin{array}{cc}& K_0=\frac{\mathcal{L}}{4\pi }{c}_V(\gamma -1)(n_{\mathrm{ad}}+1){\rho }_0\sqrt{\mathit{GMR}},\hfill \\ & \mathcal{L}=\frac{L_0}{{\rho }_0{\left(GM\right)}^{3/2}{R}^{1/2}},\hfill \end{array}$$(5)

where ℒ is the normalized luminosity, c_V is the specific heat at constant volume, γ = c_P/c_V is the ratio between the specific heat at constant pressure and volume, and n_ad = 1.5 is the adiabatic index.

The velocity field is impenetrable and stress free at all boundaries, while entropy derivatives are set to zero at θ = θ₀ and θ = π − θ₀. The magnetic field is radial at r = R and a perfect conductor boundary condition is applied at the bottom boundary. At the latitudinal boundaries, B is tangential, which means, in terms of the vector potential:

$$\begin{array}{c}\hfill A_{\mathrm{r}}=A_{\varphi }=\frac{\partial A_{\theta }}{\partial \theta }=0\mathrm{at}\theta ={\theta }_0\mathrm{and}\theta =\pi -{\theta }_0\end{array}$$(6)

Käpylä et al. (2020) showed that this latitudinal boundary condition does not generate major differences with respect to the perfect conductor boundary condition used in previous works (e.g., in Käpylä et al. 2013; Cole et al. 2014; Warnecke et al. 2014; Viviani et al. 2018).

The simulations are defined by the parameters Ω₀, ν, η, χ_SGS, K₀, ρ₀, T₀, and the energy flux at the bottom, F_bot = −K∂_rT|_r = 0.7R. As we have discussed at length in our previous papers (for an extended discussion, see Käpylä et al. 2013), the Rayleigh number in our simulations is much lower than in real stars, which means that the energy fluxes and Mach number are higher than in reality. Also the temperature and density fluctuations are several orders of magnitude larger than in the Sun. For these reasons, the pre-factor in the Kramers opacity law and how it acts on the heat conduction are all overestimated in the model. Furthermore, all the other values of diffusive quantities that are used here are enhanced from their real counterparts by orders of magnitude, to be able to run the numerical model on the finite-sized grid.

Moreover, important nondimensional input parameters are the magnetic and SGS Prandtl numbers:

$$\begin{array}{c}\hfill {\mathrm{Pr}}_{\mathrm{M}}=\frac{\nu }{\eta },{\mathrm{Pr}}_{\mathrm{SGS}}=\frac{\nu }{{\chi }_{\mathrm{SGS}}}·\end{array}$$(7)

Output parameters of the simulations are the fluid and magnetic Reynolds numbers:

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_{\mathrm{f}}},\mathrm{Rm}=\frac{u_{\mathrm{rms}}}{\eta k_{\mathrm{f}}}\end{array}$$(8)

with $u_{\text{rms}}=\sqrt{3/2{\langle U_{\mathrm{r}}^2+U_{\theta }^2\rangle }_{r\theta \varphi t}}$ the rms velocity and k_f = 2π/0.3R the wave number of the largest eddy, corresponding to the radial extent, and the Coriolis number:

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

quantifying the relative importance of rotation and convection.

Physical units are chosen using the solar radius R = 7 × 10⁸ m, the solar angular velocity Ω_⊙ = 2.7 × 10⁻⁶ s⁻¹, the density at the bottom of the solar convection zone ρ_bot = 200 Kg m⁻³, and the magnetic permeability in vacuum μ₀ = 4π × 10⁻⁷ H m⁻¹. We performed our simulations using the PENCIL CODE (Brandenburg et al. 2020)¹, a high-order, finite-difference, open source code for solving the MHD equations.

3. Results

The simulations and their defining parameters are summarized in Table 1. Runs R1 and R2 correspond to Run C3 and D of Viviani et al. (2018), where the radiative heat conductivity K was only a function of depth, as described in Käpylä et al. (2013). Here, our aim is to study the effect of the more physical treatment of heat conduction on the anti-solar-to-solar-like differential rotation transition and the transition to nonaxisymmetric magnetic fields. Run C3 was the simulation with the slowest rotation rate showing both accelerated equator and nonaxisymmetric magnetic field, and is therefore a good choice for this study. Run D had a rotation rate 2.1 times the solar value and exhibited a nonoscillatory dynamo solution with dominance of the m = 1 Fourier mode with a retrograde ADW. Run C3 behaved otherwise similarly, but the dynamo solution was oscillatory. This difference between the runs was most likely connected to the higher amount of differential rotation in C3 than in D; see also Appendix A, where we reproduce the rotation profiles of these runs. Run R3 is the extension of the three-times-solar rotation rate Run MHD2 of Käpylä et al. (2019) over the full longitudinal extent. Run MHD2 was a wedge simulation covering one-quarter of the full longitude, thus not allowing for nonaxisymmetric solutions to develop. We repeat this run in an extended azimuthal domain to study the possible topological changes of the magnetic field. Run R4 has the same setup as Run R3, but twice the rotation rate. Simulations in the same rotation range, but with fixed heat conduction profiles (e.g., Viviani et al. 2018), all showed a clear dominance of the nonaxisymmetric mode over the axisymmetric one, and ADWs with retrograde migration. For direct comparison, we keep the same numerical resolution as in the studies mentioned above. Such setups were shown to adequately describe magnetoconvection in the rotation-rate range explored here (see Viviani et al. 2018).

3.1. Convection zone structure

We define the convection zone according to the revised structure proposed by Käpylä et al. (2019), and indicate the bottom of the different layers in Figs. 1, 2 and A.1. The radial enthalpy flux is defined as $F_{\mathrm{r}}^{\mathrm{enth}}=\overline{c_{\mathrm{p}}\left(\rho u_{\mathrm{r}}\right)\prime T}$. The bottom of the buoyancy zone (BZ), in which the radial enthalpy flux is greater than zero, $F_{\mathrm{r}}^{\mathrm{enth}}>0$, and the radial entropy gradient is negative, ∂_rs <  0, is indicated with a continuous green line; we note that our BZ would be the convection zone if defined based on the Schwarzchild criterion. We denote the bottom of the Deardorff layer (DZ), in which $F_{\mathrm{r}}^{\mathrm{enth}}>0$ and ∂_rs >  0, by a dashed line, and the bottom of the overshoot zone (OZ), for which $F_{\mathrm{r}}^{\mathrm{enth}}<0$ and ∂_rs >  0, with a dash-dotted line. What we denote the convection zone is the combination of the BZ and DZ, where enthalpy flux is positive, but entropy gradient can be also positive, meaning that the DZ part of our convection zone is sub-adiabatic. In the radiative zone (RZ), $F_{\mathrm{r}}^{\mathrm{enth}}\approx 0$ and ∂_rs >  0. The values averaged over latitude and longitude for the depth of the layers are also shown in Table 1. We quote two Coriolis numbers for each simulation. The first one is obtained from Eq. (9); the second one, denoted as Co_rev in Table 1, takes into consideration the wave number of the revised convection zone (BZ and DZ); therefore, we use k_frev = 2π/(R−r_DZ) (see, also, Käpylä et al. 2019), where r_DZ is the latitudinally averaged radius of the Deardorff layer reported in Table 1.

Fig. 1. Radial Lenth normalized by L0. The arrows show the direction of $L_r^{\text{enth}}$. The continuous, dashed, and dash-dotted green lines represent, respectively, the bottom boundaries of BZ, DZ, and OZ.

Fig. 2. Differential rotation profiles. Continuous, dashed, and dash-dotted green lines are as in Fig. 1.

Run R1 has the deepest BZ of all runs, a very thin DZ and a considerable OZ. A thin RZ develops at the bottom. Run R2 has the thinnest convection zone of all runs. The OZ is also thin, and therefore the run develops a very thick RZ. In Runs R1 and R2, the thickness of the layers does not change considerably as a function of latitude, except for a slight tendency of the DZ to become thicker near the equatorial region for Run R2. For Run R3, the convection zone structure at higher latitudes again resembles that of Run R1. However, in the equatorial region the convection zone becomes very deep. Close to the tangent cylinder, the BZ becomes considerably shallower and the DZ develops a “bulge” in that region. Run R4 also exhibits a convection zone structure that varies strongly with latitude, and closely resembles the one seen in Run R3. Also, a hemispheric asymmetry develops in Run R4: the DZ “bulge” is larger and the BZ is deeper in the lower hemisphere than in the upper one.

We indicate the value of the density contrast Γ = ρ(R)/ρ(0.7R) in the last column of Table 1. Run R2 has the largest stratification, Γ = 52. These values are very different from those of the real density stratification in the Sun and stars (for the Sun, Γ ≈ 10⁶, see, e.g., Christensen-Dalsgaard et al. 1996).

3.1.1. Enthalpy flux

We inspect the radial enthalpy flux, $F_{\mathrm{r}}^{\mathrm{enth}}$, by representing the enthalpy luminosity, $L_r^{\mathrm{enth}}=4\pi r^2{F}_{\mathrm{r}}^{\mathrm{enth}}$, in Fig. 1 with black arrows. The enthalpy flux in Runs R1 and R2 is isotropic in latitude and rather radial everywhere in the BZ. There is a slight tendency of the flux being enhanced or diminished in the equatorial region for Runs R1 and R2, respectively. A weak negative flux in the equatorial region is present in the OZ for Run R1. A different situation arises for the two more rapidly rotating runs, where the convective transport of energy is stronger at low latitudes. Especially for Run R3, there is a decrease in the enthalpy flux in the regions of the tangent cylinder. A clear equatorial asymmetry is present in $L_{\mathrm{r}}^{\mathrm{enth}}$ for Run R4. This asymmetry is also reflected in the convection zone structure, as discussed above.

3.1.2. Differential rotation

The last two columns in Table 1 quantify the relative radial and latitudinal differential rotation, defined as:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{r}}=\frac{{\mathrm{\Omega }}_{\mathrm{eq}}-{\mathrm{\Omega }}_{\mathrm{bot}}}{{\mathrm{\Omega }}_{\mathrm{eq}}}\mathrm{and}{\mathrm{\Delta }}_{\theta }=\frac{{\mathrm{\Omega }}_{\mathrm{eq}}-{\mathrm{\Omega }}_{\mathrm{pole}}}{{\mathrm{\Omega }}_{\mathrm{eq}}}·\end{array}$$(10)

Here, Ω_eq is the surface rotation rate at the equator, Ω_bot the equatorial rotation rate at the bottom of the simulation domain, and Ω_pole = (Ω(R,θ₀)+Ω(R,π−θ₀))/2. We show the differential rotation profiles in Fig. 2 and compare with the corresponding simulations from other works, showing their profiles in Appendix A.

Run R1 corresponds to the simulation with the lowest rotation rate, showing an accelerated equator in Viviani et al. (2018), the rotation profile in that case appearing quite solar-like (see Fig. A.1, first panel). With an adaptable heat conduction prescription, the rotation profile is less solar-like (see, Fig. 2, left panel and Table 1): the equatorial acceleration becomes less pronounced, the angular velocity contours are more cylindrical, and additional regions of negative shear appear at mid-latitudes and in the equatorial region close to the surface. Such regions of negative shear at mid-latitudes, which appear in many simulations (e.g., Käpylä et al. 2012; Warnecke et al. 2014), have been identified as responsible for the equatorward propagation of the magnetic field at the surface. However, in this case we measure weaker relative differential rotation and also the opposite sign in comparison to the results of Viviani et al. (2018), where Δ_r = 0.07 and Δ_θ = 0.17 for Run C3. Hence, our model is even closer to the anti-solar-to-solar-like differential rotation transition than in Run C3 of these latter authors. Therefore, the contribution of the differential rotation to the large-scale dynamo should be negligible. Nevertheless, to confirm this hypothesis, a more thorough analysis would be required.

In comparison to Run D in Viviani et al. (2018) (Fig. A.1, second panel), which had weak values for the relative differential rotation (Δ_r = 0.003 and Δ_θ = 0.007), Run R2 presents stronger DR in absolute value, although with the opposite sign. The rotation profile is anti-solar with a retrograde flow at the equator. At mid-latitudes, a region of accelerated flow develops and the isorotation contours here and at higher latitudes are radially inclined. In the thick RZ, the rotation does not vary significantly in latitude or depth.

The rotation profile of Run R3 (Fig. 2) is very similar to the one from its wedge counterpart shown in the third panel of Fig. A.1; it is solar-like, showing an accelerated equator, and has a rather weak relative differential rotation in terms of Δ_r and Δ_θ. The minimum at mid-latitudes is present, and its location corresponds to the sub-adiabatic region at the top boundary, which is probably numerical in nature. A near-surface shear layer, with a negative radial gradient, is present from mid to low latitudes.

The rotation rate of Run R4 is similar to that of Run H in Viviani et al. (2018), while its Coriolis number is close to Run G^a in the same study, see also the last two panels in Fig. A.1. Hence, the use of the Kramers opacity law produces higher convective velocities and therefore smaller Co. The values for Δ_r and Δ_θ coincide with those from G^a. The rotation profile shown in the rightmost panel of Fig. 2 closely resembles that of Run H, with a deep minimum at mid-latitudes.

3.2. Dynamo solutions

All the presented runs develop large-scale dynamo (LSD) action, and thereby sustain a magnetic field. However, the magnetic Reynolds numbers are too low to allow for small-scale dynamo action (SSD). Our dynamo solutions therefore not only exhibit magnetic fields on the largest scales, but also a strong fluctuating component, that is generated by tangling of the LSD-generated magnetic field by the turbulent motions rather than from an SSD. We present the results of the decomposition of the magnetic field in the first 11 spherical harmonics (0 ≤ l, m ≤ 10) in Table 2. The decomposition was performed on the radial magnetic field component near the surface of the simulation (r = 0.98R). Although not always showing polarity reversals, all the runs show variation of the magnetic energy over time. Such variations would be akin to intensity variations in stellar light curves. Tracing these variations over time for our runs could be a possible measure for the magnetic activity. For each of the runs, we calculate the characteristic time of variation of the magnetic field, τ_cyc, from the time evolution of the dominant dynamo mode. The results are shown in the last column of Table 2, and the mode from which τ_cyc is calculated is indicated as a subscript. As described in Viviani et al. (2018), these cycles are at most quasi-periodic, and therefore Fourier analysis is not suitable here. Instead, we use a syntactic method, which means that we count how many times the dominant mode of the magnetic field peaks above its mean value, and τ_cyc is obtained by dividing the length of the full time span of measurement by the number of peak times.

Runs R1 and R2 have a dominant axisymmetric large-scale magnetic field, but also a significant contribution from the small-scale field (l, m >  5). The energy in the first nonaxisymmetric large-scale mode is less than in the mode m = 0. This is opposite to the case in Viviani et al. (2018), where simulations with the same rotation rates, but a fixed heat-conduction profile, showed a substantial m = 1 component. Run R3 is in a regime where the axisymmetric and the first nonaxisymmetric mode are of comparable strength, and therefore we characterize this run as being nonaxisymmetric. Upon closer inspection, R3 shows a weak ADW (see, also, Sect. 3.2.2). However, in order to calculate τ_cyc, we follow the same convention as in Viviani et al. (2018) for runs in this regime, and use m = 0 to obtain it. Run R4 has a strong m = 1 component, which is reflected by the presence of the ADW in the lower panel of Fig. 6.

We show the magnetic energy spectra for the radial magnetic field of R3 at three different depths in Fig. 3. The black and red curves correspond, respectively, to the top and middle of the BZ, while the blue curve is the magnetic energy in the DZ. From this figure we can see that larger scales have more energy close to the surface, while smaller scales are excited deeper in the BZ. Large-scale magnetic fields are also present in the DZ.

Fig. 3. Normalized power spectrum of the radial magnetic field at three different depths for R3.

3.2.1. Axisymmetric magnetic field

We show the $\overline{B_{\varphi }}$ near the surface as a function of time (the so-called butterfly diagram) in Fig. 4. We note that other modelling groups plot their butterfly diagrams from the base of the CZ (e.g., Fan & Fang 2014), but in our case, no migration of the magnetic field occurs there. We show the radial magnetic field in time at the bottom of the CZ in Appendix B. Run R1 is characterized by an equatorially symmetric magnetic field with nonmigrating negative polarities at low latitudes, and a poleward migrating positive field at higher latitudes. A stationary negative field is present at all times close to the latitudinal boundary. A similar, oscillatory dynamo solution was reported and analyzed in detail in Viviani et al. (2019). There, it was concluded that two dynamo modes were competing in the model, a stationary and an oscillatory one, the latter with polarity reversals. These latter authors came to the conclusion that this dynamo is driven mostly by turbulent effects, as the differential rotation was found to be weak in the model. Run R1 appears to be another incarnation of such a dynamo in the transition regime from solar-like to anti-solar differential rotation.

Fig. 4. Azimuthally averaged azimuthal component of the magnetic field, $\overline{B_{\varphi }}$, at r = 0.98R for Runs R1 (top), R2 (middle), and R3 (bottom).

$\overline{B_{\varphi }}$ in Run R2 is dipolar (equatorially antisymmetric) at the surface. The polarity is positive in the upper hemisphere and there are no signs of polarity reversals, if the weak ones at the latitudinal boundaries are not counted for. Also, no equatorward migration can be seen, although there is again a tendency for weak poleward migration with a very high frequency.

Run R3 exhibits equatorial propagation of the azimuthal magnetic field, and the pattern in the butterfly diagram (lowest panel in Fig. 4) is very similar to that of Run MHD2 in Käpylä et al. (2019), also showing a similar periodicity of ∼2 yr. In contrast to Run MHD2, the solution shows a pronounced hemispheric asymmetry, with a regular cycle in the upper hemisphere and an irregular periodicity in the lower one. The latter cycle seems to be longer than the one in the upper hemisphere. Such hemispheric dynamos have already been reported in global magnetoconvection models with various different numerical methods (see e.g., Grote & Busse 2001; Busse & Simitev 2006; Käpylä et al. 2010; Gastine et al. 2012). In comparison to our earlier work (Viviani et al. 2018), where no hemispheric dynamos were found, the additional non-linearity resulting from the Kramers opacity law increases the likelihood of obtaining such solutions.

In Fig. 5 we show butterfly diagrams at three different depths for Run R4. Two dynamo modes can be seen at the surface: a high-frequency one in the lower hemisphere, and a lower frequency one in the upper hemisphere with a periodicity similar to that seen in Run R3. Deeper down into the CZ, the high-frequency mode disappears, and we can trace its origin to the depth of 0.80 ≤ r ≤ 0.85, that is, to the bottom of the BZ. However, the lower-frequency mode persists until the OZ, and therefore we infer that it is generated there. The existence of different dynamo modes at different depths has already been reported in other studies (such as Käpylä et al. 2016, 2019).

Fig. 5. Azimuthally averaged azimuthal component of the magnetic field, $\overline{B_{\varphi }}$, for Run R4 at different depths: r = 0.98R (top), r = 0.86R (middle), and r = 0.75R (bottom).

3.2.2. Nonaxisymmetric magnetic field

In Run R3 a weak azimuthal dynamo wave is present. The upper panel of Fig. 6 shows the reconstructed m = 1 mode at 45° above the equator close to the surface as a function of time and longitude. The magnetic field would follow the inclination of this line if it was advected by the surface flow. Deviations indicate drift-independent surface flow. In Fig. 6, the magnetic field does not fall on the line for most of the time, and therefore it has its own motion as a wave, traveling in the prograde direction. Weak ADWs such as the one seen in the case of Run R3 were found to be typical in simulations that are close to the axi- to nonaxisymmetric transition (Viviani et al. 2018). In this latter study it was observed that, when the energy in the modes m = 0 and m = 1 is comparable, the ADW can be affected by the differential rotation, in which case the ADW becomes advected by it for some time intervals (see Fig. 6, upper panel, 0 yr ≤ t ≤ 15 yr). Azimuthal dynamo waves were already found in other numerical studies (Käpylä et al. 2013; Cole et al. 2014; Viviani et al. 2018), but their direction was mostly retrograde, in contrast with observational results (see, e.g., Lehtinen et al. 2016).

Fig. 6. Azimuthal dynamo wave for Runs R3 and R4 as a function of longitude and time, at latitude θ = +45° and depth r = 0.98R. The black and white dashed line represents the differential rotation at the same latitude.

A stronger prograde ADW is also present in Run R4. The wave does not persist at all times, but there are periods when it disappears. Similar behavior was also observed in the temperature maps of the active star II Peg in Lindborg et al. (2011) for example, where a very clearly defined prograde ADW persisted over ten years, but then vanished. Also, Lehtinen et al. (2012) and Olspert et al. (2015) reported ADWs for the young solar analogue star LQ Hya, but these lasted for an even shorter period of time, that is a couple of years.

We attribute the change in the ADW direction to the different heat-conduction prescription in these runs. Moreover, in Run R3 the ADW can even change direction from prograde to retrograde during some short epochs, which is most likely related to the stronger influence from the differential rotation. Such a change of direction in the migration of active longitudes was also observed in the study of, for example, Korhonen et al. (2004) in the case of the intensively studied single active star FK Coma Berenices.

We calculate the period of the ADW, P_ADW, as in Viviani et al. (2018), taking the latitudinal and temporal average of the slope of the reconstructed time evolution of the m = 1 mode. For R3, we obtain P_ADW = −24.5 yr, the minus sign indicating a negative slope, which means a retrograde direction for the dynamo wave, in contrast with the mostly prograde appearance of the migration in the upper panel of Fig. 6. As discussed above, the ADW in R3 is chaotic, even changing its direction, and also disappearing, which means that the average period is not an accurate measure of the period of the wave. P_ADW = 44.9 yr for Run R4, the positive sign indicating a prograde wave, as expected. However from the lower panel of Fig. 6 we would infer a shorter period of ∼30 yr.

4. Conclusions

This paper presents the results of our study of the effect of a dynamically adapting heat-conduction prescription, based on Kramers opacity law, in conjunction with semi-global MHD simulations. The main aim is to determine the effect of this prescription on the two major transitions reported in numerical studies (e.g., Gastine et al. 2014; Viviani et al. 2018). One concerns the rotation profiles, and is the transition from accelerated poles and decelerated equator to a solar-like profile, with a faster equator. The other involves the large-scale magnetic field, and is the transition from an axisymmetric magnetic field, as in the Sun, to a nonaxisymmetric magnetic field found in more rapid rotators. Previous studies (Viviani et al. 2018) found these transitions to occur at the same rotation rate, in contrast with the current interpretation of observations. The fact that simulations usually produce anti-solar differential rotation for the solar rotation rate could indicate that the Sun is in a transitional regime (e.g., Käpylä et al. 2014; Metcalfe et al. 2016), or could also mean that simulations still cannot fully capture the true rotational influence on turbulent convection in the Sun. Lehtinen et al. (2016) reported on the existence of nonaxisymmetric structures in stars with varying rotation rates, and were therefore able to determine quite a sharp transition point in terms of the rotation period, when fields turn from axi- to nonaxisymmetric configurations. According to dynamo theory, these two modes can compete, and there can be a transition region, where both dynamo modes co-exist, as is also clearly demonstrated by the models presented in this paper and those of Viviani et al. (2018). Therefore, the observational transition point must be regarded as a lower limit for the transition in terms of the rotation period, as it could be that the sensitivity of the current instruments is insufficient to detect the very weak nonaxisymmetric components. However, since active longitudes have not been detected on the Sun (Pelt et al. 2006), these two transitions should not be located at the same, nearly solar rotation rate.

In runs with slow rotation, the differential rotation profile is significantly affected by the Kramers opacity law and, as a result, solutions with less solar-like characteristics develop, such as for example an almost rigid body rotation and a minimum at mid-latitudes. The different heat conduction prescription also promotes the formation of a stably stratified layer in the lower quarter of the domain that is rather isotropic in latitude. For faster-rotating runs, the rotation profile is solar-like, but still maintains the minimum at mid-latitudes, and a latitudinally changing subadiabatic region forms near the equator. Also, the Coriolis number is lower than in the corresponding cases using fixed profiles for heat conduction, which is most likely the largest contributing factor to pushing the transition from anti-solar to solar differential rotation to an unwanted direction of more rapid rotation rates.

The convective transport is efficient, isotropic and almost radial everywhere in the convective region in models with slow rotation (Runs R1 and R2), while it becomes strongly concentrated to the equatorial region in runs with more rapid rotation (Runs R3 and R4). Also, the BZ becomes shallower close to the tangent cylinder in the rapid-rotation regime. Moreover, hemispheric asymmetries in the convection zone structure are seen in the run with the fastest rotation (Run R4).

The large-scale magnetic field is axisymmetric in Runs R1 and R2, while for Run R3 and R4 the first nonaxisymetric mode is dynamically more important. Both the fast rotating runs have a hemispherically asymmetric oscillating magnetic field with a periodicity of ∼2 yr. As in Viviani et al. (2018), the magnetic cycle lengths do not strongly depend on the rotation period overall. The strong magnetic field in all the runs originates from the subadiabatic layer. In Run R4 a high-frequency mode is present in the lower hemisphere. This component is generated at the bottom boundary of the BZ. The co-existence of multiple dynamo modes at different depths of the convection zone is consistent with previous studies (e.g., Käpylä et al. 2016) using prescribed profiles for heat conduction. In this latter study, the high-frequency mode was generated near the surface, while the low-frequency one formed in the middle of the CZ.

In the nonaxisymmetric runs, ADWs are present: a weak one for Run R3 and a stronger one for Run R4. In both cases, the direction is prograde, in agreement with photometric observations (Lehtinen et al. 2016). In a previous numerical study using a prescribed heat-conduction profile (Viviani et al. 2018), a preference for retrograde ADWs was found. The ADWs also show time variations. For Run R3, the ADW is rather weak and the differential rotation can advect it for some time, changing the direction of the wave. This could be caused by the comparable relative energies in the m = 0 and m = 1 modes. In Run R4 the stronger ADW disappears at certain times. Such behavior is also observed for active stars (e.g., Korhonen et al. 2004; Lindborg et al. 2011; Lehtinen et al. 2012; Olspert et al. 2015), where the active longitudes disappear or have the same velocity as the surface rotation.

In summary, in this study we show that both the major transitions related to stellar dynamos are affected by the use of a more physical description of heat conduction in global magnetoconvection simulations. The differential rotation profiles undergo a significant change near the transition from anti-solar to solar-like differential rotation, but all the runs are still in Taylor–Proudman balance, with almost cylindrical isocontours. For the same rotation rates, the convective velocities are higher, and therefore Coriolis numbers are lowered, resulting in an anticipated transition to anti-solar differential rotation, in contrast with observations. The transition from axi- to nonaxisymmetric magnetic fields is shifted towards higher rotation rates. The direction of the ADW is reverted with respect to previous studies, producing a better agreement with observations.

Acknowledgments

MV and MJK acknowledge support from the Independent Max Planck Research Group “SOLSTAR”. MV acknowledges the support from the International Max Planck Research School for Solar System Science at the University of Göttingen. MJK acknowledges the support of the Academy of Finland ReSoLVE Centre of Excellence (grant No. 307411). This project has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (project “UniSDyn”, grant agreement No. 818665). The authors would like to thank Petri Käpylä, Jörn Warnecke and Matthias Rheinhardt for constructive comments and conversations.

References

Appendix A: Comparison with rotation profiles in previous works

For comparison, we reproduce in Fig. A.1 the differential rotation profiles for the simulations with fixed heat conduction profiles (Runs C3, D, G^a and H) from Viviani et al. (2018), corresponding, respectively, to R1, R2, and R4 (the latter for G^a and H) presented in this paper. We also show again the wedge simulation MHD2 from Käpylä et al. (2019), this time showing both the hemispheres.

Fig. A.1. Differential rotation profiles for runs C3, D, Ga, and H from Viviani et al. (2018), and MHD2 from Käpylä et al. (2019). The green lines for MHD2 indicate the internal structure as for the simulations in the main text.

The rotational profile of Run D was not presented in our previous publications. An equatorial prograde flow and mid-latitude minima are present, while in its Kramers counterpart, Run R2, the profile is swapped to anti-solar with a decelerated equator and accelerated mid-latitude regions.

Appendix B: Azimuthal magnetic field at the bottom of the CZ

Figure B.1 shows the time evolution of the radial magnetic field at the bottom of the convection zone. As an indicative depth of the CZ we used rDZ from Table 1. However, this is an averaged value and, especially when anisotropies are present, as in the case of R4, it is an imprecise measure of the depth of the convective layer.

Fig. B.1. Radial magnetic field calculated at rDZ from Table 1. (a) R1[rDZ = 0.767], (b) R2[rDZ = 0.809], (c) R3[rDZ = 0.738], (d) R4[rDZ = 0.740].

In R1, the bands of negative time-varying surface magnetic field at latitudes ±15 ° −30° are not present at depth rDZ. The radial magnetic field in R2 shows multiple polarity reversals in latitude, and a poleward propagation in time at high latitudes. No trace of the equatorial propagation at mid-latitudes is present at the bottom of the CZ in R3. As already shown in the main text, the high-frequency mode in R4 is not visible, originating higher up in the CZ, at depth r ≃ rBZ.

All Tables

All Figures

	Fig. 1. Radial Lenth normalized by L0. The arrows show the direction of $L_r^{\text{enth}}$. The continuous, dashed, and dash-dotted green lines represent, respectively, the bottom boundaries of BZ, DZ, and OZ.
In the text

	Fig. 2. Differential rotation profiles. Continuous, dashed, and dash-dotted green lines are as in Fig. 1.
In the text

	Fig. 3. Normalized power spectrum of the radial magnetic field at three different depths for R3.
In the text

	Fig. 4. Azimuthally averaged azimuthal component of the magnetic field, $\overline{B_{\varphi }}$, at r = 0.98R for Runs R1 (top), R2 (middle), and R3 (bottom).
In the text

	Fig. 5. Azimuthally averaged azimuthal component of the magnetic field, $\overline{B_{\varphi }}$, for Run R4 at different depths: r = 0.98R (top), r = 0.86R (middle), and r = 0.75R (bottom).
In the text

	Fig. 6. Azimuthal dynamo wave for Runs R3 and R4 as a function of longitude and time, at latitude θ = +45° and depth r = 0.98R. The black and white dashed line represents the differential rotation at the same latitude.
In the text

	Fig. A.1. Differential rotation profiles for runs C3, D, Ga, and H from Viviani et al. (2018), and MHD2 from Käpylä et al. (2019). The green lines for MHD2 indicate the internal structure as for the simulations in the main text.
In the text

	Fig. B.1. Radial magnetic field calculated at rDZ from Table 1. (a) R1[rDZ = 0.767], (b) R2[rDZ = 0.809], (c) R3[rDZ = 0.738], (d) R4[rDZ = 0.740].
In the text