© The Authors 2025

1. Introduction

The Sun’s convection zone is differentially rotating (e.g., Thompson et al. 1996). This differential rotation is believed to be generated by the nonlinear interaction between rotation and magneto-convection. The solar-like differential rotation, with a faster equator and slower poles, has been repeatedly produced in three-dimensional (3D) numerical simulations of rotating convection in spherical shells (e.g., Brun & Toomre 2002; Miesch et al. 2008; Käpylä et al. 2011; Guerrero et al. 2013; Hotta et al. 2015a). In these numerical models, a key parameter controlling the differential rotation regime is the Rossby number, Ro = v/(2Ωℓ), or Coriolis number, Co = Ro⁻¹, where v and ℓ are the typical velocity and length scales of convection, and Ω is the rotation rate. In a high-Ro (or equivalently low-Co) regime, the differential rotation tends to become antisolar, with faster poles and a slower equator. On the other hand, in a low-Ro (high-Co) regime, the differential rotation tends to become solar-like (e.g., Gastine et al. 2013; Mabuchi et al. 2015; Featherstone & Miesch 2015; Viviani et al. 2018; Brun et al. 2022). The transition between solar-like and antisolar rotational profiles occurs at around Ro ≈ 1 (Gastine et al. 2013; Mabuchi et al. 2015; Viviani et al. 2018). The problem of numerical modeling of the Sun’s global-scale convection is that the typical Ro of the Sun is likely very close to this transition (e.g., Strugarek et al. 2017; Vasil et al. 2021), and therefore numerical results are highly dependent on subtle changes in the model parameters. In particular, the antisolar differential rotation is preferentially obtained in recent highly turbulent models when the solar values of the rotation rate, Ω_⊙, and luminosity, L_⊙, are used. This is, in fact, a striking aspect of the Sun’s convective conundrum (Hanasoge et al. 2012; O’Mara et al. 2016; Hotta et al. 2023; Käpylä et al. 2023).

A physical origin of the solar-like differential rotation in the low-Ro regime can be explained as follows. Under strong rotational influence, convection tends to be aligned with the rotational axis and takes the form of columnar rolls outside the tangential cylinder (Busse 1970). They are called Busse columns or banana cells and can be readily seen in simulations as north–south-aligned lanes of downflows across the equator (Busse 1970; Gilman & Miller 1986; Miesch et al. 2000). These banana cells have an associated Reynolds stress that transports the angular momentum equatorward and cylindrically outward to accelerate the equator (Miesch 2005; Käpylä et al. 2011; Hotta et al. 2015a; Matilsky et al. 2020).

It is well known that these convective columns propagate in a prograde direction faster than the local differential rotation speed (Miesch et al. 2008; Bessolaz & Brun 2011). This feature can be understood in terms of thermal Rossby waves (e.g., Busse 2002). They are prograde-propagating waves of z vorticity (where z denotes the rotational axis), originating from the conservation law of potential vorticity (Miesch 2005). They owe their existence to both the topographic β effect arising from the spherical curvature (e.g., Busse 2002) and the compressional β effect arising from the background density stratification (Glatzmaier & Gilman 1981; Evonuk 2008; Glatzmaier et al. 2009; Verhoeven & Stellmach 2014; Ong & Roundy 2020). In the Sun’s convection zone, the compressional β effect is dominant. Linear analyses of thermal Rossby modes were first carried out by Glatzmaier & Gilman (1981) in the solar context, and later by Hindman & Jain (2022), Jain & Hindman (2023) using a simple one-dimensional (1D) model. A more sophisticated two-dimensional (2D) spherical shell model of the Sun’s differentially rotating convection zone was recently presented by Bekki et al. (2022a). Using fully nonlinear 3D convection simulations, Bekki et al. (2022b) further demonstrated that the thermal Rossby modes are the main transporters of enthalpy (≈60% of what is required) and angular momentum (≈40% of the net Reynolds stress) in the simulation. Despite their dominant presence and significant dynamical importance in the numerical models, these columnar convective patterns have never been successfully observed on the Sun. In the solar photospheric observations, the convective power declines monotonically above the supergranular scale (e.g., Hathaway et al. 2015), implying the absence of large-scale convective columns. This puzzle is another manifestation of the convective conundrum.

A new promising explanation for how to generate the solar-like differential rotation was recently provided by Hotta et al. (2022), who showed that the small-scale Maxwell stress can transport the angular momentum radially upward against the downward transport by the Reynolds stress to help accelerate the equator. Although the role of the Maxwell stress in driving the differential rotation has been discussed by several authors (e.g., Fan & Fang 2014; Käpylä et al. 2017a; Warnecke et al. 2025), Hotta et al. (2022)’s simulation was the first successful reproduction of solar-like differential rotation in a highly turbulent regime with the solar rotation rate, Ω_⊙, and luminosity, L_⊙. It should be pointed out that, in their simulation, a low Ro is no longer needed to make the differential rotation solar-like. Complementary evidence was recently reported by Soderlund et al. (2025), who found a transition in the rotational regime from antisolar to solar-like at nearly fixed Ro by increasing the magnetic Prandtl number Pm, thereby enhancing Lorentz forces relative to inertia. This work, albeit in a Boussinesq framework, supports the view that magnetically dominated convection can favor solar-like differential rotation.

A key physical ingredient in Hotta et al. (2022)’s simulation is an efficient small-scale dynamo (SSD), which produces super-equipartition magnetic fields at small scales (see, Rempel et al. 2023, for a review on SSD). Although it has sometimes been doubted whether such an efficient SSD is possible in the Sun’s convection zone where Pm is very small (e.g., Käpylä et al. 2018), Warnecke et al. (2023) recently showed that a SSD is possible in such a low-Pm regime. Since no current simulations have yet achieved the asymptotic regime of SSD with low Pm, the actual efficiency of the SSD in the Sun is still unclear. Nonetheless, given the enormous magnetic Reynolds number in the Sun’s convection zone, Rm ≈ 𝒪(10¹⁰) (Ossendrijver 2003), it is very likely that the SSD operates more vigorously in the Sun than in any numerical models. So far, most of the numerical investigations of SSD have been carried out in nonrotating systems (e.g., Rempel 2014; Hotta et al. 2015b; Käpylä et al. 2018; Yan et al. 2021; Bhatia et al. 2022). Some numerical studies of rotating magneto-convection have reported the existence of SSD (Favier & Bushby 2012; Hotta et al. 2015a, 2016; Warnecke et al. 2025) and the associated Maxwell stress that opposes the Reynolds stress (Nelson et al. 2013; Fan & Fang 2014; Augustson et al. 2015; Käpylä et al. 2017a; Warnecke et al. 2025). However, few studies have focused on how SSD affects the properties of rotating columnar convection in the solar and stellar convection zones.

In this paper, we use a local Cartesian f-plane box model of rotating compressible convection at the equator to investigate the interaction between SSD and the columnar convective (thermal Rossby) modes. A similar local box model has been employed in previous studies of rotating magneto-convection (Käpylä et al. 2009; Favier & Bushby 2012; Masada & Sano 2014; Bushby et al. 2018). In these studies, the box was located at the pole, with the rotational axis being parallel to the gravity. Such a configuration provides helical turbulence, which can give rise to α²-type large-scale dynamo. In this study, however, we exclude the α effect from our simulations by placing the f-plane box at the equator, which enables us to focus on the influence of the SSD on the columnar convective modes arising from the background density stratification. Nevertheless, we note that this setup does not entirely rule out the possibility of large-scale dynamo action, because non-helical shear dynamos are also possible (e.g., Yousef et al. 2008; Rogachevskii & Kleeorin 2003).

The organization of the paper is as follows. In Sect. 2, our numerical model is explained. The results are presented in Sect. 3. The conclusions are summarized in Sect. 4.

2. Numerical model

2.1. Basic equations

In this study, we numerically solved 3D magnetohydrodynamic (MHD) equations in a local Cartesian f-plane box. In our definition, the x axis is directed from west to east (longitudinal) and the y axis from south to north (latitudinal), and the z axis is directed radially upward. The basic equations consist of the equation of continuity, the equation of motion, the induction equation, and the entropy equation

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_1}{\partial t}& =-\mathrm{\nabla }·[({\rho }_0+{\rho }_1)\mathit{v}],\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\rho }_0\frac{\partial \mathit{v}}{\partial t}& =-{\rho }_0\mathit{v}·\mathrm{\nabla }\mathit{v}-\mathrm{\nabla }{p}_1+{\rho }_1\mathit{g}+2{\rho }_0\mathit{v}\times {\mathbf{\Omega }}_0\hfill \\ \hfill & +\frac{1}{4\pi }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =\mathrm{\nabla }\times (\mathit{v}\times \mathit{B}),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\rho }_0{T}_0\frac{\partial s_1}{\partial t}& =-{\rho }_0{T}_0\mathit{v}·\mathrm{\nabla }{s}_1+Q_{\mathrm{rad}}.\hfill \end{array}$$(4)

Here, ρ₀, p₀, and T₀ denote values of density, pressure, and temperature of a time-independent reference background stratification, respectively. The reference background is assumed to be in an adiabatically stratified hydrostatic equilibrium with the gravitational acceleration g = −ge_z, where g is assumed to be spatially constant. The thermodynamic variables with subscript 1 (ρ₁, p₁, and T₁) represent perturbations that are small compared with the background values so that the equation of state can be linearized as

$$\begin{array}{c}\hfill p_1=p_0(\gamma \frac{{\rho }_1}{{\rho }_0}+\frac{s_1}{c_{\mathrm{v}}}).\end{array}$$(5)

For simplicity, we assume an ideal gas law with the specific heat ratio γ = c_p/c_v = 5/3, where c_p and c_v denote the specific heat at constant pressure and at constant volume, respectively. In this paper, we locate the numerical box exactly at the equator so that the rotational axis is parallel to the y axis, Ω₀ = Ω₀e_y, as is shown in Fig. 1.

Fig. 1. Sketch of the local Cartesian model of the stellar convection zone. The local f-plane box is located at the equator, where the x, y, and z axes point to the longitudinal, latitudinal, and radial directions, respectively.

2.2. Background stratification and radiative heating

The time-independent background quantities ρ₀, p₀, and T₀ are given by an adiabatically stratified polytrope as (e.g., Fan et al. 1999)

$$\begin{array}{cc}\hfill {\rho }_0(z)& ={\rho }_{\mathrm{r}}{[1-\frac{z}{(1+m)H_{\mathrm{r}}}]}^m,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill p_0(z)& =p_{\mathrm{r}}{[1-\frac{z}{(1+m)H_{\mathrm{r}}}]}^{1+m},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill T_0(z)& =T_{\mathrm{r}}[1-\frac{z}{(1+m)H_{\mathrm{r}}}],\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill H_0(z)& =\frac{p_0}{{\rho }_{0}g},\hfill \end{array}$$(9)

where m = 1/(γ − 1) is a polytropic index and ρ_r, p_r, T_r, and H_r denote the values of ρ₀, p₀, T₀, and the pressure scale height, H₀, at the bottom of our numerical domain, z = 0.

In our model, the convective instability is driven by a prescribed radiative heating and cooling term through which the radiative energy flux is injected from the bottom boundary and extracted through the upper boundary:

$$\begin{array}{c}\hfill Q_{\mathrm{rad}}(z)=-\frac{\partial }{\partial z}[F_{\mathrm{heat}}(z)+F_{\mathrm{cool}}(z)].\end{array}$$(10)

We use the same functional forms for F_heat(z) and F_cool(z) as those presented in Bekki et al. (2017). The radiative heating part is assumed to be proportional to the background pressure as (Featherstone & Hindman 2016)

$$\begin{array}{c}\hfill -\frac{\partial F_{\mathrm{heat}}}{\partial z}=f_{\mathrm{heat}}[p_0(z)-p_0(z_{\max })],\end{array}$$(11)

where the normalization factor, f_heat, is determined by the input energy flux, F_*, as

$$\begin{array}{c}\hfill f_{\mathrm{heat}}=F_{\ast }{\left[{\displaystyle {\int }_0^{z_{\max }}}(p_0(z)-p_0(z_{\max }))dz\right]}^{-1}.\end{array}$$(12)

The radiative cooling flux, F_cool, is assumed to be a Gaussian function localized near the surface

$$\begin{array}{c}\hfill F_{\mathrm{cool}}=F_{\ast }exp[-{\left(\frac{z-z_{\max }}{d_{\mathrm{cool}}}\right)}^2],\end{array}$$(13)

where the thickness of the surface cooling layer is set to be d_cool = 0.3 H_r, which is roughly the local pressure scale height at the top boundary.

2.3. Dimensionless parameters

The typical convective velocity, v_*, is estimated as

$$\begin{array}{c}\hfill v_{\ast }={\left(\frac{F_{\ast }}{{\rho }_r}\right)}^{1/3}.\end{array}$$(14)

The typical convective turnover timescale, τ_*, and the typical magnetic field strength, B_*, are then estimated, respectively, as

$$\begin{array}{c}\hfill {\tau }_{\ast }=\frac{H_r}{v_{\ast }},B_{\ast }=v_{\ast }\sqrt{4\pi {\rho }_r}.\end{array}$$(15)

In our model, the radiative energy flux, F_*, is set by specifying the modified Mach number, M_* (square root of the Euler number),

$$\begin{array}{c}\hfill {\mathrm{M}}_{\ast }=\frac{v_{\ast }}{\sqrt{p_r/{\rho }_r}}.\end{array}$$(16)

Although M_* is estimated to be ≈𝒪(10⁻³ − 10⁻⁴) in the Sun’s convection zone (Ossendrijver 2003), we use the value of M_* = 10⁻² throughout this paper in order to avoid the severe CFL condition, while keeping thermal convection sufficiently subsonic. The rotational influence on convection is controlled by

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\ast }=\frac{2{\mathrm{\Omega }}_0{H}_r}{v_{\ast }},\end{array}$$(17)

which we vary as a free input parameter from 1 to 20 as shown in Table 1. We note that the parameter Co_* is defined a priori on measurable quantities of the Sun (such as the luminosity, L_⊙, and rotation rate, Ω_⊙) and is called a flux Coriolis number (e.g., Aurnou et al. 2020; Käpylä 2024). For reference, the solar value of Co_* is estimated to be Co_⊙ ≈ 3.8 (with F_⊙ = 1.2 × 10¹¹ g s⁻³, Ω_⊙ = 2.8 × 10⁻⁶ s⁻¹, H_r = 57 Mm, and ρ_r = 0.2 g cm⁻³).

2.4. Numerical methods

We numerically solved Eqs. (1)–(4) using the fourth-order central-differencing method and four-step Runge-Kutta scheme (e.g., Vögler et al. 2005). In order to keep the divergence-free condition of the magnetic field, the nine-wave method is used (Dedner et al. 2002). In this study, we do not use any explicit diffusivities but only use the slope-limited artificial diffusion proposed by Rempel (2014) to stabilize the computations (see Appendix A for details). The kinetic and magnetic energy dissipated through this artificial diffusion is returned to the internal energy. We note that this approach leads to effective thermal and magnetic Prandtl numbers in our simulations of order unity, both of which are much higher than those estimated in the real Sun.

The size of the numerical domain is (L_x, L_y, L_z)/H_r = (19.8, 6.6, 2.2), which is sufficiently large in the longitudinal (x) direction to capture the columnar convective pattern. The density contrast between the bottom and top of the numerical domain is about 24, with the number of density scale heights N_ρ ≈ 3.2. We use a periodic boundary condition for horizontal directions and an impenetrable, stress-free boundary condition at the vertical boundaries. The magnetic field is assumed to be horizontal at the bottom (z = 0) and vertical at the top boundary (z = L_z). For all our simulations, we use the grid resolution of (N_x, N_y, N_z) = (1560, 520, 192). Owing to the nature of our diffusion scheme, this resolution determines the effective viscous, thermal, and magnetic diffusivities in the simulations. Each simulation was initiated by giving a small random perturbation to the vertical velocity, v_z, while setting the other hydrodynamic (HD) variables ρ₁, v_x, v_y, and s₁ to zero. A weak seed field of (B_x, B_y, B_z)/B_* = (0, 10⁻⁵, 0) was imposed to initiate the SSD at t = 0 in the MHD cases.

Fig. 4. Same as Fig. 3 but from the simulation cases HD-Co2 and MHD-Co2. An animation of this figure is available online.

3. Results

3.1. Simulation overview

As is summarized in Table 1, we carried out sets of nonmagnetic (denoted by “HD”) and magnetic (denoted by “MHD”) simulations with five different rotational influences (Co_* = 1, 2, 5, 10, 20). It is generally known that increasing Co_* leads to a reduction in the supercriticality of convection (e.g., Chandrasekhar 1961). Exploring the effects of varying supercriticality requires huge numerical resources and is thus beyond the scope of this paper. Since the f-plane box is located exactly at the equator, the mean flow is only established in the x (longitudinal) direction. Furthermore, we note that our MHD simulations show no clear evidence of large-scale dynamo action.

Figure 2 shows the temporal evolution of the volume-integrated kinetic and magnetic energies, E_kin and E_mag, from the MHD cases. They are defined by

Fig. 2. Temporal evolution of the volume-integrated kinetic and magnetic energies, Ekin (dashed curves) and Emag (solid curves), from the MHD runs. Different colors represent the results with different Co*.

$$\begin{array}{cc}\hfill E_{\mathrm{kin}}& ={\displaystyle {\int }_V}\frac{{\rho }_0}{2}(v_x^2+v_y^2+v_z^2)dV,\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill E_{\mathrm{mag}}& ={\displaystyle {\int }_V}\frac{B_x^2+B_y^2+B_z^2}{8\pi }dV,\hfill \end{array}$$(19)

where V = L_xL_yL_z is the total volume of the numerical domain. All our simulations operate in a highly turbulent regime with effective Reynolds and magnetic Reynolds numbers, Re_*,  Rm_* ≈ 10³ (see Appendix A for details), which are sufficiently large to ensure vigorous convection and efficient SSD action. All the simulations were run for about 80 τ_*. The SSDs are saturated and the statistically stationary states are reached by t ≈ 20 τ_*. In what follows, we focus on these statistically stationary states.

Figure 3 shows the overall spatial patterns of the vertical velocity, v_z, and vertical magnetic field, B_z, in the HD-Co1 and MHD-Co1 cases in which the rotational influence is weak (Co_* = 1). In the MHD case (Fig. 3c), the results are very similar to those of nonrotating SSD simulations (Cattaneo 1999; Vögler et al. 2005; Rempel 2014; Hotta et al. 2015b). Strong mixed-polarity magnetic fields are predominantly concentrated in the narrow downflow regions that are surrounded by slower and broader upflows. A comparison of the HD and MHD cases (Figs. 3a and b) also reveals that the convective pattern in MHD becomes much smoother than that of HD simulation. This can be attributed to the back-reaction from the SSD, i.e., the small-scale Lorentz force suppresses the small-scale turbulence (Hotta et al. 2015b).

Fig. 3. Temporal snapshots of (a) the vertical velocity, vz, from case HD-Co1, (b) vz from case MHD-Co1, and (c) the vertical magnetic field, Bz, from MHD-Co1 in the statistically stationary states at t ≈ 80 τ*. The top and middle panels show the horizontal cuts near the top surface, z = 2.2 Hr, and in the middle convection zone, z = 1.2 Hr. The bottom panels show the vertical cuts at y = 0. An animation of this figure is available online.

Figures 4, 5, 6, and 7 show the temporal snapshots of v_z and B_z in the same way as in Fig. 3 but with increasing Co_*. As Co_* increases, the columnar pattern elongated in the y direction becomes more and more apparent. The size (longitudinal extent) of the convective columns becomes smaller with the increase in Co_*. In other words, the total number of columns in the x direction increases as convection becomes more and more rotationally constrained. Furthermore, a comparison between the HD and MHD cases again shows that the columnar pattern tends to be much more coherent and self-organized in the magnetized convection (e.g., Figs. 7a and b).

Fig. 5. Same as Fig. 3 but from the simulation cases HD-Co5 and MHD-Co5. An animation of this figure is available online.

Fig. 6. Same as Fig. 3 but from the simulation cases HD-Co10 and MHD-Co10. An animation of this figure is available online.

Fig. 7. Same as Fig. 3 but from the simulation cases HD-Co20 and MHD-Co20. An animation of this figure is available online.

3.2. Notation of averages, Fourier spectra, and error estimates in this paper

In the following sections, we discuss several statistical properties of rotating magneto-convection. To this end, we define the following averaging procedures. The horizontal (x and y) average ⟨q⟩, the vertical (z) average $\overline{q}$, and the latitudinal (y) average $\stackrel{\sim }{q}$ of a quantity, q, are defined as

$$\begin{array}{cc}\hfill ⟨q⟩(z)& =\frac{1}{L_x{L}_y}{\displaystyle {\int }_0^{L_x}{\int }_0^{L_y}}q(x,y,z)dxdy,\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill \overline{q}(x,y)& =\frac{1}{L_z}{\displaystyle {\int }_0^{L_z}}q(x,y,z)dz,\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill \stackrel{\sim }{q}(x,z)& =\frac{1}{L_y}{\displaystyle {\int }_0^{L_y}}q(x,y,z)dy.\hfill \end{array}$$(22)

Fluctuations with respect to the horizontally averaged values are quoted by primes as

$$\begin{array}{c}\hfill q^{\prime }(x,y,z)=q(x,y,z)-⟨q⟩(z).\end{array}$$(23)

In addition to this, we explicitly denote the fluctuations with respect to the y-averaged quantities by y primes as

$$\begin{array}{c}\hfill q^{\prime (y)}(x,y,z)=q(x,y,z)-\stackrel{\sim }{q}(x,z).\end{array}$$(24)

We also define a Fourier transform in the horizontal directions as

$$\begin{array}{c}\hfill q(x,y,z)={\displaystyle \sum _{k_x,k_y}}\widehat{q}(k_x,k_y,z)e^{\mathrm{i}(k_{x}x+k_{y}y)},\end{array}$$(25)

where k_x and k_y are wavenumbers in the x and y directions. We note that discretized wavenumbers (k_x, k_y) = (2πn_x/L_x, 2πn_y/L_y) with integers n_x, n_y (≥0) are used.

Hereafter, we discuss the statistical properties averaged over time between 60 ≤ t/τ_* ≤ 80, with the standard errors given by $\sigma /\sqrt{N_t-1}$, where σ denotes the standard deviation and N_t is the number of temporal snapshots used. To ensure that the snapshots are statistically independent, we used the data at intervals of 1.0 τ_*, i.e., N_t = 20.

3.3. Energy spectra

Figure 8 exemplifies the 2D kinetic energy spectra in the middle convection zone (z = 1.1 H_r) for weak rotation (Co_* = 1) and strong rotation (Co_* = 20) cases. Here, the kinetic energy spectrum, ${\widehat{e}}_{\mathrm{kin}}$, is defined as

Fig. 8. Two-dimensional kinetic energy spectra shown in a horizontal wavenumber space (kx, ky) in the middle convection zone (z = 1.1 Hr) for (a) weakly rotating and (b) rapidly rotating cases. The top and bottom rows show the results of HD and MHD simulations, respectively. The insets show cuts of the same spectra at fixed ky, where the red, blue, green, and gray curves correspond to cuts through kyLy/2π = 0,  1,  2, and 3.

$$\begin{array}{c}\hfill {\widehat{e}}_{\mathrm{kin}}(\mathit{k},z)=\frac{{\rho }_0}{2}\widehat{\mathit{v}}(\mathit{k},z)·{\widehat{\mathit{v}}}^{\ast }(\mathit{k},z).\end{array}$$(26)

By definition, this satisfies

$$\begin{array}{c}\hfill \langle \frac{{\rho }_0}{2}{\mathit{v}}^2\rangle (z)={\displaystyle \sum _{\mathit{k}}}{\widehat{e}}_{\mathrm{kin}}(\mathit{k},z).\end{array}$$(27)

Although the distribution of convective power in wavenumber space (k_x, k_y) is almost isotropic in Co_* = 1 cases, the spectra become strongly anisotropic in Co_* = 20 cases wherein most power is concentrated in the Fourier domain where the latitudinal wavenumber, k_y, is small. In fact, as is indicated by the inset panels in Fig. 8b, the k_y = 0 component of the convective motions is shown to be dominant over the k_y ≠ 0 components. This corresponds to a formation of coherent columnar patterns aligned in the y direction (Figs. 6 and 7).

Figures 9a and b show the kinetic energy spectra summed over k_y for HD and MHD cases for different values of Co_*. When Co_* is small, the kinetic energy spectra roughly follow the Kolmogorov power law $\propto k_x^{-5/3}$ in the inertial spectral range. However, when Co_* is increased, ${\widehat{e}}_{\mathrm{kin}}$ begins to show spectral peaks at higher longitudinal wavenumbers. In MHD simulations, we also show the magnetic energy spectra, ${\widehat{e}}_{\mathrm{mag}}$, defined by

Fig. 9. (a) Kinetic energy spectra, ${\displaystyle \sum _{k_y}}{\widehat{e}}_{\mathrm{kin}}$, as functions of longitudinal wavenumber kx in the middle convection zone (z = 1.1 Hr) for HD cases. Different colors correspond to different Co*. The dotted black line shows a $k_x^{-5/3}$ power law. (b) Same kinetic energy spectra as panel (a) but for MHD cases. Dotted curves show the magnetic energy spectra, ${\displaystyle \sum _{k_y}}{\widehat{e}}_{\mathrm{mag}}$, at the same depth. (c) Mean wavenumber, ${\overline{k}}_{x,\mathrm{mean}}$, defined by Eq. (29) as functions of Co*. The blue and red colors correspond to the HD and MHD cases, respectively. The dotted gray lines show $\propto {\mathrm{Co}}_{\ast }^{3/5}$ and $\propto {\mathrm{Co}}_{\ast }^{1/3}$ dependences theoretically expected from the CIA and MAC balances. The inset shows the wavenumbers where the kinetic energy spectra peak, ${\overline{k}}_{x,\mathrm{peak}}$.

$$\begin{array}{c}\hfill {\widehat{e}}_{\mathrm{mag}}(\mathit{k},z)=\frac{1}{8\pi }\widehat{\mathit{B}}(\mathit{k},z)·{\widehat{\mathit{B}}}^{\ast }(\mathit{k},z).\end{array}$$(28)

We find that the SSD generates super-equipartition (i.e., ${\widehat{e}}_{\mathrm{kin}}\lesssim {\widehat{e}}_{\mathrm{mag}}$) magnetic fields at small scales. Consequently, the kinetic energy, ${\widehat{e}}_{\mathrm{kin}}$, is suppressed due to the Lorentz force feedback at these scales.

It is important to note that, in our simulations, there are two distinct characteristic length scales. One is the mean longitudinal wavenumber, k_x, mean, weighted by the kinetic energy spectrum (e.g., Christensen & Aubert 2006; Käpylä 2024),

$$\begin{array}{c}\hfill k_{x,\mathrm{mean}}(z)={\displaystyle \frac{{\sum }_{\mathit{k}}{k}_x{\widehat{e}}_{\mathrm{kin}}(\mathit{k},z)}{{\sum }_{\mathit{k}}{\widehat{e}}_{\mathrm{kin}}(\mathit{k},z)}},\end{array}$$(29)

and the other is the peak longitudinal wavenumber, k_x, peak, at which the kinetic energy spectrum reaches its maximum. In previous studies of rotating f-plane simulations in which gravity is aligned with the rotation axis, these two wavenumbers become asymptotically equivalent in the rapidly rotating limit (e.g., Käpylä 2024). However, in our model setup, in which both small-scale convective motions and large-scale columnar convective modes coexist, they do not coincide. Figure 9c shows the height-averaged mean wavenumber, ${\overline{k}}_{x,\mathrm{mean}}$, and the peak wavenumber, ${\overline{k}}_{x,\mathrm{peak}}$, as functions of Co_*. In general, the inclusion of the SSD leads to smaller ${\overline{k}}_{x,\mathrm{mean}}$ in MHD (i.e., the typical convective scale ${\overline{\ell }}_x=2\pi /{\overline{k}}_{x,\mathrm{mean}}$ becomes larger) than in HD because the SSD Lorentz force suppresses small-scale turbulence. In contrast, the peak wavenumber, ${\overline{k}}_{x,\mathrm{peak}}$, is shown to be systematically larger in MHD at Co_* ≥ 5. We find that the peak wavenumber, ${\overline{k}}_{x,\mathrm{peak}}$, largely represents the dominant scale of the k_y = 0 columnar convective (thermal Rossby) modes. The effects of the SSD on ${\overline{k}}_{x,\mathrm{peak}}$ are discussed later in Sect. 3.6.

3.4. Convective heat transport

Figures 10a–d show vertical profiles of the root-mean-square (rms) convective velocity

Fig. 10. Profiles of the rms convective velocity, vrms, and the Alfvén velocity of the rms magnetic field, $v_{\mathrm{A},\mathrm{rms}}=B_{\mathrm{rms}}/\sqrt{4\pi {\rho }_0}$. (a–d) Vertical profiles of vrms for different Co*, where blue and red represent the results from HD and MHD cases. The dashed green curves show the profiles of vA, rms from MHD cases. (e) The values of volume-averaged rms velocity and magnetic field, ${\overline{v}}_{\mathrm{rms}}$ and ${\overline{v}}_{\mathrm{A},\mathrm{rms}}$, as functions of Co*. The dotted gray lines show $\propto {\mathrm{Co}}_{\ast }^{-1/5}$ and $\propto {\mathrm{Co}}_{\ast }^{-1/3}$ dependences theoretically expected from the CIA and MAC balances.

$$\begin{array}{c}\hfill v_{\mathrm{rms}}(z)=\sqrt{⟨{(\mathit{v}(x,y,z)-U_x(z){\mathit{e}}_x)}^2⟩},\end{array}$$(30)

where U_x = ⟨v_x⟩ denotes the mean flow in the x direction. The generation mechanisms of the mean flows are discussed in a later section (Sect. 3.7). From MHD simulations, we also show the Alfvén velocity associated with the rms magnetic field,

$$\begin{array}{c}\hfill v_{\mathrm{A},\mathrm{rms}}(z)=\sqrt{\frac{⟨\mathit{B}{(x,y,z)}^2⟩}{4\pi {\rho }_0(z)}}·\end{array}$$(31)

It is evident that, due to the SSD Lorentz force feedback, v_rms in the MHD cases is suppressed by approximately 15–25% compared to the HD cases. In most of our simulations, the rms field strength, B_rms, is sub-equipartition, i.e., v_A, rms < v_rms. Although this is in accordance with the simulations by Favier & Bushby (2012) and Warnecke et al. (2025), we note that the higher-resolution simulation by Hotta et al. (2022) produces super-equipartition SSD fields throughout the whole convection zone.

It is generally known that strong rotation tends to suppress convective velocities (e.g., Aurnou et al. 2020; Vasil et al. 2021; Käpylä 2024). As is shown in Fig. 10e, our simulations also demonstrate that the volume-averaged rms convective velocity, ${\overline{v}}_{\mathrm{rms}}$, decreases with increasing Co_* in the rapidly rotating regime (Co_* ≥ 5) in both HD and MHD cases. Moreover, our study reveals that this suppression is further enhanced by the presence of SSD. In HD, ${\overline{v}}_{\mathrm{rms}}$ decreases gradually with Co_*, roughly following a power law of Co_*^−1/5, while in MHD, the decline is steeper, following Co_*^−1/3, as indicated by the dotted gray lines in Fig. 10e. The underlying mechanism by which the SSD modifies the Co_* dependence is explained in Sect. 3.5.

Next, we discuss thermal energy transport properties in our simulations. Figures 11a–d present snapshots of the entropy perturbations, ${\rho }_{0}s{\prime }_1$, for Co_* = 2 and 20. Figure 11e further shows the vertically averaged rms entropy perturbations, $\overline{{\rho }_0{s}_{\mathrm{rms}}}$, as a function of Co_*. Here, s_rms is defined by

Fig. 11. Entropy perturbations, s′1, in our simulations. (a,b) Temporal snapshots of s′1 weighted by the background density, ρ0, at fixed height (z = 1.2 Hr) and latitude (y = 0) from weakly rotating simulations HD-Co2 and MHD-Co2. (c,d) Same as panels a and b but from rapidly rotating cases HD-Co20 and MHD-Co20. (e) Volume-averaged rms entropy perturbation, $\overline{{\rho }_0{s}_{\mathrm{rms}}}$, as a function of Co*. Blue and red denote the results from HD and MHD simulations. The dotted gray lines show $\propto {\mathrm{Co}}_{\ast }^{1/5}$ and $\propto {\mathrm{Co}}_{\ast }^{1/3}$ dependences theoretically expected from the CIA and MAC balances.

$$\begin{array}{c}\hfill s_{\mathrm{rms}}(z)=\sqrt{⟨s_1^{\prime }{(x,y,z)}^2⟩}.\end{array}$$(32)

We note that the entropy perturbations shown are weighted by the background density, ρ₀, which allows us to focus on the entropy variations within the bulk convection zone, rather than those generated by strong near-surface radiative cooling. Two main effects are observed. First, we find an enhancement of the entropy perturbations, s′₁, due to the SSD action. This behavior was reported in earlier nonrotating studies (e.g., Hotta et al. 2015b), in which it was shown that the SSD magnetic fields suppress horizontal turbulent mixing of entropy between warm upflows and cold downflows, thereby enhancing entropy contrast. At small Co_*, the MHD simulations show an increase in s_rms by approximately 25% compared to the HD cases. The second effect is an enhancement of s′₁ by strong rotation at Co_* ≥ 5. It is well known that the Coriolis force makes the convective heat transport inefficient as it deflects vertical motions and reduces the entropy mixing in a vertical direction. To sustain a fixed enthalpy flux in the domain, the entropy perturbations, s′₁, are enhanced. As is shown in Fig. 11e, the volume-averaged rms entropy perturbation, $\overline{{\rho }_0{s}_{\mathrm{rms}}}$, increases with Co_* in both the HD and MHD cases. However, the scaling differs: the MHD cases follow an approximate scaling of Co_*^1/3, while the HD cases follow Co_*^1/5, indicating a stronger Co_* dependence with the presence of SSD. Obviously, these scalings are the inverse of those found in the rms convective velocity, v_rms (Fig. 10e), implying that the entropy perturbations are enhanced to compensate for the reduced convective velocities and to transport the fixed amount of enthalpy upward.

The efficiency of convective energy transport is known to be associated with the mean entropy stratification in the bulk convection zone. Figures 12a–d show vertical profiles of the superadiabaticity, δ, defined as

Fig. 12. Profiles of the superadiabaticity, δ. (a–d) Vertical profiles of δ for different Co*. Blue and red curves represent the results from HD and MHD simulations. (e) Volume-averaged values of superadiabaticity, $\overline{\delta }$, as a function of Co*. The inset shows the magnitudes of the subadiabatic layers near the bottom of the convection zone, Dsub = ∫0dsub|δ|dz (with dsub being the height below which the mean stratification is subadiabatic), denoted by the shaded areas in panels a–d.

$$\begin{array}{c}\hfill \delta =-\frac{H_0}{c_{\mathrm{p}}}\frac{\partial ⟨s_1⟩}{\partial z}·\end{array}$$(33)

When the rotational influence is weak, a weakly subadiabatic layer is formed in the lower half of the convection zone. This is because of the nonlocal heat transport by strong downflow plumes (e.g., Brandenburg 2016) and has been repeatedly reported in previous simulations of nonrotating convection (e.g., Käpylä et al. 2017b; Bekki et al. 2017; Karak et al. 2018; Käpylä 2025). It is also seen that the magnitude and width of this subadiabatic layer are enhanced and extended with the inclusion of SSD. When Co_* is increased, the mean stratification becomes more superadiabatic in the bulk convection zone because the rotation suppresses vertical entropy mixing. These trends are clearly seen in Fig. 12e. Interestingly, we find that the SSD makes the subadiabatic layer near the base of the convection zone persistent even under strong rotational influence (see the inset of Fig. 12e). This result is in striking contrast to the HD case in which the subadiabatic layer vanishes in the rapidly rotating regime. A disappearance of the subadiabatic layer in the rapidly rotating HD regime has also been reported by Käpylä (2024), who employed a radiative diffusivity based on Kramers opacity (as a function of temperature and density) instead of the spatially fixed radiative heating used in this study.

Our parameter study suggests that the weakly subadiabatic layer likely exists in the lower 40% of the Sun’s convection zone near the equator (given that Co_⊙ ≈ 3.8). A similar result is obtained by Hotta et al. (2022). We also note that the existence of this weakly subadiabatic layer in the Sun’s convection zone has been inferred by recent observations and linear modeling of the non-toroidal class of the solar inertial modes (Hanson et al. 2022; Bekki 2024).

3.5. Force balance

To understand why the convective velocity, v_rms, and entropy fluctuation, s_rms, show different Co_* dependences with and without the SSD, we investigate the dominant force balance in our simulations. We follow the method presented in Yadav et al. (2016), Aguirre Guzmán et al. (2021) and compute the rms values of the following forces:

$$\begin{array}{cc}\hfill F_{\mathrm{Adv}}& =-{\rho }_0\mathit{v}·\mathrm{\nabla }\mathit{v},\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill F_{\mathrm{Pre}}& =-{\rho }_0\mathrm{\nabla }\left(\frac{p_1}{{\rho }_0}\right),\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill F_{\mathrm{Buo}}& ={\rho }_0\frac{s_1}{c_{\mathrm{p}}}g,\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill F_{\mathrm{Cor}}& =2{\rho }_0\mathit{v}\times {\mathbf{\Omega }}_0,\hfill \end{array}$$(37)

$$\begin{array}{cc}\hfill F_{\mathrm{Mag}}& =\frac{1}{4\pi }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}.\hfill \end{array}$$(38)

Here, we note that temporally averaged component of the pressure and buoyancy forces (which represents the background hydrostatic balance and thus is irrelevant to our analysis) are subtracted from F_Pre and F_Buo. Figure 13a shows the rms forces in HD (top) and MHD (bottom) cases as functions of Co_*. In HD simulations, the advection (or inertia) term F_Adv always dominates the force balance. As Co_* is increased, a balance is achieved among advection, F_Adv, pressure, F_Pre, and Coriolis forces, F_Cor. On the other hand, in MHD simulations, the Lorentz force, F_Mag, comes into play in the leading-order force balance, where F_Mag has amplitudes very close to those of the pressure gradient force, F_Pre. This likely reflects the fact that the strong SSD fields are generated by compression of downflows where the gas pressure is locally balanced by the magnetic pressure (Hotta et al. 2022). At sufficiently high Co_*, F_Mag is even stronger than F_Adv, leading to a balance between the pressure, F_Pre, Coriolis, F_Cor, and Lorentz forces, F_Mag. This balance is called magnetostrophic (e.g., Roberts 1978). We can attribute the Co_* dependences of v_rms and s_rms discussed in Sect. 3.4 to this transition of the dominant force balance regime.

Fig. 13. Force balances in our simulations. (a) RMS amplitudes of the forces, Frms. Here, Frms is calculated such that Frms2 = V−1∫V(Fx′2 + Fy′2 + Fz′2)dV. Yellow, red, green, blue, and purple colors represent the nonlinear advection, pressure gradient, buoyancy, Coriolis, and Lorentz forces. The top and bottom panels show the HD and MHD simulations, respectively. (b) RMS amplitudes of the y-averaged forces, ${\stackrel{\sim }{F}}_{\mathrm{rms}}$. They are defined such that ${\stackrel{\sim }{F}}_{\mathrm{rms}}^2=V^{-1}{\displaystyle {\int }_V}({\stackrel{\sim }{F}}_x^{\prime 2}+{\stackrel{\sim }{F}}_z^{\prime 2})dV$. (c) Power spectra of the y-averaged forces in the middle height, z = 1.1 Hr, from HD-Co20 and from MHD-Co20 simulations. The dashed gray curves denote the ageostrophic Coriolis forces, ${\stackrel{\sim }{F}}_{\mathrm{Pre}}+{\stackrel{\sim }{F}}_{\mathrm{Cor}}$.

It is worth noting that the relative importance of the Coriolis force, F_Cor (which does not have a y component), tends to be underestimated in Fig. 13a. To properly assess the significance of the Coriolis force in the system, we also compute the same rms forces but for the y-averaged ones, as shown in Fig. 13b. We can see that there is a clear distinction between the weakly rotating regime, Co_* < 5, where the Coriolis force, ${\stackrel{\sim }{F}}_{\mathrm{Cor}}$, is not in the leading-order force balance and the strongly rotating regime, Co_* ≳ 5, where the dominant force balance is achieved between the Coriolis force, ${\stackrel{\sim }{F}}_{\mathrm{Cor}}$, and the pressure gradient force, ${\stackrel{\sim }{F}}_{\mathrm{Pre}}$. When this geostrophic balance is established at a leading order, it is necessary to look at the residual balance (e.g., Yadav et al. 2016; Aubert et al. 2017). Figure 13c displays the power spectra of the y-averaged forces in the middle convection zone (z = 1.1 H_r) from rapidly rotating HD-Co20 and MHD-Co20 cases. Here, the imbalance between Coriolis and pressure gradient forces, ${\stackrel{\sim }{F}}_{\mathrm{Cor}}+{\stackrel{\sim }{F}}_{\mathrm{Pre}}$, is denoted by dashed gray curves. In the HD case, this so-called ageostrophic Coriolis force is balanced by buoyancy, ${\stackrel{\sim }{F}}_{\mathrm{Buo}}$, and advection, ${\stackrel{\sim }{F}}_{\mathrm{Adv}}$, at large scales and predominantly by advection at small scales, manifesting the Coriolis-Inertia-Archimedean (CIA) balance (see, e.g., Aurnou et al. 2020; Käpylä 2024). In the MHD case, on the other hand, the ageostrophic Coriolis force is primarily balanced by buoyancy, ${\stackrel{\sim }{F}}_{\mathrm{Buo}}$, at large scales but is balanced by the Lorentz force, ${\stackrel{\sim }{F}}_{\mathrm{Mag}}$, at small scales, manifesting the magneto-Archimedean-Coriolis (MAC) balance (see, e.g., Davidson 2013; Augustson et al. 2019; Schwaiger et al. 2021).

Assuming the CIA and MAC balances, the following scaling relations can be derived in terms of the flux Coriolis number, Co_* (see Appendix B for detailed derivations):

$$\begin{array}{cc}\hfill \mathrm{Length}\mathrm{scale}& :\ell \propto \{\begin{array}{cc}{{\mathrm{Co}}_{\ast }}^{-3/5}\hfill & (\mathrm{CIA})\hfill \\ {{\mathrm{Co}}_{\ast }}^{-1/3}\hfill & (\mathrm{MAC}),\hfill \end{array}\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill \mathrm{Velocity}& :v\propto \{\begin{array}{cc}{{\mathrm{Co}}_{\ast }}^{-1/5}\hfill & (\mathrm{CIA})\hfill \\ {{\mathrm{Co}}_{\ast }}^{-1/3}\hfill & (\mathrm{MAC}),\hfill \end{array}\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill \mathrm{Entropy}\mathrm{fluctuation}& :s\propto \{\begin{array}{cc}{{\mathrm{Co}}_{\ast }}^{1/5}\hfill & (\mathrm{CIA})\hfill \\ {{\mathrm{Co}}_{\ast }}^{1/3}\hfill & (\mathrm{MAC}).\hfill \end{array}\hfill \end{array}$$(41)

In Figs. 9c, 10e, and 11e, our simulated Co_* dependences are compared with these theoretical scalings. It is clearly seen that our results are consistent with the predictions by CIA and MAC balances at Co_* ≥ 5.

One might question the robustness of the above results, particularly when the numerical resolution is further increased. Hints are provided by Käpylä (2024), who conducted HD simulations at varying Reynolds number Re and demonstrated that the scaling of convective velocity, v_rms, with Coriolis number remains unaffected, although the overall amplitude of v_rms increases with Re. This is because viscous force does not influence the dominant force balance. Given that the CIA and MAC balances are well established in our simulations, the above scaling relations at the high-Co_* regime are expected to remain robust even at higher resolutions, although the absolute values of v_rms (or s_rms) and the degree of suppression (enhancement) from HD to MHD will be affected.

3.6. Dominant scale of columnar convection

In this subsection, we discuss how SSD affects the dominant scale of columnar convective modes. From Figs. 7a and b, it is seen that the SSD not only makes the columnar patterns more self-organized but also shifts the dominant longitudinal scale of convection narrower. To quantify the difference in the length scales of columnar convection, we compare in Fig. 14a the kinetic energy spectra between HD and MHD cases in the middle convection zone (z = 1.1 H_r). Figure 14b shows the peak longitudinal wavenumber, ${\overline{k}}_{x,\mathrm{peak}}$, as a function of Co_*. In both the HD and MHD cases, there is a general trend that ${\overline{k}}_{x,\mathrm{peak}}$ increases with Co_*. However, it is obvious that the MHD simulations show a much more rapid increase in ${\overline{k}}_{x,\mathrm{peak}}$ with Co_* compared to the HD counterparts. As a consequence, in the high-Co_* regime, the convective columns have a smaller length scale in longitude with the presence of the SSD. This behavior cannot be explained by the CIA/MAC scaling relations (Eq. 39).

Fig. 14. (a) Normalized kinetic energy spectra, ${\displaystyle \sum _{k_y}}{\widehat{e}}_{\mathrm{kin}}$, in the middle convection zone (z = 1.1 Hr). Different colors represent different Co*. Dashed and solid lines show the results from HD and MHD simulations, respectively. (b) Longitudinal wavenumbers where the kinetic energy spectra peak, ${\overline{k}}_{x,\mathrm{peak}}$, is plotted against the flux Coriolis number, Co*. The blue and red circles are the results for the HD and MHD cases. (c) Same as panel b but as a function of the dynamical Coriolis number, Co${}_{\ell }=2{\mathrm{\Omega }}_0{\overline{\ell }}_x/{\overline{v}}_{\mathrm{rms}}$.

We find that, through the suppression in the convective velocity, v_rms, and the increase in the (mean) convective length scale, ℓ_x = 2π/k_x, mean, the SSD can change the effective rotational influence in the simulations. This is typically measured by the dynamical Coriolis number,

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\ell }=\frac{2{\mathrm{\Omega }}_0{\overline{\ell }}_x}{{\overline{v}}_{\mathrm{rms}}}·\end{array}$$(42)

As is reported in Table 1, Co_ℓ is systematically larger in MHD than in HD, implying that the MHD simulations are operating in more rotationally constrained regimes compared to the HD counterparts. Consequently, the establishment of the columnar patterns is more enhanced and the longitudinal size of these columns becomes smaller in our MHD simulations. In fact, we find that the peak wavenumbers, ${\overline{k}}_{x,\mathrm{peak}}$, in both the HD and MHD cases manifest very similar dependences on the dynamical Coriolis number, Co_ℓ (Fig. 14c). This provides supporting evidence that the dominant scale of columnar convection is largely determined by Co_ℓ, a parameter that is sensitively affected by the SSD.

We also note that there are other underlying mechanisms that can potentially explain the difference in ${\overline{k}}_{x,\mathrm{peak}}$ between HD and MHD. The first possibility is direct Lorentz force feedback. As is shown in Sect. 3.8, the large-scale shear motions in the columnar convection can generate the Maxwell stress to redistribute the angular momentum in the convection zone. As a back-reaction, the kinetic energy of columnar convection is suppressed at large scales, which can also cause a shift in ${\overline{k}}_{x,\mathrm{peak}}$ toward higher wavenumbers. Another possible mechanism comes through the change in the mean background stratification. When the SSD is present, the buoyant driving can be selectively suppressed at large scales due to the reduced superadiabaticity and the existence of the weakly subadiabatic layer near the base of the convection zone (Fig. 12). Yet, it is not easy to disentangle these mechanisms and to pinpoint the primary cause. Such an analysis is beyond the scope of this paper.

3.7. Mean flow generation

In our simulations, large-scale mean flows are generated only in the longitudinal direction, U_x(z) = ⟨v_x⟩, corresponding to the radial differential rotation at the equator. Figures 15a–d show profiles of U_x(z) established in our HD and MHD simulations for different values of Co_*. A mean flow regime can be measured by the vertical shear, ΔU_x = U_x(z_max)−U_x(z_min), as is shown in Fig. 15e. For Co_* = 1,  2, and 5, the mean flows have negative vertical shear, ΔU_x < 0, with a faster base and a slower surface. This might correspond to an antisolar differential rotation profile in a spherical shell geometry. On the other hand, for Co_* = 10 and 20, the mean flows show a positive vertical shear, ΔU_x > 0, with a faster surface and a slower base, likely corresponding to a solar-like rotational profile in the spherical case. We find that the inclusion of the SSD has the primary impact of quenching the mean flow amplitudes at high Co_*, as has been reported in the literature (e.g., Brun et al. 2004; Käpylä et al. 2017a; Matilsky et al. 2022; Warnecke et al. 2025).

Fig. 15. Profiles of the mean flows, Ux(z). (a–d) Vertical profiles of Ux(z) for different Co*. Blue and red colors represent the results from the HD and MHD simulations. (e) Vertical shear amplitudes of the mean flow, ΔUx = Ux(zmax)−Ux(zmin), as a function of Co*.

Another important effect of the SSD on the mean flow, U_x, is the excitation of torsional oscillations, i.e., temporal variations in the mean flow. In our model, these oscillations arise as a by-product of the establishment of the MAC balance, whereby torsional Alfvén waves propagate vertically as disturbances to the so-called Taylor state (Taylor 1963; Braginskiy 1970; Wicht & Christensen 2010). Our analysis indicates that their amplitudes are too small (a few percent of the rms velocity fluctuation) to have a dynamical impact in the simulation. Therefore, a more detailed discussion is deferred to the Appendix C.

We now discuss the generation mechanism of the mean flows. Taking a horizontal average of the x component of Eq. (2), we have

$$\begin{array}{c}\hfill \frac{\partial U_x}{\partial t}=-\frac{\partial }{\partial z}({\mathcal{R}}_{\mathit{xz}}+{\mathcal{M}}_{\mathit{xz}}),\end{array}$$(43)

where

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathit{xz}}={\rho }_0⟨v_x^{\prime }{v}_z^{\prime }⟩,{\mathcal{M}}_{\mathit{xz}}=-\frac{⟨B_x{B}_z⟩}{4\pi },\end{array}$$(44)

are the turbulent Reynolds and Maxwell stresses, respectively. Positive and negative values of ℛ_xz and ℳ_xz correspond to vertically upward and downward fluxes of the x (angular) momentum. Figures 16a–d show vertical profiles of ℛ_xz and ℳ_xz for different Co_*. Their height-averaged values are shown as functions of Co_* in Fig. 16e. Note that, in HD, ℛ_xz must vanish in a statistically stationary state (∂U_x/∂t ≈ 0). This reflects a balance between the diffusive and non-diffusive (the Λ effect) parts of the Reynolds stress, ℛ_xz, which transport angular momentum in opposite directions (Kitchatinov 2013; Barekat et al. 2021). In MHD, ℛ_xz and ℳ_xz must balance with each other such that their sum vanishes statistically. At low Co_* (≤5) where the mean flow is antisolar (ΔU_x < 0), the Reynolds stress is negative (${\overline{\mathcal{R}}}_{\mathit{xz}}<0)$, and the Maxwell stress is positive (${\overline{\mathcal{M}}}_{\mathit{xz}}>0$). Conversely, at high Co_* (≥10) where the mean flow is solar-like (ΔU_x > 0), the signs of the Reynolds and Maxwell stresses are flipped; ${\overline{\mathcal{R}}}_{\mathit{xz}}>0$ and ${\overline{\mathcal{M}}}_{\mathit{xz}}<0$. Thus, we conclude that the mean flows are primarily generated by the angular momentum transport (AMT) through the Reynolds stress, ℛ_xz, with the Maxwell stress, ℳ_xz, only playing a suppressive role. At low Co_*, the antisolar mean flows are established by the downward AMT by weakly rotationally constrained convection (Foukal & Jokipii 1975; Guerrero et al. 2013; Featherstone & Miesch 2015), whereas at high Co_*, the upward AMT by rotationally constrained convective columns leads to the solar-like mean flows (Hotta et al. 2015a; Matilsky et al. 2019; Bekki et al. 2022b) (see also Appendix D for a schematic explanation).

Fig. 16. Profiles of the Reynolds and Maxwell stresses, defined by Eq. (44). Positive and negative values correspond to the upward and downward AMT. (a–d) Vertical profiles of the Reynolds stress, ℛxz, for different values of Co*. Blue and red colors represent the results from the HD and MHD cases. Green curves show the vertical profiles of the Maxwell stress, ℳxz, and the dashed black curves indicate the sum ℛxz + ℳxz in the MHD simulations. (e) Volume-averaged Reynolds and Maxwell stresses, ${\overline{\mathcal{R}}}_{\mathit{xz}}$ and ${\overline{\mathcal{M}}}_{\mathit{xz}}$, plotted as functions of Co*.

To assess the impact of the columnar convective modes on AMT, we decompose the Reynolds stress into contributions from the y-averaged flows (representing the k_y = 0 columnar modes), ℛ_xz^k_y = 0, and those from the velocity fluctuations, ℛ_xz^{k_y ≠ 0}, as

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathit{xz}}={\displaystyle \underset{}{\underbrace{{\rho }_0⟨{\stackrel{\sim }{v}}_x{\stackrel{\sim }{v}}_z⟩}}{\mathcal{R}}_{\mathit{xz}}^{k_y=0}}+{\displaystyle \underset{}{\underbrace{{\rho }_0⟨v_x^{\prime (y)}{v}_z^{\prime (y)}⟩}}{\mathcal{R}}_{\mathit{xz}}^{k_y\ne 0}}.\end{array}$$(45)

Figures 17a and b show the volume-averaged Reynolds stress contributions, ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y=0}$ and ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y\ne 0}$. There is a general tendency for the k_y = 0 columnar modes to transport angular momentum upward (${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y=0}>0$) and the small-scale fluctuating velocities transport angular momentum downward (${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y\ne 0}<0$). In HD cases, these two contributions cancel each other. In particular, at high Co_*, ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y=0}>0$ corresponds to a driving of the mean flow, ΔU_x > 0, and ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y\ne 0}<0$ corresponds to a turbulent diffusion. These tendencies are consistent with scale-by-scale analyses by Mori & Hotta (2023) and Käpylä (2023), who found downward AMT from small-scale turbulence and upward AMT from large-scale columnar modes near the equator. We note, however, that their analyses are based on longitudinal scale decompositions, whereas we distinguish the k_y = 0 columnar modes from the k_y ≠ 0 fluctuations; therefore, a one-to-one correspondence should not be assumed. In MHD cases, the cancellation between ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y=0}$ and ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y\ne 0}$ breaks down. At low Co_*, the k_y ≠ 0 convective motions that are weakly influenced by rotation dominantly transport angular momentum downward. As Co_* increases, the rotationally constrained columnar modes become increasingly dominant, resulting in net upward AMT (see Appendix D). Here, the downward AMT by k_y ≠ 0 flows is substantially reduced because the effective eddy diffusion due to small-scale turbulence is suppressed by the SSD Lorentz force. The resulting net downward (upward) angular momentum flux at low (high) Co_* is then compensated for by the upward (downward) flux of the Maxwell stress, where all the contributions are from k_y ≠ 0 magnetic fields (Fig. 17c).

Fig. 17. Breakdown of the contributions from the ky = 0 and ky ≠ 0 components to the Reynolds and Maxwell stresses. (a,b) Volume-averaged Reynolds stress, ${\overline{\mathcal{R}}}_{\mathit{xz}}$, as a function of Co* from the HD and MHD simulations. Yellow and purple colors correspond to ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y=0}$ and ${\overline{\mathcal{R}}}_{\mathit{xz}}^{k_y\ne 0}$, respectively. (c) Volume-averaged Maxwell stress, ${\overline{\mathcal{M}}}_{\mathit{xz}}$, from MHD simulations.

3.8. Origin of Maxwell stress: Interplay with convective columns

In this subsection, we discuss the generation mechanisms of the Maxwell stress, ℳ_xz. To this end, we compare the spatial distribution of the magnetic field correlation $\stackrel{\sim }{B_x{B}_z}$ (averaged over y) with various y-averaged quantities, as is shown in Fig. 18. Let us first focus on the low-Co_* regime in which the Maxwell stress transports angular momentum upward (${\overline{\mathcal{M}}}_{\mathit{xz}}>0$). It is seen from Figs. 18a and b that the magnetic field correlation, $\stackrel{\sim }{B_x{B}_z}$, is not necessarily generated in downflow regions where strong magnetic energy is concentrated and where the velocity correlation, $\stackrel{\sim }{v{\prime }_{x}v{\prime }_z}$, is dominantly generated. As a result, the spatial distribution of $\stackrel{\sim }{B_x{B}_z}$ is poorly correlated with that of $\stackrel{\sim }{v{\prime }_{x}v{\prime }_z}$ (Fig. 19a). This implies that the generation of the Maxwell stress is not a consequence of a parallelization of the small-scale magnetic fields with the small-scale flows, as has been proposed by Hotta et al. (2022). Rather, in our simulations, the spatial distribution of $\stackrel{\sim }{B_x{B}_z}$ is very well correlated with those of the y vorticity, $\stackrel{\sim }{\zeta }{\prime }_y$, and also of the pressure perturbation, $\stackrel{\sim }{p}{\prime }_1$. The latter is in fact a side effect required by the geostrophic balance (anticorrelation between $\stackrel{\sim }{\zeta }{\prime }_y$ and $\stackrel{\sim }{p}{\prime }_1$).

Fig. 18. Snapshots of the y-averaged quantities from MHD simulations with Co* = 2,  5,  10, and 20. Each panel, from top to bottom, displays the y-averaged profiles of the x velocity, $\stackrel{\sim }{v}{\prime }_x$, z velocity, $\stackrel{\sim }{v}{\prime }_z$, y vorticity, $\stackrel{\sim }{\zeta }{\prime }_y$, pressure perturbation, $\stackrel{\sim }{p}{\prime }_1$, magnetic energy density, ${\stackrel{\sim }{p}}_{\mathrm{mag}}$, velocity correlation, ${\rho }_0\stackrel{\sim }{v_x\prime v_z\prime }$, and magnetic field correlation, $\stackrel{\sim }{B_x{B}_z}/4\pi$. In rapidly rotating cases MHD-Co10 and MHD-Co20 (panels c–d), black curves indicate contour lines of the y vorticity at $\stackrel{\sim }{\zeta }{\prime }_y=2{\tau }_{\ast }^{-1}$.

A good correlation between the spatial distributions of $\stackrel{\sim }{B_x{B}_z}$ and $\stackrel{\sim }{\zeta }{\prime }_y$ suggests that the shear in the columnar convection is the main origin of the Maxwell stress in our simulations. From the induction Eq. (3), we obtain

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\stackrel{\sim }{B_x{B}_z}=[\dots ]+\stackrel{\sim }{B_x^2}\frac{\partial {\stackrel{\sim }{v}}_z^{\prime }}{\partial x}+\stackrel{\sim }{B_z^2}\frac{\partial {\stackrel{\sim }{v}}_x^{\prime }}{\partial z},\end{array}$$(46)

indicating that $\stackrel{\sim }{B_x{B}_z}>0(<0)$ can be generated by shear motions in the y-averaged flows, $\partial \stackrel{\sim }{v}{\prime }_x/\partial z>0(<0)$ and $\partial \stackrel{\sim }{v}{\prime }_z/\partial x>0(<0)$. Therefore, if the main generation mechanism of the Maxwell stress is the columnar shear, $\stackrel{\sim }{B_x{B}_z}$ must be strongly correlated with the velocity strain tensor, ${\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$, defined as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}=\frac{\partial {\stackrel{\sim }{v}}_x^{\prime }}{\partial z}+\frac{\partial {\stackrel{\sim }{v}}_z^{\prime }}{\partial x}·\end{array}$$(47)

This is clearly shown to be the case from Figs. 19d–f. Since the longitudinal separation between columnar downflows is large at low Co_*, the vertical shear of the zonal velocity dominates over the zonal shear of the vertical flow, $|\partial \stackrel{\sim }{v}{\prime }_x/\partial z|\gg |\partial \stackrel{\sim }{v}{\prime }_z/\partial x|$. In such a case, both ${\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$ and $\stackrel{\sim }{\zeta }{\prime }_y$ are primarily determined by $\partial \stackrel{\sim }{v}{\prime }_x/\partial z$, which explains why $\stackrel{\sim }{B_x{B}_z}$ and $\stackrel{\sim }{\zeta }{\prime }_y$ are so well correlated in Fig. 18a and b. When averaged over the whole volume, the net correlation ⟨B_xB_z⟩ becomes negative because the y vorticity consists of spatially $\stackrel{\sim }{\zeta }{\prime }_y>0$ and broadened $\stackrel{\sim }{\zeta }{\prime }_y<0$ due to the nonlinear effect (see Appendix E for its physical explanation).

Fig. 19. (Top) Correlations between the Maxwell stress, $\stackrel{\sim }{B_x{B}_z}/4\pi$, and the Reynolds stress, ${\rho }_0\stackrel{\sim }{v{\prime }_{x}v{\prime }_z}$, and (bottom) correlations between $\stackrel{\sim }{B_x{B}_z}/4\pi$ and the velocity strain tensor, ${\rho }_0{\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}={\rho }_0({\partial }_z\stackrel{\sim }{v}{\prime }_x+{\partial }_x\stackrel{\sim }{v}{\prime }_z)$, in MHD simulations. (a,b) Probability density functions (PDFs) between $\stackrel{\sim }{B_x{B}_z}/4\pi$ and ${\rho }_0\stackrel{\sim }{v{\prime }_{x}v{\prime }_z}$ in MHD-Co2 and MHD-Co20 cases. (c) Correlation coefficient as a function of Co*. (d,e) PDFs between $\stackrel{\sim }{B_x{B}_z}/4\pi$ and ${\rho }_0{\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$ in MHD-Co2 and MHD-Co20 cases. (f) Correlation coefficients $\stackrel{\sim }{B_x{B}_z}/4\pi$ vs. ${\rho }_0{\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$ (purple) and ${\rho }_0\stackrel{\sim }{v{\prime }_{x}v{\prime }_z}$ vs. $-{\rho }_0{\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$ (green) as functions of Co*. The dashed and solid green lines show the results from HD and MHD simulations, respectively.

Next, we discuss the results in the high-Co_* regime (Figs. 18c and d). Under strong rotational influences, the longitudinal separation of the convective columns decreases and the zonal shear of the vertical flow becomes non-negligible. Consequently, $\stackrel{\sim }{B_x{B}_z}$ is no longer well correlated with $\stackrel{\sim }{\zeta }{\prime }_y$. Instead, $\stackrel{\sim }{B_x{B}_z}$ is now positively correlated with $\stackrel{\sim }{v_x\prime v_z\prime }$ (Figs. 19b and c). This is not surprising because $\stackrel{\sim }{v_x\prime v_z\prime }$ and $\stackrel{\sim }{B_x{B}_z}$ are both rooted in the y-coherent columnar pattern. We note, however, that this positive correlation between $\stackrel{\sim }{v_x\prime v_z\prime }$ and $\stackrel{\sim }{B_x{B}_z}$ is not likely due to the local alignment of small-scale velocity to small-scale magnetic fields. As is clearly manifested by the strong correlation between $\stackrel{\sim }{B_x{B}_z}$ and ${\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$ in Fig. 19e, a main generation mechanism of a positive magnetic field correlation, $\stackrel{\sim }{B_x{B}_z}>0$, remains the shear motions in the columnar convection (both $\partial \stackrel{\sim }{v}{\prime }_x/\partial z$ and $\partial \stackrel{\sim }{v}{\prime }_z/\partial x$). On the other hand, a positive velocity correlation, $\stackrel{\sim }{v_x\prime v_z\prime }>0$, is established as a consequence of the tilted pattern of the y-averaged columns (see Appendix D).

Our simulations clearly show that the Maxwell stress is primarily generated by the shear in the columnar motions, regardless of the rotational regime. This is in contrast to the results of Hotta et al. (2022), who argue that the Maxwell stress originates mainly from the Reynolds stress. One possible reason for this discrepancy lies in the differences in the simulation setups: their model covers a full spherical shell, in contrast to our equatorial local model, and the interaction between SSD and rotating convection at middle to high latitudes (where the columnar convection is absent) may play a dominant role. Moreover, in their simulations, the presence of super-equipartition magnetic fields can significantly suppress the velocity shear, which may inhibit the generation of Maxwell stress through shear.

3.9. Amplitudes of columnar modes

In all our simulations, we find that the velocity strain tensor, ${\stackrel{\sim }{\mathcal{S}}}_{\mathit{xz}}$, is not correlated with $\stackrel{\sim }{v_x\prime v_z\prime }$ but is closely correlated with $\stackrel{\sim }{B_x{B}_z}$ (Fig. 19f). This suggests that the Lorentz force associated with $\stackrel{\sim }{B_x{B}_z}$ plays a significant role as an effective viscosity on the columnar convection, leading to a suppression of the columnar convective modes. Figure 20a shows the volume-integrated kinetic energy of the y-averaged flows,

$$\begin{array}{c}\hfill {\stackrel{\sim }{E}}_{\mathrm{kin}}={\displaystyle {\int }_V}\frac{{\rho }_0}{2}({\stackrel{\sim }{v^{\prime }}}_x^2+{\stackrel{\sim }{v^{\prime }}}_z^2)dV,\end{array}$$(48)

Fig. 20. Volume-integrated (a) kinetic energy, ${\stackrel{\sim }{E}}_{\mathrm{kin}}$, and (b) y enstrophy, ${\stackrel{\sim }{Z}}_y$, of the y-averaged flow as functions of Co*. The definitions of ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ and ${\stackrel{\sim }{Z}}_y$ are given by Eqs. (48) and (49).

as a function of Co_*. There is a general tendency seen in both HD and MHD cases that ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ first increases and then decreases as Co_* increases. The initial increase in ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ is simply due to the establishment of a coherent columnar pattern. The decrease in ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ at Co_* > 5 can be attributed to the suppression of convective velocity in a strongly rotationally constrained regime (Eq. 40). It would also be instructive to use the volume-integrated y enstrophy,

$$\begin{array}{cc}& {\stackrel{\sim }{Z}}_y={\displaystyle {\int }_V}{\stackrel{\sim }{{\zeta }_y^{\prime }}}^{2}dV,\hfill \end{array}$$(49)

as a measure of the amplitudes of the columnar modes. As shown in Fig. 20b, ${\stackrel{\sim }{Z}}_y$ increases monotonically with Co_* in both the HD and MHD cases. Suppression of the columnar convection by the SSD can be seen in both ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ and in ${\stackrel{\sim }{Z}}_y$. The suppression is the most drastic in the intermediate regime Co_* = 5, in which both ${\stackrel{\sim }{E}}_{\mathrm{kin}}$ and ${\stackrel{\sim }{Z}}_y$ are reduced by about 40% from HD to MHD.

4. Summary and conclusions

In this paper, we have carried out a series of HD and MHD simulations of highly turbulent rotating compressible convection using a local f-plane box model at the equator. This model setup enabled us to focus on the effects of SSD on rotating columnar convection. We performed a parameter study in which the rotation rate was varied while the energy flux was fixed, covering the flux Coriolis number, Co_*, from 1 to 20.

We find that the SSD has a significant impact on the properties of convective heat transport. As has been reported in previous nonrotating studies, the convective velocity, v_rms, tends to be suppressed due to the Lorentz force feedback (Rempel 2014; Hotta et al. 2015b). Under strong rotational influences, we find that v_rms decreases with Co_* more rapidly in MHD cases than in their HD counterparts. That is, convection tends to be suppressed much more significantly by rotation when the SSD exists. This rapid decrease in v_rms with Co_* is compensated for by a rapid increase in the entropy fluctuations, s_rms, with increasing Co_*, which is required for rotating convection to transport the fixed amount of enthalpy upward. One of the most striking outcomes of this entropy enhancement by the SSD is the existence of a weakly subadiabatic layer near the base of the convection zone, which persists even in a strongly rotationally constrained regime (Co_* = 20). This is in significant contrast to rotating HD convection, in which the weakly subadiabatic layer vanishes when the rotational influence becomes strong enough (e.g., Käpylä 2024). The potential existence of this subadiabatic layer in rapidly rotating stars has profound implications for dynamo processes in young Sun-like stars (e.g., Brown et al. 2011; Viviani et al. 2018; Käpylä et al. 2019).

These results can be understood by investigating how SSD alters the dominant force balance. Our analysis reveals that, at small scales, the pressure gradient force is balanced by the Lorentz force rather than by inertia when the SSD is present, establishing the so-called magnetostrophic balance. On the other hand, at large scales, a geostrophic balance is always maintained. Consequently, the SSD causes a transformation of the regime of rotating convection from a quasi-geostrophic (QG) CIA balance to a QG MAC balance. In fact, it is shown that the Co_* dependences of v_rms and s_rms in our HD and MHD simulations are consistent with the theoretical scaling relations of the CIA and MAC balances in the rapidly rotating limit (Co_* ≳ 5).

In rotating columnar convection, there are two distinct characteristic length scales. When a typical convective length scale is measured by a mean longitudinal wavenumber, ${\overline{k}}_{x,\mathrm{mean}}$ (weighted by the kinetic energy spectrum), we find that ${\overline{k}}_{x,\mathrm{mean}}$ becomes smaller in MHD cases, i.e., the typical size of convective cells becomes larger, because the kinetic energy is suppressed at small scales due to the SSD Lorentz force feedback. In the rapidly rotating regime (Co_* ≳ 5), ${\overline{k}}_{x,\mathrm{mean}}$ increases with Co_* in both HD and MHD cases, roughly consistent with the CIA and MAC scalings, respectively. On the other hand, the peak wavenumber, ${\overline{k}}_{x,\mathrm{peak}}$ (at which the kinetic energy spectrum has a local maximum), represents the dominant scale of the prograde-propagating columnar convective (thermal Rossby) modes, and it does not follow the CIA and MAC scalings. Our analysis reveals that the SSD enhances the effective rotational influence on columnar convection by increasing the dynamical Coriolis number, Co${}_{\ell }\propto 1/({\overline{v}}_{\mathrm{rms}}{\overline{k}}_{x,\mathrm{mean}})$. As a result, the dominant scale of columnar convection becomes smaller when the SSD is present.

In addition to the convective heat transport, we also find that the SSD has a significant impact on the convective AMT. As has already been reported in previous studies (e.g., Gastine et al. 2013; Guerrero et al. 2013; Mabuchi et al. 2015; Featherstone & Miesch 2015), we see a transition from the “antisolar” mean flow profile (with a faster base and a slower surface) to the “solar-like” mean flow profile (with a faster surface and a slower base) as Co_* increases. We find that the SSD acts to suppress the amplitudes of the mean shear flows, as it provides a Maxwell stress that counteracts the Reynolds stress. This result is consistent with many previous studies (Brun et al. 2004; Nelson et al. 2013; Augustson et al. 2015; Käpylä et al. 2017a; Warnecke et al. 2025). Some spherical dynamo simulations (e.g., Fan & Fang 2014; Karak et al. 2015; Hotta et al. 2022; Soderlund et al. 2025) further reported that the magnetic field can change the rotational profile from antisolar to solar-like near the solar parameter regime (Co_* ≈ 2 − 5). Such a drastic transition was not seen in this study. This can be primarily attributed to the lack of AMT via meridional circulation – an essential ingredient for the generation of differential rotation in spherical models – in our simplified f-plane box model.

In Hotta et al. (2022)’s simulation, a key to achieving solar-like differential rotation was the radially outward AMT by the SSD Maxwell stress. They proposed that this Maxwell stress arises as a back-reaction to the Reynolds stress because the small-scale magnetic fields tend to be parallel to the small-scale flows. In our simulations, we do not find clear evidence of this so-called punching-ball effect. In fact, we find that the y-averaged profile of the Maxwell stress, $\stackrel{\sim }{B_x{B}_z}$, is poorly correlated with that of the Reynolds stress, $\stackrel{\sim }{v_x\prime v_z\prime }$, but is well correlated with those of the y vorticity, $\stackrel{\sim }{\zeta }{\prime }_y$, or the vertical velocity shear, $\partial \stackrel{\sim }{v}{\prime }_x/\partial z$, of the columnar convection. This implies that the Maxwell stress in our simulations is mainly generated by the velocity shear associated with the large-scale convective columns, rather than by the Reynolds stress.

Our study highlights various impacts that SSD has on the properties of rotating columnar convection. These effects need to be taken into account to properly model rotating magneto-convection, particularly in rapidly rotating stars where SSD significantly alters the dominant force balance. Given that fully resolving SSD in global models of solar and stellar convection requires a substantial amount of computational resources, a promising alternative would be to develop subgrid-scale models that encapsulate the essential effects of SSD, as has been attempted in some earlier studies (Bekki et al. 2017; Hotta 2017; Karak et al. 2018).

In relation to the possible implications for the solar convective conundrum, one of the most striking issues is the absence of columnar patterns in the solar surface observations. In this work, we have shown that the SSD acts as an effective viscosity on columnar convection, contributing to the suppression of their mode power. In fact, the kinetic energy and enstrophy associated with columnar modes are reduced by approximately 30–50% due to the SSD effect (see Fig. 20). However, even in the presence of vigorous SSD, our simulations exhibit highly coherent columnar structures near the solar regime (Figs. 4–5). Therefore, the SSD effect alone is likely insufficient to fully account for the absence of observational evidence for columnar convection in the Sun. Some authors (e.g., Guerrero et al. 2013; Noraz et al. 2025) have suggested that the columnar patterns present in the deep convection zone may be obscured by small-scale granular convection near the surface. The density stratification in our simulations (N_ρ ≈ 3.2) is not strong enough to produce such a concealing effect. Future investigations are required to determine how the results presented would change under conditions of stronger density stratification.

Data availability

Movies associated with Figs. 3–7 are available at https://www.aanda.org

Acknowledgments

The author appreciates an anonymous referee for the valuable remarks and suggestions. The author further thanks J. Wicht, T. Bhatia, and H. Hotta for their constructive comments. Y. B. acknowledges a support from ERC Synergy Grant WHOLESUN 810218. All the numerical computations were performed at GWDG and the Max-Planck supercomputer RZG in Garching. This work benefited from discussions within the NORDITA workshop “Stellar Convection: Modeling, Theory, and Observations”.

References

Appendix A: Numerical diffusion scheme

In this study, we do not use any explicit diffusivities but instead employ a slope-limited artificial diffusion scheme (Rempel 2014). This Appendix outlines how this artificial diffusion is implemented and how the effective diffusivities are estimated in our simulations.

The numerical diffusive terms are added in Eqs. (1)–(4) as

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_1}{\partial t}& =[...]-\mathrm{\nabla }·{\mathit{f}}^{\mathrm{diff}}({\rho }_1),\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill {\rho }_0\frac{\partial \mathit{v}}{\partial t}& =[...]-\mathrm{\nabla }·[{\rho }_0{\mathit{f}}^{\mathrm{diff}}(\mathit{v})],\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =[...]-\mathrm{\nabla }·{\mathit{f}}^{\mathrm{diff}}(\mathit{B}),\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill {\rho }_0{T}_0\frac{\partial s_1}{\partial t}& =[...]-\mathrm{\nabla }·[{\rho }_0{T}_0{\mathit{f}}^{\mathrm{diff}}(s_1)]+Q_{\mathrm{visc}}^{\mathrm{diff}}+Q_{\mathrm{joul}}^{\mathrm{diff}},\hfill \end{array}$$(A.4)

where f^diff(q) represents the numerical diffusive flux of a variable q which is non-linearly extrapolated at the grid interfaces (see Sect. 2.1 in Rempel 2014, for its detailed form). For a parameter controlling the strength of diffusion (expressed in Eq.10 in Rempel 2014), we set h = 1.5. In the entropy equation, we add the heating terms arising from the kinetic and magnetic energy dissipation, $Q_{\mathrm{visc}}^{\mathrm{diff}}$ and $Q_{\mathrm{joul}}^{\mathrm{diff}}$. They are estimated as

$$\begin{array}{cc}\hfill Q_{\mathrm{visc}}^{\mathrm{diff}}& =-{\rho }_0{\displaystyle \sum _{i,k}}\frac{\partial v_i}{\partial x_k}{f}_k^{\mathrm{diff}}(v_i),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill Q_{\mathrm{joul}}^{\mathrm{diff}}& =-\frac{1}{4\pi }{\displaystyle \sum _{i,k}}\frac{\partial B_i}{\partial x_k}{f}_k^{\mathrm{diff}}(B_i).\hfill \end{array}$$(A.6)

To infer the effective viscous and magnetic diffusivities, ν_eff and η_eff in our simulations, we compared the above numerical energy dissipation rates with those inferred from explicit diffusivities as (Hotta 2017; Rempel 2018)

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{eff}}& =\frac{⟨Q_{\mathrm{visc}}^{\mathrm{diff}}⟩}{{\rho }_0\langle {\displaystyle {\sum }_{i,k}}\frac{\partial v_i}{\partial x_k}[\frac{\partial v_i}{\partial x_k}+\frac{\partial v_k}{\partial x_i}-\frac{2}{3}{\delta }_{\mathit{ik}}(\mathrm{\nabla }·\mathit{v})]\rangle },\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{eff}}& =\frac{⟨Q_{\mathrm{joul}}^{\mathrm{diff}}⟩}{{\langle |\mathrm{\nabla }\times \mathit{B}|}^2\rangle /4\pi }.\hfill \end{array}$$(A.8)

Using the volume-averaged values of these diffusivities, we can estimate the effective Reynolds and magnetic Reynolds numbers, Re_* and Rm_*, as

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\ast }=\frac{v_{\ast }{H}_r}{{\overline{\nu }}_{\mathrm{eff}}},\mathrm{and}{\mathrm{Rm}}_{\ast }=\frac{v_{\ast }{H}_r}{{\overline{\eta }}_{\mathrm{eff}}}.\end{array}$$(A.9)

The computed values are reported in Table 1.