© The Authors 2024

1. Introduction

Understanding the spatial distribution of the light emission at the surface of stars is of prime importance to characterize their properties. This is not only necessary to convert the observed magnitude of a star into its actual luminosity, but also essential to help analyze the detailed mapping of surface features obtained by interferometry or Zeeman-Doppler reconstruction techniques (e.g., Roettenbacher et al. 2016; Finociety et al. 2021; Petit et al. 2022), as well as interpret the photometric light curves of eclipsing binaries or transiting exoplanets (e.g., Valio 2017; Özavcı et al. 2018).

In addition to the limb-darkening effect, which refers to the decrease in the light flux emitted onto the line of sight from the center of the stellar disk to its limb (e.g., Morello et al. 2017; Maxted 2023), internal dynamics itself can affect the brightness distribution on the surface of stars. Stellar rotation is known to be a key ingredient in this regard. In rapidly rotating, hot (early-type) stars with a stably stratified radiative envelope, the centrifugal structural flattening along the rotation axis makes the poles hotter and brighter than the equator; this is the so-called gravity-darkening effect, which was first formulated by von Zeipel (1924) and well constrained by interferometric data and advanced 2D stellar models (e.g., Espinosa Lara & Rieutord 2011; Domiciano de Souza et al. 2014; Bouchaud et al. 2020). In contrast, in cooler (late-type) stars with a convective envelope, the picture remains unclear. Modeling the outgoing heat flux through turbulent rotating envelopes is indeed a complex task and requires one to account for the impact of three main physical ingredients on the convective flows. The first one is naturally the effect of the Coriolis acceleration that tends to organize the convection in a columnar shape; the second one is the effect of magnetic fields; and the third one is the effect of the centrifugal acceleration. The latter is expected to make the star oblate as well as decrease the convective buoyancy driving in the regions close to the equator. Understanding the impact of the centrifugal acceleration on the surface brightness distribution of stars is difficult and is deferred to future works.

Regarding the Coriolis force, its effect has already been investigated in previous work (Raynaud et al. 2018, hereafter Paper I). Using numerical simulations of highly stratified rotating convective envelopes, we have demonstrated that in turbulent regimes (i.e., far from the onset of convection), the latitudinal distribution of brightness is determined by the competition between the inertia of convective eddies and the Coriolis acceleration. When the Coriolis acceleration dominates inertial forces at the surface (i.e., low surface Rossby number), a strong prograde azimuthal flow is sustained close to the equator owing to Reynolds stresses induced by Rossby-like eddies. This tends to inhibit the outgoing heat flux in the equatorial region, making the poles brighter. On the contrary, when inertia tends to dominate the Coriolis acceleration (i.e., large surface Rossby number), the heat flux tends to become uniform in latitude and decouple from the surface differential rotation profile. To avoid any misleading confusion with the original concept of gravity darkening, we decide in this second paper of the series to term Coriolis darkening this modulation of the surface brightness. In this work, we aim to extend the results of Paper I by taking magnetic effects generated by dynamo action into account.

As shown by spectropolarimetric observations, magnetic fields are ubiquitous in low-mass stars and exhibit diverse topologies, ranging from small-scale fields to large-scale dipolar fields (Donati & Landstreet 2009; Morin et al. 2010; Folsom et al. 2018; Kochukhov 2021; Donati et al. 2023). Such magnetic fields result from dynamo mechanisms that are ruled by the turbulent motions in the convective envelopes of these stars. These mechanisms permit to transfer a part of the kinetic energy into magnetic energy through nonlinear amplification processes of electrical currents, allowing magnetic fields to be sustained against Ohmic dissipation. Since the pioneer works of Gilman (1983), Glatzmaier (1985), Glatzmaiers & Roberts (1995), numerical modeling of this multiscale mechanism has been developing considerably. However, direct numerical simulations still cannot meet realistic stellar conditions and have to rely on some convective approximations to simplify the problem, such as the Boussinesq or the anelastic approximations (e.g., Ogura & Phillips 1962; Gough 1969; Braginsky & Roberts 1995; Lantz & Fan 1999). Both are sound-proof approximations of the magnetohydrodynamic (MHD) equations that are well adapted to model low-Mach turbulence. Either in the incompressible Boussinesq paradigm (e.g., Christensen & Aubert 2006; Schrinner et al. 2012; Yadav et al. 2013; Petitdemange 2018; Menu et al. 2020) or in the anelastic framework that allows for variations in the medium density (e.g., Gastine et al. 2012; Duarte et al. 2013; Schrinner et al. 2014; Raynaud et al. 2015; Schwaiger et al. 2021; Zaire et al. 2022), magnetic fields obtained in numerical simulations can be roughly classified as either dominated by the large-scale poloidal dipolar component or not, even though a more sophisticated description in terms of force balance and temporal evolution is clearly suggested by the numerous results on the subject. In this framework, anelastic dynamos were shown to exhibit some specific behaviors compared to Boussinesq models. For instance, an increase in the density stratification seems to disadvantage dipolar configurations (Gastine et al. 2012; Jones 2014), unless the magnetic Prandtl number, Pm, is large enough (Schrinner et al. 2014). Such differences between Boussinesq and anelastic configurations result from the difference of mass distribution and gravity profile in both cases, which directly affect the properties of the convective cells with depth (Raynaud et al. 2014). For instance, Raynaud et al. (2015) showed that a large density drop toward the outer regions is associated with an increase in the convective inertia, which can destabilize at some point a dipolar magnetic configuration. Besides, when the Lorentz force plays a major role at large scales, that is, above the dynamo threshold and with large values of the magnetic Prandtl number, strong dipole magnetic fields can be maintained even in highly stratified fluids. This regime (or strong branch) is characterized by a balance between the Coriolis acceleration and the Lorentz force although the convective inertia remains significant close to the surface (Dormy 2016; Dormy et al. 2018; Menu et al. 2020; Schwaiger et al. 2021; Zaire et al. 2022).

These previous studies constitute an appropriate reference basis to investigate the combined effects of the Coriolis acceleration and self-sustained magnetic fields on the brightness distribution at the surface of convective stellar envelopes. Actually, Yadav et al. (2016) have already partially tackled the question. In the regime where the Coriolis acceleration dominates turbulent inertia at the surface, they have shown that self-sustained magnetic fields can globally quench the azimuthal mean flow close to the equator, relaxing its blocking action on the outgoing heat flux in this region. However, this last study is restricted to Boussinesq models with a constant magnetic Prandtl number (namely, equal to unity). In this work, we aim to question this result exploring a wider parameter space (e.g., large density stratification, varying magnetic Prandtl numbers), and to extend at the same time the results of Paper I to the magnetic case. To do so, we use about 100 three-dimensional, self-consistent dynamo models computed by direct numerical simulations. Our goal is to understand general physical trends as a function of the main control parameters. Hereafter, we specifically focus on stars cooler than the Sun with thick convective envelopes. For the sake of simplicity, the aspect ratio of the simulated convective shell is set to χ = 0.35, where χ = r_i/r_o with r_i and r_o the inner and outer radii of the shell, respectively. We recall that the Sun is characterized by a relatively thin convective envelope, with χ_⊙ ≈ 0.7. In comparison, our setup can be regarded as representative of different evolutionary stages: M-dwarf stars of about 0.35 M_⊙ (e.g., Chabrier & Baraffe 1997); pre-main sequence stars leaving the Hayashi phase while developing a radiative core (e.g., Siess et al. 2000; Gregory et al. 2016; D’Antona 2017); low- and intermediate-mass red giant stars at the beginning of the red giant branch (e.g., Pinçon et al. 2020). On the one hand, cool M dwarfs are currently one of the most interesting targets to characterize rocky Earth-like exoplanets by the method of the transits owing to the accessible large signal-to-noise ratio (e.g., Garcia et al. 2022). On the other hand, pre-main sequence evolution is still subject to huge uncertainties related for instance to disk accretion phenomena and stellar magnetism (e.g., Alecian & Lebreton 2020). Finally, red giant stars, with their wide mass range (and thus age range) and their large luminosity, stand for a mine of information for galactic archaeology research (e.g., Mosser et al. 2016; Miglio et al. 2021). All these types of stars are therefore important targets and deserve special attention, motivating our choice. The numerical simulations used in this work are introduced in Sect. 2. The different magnetic morphologies encountered in such simulations are briefly recalled in Sect. 3. The effect of the magnetic fields on the mean latitudinal luminosity distribution is then addressed in Sect. 4. The results are discussed in Sect. 5 and the conclusions are formulated in Sect. 6.

2. Anelastic simulations of convective envelopes

2.1. Physical setup

Our setup is similar to Raynaud et al. (2015). We consider a convective spherical shell of width d that rotates about the z-axis at angular velocity Ω; it is filled with a perfect, electrically conducting gas of kinematic viscosity ν, thermal diffusivity κ, specific heat c_p, and magnetic diffusivity η (all supposed to be constant). We assume that the mass is concentrated inside the inner sphere and neglect the centrifugal force. The gravitational acceleration is thus approximated by g = −(GM/r²) e_r, where G is the gravitational constant, M is the central mass, and e_r is the radial unit vector. For the reference state, we assume an adiabatic hydrostatic equilibrium (i.e., at marginal stability with respect to convection). The reference state is characterized by a constant polytropic index n (hereafter, n = 2 by hypothesis¹) and a given number of density scale heights between the bottom boundary and the surface, that is,

$$\begin{array}{c}\hfill N_{\rho }=ln\left[\frac{\rho (r_{\mathrm{i}})}{\rho (r_{\mathrm{o}})}\right],\end{array}$$(1)

where ρ(r) is the reference state density. Convection is then driven by imposing a fixed entropy drop ΔS between the inner and outer spheres. We apply stress-free boundary conditions for the velocity field, while the magnetic field matches a potential field inside and outside the fluid shell.

The governing equations rely on the LBR formulation of the anelastic approximation (Jones et al. 2011) and are recalled in Appendix A for completeness. The magnetohydrodynamic system in Eqs. (A.4)–(A.8) involves four other control parameters, namely the Rayleigh number, the Ekman number, and the thermal and magnetic Prandtl numbers, which are defined respectively as

$$\begin{array}{c}\hfill \mathrm{Ra}=\frac{GMd\mathrm{\Delta }S}{\nu \kappa c_{\mathrm{p}}},E=\frac{\nu }{\mathrm{\Omega }{d}^2},\mathrm{Pr}=\frac{\nu }{\kappa }\mathrm{and}\mathrm{Pm}=\frac{\nu }{\eta }·\end{array}$$(2)

In the convective envelope of stars, these dimensionless parameters reach extreme values owing to the very low values of the viscosity and diffusivity; for instance, in the solar convective zone, we expect Ra ∼ 10²⁰, E ∼ 10⁻¹⁵, Pr ∼ 10⁻⁷ and Pm ∼ 10⁻⁵. Considering such regimes in direct numerical simulations of extended stellar magnetized convective zones is currently impossible, as it will demand a much too long computing time to resolve the smallest dissipation length scale (e.g., Kupka & Muthsam 2017). Given these numerical limitations, we can instead rely on systematic parameter studies in more moderate and accessible regimes; the aim is to obtain valuable insights on the global dynamics, applicable to some extents to the stellar case (e.g., see the discussion in Sect. 5.1). The unrealistic large values of the viscosity and diffusivity can therefore be regarded as effective transport coefficients mimicking the effect of turbulent small-scale motions on the large-scale flow. Keeping this in mind for our set of simulations, the Ekman number, E, is set to low accessible values between 3 × 10⁻⁵ and 3 × 10⁻⁴, and the thermal Prandtl number, Pr, remains close to unity. The parameter N_ρ covers values between 3 and 6 (i.e., a density contrast between the inner and outer boundaries going from about 20 to 400). We notice that even setting N_ρ = 6 does not permit to model the thin surface layer of stars where the density is expected to drop by several orders of magnitude. In this very narrow shell, the anelastic approximation fails as the flow becomes supersonic and is dominated by small-scale eddies. Nevertheless, this setup is expected to grasp the general large-scale dynamics in a substantial part of the convective bulk, independently of these very near-surface layers. We also explore increasing levels of turbulence by varying the ratio Ra/Ra_c from unity up to about 30, with Ra_c the critical Rayleigh number at the linear onset of convection. This latter has been calculated solving the boundary value problem of the linearized hydrodynamic equations as in Jones et al. (2009). They are reported in Table B.1 for the control parameters considered in this study, with typical values on the order of 10⁶. The magnetic Prandtl number, Pm, is then systematically varied between 0.1 and 5, which enables us to explore different magnetic field intensities and morphologies. The covered parameter space is briefly summarized in Table B.2. Finally, as initial conditions, both strong dipoles or random seeds have generally been used to account as much as possible for the two dynamo branches at low Rossby numbers (e.g., Raynaud et al. 2015).

2.2. Numerical setup and output diagnostics

For all the considered models, the equations are integrated at least over about one magnetic diffusion timescale τ_η = d²/η with the 3D-MHD codes PARODY (Dormy et al. 1998; Schrinner et al. 2014) or MAGIC (Gastine et al. 2012; Schaeffer 2013). In both codes, the vector fields are transformed into scalars using a poloidal-toroidal decomposition. The equations are discretized in the radial direction with a finite-difference scheme for PARODY and an expansion onto the Chebyshev polynomials for MAGIC; on each concentric sphere, variables are expanded onto the spherical harmonic basis, with angular degrees ℓ and azimuthal orders m. Typical resolutions in the radial direction are up to 240 mesh points and the horizontal spectral decomposition is truncated at ℓ_max ∼ m_max ≤ 400 for the more turbulent and stratified models. Both solvers give comparable results for the set of models considered here. For each model, the total kinetic energy spectrum has been checked to drop monotonically by at least two orders of magnitude from the scale at maximum power spectrum toward the smallest considered scale. This represents our empirical criterion to ensure numerical convergence.

By construction, the entropy is set to a constant value and the radial velocity vanishes at the surface, so that the surface heat flux results only from the radial entropy gradient. It is commonly represented by the surface Nusselt number Nu(r_o, θ, φ)=(e_r ⋅ ∇S)/(e_r ⋅ ∇S​_c), where S​_c is the conductive entropy profile, that is, the profile of entropy that we would get if the energy was transported by thermal diffusion only. The Nusselt number hence measures the efficiency of the heat transport by convection. In this work, we are interested in studying the combined effects of the magnetic fields and the Coriolis acceleration on the global surface brightness distribution. For the sake of convenience, we consider the mean surface Nusselt number averaged over the duration of the simulation T_run (i.e., the time spent after reaching a stationary state and larger at least than about τ_η) and over the azimuthal direction, that is (see Paper I for details),

$$\begin{array}{c}\hfill {\overline{\mathrm{Nu}}}_{r_{\mathrm{o}}}(\theta )=-\frac{(1-e^{-N_{\rho }})r_{\mathrm{o}}^2{w}_{\mathrm{o}}}{n{c}_1}\frac{1}{2\pi T_{\mathrm{run}}}{\displaystyle \int {\int }_0^{2\pi ,T_{\mathrm{run}}}}{\frac{\partial S}{\partial r}|}_{r_{\mathrm{o}}}\mathrm{d}\phi \mathrm{d}t.\end{array}$$(3)

This diagnostic is sufficient in a first step to investigate the effect of the magnetic field on the mean heat flux distribution observed in the hydrodynamic case (Paper I). It will nevertheless not permit to account for the rapid temporal variations generated by nonaxisymmetric features (e.g., nonpolar stellar spots), which deserve a dedicated study out of the scope of this paper.

3. Magnetic morphologies in stratified spherical dynamo models

In this section, we briefly recall the general magnetic morphologies that are encountered in numerical simulations of thick convective envelopes and that have already been observed in previous studies. For later purpose, we also emphasize the global impact of magnetic fields on the large-scale differential rotation.

3.1. A diversity of dynamos

Based on the previous studies, we choose to classify the magnetic fields observed in spherical dynamo simulations into three main categories: axial dipoles, equatorial dipoles, and multipolar configurations, among which we also distinguish hemispherical and coherently oscillating dynamos.

Axial dipoles. Dipolar magnetic fields almost aligned with the rotation axis were shown to be stable, firstly, when the differential rotation between the poles and the equator is moderate and do not enter in conflict with the growth of large-scale axisymmetric components; and secondly, when the Coriolis force dominates the influence of turbulent inertia in a subsequent part of the convective bulk (e.g., Menu et al. 2020; Zaire et al. 2022, and references therein). Raynaud et al. (2015) showed that this second criterion is met when Ro_surf ≲ 0.1, with the surface (convective) Rossby number

$$\begin{array}{c}\hfill {\mathrm{Ro}}_{\mathrm{surf}}=\frac{v_{\mathrm{rms}}}{{\ell }_v\mathrm{\Omega }},\end{array}$$(4)

where v_rms and ℓ_v are the root mean square and the characteristic length scale of the nonaxisymmetric component of the surface convective velocity, respectively. In the covered parameter space, Ro_surf <  0.1 is achieved for weakly-magnetized models with moderate stratification (typically, up to N_ρ ≈ 3.5 and Pm ≈ 3). In the strong branch (i.e., Pm ≳ 10), Menu et al. (2020) and Zaire et al. (2022) showed that the dipolar branch can survive at higher values of Ro_surf since the Lorentz force starts playing an important role and competing turbulent inertia. Raynaud et al. (2015) also demonstrated that the stability domain of the dipolar branch is associated with a narrower and narrower range of values of Ra/Ra_c as the density stratification increases, so that it seems more and more difficult to generate magnetic dipoles. We show the surface map of the radial component of the magnetic field for such a dipolar model in Fig. 1a, with N_ρ = 3.5 and Ra/Ra_c ≈ 2.5. This snapshot exhibits relatively small-scale structures close to the surface, even if the dipole character remains obvious. The butterfly diagram in Fig. 2a clearly confirms the mean dipole behavior of the model over a magnetic diffusion timescale. At the very beginning of the time series, we can see the occurrence of a polarity reversal, which is common in such a model (e.g., Menu et al. 2020). As another diagnosis, we can also consider the time series of the $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ indicator, which is defined as the surface ratio of the dipole magnetic energy to the magnetic energy of the angular degrees between ℓ = 1 and a certain cut-off degree ℓ_cut, that is,

Fig. 1. Hammer projection of the radial component of the magnetic field at r = 0.9 ro and at a given time, for four typical dynamo models observed in direct numerical simulations. (a) Axial dipole with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac and Pm = 5. (b) Equatorial dipole with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. (c) Hemispherical dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 4 Rac, and Pm = 2. (d) Coherently oscillating dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 8 Rac, and Pm = 2.

Fig. 2. Surface butterfly diagram of the four typical examples shown in Fig. 1, i.e., the azimuthal average of the magnetic field radial component as a function of time and colatitude at the surface. (a) Axial dipole with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac, and Pm = 5. (b) Equatorial dipole with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. (c) Hemispherical dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 4 Rac, and Pm = 2. (d) Coherently oscillating dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 8 Rac, and Pm = 2.

$$\begin{array}{c}\hfill f_{\mathrm{dip}}^{{\ell }_{\mathrm{cut}}}={\left(\frac{{\sum }_{m=-1}^1\int {\mathit{B}}_{\ell =1,m}^2(r_{\mathrm{o}},\theta ,\phi )sin\theta \mathrm{d}\theta \mathrm{d}\phi }{{\sum }_{\ell =1}^{{\ell }_{\mathrm{cut}}}{\sum }_{m=-\ell }^{m=+\ell }\int {\mathit{B}}_{\ell ,m}^2(r_{\mathrm{o}}\theta ,\phi )sin\theta \mathrm{d}\theta \mathrm{d}\phi }\right)}^{1/2}.\end{array}$$(5)

This quantity is plotted for the present axial dipole in Fig. 3. Aside from the short reversal event at the very beginning of the time series, we can see that $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ remains higher than 0.6 for any value of ℓ_cut, as expected. A value of ℓ_cut = 12 is generally considered when tackling the question of the Earth’s core magnetism since it corresponds to the maximal spatial resolution accessible by observations (e.g., Christensen & Aubert 2006). Using the same value, we can see that $f_{\mathrm{dip}}^{12}$ drops on average around 0.5 during the reversal event, showing the partial transfer of energy from the ℓ = 1 components toward higher degrees. It is also worth considering smaller values for ℓ_cut. For instance, by choosing ℓ_cut = 5, we can see that $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^5$ remains above 0.5 even during the reversal event. The choice of ℓ_cut naturally appears important if one wants to define relevant diagnostics of the dipolar field strength in stars. This is true in particular when confronting theoretical models to observations, the latter being limited regarding the spatial resolution (see also discussion in Sect. 5.3).

Fig. 3. Dipolarity indicator as defined in Eq. (5) for the axial dipole considered in the panel Fig. 1a, with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac, and Pm = 5. The relative contribution of the dipole is computed for the cut-off degrees ℓcut = ℓmax, 12, and 5 (black, blue, and red lines, respectively).

Equatorial dipoles. In some cases, the generated magnetic field can also be structured as a large-scale dipole orthogonal to the rotation axis, that is, with the energy mainly contained in the ℓ = 1 and |m|=1 components; this is commonly called an equatorial dipole (e.g., Schrinner et al. 2014; Raynaud et al. 2014). The surface radial component of the magnetic field for a model exhibiting such a topology is plotted in Fig. 1b, with Ra = 2.48 Ra_c and N_ρ = 3. This configuration is generally observed for models close to the onset of convection and dynamo threshold, and at low values of the Rossby numbers, where it is favored by the columnar shape of the convective flow. This columnar structure can be discerned near the equator in Fig. 1b. While the magnetic axis can remain near the equator over time at low values of Ra/Ra_c, it may start oscillating around the equator over time as the value of Ra/Ra_c is progressively increased. This is actually the case for this model, as can be seen in the butterfly diagram plotted in Fig. 2b. The mean polarity alternates between positive and negative values and are always opposite in each hemisphere. To better reveal this oscillation, it is instructive to show the $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ indicator and the contribution of the axial component alone (i.e., m = 0), denoted with $f_{\mathrm{dip}}^{\mathrm{ax}}$; this is plotted in Fig. 4. As we can see, the value of $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ remains on average around 0.5. In contrast, the axial contribution $f_{\mathrm{dip}}^{\mathrm{ax}}$ oscillates between 0 and a maximum of about 0.4, the maximum being reached when the dipole axis tends to align with the rotation axis. Otherwise, $f_{\mathrm{dip}}^{\mathrm{ax}}$ is very low. Indeed, as Ra/Ra_c is increased from the onset of convection, an equatorial prograde jet progressively forms and enters at some point in conflict with the nonaxisymmetric part of the magnetic field. The system then leads to the generation of an axisymmetric dipole, which in contrast tends to kill the differential rotation. A nonaxisymmetric equatorial dipole is then rebuilt, whence the oscillations observed in the butterfly diagram. Increasing slightly again the value of Ra/Ra_c, the (axial) dipole branch discussed in the previous paragraph may be retrieved with adequate initial conditions, especially for large values of Pm. Otherwise, the reinforcement of the differential rotation level as Ra/Ra_c increases tends to disrupt the emergence of equatorial dipoles and favors multipolar morphology, which we now discuss (e.g., Raynaud et al. 2015).

Fig. 4. Same as in Fig. 3, but for the equatorial dipole considered in Fig. 1b, with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. The relative contribution of the axial dipole alone (m = 0) for the cut-off degrees ℓcut = ℓmax, denoted with $f_{\mathrm{dip}}^{\mathrm{ax}}$, is also represented in green.

Multipolar dynamos. When the conditions to sustain a dominant large-scale dipolar magnetic field are not met, the possible magnetic topologies are generally encapsulated into the complementary branch, that is, the so-called multipolar branch (e.g., Kutzner & Christensen 2002; Christensen & Aubert 2006; Goudard & Dormy 2008; Gastine et al. 2012; Schrinner et al. 2014, and references therein). We distinguish here two sub-classes in this branch. First, we have hemispherical dynamos for which a given hemisphere is significantly more magnetized than the other one. This kind of models is generally observed for moderate values of Ra/Ra_c (Schrinner et al. 2014; Raynaud et al. 2014; Raynaud & Tobias 2016). As an illustration, we plot in Fig. 1c the surface radial component of the magnetic field for such an hemispherical model, with Ra/Ra_c = 4 and N_ρ = 6. As this figure suggests, this configuration is the result of the superposition of a symmetric and an antisymmetric dynamo modes with respect to the equator that share similar amplitudes, leading to the global canceling of the magnetic field in one given hemisphere. We can see in the butterfly diagram plotted in Fig. 2c that this cancellation is maintained over about one magnetic diffusion timescale. In the magnetized hemisphere, the dynamics seems noisy close to the pole but oscillating features tend to emerge and survive very close to the equator, somehow overshadowed by the antisymmetric mode. Such equatorial dynamo waves remarkably look like the structures observed in the second classes of the multipolar branch, that is, the so-called coherently oscillating dynamos. An example of models exhibiting such a behavior is plotted in Fig. 1d, for N_ρ = 6 and Ra/Ra_c = 8. This kind of dynamo is generally obtained for large values of Ra/Ra_c. In this case, the convective velocity fluctuations are large enough and the inertia dominates the Coriolis force in a substantial part of the convective bulk so that any dipolar component would collapse. As shown by Raynaud et al. (2015), this seems to be met at lower and lower values of Ra/Ra_c as N_ρ increases. As seen in Fig. 1d, the dynamics is dominated by small-scale structures. Nevertheless, the butterfly diagram in Fig. 2d shows that the magnetic structures migrate coherently on average from the equator to the poles in both hemispheres. We note that the migration tends to be quasi symmetrical with respect to the equator, even if some asymmetries persist. These oscillatory dynamo modes have been shown to be well explained by the Parker’s wave formalism to a good approximation (Parker 1955; Yoshimura 1975; Gastine et al. 2012; Warnecke et al. 2014; Duarte et al. 2016; Warnecke 2018). Nevertheless, even if this attractive scenario has long been proposed to explain the origin of the solar cycle, it may not explain all the complexity of stellar magnetic cycles (e.g., Beaudoin et al. 2016; Raynaud & Tobias 2016; Strugarek et al. 2017, 2018; Brandenburg et al. 2017; Beer et al. 2018; Brun et al. 2022; Ortiz-Rodríguez et al. 2023).

3.2. Magnetic quenching of the differential rotation

One of the most important properties of thermal convection in rotating shells is the generation of a large-scale differential rotation (Christensen 2002; Aurnou et al. 2007; Gastine et al. 2012, 2013; Guerrero et al. 2013; Brun et al. 2017) which can play an important role in the dynamo process (e.g., Schrinner et al. 2011, 2012; Brun & Browning 2017; Käpylä et al. 2023). In Paper I, the outgoing surface heat flux has been shown to remain uniform when the influence of turbulent inertia on the convective motions dominates the Coriolis acceleration close to the surface, that is, when Ro_surf ≳ 1, which is usually met for large values of Ra. In contrast, when Ro_surf ≲ 1, the Coriolis acceleration can significantly affect the surface turbulent convective rolls and couple the heat flux with the surface differential rotation. In this case, the presence of a strong prograde equatorial jet at the surface is associated with a darker equatorial band.

Before addressing the impact of the magnetic fields on the surface heat flux distribution, it is thus worth questioning its impact on the differential rotation. To illustrate this point, we consider models with N_ρ = 6, Ra/Ra_c = 8, Pr = 1, and E = 3 × 10⁻⁴, and plot their azimuthally averaged zonal velocity in Fig. 5. First, in Fig. 5a, we consider the hydrodynamic case without magnetic field. We see that the mean zonal flow exhibits a cylindrical geometry, as expected from the Taylor-Proudman constraint. The outer column is prograde and the inner cylinder embedding the inner core is retrograde. Indeed, the convective Rossby number remains smaller than unity in a substantial part of the convective bulk and the turbulent-mixed superficial layers are not able to evince the resulting positive angular momentum deposition induced by Reynolds stress in the external layers (Aurnou et al. 2007; Gastine et al. 2013). In turn, the inner part has to rotate in the opposite sense in order to conserve the total angular momentum. Then, accounting for dynamo magnetic fields, we are left with a new control parameter, which is the magnetic Prandtl number Pm. In Fig. 5b, we show the result obtained with Pm = 1. This model exhibits a coherently oscillating multipolar topology, similar to that illustrated in Fig. 1d. We see that the magnetic Lorentz force has two main effects on the mean zonal flow. First, it breaks the cylindrical geometry of the flow: the prograde region extends in deeper layers close to the equator while its latitudinal extension is reduced compared to the case without magnetic fields, and the retrograde region takes a conical shape. At the surface, the development of small-scale surface zonal jets at high latitudes is favored; these are very similar to the alternation of retrograde and prograde bands observed in Jupiter or Saturn, which have already been reproduced in numerical simulations (e.g., Gastine et al. 2014a; Heimpel et al. 2016, 2022; Christensen et al. 2020; Yadav et al. 2022). Second, the magnetic field clearly reduces the intensity of the zonal mean flow (by a factor of about 3 for Pm = 1, see the colorbar). In Figs. 5c and d, we plot the results for two other models computed by Pm = 2 and 3. We note that the model with Pm = 2 is the same as the one considered in Figs. 1d and 2d. When increasing the value of Pm, we actually increase the intensity of the magnetic field. Indeed, at a given value of Ra/Ra_c, we expect the Reynolds number Re = v_rmsd/ν, with v_rms the root mean square of the convective fluctuations of velocity, to remain about constant when varying Pm. In turn, when increasing Pm, the magnetic Reynolds number, Rm = v_rmsd/η = Pm Re, increases. In other words, the effect of the Ohmic dissipation decreases and the dynamo-induced magnetic fields have larger magnitudes. We can see in Figs. 5c and d that increasing the magnetic intensity leads to a decrease in the differential rotation. This clearly indicates that magnetic fields globally quench the differential rotation compared to the case without magnetic fields. We check that this trend holds for the set of models we have considered in this work, as we subsequently show in detail in Sect. 4. This result extends the conclusions of Yadav et al. (2016) and Käpylä et al. (2017) obtained for Boussinesq and moderately stratified models, respectively, to highly-stratified models and a wider parameter space: whatever their morphology, the self-sustained magnetic fields tend to globally quench the large-scale differential rotation in typical thick convective envelopes of cool stars. Since the dynamo is triggered at higher values of Pm when N_ρ is increased (Raynaud et al. 2015), we also expect that the higher the density stratification, the larger the value of Pm for magnetic fields to start quenching the differential rotation.

Fig. 5. Meridional slices of the azimuthally averaged zonal velocity, for different values of the magnetic Prandtl number Pm, in unit of the viscous velocity scale ν/d. Models with Pr = 1, E = 3 × 10−4, Nρ = 6, and Ra = 8 Rac. (a) Hydro, (b) Pm = 1, (c) Pm = 2, (d) Pm = 3.

4. Effect of self-sustained magnetic fields on the global surface heat flux distribution

After having reviewed the different magnetic topologies expected in stratified dynamo models, we now investigate the global effect of magnetic fields on the mean latitudinal distribution of luminosity at the surface of cool stars, with thick convective envelopes. To that end, we compute the mean latitudinal Nusselt number in Eq. (3) for each simulation. We go step-by-step and first tackle the case of models close to the onset of convection before dealing with more turbulent models, which are more representative of stellar conditions. The results obtained from our set of simulations are summarized through the figures presented in Appendices C and D, showing the mean Nusselt and large-scale zonal flow profiles, respectively. We note that for most of our models, the flow is rather turbulent close to the surface especially for N_ρ >  3, so that the mean Nusselt profile may fluctuate over time and appear noisy in case the time average is performed on only few snapshots (e.g., see time variations in Fig. 2 of Paper I). We checked that averaging on few dozens of snapshots taken at fixed intervals over a magnetic diffusion timescale is sufficient to grasp the mean brightness and velocity profiles. This is systematically performed for each model. Typical standard deviations around the mean value of the Nusselt profile depends on the latitude but we expect through a visual inspection in few models that it is generally lower than 10−15%; the mean brightness distribution as computed in this work thus appears relevant for a first step investigation.

4.1. Close to the onset of convection

At the onset of convection, the first unstable modes take the form of columns aligned with the rotation axis. In the hydrodynamic case, it was shown in Paper I that the surface heat flux is preferentially conveyed close the equator. As the value of N_ρ increases at a given value of Ra/Ra_c, the convective cells get more and more concentrated near the surface layers, and the heat flux tends to be maximum near the tangent cylinder² while keeping an equator brighter than the poles. A typical example of the surface Nusselt profile in this case is shown by the solid black line in Fig. 6a, for a model with N_ρ = 3, E = 3 × 10⁻⁵, Pr = 1, and Ra = 1.7 Ra_c; the value of Pm is equal to 0.5, which is too low to trigger a magnetic dynamo (see Sect. 3).

Fig. 6. Mean latitudinal profile of the surface Nusselt number defined in Eq. (3), for three typical models. (a) E = 3 × 10−5, Pr = 1, Nρ = 3, and Ra ≈ 1.7 Rac close to the onset of convection (b) E = 3 × 10−5, Pr = 0.3, Nρ = 4, and Ra ≈ 3 Rac at the transition toward turbulent flows (c) E = 3 × 10−4, Pr = 1, Nρ = 6, and Ra ≈ 8 Rac. Left: variations with the magnetic Prandtl number Pm. Solid black lines represent either pure hydrodynamic models or MHD models without magnetic dynamo (i.e., too low values of Pm). Right: spherical representation for the nonmagnetized model (i.e., the lowest value of Pm) and the most magnetized dynamo model (i.e., the highest value of Pm) in each case. Note that we use different colorbars in each row to enhance the visualization contrast.

We now investigate the effect of magnetic fields on the mean surface heat flux distribution, remaining first in a regime close to the onset of convection. In this case, the magnetic Prandtl number, Pm, is varied above the minimal critical value needed to trigger a self-sustained dynamo while attempting to cover both the dipolar (axial and equatorial) and multipolar branches (for given E, Ra/Ra_c, N_ρ, and Pr). As explained in Sect. 3.2, this in addition allows us to quantify the impact of the intensity of the magnetic field on the surface Nusselt profile; the larger Pm for given parameters, the larger the magnetic intensity. For the previous illustrative model in Fig. 6a, we can clearly see that close to the onset of convection, magnetic fields generally tend to increase, on average, the integrated surface Nusselt number. In particular, the poles and the regions near the tangent cylinder can become brighter than in the hydrodynamic case as Pm increases. This is a systematic trend observed in several models close to the onset of convection, as we can see in Figs. C.1a,b, C.2a–d, and C.3a, for values of N_ρ between 3 and 4 and Ra/Ra_c varying up to a value of 4. In all these figures and in the subsequent sections of the paper, we note that the models represented by solid black lines are either purely hydrodynamic or without dynamo (i.e., with a value of Pm too low to sustain a magnetic field). In this regime, a surface equatorial prograde jet starts developing but it is inhibited by magnetic fields, as expected from Sect. 3.2. For example, this can be seen in Figs. D.1a and D.2a for two models with N_ρ = 3.5 and 4, Ra ∼ 2.5 Ra_c. Similarly to the hydrodynamic case (see Paper I), the surface heat flux close to the onset of convection thus appears not to be anticorrelated with the surface differential rotation, but to be a property of the first unstable convective modes, in contrast with more turbulent models as we see in the next sections.

4.2. Transition toward turbulent flows

Continuously increasing the Rayleigh number from the onset of convection progressively reinforces the equatorial prograde jet. As already explained, this jet results from the Reynolds stress induced by coherent convective rolls, and that can redistribute positive angular momentum to the mean flow in the outer layers. This typical enhancement of the surface differential rotation as the Rayleigh number is increased from the onset of convection toward turbulent flows is observed in most of our models with no dynamo (e.g., see Figs. D.1–D.3a–c, solid black lines). When the equatorial jet is strong enough, it progressively inhibits the outgoing heat flux close to the equator, making at some point the equator as bright as or even darker than the poles, as expected from Paper I. A typical example of the surface Nusselt profile in this situation is shown in Fig. 6b (solid black line) for a model with E = 3 × 10⁻⁵, Pr = 0.3, N_ρ = 4, and Ra ≈ 3 Ra_c. When adding magnetic fields, we expect according to Sect. 3.2 the surface differential rotation to be quenched by the Lorentz force. For this example model, Fig. 7 indeed shows that the prograde equatorial jet slows down as the magnitude of the magnetic field (i.e., the value of Pm) increases. As seen in Fig. 6b, the equator concomitantly becomes brighter than the poles when going from Pm = 0.3 to Pm = 0.5 (black and blue lines). Increasing Pm to a value of two in a second step (red line), we can see that the heat flux tends to increase in the region inside the tangent cylinder, producing a slight dip in the vicinity of the equator. The change in the Nusselt profile resulting from this second increase in Pm is very similar to the trends observed for models close to the onset of convection and depicted in Sect. 4.1. As a matter of fact, in weakly turbulent convection, the emergence of dynamo magnetic fields can retroactively lead to the laminarization of the flow (e.g., Petitdemange 2018). Therefore, the properties of the convective shell are then expected to be very similar to those of a model with a slightly lower value of the Rayleigh number, that is, close to the linear onset of convection. For the models with Pm = 0.5 and Pm = 2 represented in Fig. 7, we can check this interpretation holds valid in a simple way by considering the model in Fig. C.2d obtained with the same parameters but with a slightly lower value of Ra/Ra_c. Although the mean value of the Nusselt number remains naturally larger for the models with the largest value of Ra, as expected from simple scaling laws, we see that the shape of the profiles are very similar in both figures, even between the nonmagnetic and the magnetic models at Pm = 0.5, hence supporting the laminarization effect.

Fig. 7. Mean surface azimuthal velocity as a function of the colatitude, for the models in Fig. 6b, i.e., with E = 3 × 10−5, Pr = 0.3, Nρ = 4, and Ra = 3 Rac. The magnetic Prandtl number Pm is varied and the solid black line represents the model with no dynamo (i.e., too low values of Pm).

Our set of simulations comprises other examples of such inversion. A clear example can be found in Fig. C.1c for E = 3 × 10⁻⁵, Pr = 1, N_ρ = 3, and Ra ≈ 3.7 Ra_c. In this figure, the model with Pm = 1 exhibits an hemispherical dynamo while the model with Pm = 3 has a dominant dipole component. The laminarizing effect of the magnetic field is again suggested when comparing this figure with the models in Fig. 6a for Ra ≈ 1.7 Ra_c. Another example can be found in Fig. C.4a, for a model with E = 3 × 10⁻⁴, Pr = 1, N_ρ = 6, and Ra ≈ 3 Ra_c, exhibiting an oscillatory dynamo. We note that this model with such a low value of Ra/Ra_c requires values of Pm larger than about 6 to trigger a magnetic dynamo. We thus conclude that, whatever the magnetic topology, magnetic fields can globally induce an inversion of the pole-equator brightness contrast, going from a darker to a brighter equator as the magnetic intensity increases. This generally happens for models at the transition between the onset of convection and more turbulent flows.

4.3. Turbulent models

We finally consider turbulent models with values of Ra much larger than the threshold value above which the equator becomes darker than the poles (i.e., typically Ra ≳ 5 Ra_c). As an illustration, we plot in Fig. 6c the Nusselt profile of such models with E = 3 × 10⁻⁴, Pr = 1, N_ρ = 6, and Ra ≈ 8 Ra_c. We remind that the mean azimuthal velocity of these models is represented in Fig. 5. The associated latitudinal profiles at the surface are also plotted in Fig. D.3c for more details. As already shown in Paper I and in Sect. 4.2 for less turbulent models, we can perceive the strong correlation between the prograde jet and the dark band close to the equator in the hydrodynamic case. Dynamo models are then compared in the same figures. We can see in Fig. 6c that increasing Pm and thus the magnetic field intensity tends to flatten the pole-equator brightness contrast. This is actually a general trend that is systematically observed in our turbulent models. This is summarized in Figs. C.2f, C.3b–j and C.4b–d, for values of N_ρ between 3.5 and 6, and values of Ra/Ra_c between about 5 and 30. For such values of Ra/Ra_c, the models generally exhibit a multipolar magnetic structure. This result suggests that far enough from the onset of convection, the magnetic fields globally tend to reduce the surface brightness contrast between the poles and the equator. This effect can be significant: if Pm is large enough, a luminosity contrast between the equator and the pole of about 150% such as depicted in Paper I in the hydrodynamic case can be totally smoothed out.

Similarly to the transition models discussed in Sect. 4.2, the global reduction of the pole-equator brightness contrast depicted in Fig. 6c is clearly correlated with the quenching of the surface differential rotation by the magnetic fields, as we can see Fig. D.3c. This is again a systematic trend observed in all the turbulent models of our sample. Several examples can be found for instance in Figs. D.1c, D.2b–i and D.3b–d, for values of N_ρ between 3.5 and 6, as well as values of Ra/Ra_c between about 5 and 20. For the models with no dynamo in these figures (solid black lines), the prograde equatorial jet reinfoces as Ra increases, as noted before. In contrast, increasing Pm in all these figures results in the reduction of the large scale differential rotation, which is directly responsible for the homogenization of the surface heat flux in this regime.

When increasing the value of Ra/Ra_c toward very large values (i.e., typically Ra/Ra_c ≳ 30 in our sample), the effect of inertia on the flow then starts dominating the Coriolis force at some point, reducing ipso facto the angular momentum transfer from the small scales to the mean flow and hence, the magnitude of the surface differential rotation. In order to conserve the total angular momentum, the system has no choice but to generate a retrograde equatorial zonal jet at the surface, with a sharp transition between both regimes (e.g., Gastine et al. 2013). In the hydrodynamic case, we have shown in Paper I that the mean surface brightness contrast tends to decouple from the differential rotation and homogenize as the surface convective envelopes become more and more turbulent. Exploring this kind of dynamo models in the present parametric study is difficult since they require long numerical integration. We census only eight models around the transition between a solar (prograde) and an antisolar (retrograde) differential rotation in our set of simulations. First, we have two models with E = 10⁻⁴, Pr = 1, N_ρ = 5, and Ra = 32 Ra_c, whose surface azimuthal velocity profiles are plotted in Fig. D.2j. A quick comparison with the nonmagnetic models at lower values of Ra/Ra_c in Figs. D.2f–i (black lines) suggests that these models are very close to the transition: indeed, we can see in Fig. D.2j that the equatorial jet starts decreasing slightly around Ra = 32 Ra_c in the hydrodynamical case (black line). We see in the same figure that the dynamo magnetic fields also quench the differential rotation in this case (blue line), without any clear differences with less turbulent models. Second, we have six other models with E = 10⁻⁴, Pr = 1, N_ρ = 4, and Ra ≈ 29 Ra_c, and with E = 3 × 10⁻⁴, Pr = 1, N_ρ = 6, and Ra ≈ 32 Ra_c, the surface azimuthal velocity profiles of which are plotted in Figs. D.2e and D.3e, respectively. In both figures, we can see that the models without magnetic fields (solid black lines) exhibit a solar differential rotation. Regarding the dynamo models computed with the same parameters, we see either that the zonal equatorial jet is quenched or that a strong antisolar jet develops when adding magnetic fields; we nevertheless do not depict any clear dependence on the value of Pm. Actually, the transition between a solar and an antisolar differential rotation is known to be hysteretic in the hydrodynamic case and thus depends on the initial conditions (e.g., Gastine et al. 2014b). In contrast, Karak et al. (2015) gave some evidence that dynamo magnetic fields may remove the flow bistability in this regime. In this context, it is thus difficult to explain our results through a simple trend as a function of Pm. In any case, the surface Nusselt profiles of these eight models displayed in Figs. C.3e,j and C.4e demonstrate that the surface brightness contrast is always reduced compared to the hydrodynamic case. This indicates that the main trend observed for slightly less turbulent models with a solar differential rotation regime also holds true for more turbulent models with an antisolar differential rotation. We nonetheless acknowledge that the transition regime toward an antisolar differential rotation deserves a dedicated study (e.g., Simitev et al. 2015; Noraz et al. 2022a,b), which is beyond the scope of this paper.

5. Discussion

In this section, we discuss the results of this work in view of the past studies and observational context.

Fig. 8. Maximum brightness contrast as a function of the mean Elsasser number, i.e., a proxy of the mean magnitude of the magnetic field, computed for models with E = 3 × 10−4, Nρ = 6 (stars), E = 10−4, Nρ = 5 (diamonds), and E = 10−4, Nρ = 4 (squares); the color codes the Rayleigh number. For all these models, Pr = 1.

5.1. Magnetic brightness flattening: Applicability and limitations

Our parametric study demonstrates that magnetic fields and the mean brightness contrast between the poles and equator in rotating cool stars (i.e., with relatively thick turbulent convective envelopes) are tightly related. This was shown to result from the systematic quenching of the large-scale differential rotation by the magnetic fields (see Sect. 3.2), which in turn affects the convective transport of energy, as previously expected from Yadav et al. (2016) and Paper I. The larger the magnetic intensity, the lower the differential rotation, and the lower the surface heat flux variations at the surface. Moreover, this turns out to hold true at leading order whatever the magnetic topology. This is illustrated in Fig. 8 for some of the most stratified models of our sample with Pr = 1, N_ρ between 4 and 6, and different values of Ra/Ra_c; the mean surface Nusselt and zonal flow of these models are plotted in Figs. C.4 and D.3 for N_ρ = 6, in Figs. C.3h and D.2h for N_ρ = 5, and Figs. C.3d and D.2d for N_ρ = 4. In Fig. 8, the relative difference between the maximum and minimum values of the surface Nusselt number (normalized by its mean value) is plotted as a function of the global Elsasser number $\mathrm{\Lambda }=\overline{B^2}/(\mathrm{\Omega }{\rho }_{\mathrm{c}}\mu \eta )$. Since the surface Nusselt profile can sometimes exhibit localized oscillating variations around the mean trend, its maximum and minimum values are taken by computing a rolling average over a latitude interval of 10° beforehand.

First, Fig. 8 shows that the brightness contrast tends to increase as N_ρ decreases whatever the value of Λ, as expected from Paper I. Indeed, as N_ρ increases, the surface convective Rossby number tends to increase so that the turbulent standardization of the surface heat flux is more efficient (e.g., see Figs. 5 and 7 in Paper I). The trend observed in Paper I in the hydrodynamical case is thus shown to hold valid when magnetic fields set in. Second, focusing on N_ρ = 6 models only, we can observe a nonmonotonic behavior as a function of Ra/Ra_c in the nonmagnetic case (i.e., Λ = 0); the brightness contrast first increases from Ra/Ra_c = 3 to Ra/Ra_c = 4, before decreasing for larger values of Ra/Ra_c. This change of regime has also already been observed in Paper I. It has been shown to result initially from the progressive reinforcement of the prograde equatorial jet and, then, the progressive homogenization of the heat transport as turbulence is enhanced. We have checked that the surface Rossby number of these models increases from about 0.1 to about 2 as Ra/Ra_c increases from 3 to 32, that is, covering the above-mentioned transition from the onset of convection toward the start of the convective homogenization (e.g., see Fig. 7 in Paper I). Besides, in the dynamo case, we do distinguish two regimes already depicted in Sect. 4:

• On the one hand, the model with N_ρ = 6 and Ra/Ra_c = 3 exhibits a typical reversal of the pole-equator brightness contrast in presence of magnetic fields, resulting in an increase in its value as Λ increases. Such a case is associated with models close to the linear onset of convection and, although being of interest from a theoretical point of view, it is likely not to be applicable to very turbulent stellar convective envelopes.

• On the other hand, for more turbulent models, we clearly see the systematic decrease of the surface brightness contrast as Λ increases. This relation between the magnetic field intensity and the global brightness contrast, although counterintuitive at first, has been demonstrated to be physically grounded and could have implications for interpreting the combination of photometric and spectropolarimetric data.

From a quantitative point of view, the surface heat flux contrast predicted by our models in Fig. 8 can reach an upper limit of about 100% of the mean value in the hydrodynamic case (see also Paper I). In order to interpret such differences, it is worth translating them in terms of temperature variations. By construction and for sake of simplicity, a constant temperature (or entropy) is imposed at the surface of our models. A simple way to circumvent this limitation is to use the Stefan-Boltzmann’s law at the outer sphere to link the relative variations in the Nusselt number with those in temperature, that is, ΔNu/Nu ∼ 4 ΔT/T. The surface temperature deviations are therefore expected to reach about 25% at most in our simulations. As an illustration, we know that the temperature contrast between the quiet photosphere and cool nonaxisymmetric stellar spots in G, K, and M dwarfs can reach up to about 20% (e.g., Solanki 2003; Afram & Berdyugina 2015; Garcia et al. 2022). In absence of additional constraints on the mean axisymmetric brightness distribution in thick stellar convective envelopes (see Sect. 5.2), this comparison, although simplistic, does not rule out the somehow large surface temperature and brightness contrasts predicted by the present set of simulations; the large relative Nusselt deviations predicted and represented in Fig. 8 thus do not appear untenable.

Nevertheless, before going further in the comparison, it is necessary to discuss in more detail the applicability and limitations of these findings in an astrophysical context. Indeed, the parameter space covered by our simulations is ineluctably limited and far from stellar regimes for numerical reasons. As a first attempt to compare simulations and actual stars, we usually rely on the predominant force balance. In all our models, the geostrophic balance (i.e., between the pressure gradient and the Coriolis acceleration) is globally met at large scales. Regarding the convective fluctuations around this global state, we have seen that the local Rossby number, Ro_surf, lies between 0.1 and about unity at the surface, and thus is even smaller in deeper layers where it decreases (e.g., see Fig. 5 in Paper I). Therefore, the Coriolis acceleration is mainly responsible for the dynamical properties of convection, even if buoyancy starts acting in the surface layers for the most turbulent models. Using such considerations, we can attempt to define a first applicability domain: cool stars with relatively thick convective envelopes and Ro_surf ≲ 1. This actually corresponds to the range of values where the contrast in the brightness distribution is expected to be the largest in the hydrodynamic case, with predicted contrasts reaching up to 150% (e.g., see Fig. 7 of Paper I). According to observations and stellar modeling, this low Rossby number regime should be met in most of the K and M dwarfs (e.g., Wright et al. 2011, 2018) as well as pre-main sequence stars (e.g., Landin et al. 2023). However, this is not met for the other usual type of cool stars, namely evolved red giant stars. Indeed, using asteroseismic data, we know for such stars that the surface granulation turnover timescale τ_conv ≲ 1 day and that the mean rotation period, P ≳ 10 days (e.g., Samadi et al. 2013; de Assis Peralta et al. 2018; Gehan et al. 2018), which permits us to estimate a lower limit for the surface local Rossby number, that is, Ro_surf ∼ P/τ_conv ≳ 10. Such values are out of the parameter space considered in the present paper. Besides, it is also worth mentioning RS CVn close binary stars in which the red-giant component can rotate very rapidly. However, the presence of strong tides adds important dynamical effects that go beyond the framework of this work. More detailed investigations are needed in order to tackle the regime of large Rossby number, which will demand to study the transition from a solar to an antisolar differential rotation and very turbulent flows. According to Paper I and this work, we expect a rather uniform mean brightness distribution as the turbulence increases at the surface and the heat flux decouples from the differential rotation. We recall that this agrees with the Sun whose surface Rossby number is quite large and photosphere is uniform at leading order, regardless of the strong prograde differential rotation. Of course, another current question is to understand the formation and persistence of such strong prograde differential rotation in such regimes, but this is out of the scope of this paper. As a consequence, to tackle the regime of large Rossby numbers, it will be necessary to go beyond the large-scale and mean diagnostic used in this work and to account for other higher-order effects such as the impact of small-scale and time-variant surface magnetic features.

5.2. Observational tests and potential implications

In order to test our theoretical results in a direct way, we need observational targets for which we are able to measure both the surface magnetic fields and brightness distribution in addition to their fundamental properties (i.e., mean density, rotation rate, luminosity). Zeeman-Doppler imaging techniques permit to map magnetic fields, but it can only give insights on nonaxisymmetric brightness structures (e.g., Kochukhov 2016). In that context, interferometric imaging seems the most adequate way to map directly the mean axisymmetric brightness distribution of stars. Recent instruments such as GRAVITY (e.g., GRAVITY Collaboration 2017) can provide an angular resolution of about 1 mas on order of magnitude. While such precision enables us to estimate the apparent radii of the closest K or M dwarfs (e.g., Kervella et al. 2017), it does not allow to map the surface of these stars directly. Indeed, assuming typical radii of about 0.5 R_⊙ (e.g., Chabrier & Baraffe 1997), resolving for instance a fifth of the stellar disk at least with this precision requires parallaxes larger than about 1″. For pre-main sequence stars, assuming typical radii of about 2 R_⊙ (e.g., D’Antona 2017), this would still require parallaxes larger than about 300 mas. Both estimated parallaxes are unrealistic, but simple checks in the available astrophysical database (e.g., the Gaia database) show that increasing the angular resolution by a factor of three would permit to have access to an interesting sample of potential targets. This sample is composed of one pre-main sequence star of about 20 Myr (i.e., V* AP Col), and of about twenty K and M dwarfs, which is quite promising for forthcoming projects. We can note especially the recent instrument CHARA/SPICA at Mount Wilson that is expected to improve the accessible angular resolution by a factor of two at least in a very near future (e.g., Mourard et al. 2017; Ligi et al. 2023).

While direct tests through interferometry remain currently challenging, indirect measurements could also help in probing the (magnetic) Coriolis darkening effect. One of this indirect method is to compare the spectral energy distribution of stars with a template assuming a uniform surface brightness. Indeed, the differences in temperature in the dark and bright surface bands should induce local differences in the spectral lines and thus deviations in the total integrated spectrum. Such method has already been used to probe the gravity darkening in fastly-rotating intermediate-mass stars (e.g., Girardi et al. 2019). While the effect is certainly subtle and the result depends on the model of darkening considered, analyzing a large sample of cool stars of similar mass with different inclination angle can be envisaged to lift degeneracies (e.g., Zorec 2023). Another indirect method to test the Coriolis darkening consists in observing a companion body passing in front of the star, and to use the photometric and spectroscopic measurements of this event to infer the presence and properties of the star’s surface axisymmetric structure. Indeed, assuming an equator much darker than the poles at the stellar surface, as expected from our simulations, the transit depth in the observed light curve is expected to vary with the transit phase, with a minimum occurring when the companion obscures the dim stellar equatorial band. For that purpose, transiting close-in exoplanets are excellent candidates (e.g., Morris et al. 2017) as these bodies are dark compared to the star and have a small relative size (in opposition to stellar companions), enabling the discrete occultation of finer details on the host star photosphere. Nonetheless, the detection of such transit events will only be possible under stringent orbital configurations: it requires the orbital plane of the transiting exoplanet to be sufficiently misaligned with the equatorial plane of the star to produce a change in the apparent brightness that differs from the change expected during the transit of a simple limb-darkened body. Such a configuration, and its resulting light curve, are illustrated in Fig. 9. In order for the detection of this structure to be significant, it will be interesting to combine photometric measurements with time-series spectroscopic measurements that are able to encapsulate Doppler effects due to the surface differential rotation of the star. While this could be achievable through the modeling of the Rossiter–McLaughlin effect (Rossiter 1924; Ohta et al. 2005), the extent to which a proper model comparison can be made in order to decipher an equatorial jet from a simple stellar spot occultation is under investigation, and will be featured in a future study.

Fig. 9. Transit light curve of an exoplanet occulting a star with an equatorial jet dimmer than the rest of its photosphere. This simulation made use of starry (Luger et al. 2019). It features a limb-darkened star with an adimensional unit radius and an occulting spherical body of relative radius 0.15. The equatorial band has a contrast of 1 (for illustration purposes) and an angular radius of 0.2 radians. The upper plot shows the system at scale for six different times. The bottom plot shows the occultation light curve (black line), the differential flux measured for the six times shown above (points), as well as the transit signal that would be observed for a uniform limb-darkened star (gray).

In case the spin-orbit angle is close to 90°, that is when the transit chord of the planet lies within the equatorial band, it is worth mentioning that a significant pole-equator brightness contrast as predicted in this work may also have indirect implications on the interpretation of exoplanets transmission spectroscopy measurements. During a transit, the light coming from the star is transmitted through the planet’s atmosphere. Owing to the wavelength dependence of its atmosphere opacity, the apparent radius of the planet inferred via the analysis of a transit can vary as a function of the observed wavelength, leading to a transmission spectrum that contains precious information about the planet composition (e.g., Kreidberg et al. 2014). In practice, this technique often assumes that the transited chord is representative of the rest of the stellar photosphere. However, any difference between this chord and the unocculted stellar disk results in the contamination of the planet’s transmission spectrum by this unocculted part. The low effective temperatures of cool stars lead to stellar atmospheres containing complex molecules (and clouds for the latest types) that are hard to separate from the features of cool exoplanets atmospheres (e.g., Rackham et al. 2018). For this reason, stellar contamination is a major concern for the study of exoplanets atmospheres around late M dwarfs, and a great obstacle to the characterization of rocky exoplanets envisioned for the next decades (Rackham et al. 2023). Given the large-scale structures resulting from the magnetic Coriolis darkening presented in this work, this effect results in a greater level of transmission spectra stellar contamination, with the worst-case scenario being when the spin-orbit angle is equal to 90°. In this case, the transit chord is everywhere different from the rest of the stellar disk, with no observable capable of identifying the presence of an equatorial band. A quantitative analysis of the impact of such structures on the interpretation of exoplanets transmission spectra is out of the scope of this paper, but certainly deserves further investigation.

5.3. On the effect of the magnetic morphology

According to our work, the magnetic morphology does not affect the leading-order magnetic flattening of the surface brightness distribution. Nevertheless, going beyond this leading-order trend, there are clues that the magnetic field structure may impact this distribution. A first example at the transition between the onset of convection and more turbulent models can be found in Fig. C.1c for N_ρ = 3 and Ra = 3.7 Ra_c. In this figure, the dynamo for Pm = 1 is hemispherical with a more magnetized south pole than the north one, while the dynamo is dipolar for Pm = 3. We can see that while the heat flux is symmetric with respect to the equator for the dipolar model, it becomes somehow asymmetric in the hemispherical model with a brighter north hemisphere. Another striking example is given in Fig. C.4b for N_ρ = 6 and Ra/Ra_c = 4. These dynamo models are also hemispherical with a more magnetized south hemisphere, as we can see in Figs. 1c and 2c for Pm = 2. Inspecting the surface Nusselt profile, we can see that although magnetic fields reduce the global brightness contrast, the north pole also becomes brighter than the south pole. The concentration of the magnetic energy in the south seems to inhibit the heat flux in this hemisphere. It is also important to note that this south-north deviation is not directly correlated with the differential rotation (see Fig. D.3b), supporting the magnetic field topology as the responsible of this trend. The magnetic impingement of the brightness contrast still holds from a global point of view, but it is clear that these higher-order effects remain to be addressed properly in the future. This point also encompasses the impact of potential polar spots, which are already known to be able to play a significant role (e.g., Yadav et al. 2015a).

In the same spirit as in the previous section, it is also worth discussing the way to test the relation between the magnetic topology and the surface heat flux distribution. To do so, a relevant indicator is necessary to compare the predictions of dynamo simulations and observations. The simplest diagnostic consists in classifying the magnetic topology into the two main families, that is, large-scale dipoles or small-scale multipoles. For this purpose, we usually rely on the mean dipole field strength indicator, which is the time average of $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ defined in Eq. (5) (e.g., Christensen & Aubert 2006). Above a value of 0.5, we merely consider that the topology is dipolar; otherwise, it is multipolar. The value of this indicator directly depends on the value of the cut-off degree, ℓ_cut. For sake of the consistency, ℓ_cut is generally set to a value corresponding to the maximum spatial resolution reachable with the considered observational techniques, that is, between about 5 and 12 for Zeeman-Doppler reconstruction (e.g., Yadav et al. 2015b; Zaire et al. 2022). As seen in Figs. 3 and 4, changing ℓ_cut from 5 to 12 can decrease the value of $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ by about 10%. For the dipole model in Fig. 3, the diagnostic does not depend on the value of ℓ_cut globally, except close to polarity reversal events where $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ can be either greater or smaller than 0.5. Without the knowledge of the long-term evolution of $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$, the diagnostic can thus remain ambiguous, especially if the measurement is performed close to a polarity reversal (e.g., Petit et al. 2009). As another example, the equatorial dipole considered in Fig. 4 exhibits $f_{\mathrm{d}\mathrm{i}\mathrm{p}}^{{\ell }_{\mathrm{c}\mathrm{u}\mathrm{t}}}$ very close to the threshold 0.5, because of the substantial contribution of the small-scale dynamics superimposed to the dipole component, especially for highly stratified models with a well-developed turbulence in the near-surface layers (because of the sharp density drop). When confronting magnetic morphologies from dynamo simulations and observations of stars, the knowledge of the time evolution of the magnetic field as well as the contribution of the intermediate- and small-scale dynamics therefore appears crucial to make a relevant diagnostic of the magnetic morphology in such stratified systems.

6. Conclusion

In Raynaud et al. (2018, Paper I), the effect of the Coriolis acceleration on the mean latitudinal surface brightness distribution of cool stars (i.e., with a relatively thick convective envelope) has been studied through direct numerical simulations. This latter study has highlighted the existence of the so-called Coriolis darkening phenomenon, that is, the emergence of a darker surface equatorial band in fastly-rotating convective shells, and whose nature differs from the usual gravity darkening effect expected in the oblate radiative envelopes of hotter highly-rotating stars. In this paper, we extend this latter work taking into account the additional effect of magnetic fields generated by a self-sustained convective dynamo. To do so, we perform a parametric study based on a large set of direct numerical simulations of convection in rotating shells with large density stratification. We focus on stars with a convective envelope thickness of about 65% of the total radius, which is representative of ∼0.35 M_⊙ cool M dwarfs, stars at the beginning of the ascent of the red giant branch and stars on the pre-main sequence.

Based on the computation of the azimuthally and time averaged latitudinal Nusselt profile, we distinguish three different leading-order regimes, depending on how far we are from the onset of convection. First, close to convection threshold, while the mean heat flux is maximum at the equator for nonmagnetic models, dynamo magnetic fields tend in contrast to increase it globally inside the tangent cylinder (i.e., close to the rotation poles). Second, for weakly turbulent models, that is, for models in which the emergence of a mean zonal surface jet starts inhibiting the heat flux close to the equator in the hydrodynamic case, the presence of dynamo magnetic fields can inverse the pole-equator brightness contrast, hence generating a brighter equator on average. Owing to the laminarization effect of the magnetic fields on the convective flow, we retrieve at some point, when the magnetic fields are strong enough, the same global behavior as that observed close to the onset of convection; the shape of the Nusselt profile is similar in both cases but the mean value consistently remains larger for larger values of the Rayleigh number. Third, for well-developed turbulent flows (i.e., for large values of the Rayleigh number Ra ≳ 5 Ra_c), magnetic fields tend to smooth out the mean latitudinal brightness distribution. In all these regimes, we show that the magnetic fields systematically reduce the surface large-scale differential rotation level. Whereas the mean azimuthal velocity exhibits a cylindrical symmetry in the nonmagnetic case according to the Taylor-Proudman constraint, magnetic fields tend to rearrange the flow into a conical-like symmetry while progressively quenching the latitudinal surface variations as its magnitude increases. The lowering of the mean equatorial jet turns out to be directly correlated with the increase of the heat flux close to the equator in the case of either a pole-equator brightness contrast inversion in the second regime or of its global reduction in the third turbulent regime. At least, this generally happens for models with a surface convective Rossby number, Ro_surf, lower than unity. Besides, for the most turbulent models (i.e., Ro_surf ≳ 1), the surface heat flux is expected to decouple from the differential rotation and to become homogeneous; this happens close to the transition from a solar to an antisolar surface differential rotation. However, the same mean behavior is observed: magnetic fields tend to reduce independently the global surface brightness contrast during this transition.

These predictions obtained about the global magnetic Coriolis darkening now have to be tested through observations. We expect low-mass nearby M dwarfs to be one of the most interesting targets in this context. Direct mapping of the surface brightness distribution of such stars remains nonetheless challenging as it would require very-high-resolution interferometric imaging, still unreachable with the current instruments. Indirect measurements exploiting the variations in the light curve of M-dwarf stars induced by the transits of close-in exoplanets can be helpful too to scan the surface features, although the whole problem is somehow degenerate with the unknown planet parameters. The potential of such indirect constraints, as well as the natural repercussions of the Coriolis darkening effect on the estimate of the exoplanet properties, certainly deserve a dedicated study in a near future. As the next step in the investigation, the systematic computation of light curves based on the mean brightness distribution profile derived from MHD simulations (e.g., Yadav et al. 2015a) will be necessary to make the most relevant possible comparison with observational constraints and their intrinsic limitations (e.g., spatial and time resolution). This will also provide practical models of the stellar disk light that could be used for applications in interferometry observations, Zeeman-Doppler imaging or exoplanet transit analyses. In addition, this will open the door to the study of the temporal dimension of the problem by means of global dynamo simulations, which has been intentionally put aside in this work for a first-step study. Focusing on the leading-order average latitudinal brightness distribution, we indeed did not take into account the impact of rotating nonaxisymmetric structures or of dynamo cycles. Going further in this direction, we will question the notion of stellar spots, their emergence, their persistence, or the relation between the surface magnetic structures and the emitted light. Regarded as higher-order effects compared to the mean behavior in this study, we can expect them to dominate the Coriolis darkening in case the mean brightness distribution tends to become uniform, especially for slowly-rotating cool stars, as for instance the Sun. Over a longer term, such a project will certainly need to account for a relevant description of the radiative transfer in the stellar atmosphere. While numerical tools already exist to tackle this issue locally (e.g., Panja et al. 2020; Witzke et al. 2022; Ludwig et al. 2023), with imposed initial equilibrium MHD properties (e.g., magnetic field flux intensity, bottom velocity boundary conditions), a strong effort will be needed in the future to couple such local atmosphere models with the global dynamo simulations performed in this study.

Acknowledgments

We are indebted to C. Catala for very fruitful discussions on the paper. During this work, C.P. was funded by postdoctoral grants from Sorbonne Université (Émergence project) and the Research Federation PLAS@PAR. C.P. and M.R. acknowledge financial support from Centre National d’Études Spatiales (CNES). M.R. also acknowledges support from the European Research Council (ERC) under the Horizon Europe programme (Synergy Grant agreement N°101071505: 4D-STAR). While funded by the European Union, views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This study was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This work was also supported by the “Programme National de Physique Stellaire” (PNPS) and the “Programme National des Hautes Énergies” (PNHE) of CNRS/INSU co-funded by CEA and CNES, by a grant from Labex OSUG@2020 (Investissements d’avenir – ANR10 LABX56), and by l’Agence Nationale de la Recherche (ANR), project ANR-22-CE31-0020. It has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium), as well as of the SIMBAD database, operated at CDS, Strasbourg, France.

References

Appendix A: MHD anelastic equations

Following Jones et al. (2011), we use the LBR formulation of the anelastic approximation (Lantz & Fan 1999; Braginsky & Roberts 1995), in which the closure relation for the heat flux is expressed as a diffusive law of the entropy, that is, −κρc_pT∇S with c_p the specific heat at constant pressure, ρ the mean density and S the entropy. The reference state is assumed to be close to the adiabatic hydrostatic equilibrium, that is, following the polytropic structure for the pressure, density, and temperature

$$\begin{array}{c}\hfill P=P_{\mathrm{c}}{w}^{n+1},\rho ={\rho }_{\mathrm{c}}{w}^n,T=T_{\mathrm{c}}w,\end{array}$$(A.1)

with

$$\begin{array}{cc}& w=c_0+\frac{c_{1}d}{r},c_0=\frac{2w_{\mathrm{o}}-\chi -1}{1-\chi },\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& c_1=\frac{(1+\chi )(1-w_{\mathrm{o}})}{{(1-\chi )}^2},w_{\mathrm{o}}=\frac{\chi +1}{\chi exp(N_{\rho }/n)+1}.\hfill \end{array}$$(A.3)

In the above expressions, n is the polytropic index and N_ρ = ln[ρ(r_i)/ρ(r_o)] is the number of density scale heights considered in the simulation domain. The values ρ_c, T_c, and P_c are the reference-state density, temperature, and pressure midway between the inner and outer boundaries. A small deviation from an adiabatic state is nevertheless needed if we want convection to start. Convection is thus triggered by imposing an entropy drop ΔS between r_i and r_o. By scaling length by the shell width d, time by the magnetic diffusion timescale d²/η, entropy by ΔS, pressure by Ωρ_cη, density by ρ_c, temperature by T_c, and magnetic field by $\sqrt{\mathrm{\Omega }{\rho }_{\mathrm{c}}\mu \eta }$, where μ is the magnetic permeability, the equations governing the anelastic dynamics in the frame co-rotating with the reference shell are given by

$$\begin{array}{cc}\hfill \frac{\partial \mathit{v}}{\partial t}+\mathit{v}·\mathbf{\nabla }\mathit{v}=& -\frac{\mathit{Pm}}{E}\mathbf{\nabla }\left(\frac{P^{\prime }}{w^n}\right)+\frac{RaP{m}^2}{\mathit{Pr}}\frac{S}{r^2}{\mathit{e}}_r-\frac{2Pm}{E}{\mathit{e}}_z\wedge \mathit{v}\hfill \\ \hfill & +\frac{\mathit{Pm}}{E{w}^n}(\mathbf{\nabla }\wedge \mathit{B})\wedge \mathit{B}+Pm{\mathit{F}}_{\nu },\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill \frac{\partial S}{\partial t}+\mathit{v}·\mathbf{\nabla }S=& \frac{\mathit{Pm}}{w^{n+1}Pr}\mathbf{\nabla }·\left(w^{n+1}\mathbf{\nabla }S\right)\hfill \\ \hfill & +\frac{\mathit{Di}}{w}[\frac{{(\mathbf{\nabla }\wedge \mathit{B})}^2}{E{w}^n}+Q_{\nu }]\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}=& \mathbf{\nabla }\wedge (\mathit{v}\wedge \mathit{B})+{\mathbf{\nabla }}^2\mathit{B}\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\mathit{B}=& 0\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\left(w^n\mathit{v}\right)=& 0\hfill \end{array}$$(A.8)

with P′ the pressure perturbation compared to the reference state. The viscous force F_ν is equal to F_ν = w⁻ⁿ∇S, where S is the rate of strain tensor provided by

$$\begin{array}{c}\hfill {\mathsf{S}}_{\mathit{ij}}=2w^n(e_{\mathit{ij}}-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }·v),e_{\mathit{ij}}=\frac{1}{2}(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}).\end{array}$$(A.9)

Moreover, the dissipation parameter and the viscous heating are respectively equal to

$$\begin{array}{c}\hfill Di=\frac{c_{1}Pr}{\mathit{PmRa}},Q_{\nu }=2[e_{\mathit{ij}}{e}_{\mathit{ij}}-\frac{1}{3}{(\mathbf{\nabla }·\mathit{v})}^2].\end{array}$$(A.10)

Appendix B: Simulation parameters

The critical values of the Rayleigh number for the set of parameters used in this work are provided in Table B.1. An overview of the parameter space covered by our simulations is given in Table B.2. Finally, the set of simulations considered for the paper is presented in more detail in Table B.3.

Appendix C: Mean surface Nusselt profiles

Fig. C.1. Mean latitudinal profile of the Nusselt number defined in Eq. (3), for models with Nρ = 3, E = 3 × 10−5, and Pr = 1. The value of the Rayleigh number Ra increases from the left to the right and the magnetic Prandtl number Pm is varied. Solid black lines represent either pure hydrodynamic models or MHD models without magnetic dynamo (i.e., too low values of Pm). Note that panel (a) is similar to Fig. 6a but is added here for the sake of completeness and comparison.

Fig. C.2. Same as in Fig. C.1 but for models with E = 3 × 10−5, Pr = 0.3, and moderate values of the Rayleigh number. The first and second rows correspond to Nρ = 3 and 4, respectively. Note that panel (e) is similar to Fig. 6b but is added here for the sake of completeness and comparison.

Fig. C.3. Same as in Fig. C.2, but for models with E = 10−4 and Pr = 1. The first and second rows correspond to Nρ = 3.5 and 4, respectively, while the last rows correspond to Nρ = 5.

Fig. C.4. Same as in Fig. C.3, but for models with E = 3 × 10−4, Pr = 1, and Nρ = 6. The model with Ra = 32 Rac requires a very high resolution and may be not fully resolved or relaxed. We checked this has no effect on the shape of the differential rotation; we nevertheless have to be cautious when considering the absolute value of the surface Nusselt. Note that panel (c) is similar to Fig. 6c but is added here for the sake of completeness and comparison.

Appendix D: Surface azimuthal velocity profiles

Fig. D.1. Same as in the second row of Fig. C.2, but for the mean azimuthal velocity profile as a function of the colatitude. Note that panel (b) is similar to Fig. 7 but is added here for the sake of completeness and comparison.

Fig. D.2. Same as in Fig. C.3, but for the mean azimuthal velocity profile as a function of the colatitude.

Fig. D.3. Same as in Fig. C.4, but for the mean azimuthal velocity profile as a function of the colatitude.

All Tables

All Figures

Fig. 1. Hammer projection of the radial component of the magnetic field at r = 0.9 ro and at a given time, for four typical dynamo models observed in direct numerical simulations. (a) Axial dipole with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac and Pm = 5. (b) Equatorial dipole with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. (c) Hemispherical dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 4 Rac, and Pm = 2. (d) Coherently oscillating dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 8 Rac, and Pm = 2.

In the text

Fig. 2. Surface butterfly diagram of the four typical examples shown in Fig. 1, i.e., the azimuthal average of the magnetic field radial component as a function of time and colatitude at the surface. (a) Axial dipole with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac, and Pm = 5. (b) Equatorial dipole with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. (c) Hemispherical dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 4 Rac, and Pm = 2. (d) Coherently oscillating dynamo with E = 3 × 10−4, Pr = 1, Nρ = 6, Ra = 8 Rac, and Pm = 2.

In the text

	Fig. 3. Dipolarity indicator as defined in Eq. (5) for the axial dipole considered in the panel Fig. 1a, with E = 10−4, Pr = 1, Nρ = 3.5, Ra = 2.5 Rac, and Pm = 5. The relative contribution of the dipole is computed for the cut-off degrees ℓcut = ℓmax, 12, and 5 (black, blue, and red lines, respectively).
In the text

	Fig. 4. Same as in Fig. 3, but for the equatorial dipole considered in Fig. 1b, with E = 3 × 10−5, Pr = 1, Nρ = 3, Ra = 2.48 Rac, and Pm = 2. The relative contribution of the axial dipole alone (m = 0) for the cut-off degrees ℓcut = ℓmax, denoted with $f_{\mathrm{dip}}^{\mathrm{ax}}$, is also represented in green.
In the text

	Fig. 5. Meridional slices of the azimuthally averaged zonal velocity, for different values of the magnetic Prandtl number Pm, in unit of the viscous velocity scale ν/d. Models with Pr = 1, E = 3 × 10−4, Nρ = 6, and Ra = 8 Rac. (a) Hydro, (b) Pm = 1, (c) Pm = 2, (d) Pm = 3.
In the text

Fig. 6. Mean latitudinal profile of the surface Nusselt number defined in Eq. (3), for three typical models. (a) E = 3 × 10−5, Pr = 1, Nρ = 3, and Ra ≈ 1.7 Rac close to the onset of convection (b) E = 3 × 10−5, Pr = 0.3, Nρ = 4, and Ra ≈ 3 Rac at the transition toward turbulent flows (c) E = 3 × 10−4, Pr = 1, Nρ = 6, and Ra ≈ 8 Rac. Left: variations with the magnetic Prandtl number Pm. Solid black lines represent either pure hydrodynamic models or MHD models without magnetic dynamo (i.e., too low values of Pm). Right: spherical representation for the nonmagnetized model (i.e., the lowest value of Pm) and the most magnetized dynamo model (i.e., the highest value of Pm) in each case. Note that we use different colorbars in each row to enhance the visualization contrast.

In the text

	Fig. 7. Mean surface azimuthal velocity as a function of the colatitude, for the models in Fig. 6b, i.e., with E = 3 × 10−5, Pr = 0.3, Nρ = 4, and Ra = 3 Rac. The magnetic Prandtl number Pm is varied and the solid black line represents the model with no dynamo (i.e., too low values of Pm).
In the text

	Fig. 8. Maximum brightness contrast as a function of the mean Elsasser number, i.e., a proxy of the mean magnitude of the magnetic field, computed for models with E = 3 × 10−4, Nρ = 6 (stars), E = 10−4, Nρ = 5 (diamonds), and E = 10−4, Nρ = 4 (squares); the color codes the Rayleigh number. For all these models, Pr = 1.
In the text

Fig. 9. Transit light curve of an exoplanet occulting a star with an equatorial jet dimmer than the rest of its photosphere. This simulation made use of starry (Luger et al. 2019). It features a limb-darkened star with an adimensional unit radius and an occulting spherical body of relative radius 0.15. The equatorial band has a contrast of 1 (for illustration purposes) and an angular radius of 0.2 radians. The upper plot shows the system at scale for six different times. The bottom plot shows the occultation light curve (black line), the differential flux measured for the six times shown above (points), as well as the transit signal that would be observed for a uniform limb-darkened star (gray).

In the text

Fig. C.1. Mean latitudinal profile of the Nusselt number defined in Eq. (3), for models with Nρ = 3, E = 3 × 10−5, and Pr = 1. The value of the Rayleigh number Ra increases from the left to the right and the magnetic Prandtl number Pm is varied. Solid black lines represent either pure hydrodynamic models or MHD models without magnetic dynamo (i.e., too low values of Pm). Note that panel (a) is similar to Fig. 6a but is added here for the sake of completeness and comparison.

In the text

	Fig. C.2. Same as in Fig. C.1 but for models with E = 3 × 10−5, Pr = 0.3, and moderate values of the Rayleigh number. The first and second rows correspond to Nρ = 3 and 4, respectively. Note that panel (e) is similar to Fig. 6b but is added here for the sake of completeness and comparison.
In the text

	Fig. C.3. Same as in Fig. C.2, but for models with E = 10−4 and Pr = 1. The first and second rows correspond to Nρ = 3.5 and 4, respectively, while the last rows correspond to Nρ = 5.
In the text

Fig. C.4. Same as in Fig. C.3, but for models with E = 3 × 10−4, Pr = 1, and Nρ = 6. The model with Ra = 32 Rac requires a very high resolution and may be not fully resolved or relaxed. We checked this has no effect on the shape of the differential rotation; we nevertheless have to be cautious when considering the absolute value of the surface Nusselt. Note that panel (c) is similar to Fig. 6c but is added here for the sake of completeness and comparison.

In the text

	Fig. D.1. Same as in the second row of Fig. C.2, but for the mean azimuthal velocity profile as a function of the colatitude. Note that panel (b) is similar to Fig. 7 but is added here for the sake of completeness and comparison.
In the text

	Fig. D.2. Same as in Fig. C.3, but for the mean azimuthal velocity profile as a function of the colatitude.
In the text

	Fig. D.3. Same as in Fig. C.4, but for the mean azimuthal velocity profile as a function of the colatitude.
In the text