© The Authors 2026

1. Introduction

Stars and planets are often characterized by magnetic fields generated by a magnetic dynamo, namely a mechanism that converts kinetic energy into magnetic energy via magnetoconvection in a rotating, electrically conducting fluid (such as a stellar plasma). Helical motions, generated by the combined effects of convection and rotation, break mirror symmetry and enable the amplification and maintenance of significant magnetic fields (Moffatt 1978). Although the specific configuration of the dynamo mechanism depends on the internal structure of the convective zone, some general features are common across various stellar types and some planets (Christensen et al. 2009).

In stars with a radiative core, differential rotation often develops within the convection zone, along with a velocity shear layer at the boundary between the radiative core and the convective envelope, similar to the solar tachocline. In contrast, fully convective stars lack a radiative zone and usually display nearly solid-body rotation (e.g., Donati et al. 2006; Morin et al. 2008; Davenport et al. 2020). This observational evidence makes it possible to distinguish between scenarios for the dynamo in stars with significant differential rotation and the dynamo in stars that show differential rotation compatible with solid-body rotations. Specifically, when differential rotation is significant, it produces the Ω effect, which, together with the helical turbulence (the α effect), generates a kind of dynamo classified as the α-Ω dynamo (Parker 1955). In contrast, when differential rotation does not play a significant role, the mechanism is classified as the α² dynamo. In an α² dynamo, the magnetic field is amplified primarily through cyclonic turbulence generated by magnetoconvection in conjunction with stellar rotation. This type of dynamo is conjectured to be responsible for magnetic field generation in fully convective stars and planets like the Earth, as observations are consistent with the absence of velocity shear and differential rotation in their convection zones.

In a simplified framework, an α² dynamo can be modeled through Rayleigh–Bénard (RB) convection in a magnetohydrodynamic (MHD) fluid that, under the effect of rotation, produces a helical flow capable of amplifying and sustaining intense magnetic fields. RB convection is a paradigm of a fluid layer heated from the bottom and cooled at the top, where the Nusselt number (Nu) is a dimensionless parameter that measures the efficiency of heat transport across the layer (Özisik 1985; White 1988).

The challenge of characterizing heat transport in astrophysical convection arises from the enormous range of scales involved. In the solar convection zone, an illustrative example, spatial scales span from giant cells ∼200 Mm near the base, with lifetimes of about 1 month, to granulation ∼1 Mm and sub-granular features down to tens of kilometers at the solar surface that evolve in minutes (Hathaway et al. 2013; Miesch 2005; Nordlund et al. 2009; Kuridze et al. 2025). Consequently, the enormous range of spatial and temporal scales implies an extremely high Reynolds number, possibly reaching Re ∼ 10¹⁵ or even higher. Rieutord & Rincon (2010) estimate Re values of 10¹⁰ to 10¹³ in the solar convection zone, based on typical convective velocities and the vertical extent of the convective layer, and Hanasoge et al. (2012) report Re ∼10¹² − 10¹⁶ and Ra ∼ 10¹⁹ − 10²⁴. We expect similar values for the convection zone in other late-type stars.

Unfortunately, the convection regimes observed in astrophysical contexts are not yet accessible to current direct numerical simulations (DNSs). To probe the asymptotic regime of magnetoconvection, where both turbulent and magnetic Reynolds numbers (Rm) are very large, one must therefore rely on reduced or parameterized models. Simplified approaches such as mean-field, low-order, or shell-model formulations provide an alternative route to explore convection and dynamo processes under extreme parameter conditions. Although shell models neglect geometric and stratification effects, they retain key nonlinear interactions and allow a wide exploration of parameter space at a relatively low computational cost (Biferale 2003). In particular, by representing turbulent cascades as a hierarchy of interacting Fourier shells, they are well suited to studying asymptotic trends and scaling laws (e.g., Frick & Sokoloff 1998; Brandenburg et al. 1996, 1997; Nigro et al. 2005; Nigro & Veltri 2011; Perrone et al. 2011; Nigro 2013; Nigro & Carbone 2015; Brandenburg & Larsson 2023).

Understanding how Nu asymptotically scales with the control parameters of magnetoconvection, such as the Rayleigh (Ra) and Prandtl (Pr) numbers, is important because this scaling can provide a quantitative measure of the efficiency of convective heat transport in regimes pertinent to naturally occurring dynamos. In stellar interiors, where convection is highly supercritical, the resulting turbulent flows might yield large Nu, indicating efficient heat transport and reflecting the dominance of convective over radiative and conductive transport (Brun & Browning 2017).

In stellar interiors, the luminosity sets the outward energy flux that must be transported through the convective envelope. In fully convective stars, where convection carries most or all of this flux, the Nusselt number provides a natural dimensionless measure of how efficiently turbulent convective motions transport the required energy flux. In a dynamo-active convective region, the self-generated magnetic field can influence the efficiency of convection in transporting the energy flux by redistributing part of the convective energy into magnetic energy. In MHD simulations of stellar convection, the same convective motions that carry the luminosity also amplify the magnetic field through nonlinear dynamo action (e.g., Miesch et al. 2000; Brun et al. 2004; Browning et al. 2006; Browning 2008; Pongkitiwanichakul et al. 2016; Brun & Browning 2017; Brun et al. 2022). The behavior of the Nusselt number in dynamo-active systems is therefore expected to depend not only on thermal forcing but also on the nonlinear feedback of the self-generated magnetic field on the flow. These factors are crucial for understanding how convective heat transfer, magnetic-field amplification, and large-scale dynamo action become dynamically interconnected in stellar and planetary interiors. In line with this view, Christensen & Aubert (2006) and Christensen et al. (2009) have emphasized the central role of convection and heat flux in the regulation of magnetic field generation through dynamo processes.

Previous studies of non-magnetized RB convection have established that Nu generally follows a power-law relationship with Ra and Pr, such as Nu ∝ Ra^γPr^β, with exponents β and γ that depend on the turbulent regime (Ahlers et al. 2009; Chillà & Schumacher 2012; Pandey & Sreenivasan 2021). In magnetized convection, the situation is more complex. The scaling properties of heat transport have been investigated in a variety of idealized configurations, both theoretically and numerically, as extensively reviewed by Proctor & Weiss (1982) and in the monograph by Weiss & Proctor (2014). These studies, although restricted to simplified 2D configurations and in regimes with lower Reynolds numbers than those considered here, have identified several asymptotic regimes and magnetic-field-dependent modifications of the classical Nu(Ra) scaling, including quasi-linear behaviors in strongly magnetized configurations. More recently, numerical simulations of RB magnetoconvection with high Ra and imposed magnetic fields have identified distinct dynamical regimes and the corresponding heat-transport scalings (e.g., Yan et al. 2019; Cresswell et al. 2023). Although such simulations were typically performed in idealized settings where, in some cases, the magnetic field is imposed rather than self-generated, their results, together with those provided in the review by Plumley & Julien (2019), highlight the wide variety of possible Nu(Ra) behaviors in magnetized flows. Nevertheless, the asymptotic scaling of convective heat transport in dynamo-active systems and in fully turbulent regimes characterized by simultaneously very large Ra and Rm, as relevant to stellar and planetary interiors, remains poorly constrained and difficult to access with DNSs. This motivates the use of reduced models to explore extreme parameter regimes.

We used a thermally driven shell model (Nigro 2022; Nigro et al. 2026) to investigate the asymptotic behavior of convective heat transport in an α² dynamo with extreme parameter values. Although the model neglects density stratification, explicit geometry, and differential rotation, it retains the essential nonlinear couplings between convection and magnetic-field generation, making it well suited to explore regimes inaccessible to present DNSs.

The article is organized as follows. Section 2 describes the model, Sect. 3 presents the numerical results and the derived scaling laws, and Sect. 4 summarizes our findings and discusses them in the context of previous studies of asymptotic convective heat transport.

2. Numerical model

The model adopted in this study is a thermally driven magneto-convective shell model, introduced by Nigro (2022), which is described by the following coupled equations:

$$\begin{array}{ccc}& (\frac{d}{\mathit{dt}}+\nu k_n^2)u_n=-\stackrel{\sim }{\alpha }{\theta }_n+i{k}_n[(u_{n+1}{u}_{n+2}-b_{n+1}{b}_{n+2})+\hfill & \hfill \\ \hfill & -\frac{\xi }{2}(u_{n-1}{u}_{n+1}-b_{n-1}{b}_{n+1})-\frac{1-\xi }{4}(u_{n-2}{u}_{n-1}-b_{n-2}{b}_{n-1}){]}^{\ast }\hfill \end{array}$$(1)

$$\begin{array}{cc}& (\frac{d}{\mathit{dt}}+\eta k_n^2)b_n=i{k}_n[(1-\xi -{\xi }_m)(u_{n+1}{b}_{n+2}-b_{n+1}{u}_{n+2})+\hfill \\ \hfill & +\frac{{\xi }_m}{2}(u_{n-1}{b}_{n+1}-b_{n-1}{u}_{n+1})+\frac{1-{\xi }_m}{4}(u_{n-2}{b}_{n-1}-b_{n-2}{u}_{n-1}){]}^{\ast }\hfill \end{array}$$(2)

$$\begin{array}{cc}& (\frac{d}{\mathit{dt}}+\chi k_n^2){\theta }_n=i{k}_n[{\alpha }_1{u}_{n+1}^{\ast }{\theta }_{n+2}^{\ast }+{\alpha }_2{u}_{n+2}^{\ast }{\theta }_{n+1}^{\ast }+\hfill \\ \hfill & +{\beta }_1{u}_{n-1}{\theta }_{n+1}-{\beta }_2{u}_{n+1}{\theta }_{n-1}+{\gamma }_1{u}_{n-1}{\theta }_{n-2}+{\gamma }_2{u}_{n-2}{\theta }_{n-1}{]}^{\ast }+f_n,\hfill \end{array}$$(3)

where u_n are the velocity field fluctuations, b_n magnetic field fluctuations, and θ_n temperature fluctuations at the shell n, where n = 1, …, N and N is the total number of shells, ν is the kinematic viscosity, η the magnetic diffusivity, and χ is the thermal diffusivity, and $\stackrel{\sim }{\alpha }$ denotes the thermal convection coefficient. The values of the model coefficients ξ = 1/2 and ξ_m = 1/3 are chosen in such a way that Eqs. (1) and (2) conserve the quadratic invariants of the MHD equations – total energy, cross helicity, and magnetic helicity – in the ideal case, namely where there is no external forcing or dissipation, and when $\stackrel{\sim }{\alpha }=0$ (Frick & Sokoloff 1998; Giuliani & Carbone 1998). The values of the coefficients in the temperature equations are α₁ = α₂ = 1, β₁ = β₂ = 1/2, and γ₁ = γ₂ = −1/4, corresponding to the coupling with the model of Jensen et al. (1992, see also Mingshun & Shida 1997). Finally, the forcing (f_n) is a random Gaussian signal with unit standard deviation, with an autocorrelation time of one time unit. The equations are normalized to the typical physical parameters, namely the velocity is normalized to the freefall velocity $U=\sqrt{\alpha gL\mathrm{\Delta }T}=\sqrt{\nu \chi \mathrm{Ra}}/L$, considering the typical magnetic field B₀ = U (the typical magnetic field obtained when the freefall velocity U is equal to the Alfvén speed), the temperature is measured in terms of ΔT; and finally, the length in units of L, that is, the characteristic length of the system and the time in units of freefall time ${(L/\widehat{\alpha }g\mathrm{\Delta }T)}^{1/2}$.

The shell model describes the evolution of velocity, temperature and magnetic field fluctuations in a hierarchy of exponentially spaced wavenumber shells k_n = k₀λⁿ, with n = 1, …, N, where λ = 2 is the shell-spacing parameter. Each shell represents a scale of turbulent eddies, and the nonlinear interactions are constructed to mimic the local cascade processes in MHD turbulence.

The large-scale magnetic field (b₁) is modified to include a supercritical pitchfork bifurcation term (Nigro & Carbone 2010):

$$\begin{array}{c}\hfill \frac{d{b}_1}{\mathit{dt}}=-\eta k_1^2{b}_1+i\frac{k_1}{6}{(u_2{b}_3-b_2{u}_3)}^{\ast }+\mu b_1(1-\frac{b_1^2}{B_0^2}).\end{array}$$(4)

The pitchfork bifurcation term μb₁(1 − b₁²/B₀²) accounts for magnetic field reversals, as b₁ = ±B₀ are two stable equilibrium states for the large-scale magnetic field b₁, while b₁ = 0 is an unstable equilibrium for μ > 0 (Hoyng et al. 2001; Benzi & Pinton 2010; Nigro et al. 2013). This term is introduced in the spirit of mean-field dynamo theory to mimic the nonlinear transfer of magnetic energy and helicity from small to large scales (for more details, see Nigro et al. 2026). In this term, the parameter μ is a positive dimensionless coefficient that controls the growth rate of the large-scale magnetic field component (b₁) and the relaxation rate toward the two saturated states. In a reduced way, μ represents the inductive effect associated with the kinetic helicity of flow motions (Pouquet et al. 1976). In the present dimensionless formulation, we set μ = 1, so that the magnetic-field growth timescale is unity, while the saturated magnetic equilibria remain fixed at b₁ = ±B₀.

We numerically evolved the system using the fourth-order Runge–Kutta numerical scheme. We considered a wide range of model parameters, focusing in particular on variations in the thermal diffusivity χ, the kinematic viscosity ν and the magnetic diffusivity η. These parameters determine the Rayleigh number, Ra $=\stackrel{\sim }{\alpha }{\theta }_0{L}^3/(\nu \chi )$, with ${\theta }_0=\sqrt{⟨}{\displaystyle \sum _{n=1}^N}({\theta }_n{\theta }_n^{\ast }){⟩}_t$, according to the shell-model description; the Reynolds number (Re =Lu₀/ν, where $u_0=\sqrt{⟨}{\displaystyle \sum _{n=1}^N}{u}_n{u}_n^{\ast }{⟩}_t$); and the magnetic Reynolds number (Rm =Lu₀/η). The Prandtl number (Pr) and the magnetic Prandtl number (Pm) are likewise determined by the ratios ν/χ and ν/η, respectively.

The objective of this model is to reproduce the essential features of an α² dynamo in a convective MHD plasma with high kinetic helicity that we assume is produced by rapid rotation, i.e., the ratio between kinetic helicity and kinetic energy is equal to 1, corresponding to the maximum kinetic helicity. An important feature of this model is its computational efficiency, which enables the exploration of turbulent regimes with Ra, Re, and Rm values well beyond the reach of current DNSs, although spatial information, such as density stratification, is neglected according to the shell model assumptions. The model captures essential aspects of convective heat transfer and dynamo action, including magnetic-field reversals and their relationship to the level of turbulence and convective transport.

3. Numerical results

The Nusselt number (Nu) is a dimensionless measure of convective heat transport expressed in terms of the ratio of total heat flux (convective and conductive) to conductive heat flux (Ahlers et al. 2009; Chillà & Schumacher 2012; Verma 2018; Pandey et al. 2021). Taking into account v_z as the velocity component along the acceleration of gravity, i.e., along the z axis, u as the vertical turbulent velocity and Θ the temperature fluctuations, following Benzi et al. (1998), we can write the following:

$$\begin{array}{c}\hfill Nu=\frac{⟨v_{z}T⟩-\chi \frac{\partial ⟨T⟩}{\partial z}}{\chi \frac{\mathrm{\Delta }T}{L}}=\frac{⟨u\mathrm{\Theta }⟩}{\chi \frac{\mathrm{\Delta }T}{L}},\end{array}$$(5)

where the brackets represent the space and time averages.

In DNSs, there are different approaches to calculating Nu, although this should converge to a consistent value if the DNS is well resolved. The definition adopted for Nu is therefore crucial to characterizing heat transport. For example, Xu et al. (2020) tested different approaches to computing Nu, showing that the relation $Nu=\sqrt{\mathit{RaPr}}{⟨v_z\mathrm{\Theta }⟩}_{V,t}$, where the brackets represent the volume V and the average time, provides the best definition, among those they tested, to better describe the heat transport in their setup. In our shell model framework, the Nusselt number is defined consistently with the model formalism as

$$\begin{array}{c}\hfill Nu=\sqrt{\text{Ra Pr}}\left|{\langle {\displaystyle \sum _n}(u_n{\theta }_n^{\ast }+u_n^{\ast }{\theta }_n)\rangle }_t\right|,\end{array}$$(6)

where the brackets ⟨⟩_t indicate a time average, the quantity 𝒞_vθ = |⟨∑_n(u_nθ_n^*+u_n^*θ_n)⟩_t| indicates the correlation between velocity and temperature fluctuations in dimensionless units.

Individually varying the dissipative coefficients samples different paths in the control-parameter space, and this must be taken into account when interpreting Fig. 1. In particular, when the thermal diffusivity χ decreases at fixed ν, both Ra ∝ (νχ)⁻¹ and Pr = ν/χ increase simultaneously. Along this path obtained for decreasing χ, we find approximately Nu ∼ Ra and Nu ∼ Pr over the explored range.

Fig. 1. Asymptotic behavior of the Nusselt number (Nu) as the Rayleigh number (Ra) and the Prandtl number (Pr) increase.

Therefore, the Nu(Pr) scaling shown in the right panel of Fig. 1 should not be interpreted as a pure Prandtl-number dependence at a fixed Rayleigh number. Rather, it reflects a combined variation of Ra and Pr, together with the explicit prefactor $\sqrt{\mathrm{Ra}\mathrm{Pr}}$ in Eq. (6). Indeed, the nearly linear increase in Nu with Pr in the right panel of Fig. 1 is not obtained with a constant Ra.

The scaling Nu ∼ Ra obtained when varying χ is steeper than the classical hydrodynamic (HD) scalings. Classical HD studies propose Nu ∼ Ra^1/3 (Malkus 1954a,b; Spiegel 1962) and Nu ∼ Ra^1/2 from Kraichnan’s rigorous upper-bound argument (Kraichnan 1962; Spiegel 1971). The precise scaling of Nu with Ra has been the subject of extensive theoretical and experimental investigation, sometimes leading to different or contradictory results (Niemela et al. 2000; Niemela & Sreenivasan 2006; Ahlers et al. 2009; Chillà & Schumacher 2012; He et al. 2012; Urban et al. 2012; Doering et al. 2019; Iyer et al. 2020). The weak dependence of Nu on Pr, subsequently revealed by high-precision measurements (Pandey & Sreenivasan 2021; Bhattacharya et al. 2021), together with the difficulty of reaching asymptotic regimes, makes the precise scaling of Nu difficult to capture (Verma 2018). In our case, the linear dependence of Nu on Pr reported in the right panel of Fig. 1 arises because Pr and Ra are varied simultaneously and therefore does not contradict previous HD studies reporting a weak Prandtl-number dependence at fixed or nearly fixed Ra.

Although the dependence of Nu on Pr has turned out to be very small in HD convection, it can still be relevant for solar and stellar magnetoconvection, as the Pr values are very small. In addition, it became evident that Nu depends not only on Ra and Pr, but also on the correlation between turbulent velocity and temperature fluctuations, which in turn depends on the particular state of the turbulence (Pandey et al. 2021).

The abovementioned studies address the asymptotic law for HD flows, whereas the presence of an external and/or a dynamo-active magnetic field adds additional complexity in revealing the exact asymptotic law for Nu. This may be a crucial point, as the presence of a magnetic field is known to change the behavior of many fluid instabilities (Chandrasekhar 1961).

Asymptotic studies of convective heat transport in MHD systems, although in an idealized description with imposed magnetic fields and lower control-parameter values than those explored here, propose quasi-linear or sublinear scalings of Nu with Ra (Proctor & Weiss 1982; Weiss & Proctor 2014; Yan et al. 2019). In particular Proctor & Weiss (1982) and Weiss & Proctor (2014) derived a quasi-linear dependence of Nu on Ra with logarithmic corrections, while recent simulations with imposed magnetic fields identified distinct dynamical regimes, including magnetically constrained states with approximately linear scalings and more turbulent weak-field regimes with sublinear behavior Nu ∼ Ra^2/7, closer to the classical HD scaling (Yan et al. 2021; Cresswell et al. 2023). Moreover, Yan et al. (2019) and Plumley & Julien (2019) suggest a progressive approach to Nu ∼ Ra in the strong-field limit. Although these results are obtained in idealized MHD systems and at lower Ra, Re, and/or Rm values than the extreme parameter ranges explored by the present shell-model calculations, they remain broadly consistent with the trends observed in our work. However, it is important to note that the asymptotic behavior of the Nusselt number in fully turbulent, dynamo-active MHD systems at extreme Ra and Pr, where the magnetic field is generated self-consistently by the flow, has not yet been systematically investigated, as it remains a very challenging task in DNSs and laboratory experiments.

Our model, owing to the shell-model formulation, allows us to validate and extend to extreme parameter regimes the indication already proposed in previous MHD studies, supporting the idea that magnetic fields may modify convective heat transport. To test this hypothesis and thus explore the effect of the magnetic field dynamics on the Nu scaling law, we compared two sets of simulations: one carried out using the MHD shell model discussed here, and the other set of simulations performed in the purely HD limit, obtained by setting the magnetic field components b_n equal to zero in the model equations. In both sets of simulations, the Rayleigh number was increased by varying the thermal expansion coefficient $\stackrel{\sim }{\alpha }$ while keeping the other parameters fixed. We find that the resulting scaling laws differ significantly between the full dynamo-active model and the corresponding HD limit: in the HD case, we obtain Nu ∝ Ra^γ with γ = 0.50 ± 0.07, whereas in the corresponding MHD simulations we find γ = 1.0 ± 0.1 (see Fig. 2). This comparison shows that the full dynamo-active model displays a steeper asymptotic scaling law than the corresponding HD limit. However, this does not necessarily imply that magnetized convection generally transports more heat at fixed Rayleigh numbers, since the critical Rayleigh number and pre-factors can differ.

Fig. 2. In the asymptotic regimes, the full dynamo-active version of the model exhibits a steeper convective heat-transport scaling from Nu ∝ Ra0.5 in a corresponding HD model to Nu ∝ Ra in our MHD shell model. The depicted scaling laws were obtained considering simulations with χ = 10−4, ν = 10−3, and (in the MHD case) η = 10−4. The increase in the Rayleigh number was obtained by increasing the thermal expansion coefficient ($\stackrel{\sim }{\alpha }$) starting from the dimensionless value of 0.05 to 5000.

Across the full set of simulations, the Nusselt number (Nu) is predominantly correlated with the Rayleigh number (Ra), which acts as the dominant control parameter for heat transport. Ra shows the strongest correlation with Nu: the Pearson correlation coefficient is ρ_ℓ = 0.99 on linear scales and ρ = 0.77 on logarithmic scales. This indicates the existence of an approximate scaling law, although with a non-negligible dispersion. This dispersion reflects the fact that the Nu-Ra scaling is not controlled by Ra alone but is conditioned by the dependence on other control parameters, most notably the Prandtl number Pr and, more weakly, the magnetic diffusivity η (or, equivalently, Rm). Accordingly, the scaling can be written in the form

$$\begin{array}{c}\hfill Nu\sim {\text{Ra}}^{\gamma }{\text{Pr}}^{\beta }f(\text{Rm}),\end{array}$$(7)

where the values of the exponents γ and β depend on the path considered in the parameter space.

No significant global dependence of Nu on Rm is observed when all simulations are considered together, as indicated by a weak logarithmic Pearson coefficient ρ ∼ −0.095. However, when subsets of approximately fixed Ra are considered (i.e., fixed χ and $\stackrel{\sim }{\alpha }$), a clear secondary trend emerges: Nu decreases systematically with increasing Rm, since the log-log Pearson correlation between Nu and Rm is always negative and remains strong, with values typically in the range ρ ≈ −0.8 to −0.7 (see Table 1). In all subsets, the log-log slope of the Nu–Rm relation is negative, with a typical value of ∼ − 0.2, suggesting an approximate scaling of the form

$$\begin{array}{c}\hfill Nu\propto {\text{Rm}}^{-0.2}\end{array}$$(8)

at fixed Ra within the range of explored parameters (see Fig. 3). This behavior is consistent with a magnetically constrained convective regime in which reduced magnetic diffusivity enhances the Lorentz feedback on the flow, thereby weakening the velocity–temperature correlation and, consequently, the efficiency of convective heat transport.

Fig. 3. Panels (a)–(e): Nu versus Rm at fixed χ, showing a decrease in convective transport with increasing Rm. Panel (f): Selected cases showing no clear monotonic trend of Nu with Re. Additional simulations (not shown) display even more irregular behavior.

Our results reveal that, although the full dynamo-active model exhibits a steeper asymptotic scaling Nu(Ra) (Fig. 2), an increase in the magnetic Reynolds number Rm leads to a systematic reduction of the Nusselt number (Nu) at a fixed Rayleigh number. This highlights an anisotropic behavior of the system in the parameter space: the full dynamo-active model exhibits a steeper asymptotic scaling of convective transport along the Rayleigh direction, whereas heat transport is suppressed along the magnetic Reynolds number direction.

The dependence of Nu on the HD Reynolds number is not robust. While the selected cases depicted in Fig. 3f do not show a clear monotonic trend in Nu with a varying Re (see also the red squares in Fig. 1), the remaining cases (not shown) display even more scattered behaviors, as they exhibit variations that do not follow a common direction. The lack of a common trend with Re suggests that a higher level of turbulence does not necessarily imply an enhancement of convective heat transport.

Therefore, these results provide a physical interpretation of the different Nu(Ra) behaviors obtained when varying χ, ν, or η. These variations do not correspond to a single path through the parameter space, but rather to different trajectories in the multidimensional space spanned by Ra, Pr, and Rm. As discussed above, when χ decreases at fixed ν, both Ra and Pr increase simultaneously, and the explicit prefactor $\sqrt{\mathit{RaPr}}$ in Eq. (6) also increases. By contrast, when ν is varied at fixed χ, Ra and Pr vary in opposite directions. If only the explicit dependence on the control parameters is considered, the dependence on ν cancels out in the product of the Rayleigh and Prandtl numbers, apart from the implicit dependence through θ₀. Along this path, the resulting value of Nu is therefore mainly controlled by the velocity–temperature correlation (C_vθ) rather than by Ra alone. Finally, varying η primarily changes the magnetic coupling, or equivalently Rm, and therefore modifies the Lorentz-force feedback on the flow rather than the thermal forcing itself.

We therefore do not find evidence for a universal collapse of all the data with a function of Ra and Pr alone. A general empirical description could be given in Eq. (7), where the exponents depend on the path in the parameter space. In this expression, the function f(Rm) accounts for the magnetic feedback and for the path dependence. This interpretation is consistent with the decrease in Nu with increasing Rm observed at approximately fixed Ra, suggesting that stronger magnetic coupling weakens the velocity–temperature correlation and reduces the efficiency of convective heat transport.

As the global Nu increases, our model shows a clear tendency toward a higher frequency of magnetic-polarity reversals, as shown in Fig. 4. This trend is consistent with previous results obtained in the same dynamo framework, where rapid increases in convective heat flux were found to trigger magnetic polarity reversals (see Nigro 2022 for more details) and larger values of global Nu were shown to be associated with stronger magnetic variability. Although the occurrence of magnetic polarity reversals is probably not controlled by Nu alone, since Nu also depends on other parameters such as Ra, Pr, and Rm, the general trend suggests a coupling between vigorous convective transport and enhanced magnetic variability (for more details, see Nigro et al. 2026).

Fig. 4. Simulations with higher Nusselt numbers (Nu) produce more magnetic polarity reversals. This trend is depicted for four simulations for which Nu increases from the top to the bottom panel.

4. Discussion and conclusions

We have investigated the asymptotic behavior of convective heat transport in a simplified magneto-convective system based on a thermally driven shell model designed to capture the essential dynamics of an α² dynamo. The main advantage of this approach lies in its ability to access extreme parameter regimes, characterized by very high values of Ra, Re, and Rm, while maintaining reasonable computational cost. This makes it particularly suited for investigating asymptotic behaviors that are otherwise inaccessible to current DNSs.

The scaling behavior of the Nusselt number (Nu) with the Rayleigh and Prandtl numbers (Ra and Pr) reveals a clear power-law dependence, with slopes steeper than those typically found in purely HD convection (e.g., Niemela & Sreenivasan 2006; Ahlers et al. 2009; Chillà & Schumacher 2012; He et al. 2012; Urban et al. 2012; Iyer et al. 2020; Pandey & Sreenivasan 2021). This difference is probably related to the turbulent regime realized in the present model, in which the correlation between the temperature and velocity field fluctuations, and therefore Nu, increases with Ra more strongly than in the HD studies cited above.

With the aim of clarifying a likely difference in asymptotic power law between an MHD fluid and a non-magnetized fluid, we performed additional simulations of the purely HD version of the model. These simulations revealed significantly shallower scaling than in the corresponding MHD simulations. The comparison between HD and MHD cases suggests that the magnetic-field dynamics can modify the correlation between velocity and temperature fluctuations, thereby affecting the asymptotic heat-transport scaling. However, this result does not necessarily imply that magnetized convection transports more heat at a fixed Ra, since the critical Rayleigh number and pre-factors can differ in different systems.

Therefore, the presence of magnetic fields introduces additional nonlinear couplings that can alter the classical HD scaling exponents and energy pathways, extending the theoretical framework of turbulent convection into the magnetized domain. The approximately linear growth of Nu with Ra observed in our extended parameter space is reminiscent of the quasi-linear scalings reported in previous studies of RB magnetoconvection with imposed magnetic fields or weakly turbulent flows (Proctor & Weiss 1982; Yan et al. 2019; Cresswell et al. 2023). However, in our model, this behavior emerges in a fully turbulent, dynamo-active regime, where the magnetic field is generated self-consistently and strongly coupled to the convective flow.

We also note that, in the MHD version of the model, the large-scale magnetic-field equation includes a pitchfork bifurcation term that acts as an α-term, introduced in the spirit of mean-field dynamo theory to mimic transfers of magnetic energy and helicity from small to large scales. No analogous nonlocal term is present in the velocity equation and, therefore, no corresponding transfer between small- and large scales is explicitly modeled in the kinetic-energy channel. This asymmetry can contribute to the quantitative differences between the MHD and HD cases. However, assessing its specific role requires a more detailed investigation and will be the subject of future work.

Considering the behavior of Nu with Re at approximately fixed Ra (Fig. 3f), we suggest that high levels of turbulence do not necessarily imply enhanced convective heat transport. What ultimately determines the efficiency of heat transfer is the degree of correlation between temperature and velocity fluctuations. In certain turbulent regimes, these two fields can become partially de-correlated, thus reducing the effective heat flux (e.g., Pandey & Sreenivasan 2021). In particular, simulations showing a nearly constant behavior of Nu as a function of Ra and Pr (simulations represented by red squares in Fig. 1 and Fig. 3f) are cases in which increasing the levels of turbulence does not substantially increase the transport of convective heat.

Another important result concerns the dependence of Nu on Rm. Although no significant global correlation between Nu and Rm is found when all simulations are considered together, a clear secondary trend emerges when subsets at nearly fixed Ra are analyzed: Nu decreases systematically with increasing Rm. This result indicates that within the isotropic MHD framework described by a shell model, increasing magnetic coupling tends to suppress convective heat transport, as the system redistributes energy toward the magnetic energy channel. This trend along the Rm direction in the parameter space is likely due to the enhancement of the term describing Lorentz-force feedback, which reduces the correlation between velocity and temperature fluctuations. This behavior, which emerges when Ra is kept approximately constant, coexists with the steepening of the asymptotic Nu(Ra) scaling observed in the dynamo-active case, revealing an anisotropic behavior in the parameter space: increasing Ra in the dynamo-active regime leads to a steepening of the asymptotic heat-transport scaling, whereas increasing Rm at fixed Ra systematically reduces the Nusselt number.

This anisotropy in the parameter space suggests that magnetic fields do not act as a purely enhancing or inhibiting mechanism for convective heat transport; rather, they restructure the nature of turbulent convection through competing effects. This interpretation is consistent with previous results showing that magnetic fields were found to enhance heat-transfer efficiency only in specific dynamical regimes, for example, by suppressing zonal shear or under conditions close to magnetostrophic balance (Yadav et al. 2016).

It is worth noting that the exact indices of the scaling laws obtained here must be considered in the limit of the approximation adopted in the model description, which neglects explicit geometry, stratification, and field anisotropy. In particular, the model cannot capture effects known to arise in geometrically constrained systems with imposed magnetic fields, where the suppression of small-scale turbulence can coexist with the emergence of large quasi-2D coherent structures aligned with the field (Sommeria & Moreau 1982). In such specific configurations, these elongated columnar vortices can promote efficient heat transport (Cioni et al. 2000; Murali et al. 2021), but are exceptions to the general trend of decrease.

The dual behavior of magnetic fields in either enhancing or inhibiting convective heat transport is particularly relevant in an astrophysical context. In stellar interiors, the luminosity sets the total outward energy flux that must be transported through the convective envelope. In fully convective stars, convection is responsible for carrying most or all of this flux, so the Nusselt number can be viewed as a dimensionless measure of the efficiency with which turbulent motions transport the energy required by the stellar luminosity. In dynamo-active systems, however, the same convective motions that carry this flux also amplify the magnetic field through nonlinear dynamo action, so that part of the convective energy is redistributed into the magnetic energy channel. In this sense, the self-generated magnetic field does not simply add an extra channel of transport, but modifies the efficiency of convective heat transport itself through its nonlinear feedback on the flow. This view is consistent with global simulations of stellar convection, where the enthalpy flux dominates the energy transport through most of the convective envelope, while magnetic energy generation and feedback alter the flow dynamics and the partition of energy among the various reservoirs (e.g., Browning 2008; Brun et al. 2022).

Our model, although it is a low-dimensional model and lacks geometric complexity, describes some of the main dynamical and nonlinear features of an MHD fluid. It conserves the main physical invariants in the ideal case (i.e., the non-dissipative and unforced case) and retains the essential nonlinear interactions among velocity, temperature, and magnetic field fluctuations without requiring excessive computational effort. This makes it a valuable theoretical tool for exploring parameter space inaccessible to global 3D simulations, guiding the interpretation of more realistic calculations and observations, and for identifying potential scaling regimes relevant to stellar and planetary interiors.

Acknowledgments

We are very grateful to the anonymous reviewer who offered constructive comments and suggestions that improved the article. This work is supported by the MELODY research project, funded by the Italian Ministry of Universities and Research (MUR) under the PNRR Young Researchers program 2022 (MELODY SoE project, grant agreement No SOE_0000119, CUP E53C22002450006).

References

All Tables

All Figures

	Fig. 1. Asymptotic behavior of the Nusselt number (Nu) as the Rayleigh number (Ra) and the Prandtl number (Pr) increase.
In the text

Fig. 2. In the asymptotic regimes, the full dynamo-active version of the model exhibits a steeper convective heat-transport scaling from Nu ∝ Ra0.5 in a corresponding HD model to Nu ∝ Ra in our MHD shell model. The depicted scaling laws were obtained considering simulations with χ = 10−4, ν = 10−3, and (in the MHD case) η = 10−4. The increase in the Rayleigh number was obtained by increasing the thermal expansion coefficient ($\stackrel{\sim }{\alpha }$) starting from the dimensionless value of 0.05 to 5000.

In the text

	Fig. 3. Panels (a)–(e): Nu versus Rm at fixed χ, showing a decrease in convective transport with increasing Rm. Panel (f): Selected cases showing no clear monotonic trend of Nu with Re. Additional simulations (not shown) display even more irregular behavior.
In the text

	Fig. 4. Simulations with higher Nusselt numbers (Nu) produce more magnetic polarity reversals. This trend is depicted for four simulations for which Nu increases from the top to the bottom panel.
In the text