Physically motivated heat conduction treatment in simulations of solar-like stars: effects on dynamo transitions

Context. Results from global magnetoconvection simulations of solar-like stars are at odds with observations in many respects: They show a surplus of energy in the kinetic power spectrum at large scales, anti-solar differential rotation profiles, with accelerated poles and a slow equator, for the solar rotation rate, and a transition from axi- to non-axisymmetric dynamos at a much lower rotation rate than what is observed. Even though the simulations reproduce the observed active longitudes in fast rotators, their motion in the rotational frame (the so-called azimuthal dynamo wave, ADW) is retrograde, in contrast to the prevalent prograde motion in observations. Aims. We study the effect of a more realistic treatment of heat conductivity in alleviating the discrepancies between observations and simulations. Methods. We use physically-motivated heat conduction, by applying Kramers opacity law, on a semi-global spherical setup describing convective envelopes of solar-like stars, instead of a prescribed heat conduction profile from mixing-length arguments. Results. We find that some aspects of the results now better correspond to observations: The axi- to non-axisymmetric transition point is shifted towards higher rotation rates. We also find a change in the propagation direction of ADWs so that also prograde waves are now found. The transition from anti-solar to solar-like rotation profile, however, is also shifted towards higher rotation rates, leaving the models into an even more unrealistic regime. Conclusions. Although a Kramers-based heat conduction does not help in reproducing the solar rotation profile, it does help in the faster rotation regime, where the dynamo solutions now match better with observations.


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 have allowed also to uncover the subsurface rotation (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 angular velocity contours constant on cylinders would be observed. Hence, the Sun is able to break the Taylor-Proudman balance with 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, 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 on an enhancement of ∼ 2.5 K at the Sun's poles. One possibility to explain such a temperature gradient theoretically is to argue for 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).
Modelling efforts of stellar convection in spherical or semispherical shells still struggle in producing 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 to imply 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 of present-day numerical setups are using 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 observations (Gizon & Birch 2012). All the above mentioned discrepancies between observations and numerical models are collectively known as "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 surfacedriven 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 identified first in the Earth's atmosphere (Deardorff 1961(Deardorff , 1966, and hence was termed as "Deardorff layer". The formation of such sub-adiabatic layers has been reported in the hydrodynamic studies of, e.g., Hotta (2017), Korre et al. (2017), and Käpylä et al. (2017). Especially the study by Käpylä et al. (2017) is relevant here, as they demonstrated the emergence of a substantial sub-adiabatic layer, and the existence of non-local surface driving of convection, by 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 buoyancy zone, plus the sub-adiabatic part, now denoted as Deardorff layer. The depth of the layers was not determined a priori, but rather was 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 observed 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 too 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 e.g., from the action of a small-scale dynamo instability (see, e.g., Kazantsev 1968) operating the the CZ. Thereby 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, the ratio of the viscosity and resistivity of the fluid. Karak et al. (2018) investigated such a situation numerically, and 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) and Käpylä et al. (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 roughly by 180 degrees in longitude (e.g., Jetsu 1996). Active longitudes usually migrate in the rotational frame of the star, forming azimuthal dynamo waves (ADW) (e.g., Berdyugina & Tuominen 1998). The direction of migration of these structures can follow the plasma rotation, and in this case we will talk about prograde ADWs; they can also drift in the opposite direction (retrograde ADWs), or they can stand still with respect to the observer's point of view (standing ADW). 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 shortlived 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 non-axisymmetric dynamo modes operating in these stars (Tuominen et al. 2002), in contrast to the axisymmetric dynamo operating in less active stars. The transition from non-axisymmetric to axisymmetric dynamos has also been studied numerically (Cole et al. 2014;Viviani et al. 2018), but these studies reported a majority of retrograde ADWs, and a transition from axi-to non-axisymmetric solutions at too low rotation rates in comparison to observations. Both studies were using prescribed and MLT-motivated profiles for heat conduction, resulting in a priori fixed depth of the convection zone.
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 to be able to study both axi-and non-axisymmetric dynamo solutions.

Setup and Model
We apply a similar setup as in Käpylä et al. (2013) and Käpylä et al. (2019), representing the outer envelope of a solar-like star, We solve numerically the system of MHD equations: where ρ and U are the density and the velocity field, g = −GM/r 3 is the gravitational acceleration, with G being the gravitational constant and M the mass of the star. Ω 0 = Ω 0 (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 µ 0 are the pressure, the viscosity, and the magnetic permeability in vacuum, while η is the magnetic diffusivity. S is the rate-of-strain tensor. F rad and F SGS are the radiative and sub-grid scale (SGS) fluxes, expressed by: (2) K is the radiative heat conductivity and χ SGS is the SGS heat diffusivity, assumed to be constant. s is the fluctuating entropy, s = s − s θ , 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. Its flux is expressed as: 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;Käpylä 2019;Käpylä et al. 2019;Käpylä et al. 2020): where ρ 0 and T 0 are reference values for density and temperature. The constant K 0 is defined via: L 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 θ = θ 0 and θ = π − θ 0 . 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: 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 Ω 0 , ν, η, χ SGS , K 0 , ρ 0 , T 0 and the energy flux at the bottom, F bot = −K∂ r T | r=0.7R .
Moreover, important non-dimensional input parameters are the magnetic and SGS Prandtl numbers: Output parameters of the simulations are the fluid and magnetic Reynolds numbers: with u rms = 3/2 U 2 r + U 2 θ rθφt rms velocity and k f = 2π/0.3R the wavenumber of the largest eddy, corresponding to the radial extent, and the Coriolis number: quantifying the relative importance of rotation and convection. Physical units are chosen using the solar radius R = 7 · 10 8 m, the solar angular velocity Ω = 2.7 · 10 −6 s −1 , the density at the bottom of the solar convection zone ρ bot = 200 Kg m −3 , and the magnetic permeability in vacuum µ 0 = 4π · 10 −7 H m −1 . We performed our simulations using the PENCIL CODE 1 , a high-order, finite-difference, open source code for solving the magnetohydrodynamic equations.

Results
The simulations and their defining parameters are summarized in Table 1. Run 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 for heat conduction on the anti-solar to solar-like differential rotation transition and the transition to non-axisymmetric magnetic fields. Run C3 was the simulation with the slowest rotation rate showing both accelerated equator and non-axisymmetric magnetic field, hence it is a good choice for this study. Run D had a rotation rate 2.1 times the solar value, exhibited a non-oscillatory 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 in between the runs was most likely connected to the higher amount of differential rotation in C3 than in D, see also Sect.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 1/4 of the full longitude, hence not allowing for non-axisymmetric 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.

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 Figures 1 and 2 and Figure A.1 in Appendix A. The radial enthalpy flux is defined as F enth r = c p (ρu r ) T . The bottom of the buoyancy zone (BZ), in which the radial enthalpy flux is greater than zero, F enth r > 0, and the radial entropy gradient is negative, ∂ r s < 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 enth r > 0 and ∂ r s > 0, by a dashed line, and the bottom of the overshoot zone (OZ), for which F enth r < 0 and ∂ r s > 0, with a dash-dotted line. What we denote as 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 enth r ≈ 0 and ∂ r s > 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 Equation (9); the second one, denoted as Co rev in Table 1, takes into consideration the wavenumber of the revised convection zone (BZ and DZ), therefore we use k f rev = 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.
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. Also the OZ is thin, hence the run develops a very thick RZ. In Runs R1 and R2 the thickness of the layers does not change considerably as function of latitude, except for a slight tendency of the DZ of becoming thicker near the equatorial region for Run R2. For Run R3 the convection zone structure at higher latitudes again resembles the one of Run R1. In the equatorial region, however, the convection zone becomes very deep. Close to the tangent cylinder the BZ becomes considerably shallower, and the DZ developes 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.

Enthalpy flux
We inspect the radial enthalpy flux, F enth r , by representing the enthalpy luminosity, L enth = 4πr 2 F enth r , in Figure 1 with black arrows. The enthalpy flux in Run R1 and R2 is isotropic in latitude and rather radial everywhere in the BZ. There is a slight tendency of the flux being enhanced/decreased in the equatorial region for Run R1/R2. 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 of the enthalpy flux in the regions of the tangent cylinder. A clear equatorial asymmetry is present in L enth r for Run R4. This asymmetry is also reflected in the convection zone structure, as discussed above.

Differential rotation
The last two columns in Table 1 quantify the relative radial and latitudinal differential rotation, defined as: 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, θ 0 ) + Ω(R, π − θ 0 )) /2. We show the differential rotation profiles in Figure 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 being quite solar-like (see Figure A.1, first panel). With an adaptable heat conduction prescription the rotation profile is less solar-like (see, Figure 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, 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. In this case, however, 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 Run C3 of theirs. Therefore, the contribution of the differential rotation to the large-scale dynamo should be negligible. To assert this, though, a more thorough analysis would be required.
In comparison to Run D in Viviani et al. (2018) (Figure 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 much in latitude and depth.
The rotation profile of Run R3 (Figure 2) is very similar to the one from its wedge counterpart shown in Figure A.1, third panel. 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 corresponds to 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 Figure A.1. Hence, the usage of the Kramers opacity law produces higher convective velocities and, therefore, smaller Co. The values for ∆ r and ∆ θ coincide with the ones from G a . The rotation profile, shown in Figure 2, rightmost panel, closely resembles that of Run H, with a deep minimum at mid-latitudes.

Dynamo solutions
All the presented runs develop large-scale dynamo (LSD) action, and thereby sustain a magnetic field. The magnetic Reynolds numbers, however, are too low to allow for small-scale dynamo action (SSD). Our dynamo solutions, therefore, 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). 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 Table 2, last column, 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, hence 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.
Run 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, were showing a substantial m = 1 component. Run R3 is in a regime where the axisymmetric and the first nonaxisymmetric mode have comparable strengths, hence we characterize this run as being non-axisymmetric. R3, in fact, shows a weak azimuthal dynamo wave (ADW, see, also 3.2.2). For calculating τ cyc , however, 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 Figure 5, lower panel.

Axisymmetric magnetic field
We show the azimuthally averaged longitudinal magnetic field near the surface as a function of time (the so-called butterfly diagram) in Figure 3. Run R1 is characterized by equatorially symmetric magnetic field, with non-migrating negative polarities at low latitudes and 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 analysed 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. This dynamo was concluded to be driven mostly by turbulent effects, as the differential rotation was found 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.
The longitudinally averaged B φ 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, albeit there is again a tendency for a weak poleward migration with a very high frequency. Fig. 2. Differential rotation profiles. Continuous, dashed and dash-dotted green lines as in Figure 1.  Table 2. Magnetic energy from the decomposition in the first 11 spherical harmonics (0 ≤ l, m ≤ 10) of the near-surface (r = 0.98R) radial magnetic field. The labels E dec m indicate the energy in the corresponding m mode, in units of 10 5 J m −3 . We define 0 ≤ l, m ≤ 5 the large-scale field. The numbers in parenthesis represent the power of ten. Last column indicates the characteristic time of the dynamo, calculated on the dominant mode, indicated in the subscript.
Run R3 exhibits equatorial propagation of the azimuthal magnetic field, and the pattern in the butterfly diagram (lowest panel in Fig. 3) 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.
In Fig. 4 we show butterfly diagrams at three different depths for Run R4. At the surface, it shows two dynamo modes: a highfrequency one in the lower hemisphere and a lower-frequency one, with a periodicity similar to Run R3, in the upper hemisphere. As we go deeper down to the convection zone, the highfrequency mode disappears, and we can trace its origin to the depth of 0.80 ≤ r ≤ 0.85, therefore to the bottom of the BZ. The lower-frequency mode, however, persists until the OZ, and therefore we infer that it is generated there. The existence of different dynamo modes at different depths has been reported already in other studies (such as Käpylä et al. 2016Käpylä et al. , 2019.

Non-axisymmetric magnetic field
In Run R3 a weak azimuthal dynamo wave is present. In Figure 5, upper panel, we plot the reconstructed m = 1 mode at 45 o above the equator close to the surface as a function of time and longitude. The black-white dashed line represents the pattern of differential rotation at the same latitude. In the absence of a dynamo wave, the magnetic field would follow the propagation speed of this pattern. Instead, in Figure 5, the magnetic field does not fall on the line for most of the time, hence, it has its own motion as a wave, travelling 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 non-axisymmetric transition (Viviani et al. 2018). In this 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 it becomes advected by it for some time intervals (see Figure 5, left panel, 0 yr ≤ t ≤ 15 yr). ADWs 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). 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 behaviour was also observed, for example, in the temperature maps of the active star II Peg in Lindborg et al. (2011), 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 on ADWs on the young solar analogue star LQ Hya, but these lasted even for a shorter period of time, of about 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, most likely related to the stronger influence from the differential rotation. Such a change of direction in the migration of active longitudes has also been observed, e.g. in the study of 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, therefore a retrograde direction for the dynamo wave, in contrast with the mostly prograde appearance of the migration in Figure 5, upper panel. As discussed before, the ADW in R3 is chaotic, even changing its direction, and also disappearing, and this causes the average period to be an inaccurate measure. P ADW = 44.9 yr for Run R4, the positive sign indicating a prograde wave, as expected. From the lower panel of Figure 5, however, we would infer a shorter period ∼ 30 yr.

Conclusions
In this paper we studied the effect of a dynamically adapting heat conduction prescription, based on Kramers opacity law, on semi-global MHD simulations. The main aim was to determine its effect 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 faster equator. The other one involves the large-scale magnetic field, and is the transition from an axisymmetric magnetic field, as in the Sun, to a non-axisymmetric one found in more rapid rotators. Previous studies (Viviani et al. 2018) reported these transitions occurring 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 it could also mean that simulations still cannot fully capture the right rotational influence on turbulent convection in the Sun. The study of Lehtinen et al. (2016) reported on the existence of non-axisymmetric structures in stars with varying rotation rates, and hence could determine quite a sharp transition point, in terms of the rotation period, when fields turn from axito non-axisymmetric 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 Viviani et al. (2018). Hence, the observational transition point must be regarded as a lower limit, in terms of the rotation pe- riod, for the transition, as it could be that the sensitivity of the current instruments is not high enough to detect the very weak non-axisymmetric components. Since active longitudes have not been detected on the Sun (Pelt et al. 2006), though, 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, like an almost rigid body rotation and a minimum at mid-latitudes, develop. The different heat conduction prescription also promotes the formation of a stably stratified layer, rather isotropic in latitude, in the lower quarter of the domain. For faster rotating runs, the ro- tation 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 push the anti-solar to solar differential rotation transition 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 gets 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 Run R1 and R2, while for Run R3 and Run R4 the first non-axisymetric mode is dynamically more important. Both the fast rotating runs have a hemispherically asymmetric oscillating magnetic field, with a periodicity of ∼ 2 years. As in Viviani et al. (2018), the magnetic cycle lengths do not overall depend strongly on the rotation period. 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 study, the highfrequency mode was generated near the surface, while the lowfrequency one in the middle of the CZ.
In the non-axisymmetric 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 the previous numerical study using a prescribed heat conduction profile (Viviani et al. 2018), we found a preference for retrograde ADWs. 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 what is 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 have shown that both the major transitions related to stellar dynamos are affected by the use of a more physical description of heat conduction in global magneto-convection simulations. The differential rotation profiles undergo a significant change near the anti-solar to solarlike differential rotation transition, but all the runs are still in Taylor-Proudman balance, with almost cylindrical isocontours. For the same rotation rates, the convective velocities are higher, hence Coriolis numbers are lowered, resulting in an anticipated transition to anti-solar differential rotation, in contrast with observations. The transition from axi-to non-axisymmetric 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.