Convective Velocity Suppression via the Enhancement of Subadiabatic Layer:
Role of the Effective Prandtl Number
Abstract
It has recently been recognized that the convective velocities achieved in the current solar convection simulations might be overestimated. The newlyrevealed effects of the prevailing smallscale magnetic field within the convection zone may offer possible solutions to this problem. The smallscale magnetic fields can reduce the convective amplitude of smallscale motions through the Lorentzforce feedback, which concurrently inhibits the turbulent mixing of entropy between upflows and downflows. As a result, the effective Prandtl number may exceed unity inside the solar convection zone. In this paper, we propose and numerically confirm a possible suppression mechanism of convective velocity in the effectively highPrandtl number regime. If the effective horizontal thermal diffusivity decreases (the Prandtl number accordingly increases), the subadiabatic layer which is formed near the base of the convection zone by continuous depositions of low entropy transported by adiabatically downflowing plumes is enhanced and extended. The global convective amplitude in the highPrandtl thermal convection is thus reduced especially in the lower part of the convection zone via the change in the mean entropy profile which becomes more subadiabatic near the base and less superadiabatic in the bulk.
=1
1 Introduction
Thermal convection plays a critical role inside the solar interior; it transports heat while powering dynamos. For investigating the convective properties and the dynamo processes occurring inside the convection zone, three dimensional (3D) fullspherical hydrodynamic (HD) or magnetohydrodynamic (MHD) simulations have been conducted (e.g., Glatzmaier 1984; Miesch et al. 2000; Brun et al. 2004). However, it has been recently recognized that the recent 3D solar convection simulations which aim to achieve more realistic regimes with higher resolutions tend to overestimate the convective amplitudes, which is commonly known as the convective conundrum (see Hanasoge et al. 2016, for details).
Recent local helioseismic observations by Hanasoge et al. (2012) inferred the largescale convective amplitude which is about two orders of magnitude smaller than the value typically predicted by mixinglength model or obtained in global HD simulations (Miesch et al. 2008). It should be noted that the helioseismic estimates depend on the inversion technique employed, for Greer et al. (2015) instead found much larger flow velocities that are consistent with the theoretical models. Indeed, realistic surface convection simulations also reported the largescale velocity power an order of magnitude larger than photospheric observations, also suggesting that the numerical simulations are overpowering the convection larger than supergranular scales (Lord et al. 2014).
Another evidence comes from the problem related with the differential rotation profile. Recently, many highresolution global simulations with solar parameters, such as the solar luminosity and the rotation rate , reported that the antisolar differential rotations with faster rotating pole and slowly rotating equator were unexpectedly achieved (Gastine et al. 2013; Käpylä et al. 2014; Fan & Fang 2014; Hotta et al. 2015b; Karak et al. 2015). It is generally believed that this antisolar differential rotation profile is attributed to a large Rossby number , which measures the effect of convection with respect to the rotation, where , and are the rotation rate, typical convective velocity and length scale, respectively. If the simulated convective velocities are overestimated, the Rossby number becomes large and then the angular momentum is transported radially inward by the turbulent Reynolds stress, which finally results in producing the antisolar differential rotation profile (Featherstone & Miesch 2015). We believe therefore that if the problem of too fast convection obtained in the recent highresolution simulations is solved, it will highly contribute to alleviating the problem of antisolar differential rotation simultaneously.
The problem described so far suggests the existence of some physical processes that the recent global (M)HD simulations neglected or just could not capture well enough. A possible solution to the problem is smallscale magnetic fields generated by smallscale dynamos (e.g., Brandenburg & Subramanian 2005). Smallscale dynamos can operate when the magnetic Reynolds number exceeds the critical value and can generate magnetic fields on scales smaller than the energycarrying scale of the turbulent flow. The effects of smallscale magnetic fields are not included in the HD simulations, and even for MHD case, it is highly likely that most global simulations have insufficient resolutions to achieve large enough to resolve the inertial range of the turbulence, which is essential for the smallscale dynamo. The importance of the smallscale dynamo near the surface has been widely recognized in explaining the observed mixedpolarity field in the photosphere (Cattaneo 1999; Vögler & Schüssler 2007; Rempel 2014) where the flows are essentially nonhelical due to a sufficiently large . On the other hand, in the deeper convection zone where largescale dynamos can be easily excited by helical turbulence, little attention has been paid so far to the effects of the smallscale dynamo.
Recently, however, Hotta et al. (2015a, 2016) conducted highresolution MHD simulations of the solar convection zone with and without rotation and finally demonstrated that the smallscale dynamos can also be efficiently excited throughout the convection zone even when the rotational effects are included. It is shown that the magnetic energy becomes close to equipartition or even superequipartition with the convective kinetic energy at small scales. The Maxwell stress related to these superequipartition magnetic fields is found to act similarly to a viscous stress for convective motions. As a result, the convective velocities are suppressed via the Lorentz force feedback by in the highestresolution MHD simulation of Hotta et al. (2015a). Although, the smallscale dynamo is still not saturated in the sense that it does not converge in resolution, the suppression of convective velocities directly through the Lorentzforce seems not enough for explaining the huge discrepancy between observations and simulations.
Another effect of the smallscale magnetic field which prevails within the convection zone was pointed out and first investigated by O’Mara et al. (2016). The smallscale magnetic field not only inhibits shearing motions but also reduces the horizontal turbulent mixing of the entropy perturbations between warm upflows and cold downflows (Hotta et al. 2015a). In other words, the effective viscosity is enhanced while the effective thermal diffusivity is reduced, leading to an increase in the effective Prandtl number which can be larger than unity. O’Mara et al. (2016) conducted a set of solar convection simulations in which is decreased while is kept fixed (and thus is increased). It was shown that the convective velocity systematically decreases as is increased. They attribute this decrease in convective velocities to the large thermal contents of downflowing plumes originating from a stronglysuperadiabatic surface layer, which is formed by their vertical thermal conductivetype upper boundary condition. In short, convection becomes weak in high regime so that the fixed solar luminosity can be transported outward with larger thermal contents possessed by warmer upflows and colder downflows.
Let us distinguish these two effects of the smallscale magnetic field by referring to the former one (suppression by the Lorentz force) as a dynamical effect and the latter (suppression via the change of ) as a thermal effect. Through the dynamical effect of the Lorentz force, the kinetic energy of the convective flows can be reduced via the energy exchanges between the dynamogenerated magnetic energy. The importance of the thermal effect via the decrease in is significant because it describes that the reduction in the kinetic energy can also result from an increase in the internal energy. Since the internal energy is far larger than both the kinetic and the magnetic energy in the deep convection zone owing to a very small Mach number () and a very large plasma beta () (Ossendrijver 2003), the thermal effect will offer another possibility to resolve the huge discrepancies between simulations and observations.
In this paper, we propose a velocity suppression mechanism that can be achieved in high regime which is distinct from what was described in the previous study of O’Mara et al. (2016). We further confirm the suppression mechanism by conducting a set of HD simulations of Cartesian box where the effects of smallscale magnetic fields are modeled as subgridscale (SGS) diffusivities, and . Previous smallscale dynamo simulations revealed that the horizontal rootmeansquare (rms) velocity is more decreased than the vertical one, which leads to a main reduction of a “horizontal” heat transport by the turbulence (Hotta et al. 2015a). We therefore introduce anisotropy in the SGS thermal diffusion and examine a dependence of the convective amplitudes on each of the horizontal and vertical SGS thermal diffusivities.
The organization of the paper is as follows. In Section 2, we describe a possible convective velocity suppression mechanism. In Section 3, our numerical model is explained. The simulation results for both isotropic and anisotropic thermal diffusion cases are presented in Section 4. We discuss implications on the solar convection zone dynamics in Section 5, and the conclusions are summarized at last in Section 6.
2 Effects of the Decrease in the Effective Thermal Diffusivity
In this section, we examine possible effects of a decrease in the effective thermal diffusivity which can be brought about by smallscale magnetism. We argue here that the convective velocity can be suppressed via the enhancement of a weakly subadiabatic layer near the base of the convection zone.
Figure 1 shows a schematic illustration of this process which is explained as follows. If thermal diffusivity is reduced, cold downflows become able to fall down without losing their large thermal contents via the thermal diffusion so that the low entropy fluids are transported more adiabatically to the bottom. Note that continuous depositions of the low entropy near the bottom may naturally form a weaklysubadiabatic (convectivelystable) layer. The formation of this subadiabatic layer near the base would be therefore enhanced and extended in low (and thus high) regime. The influence of the mean entropy stratification on convection is evaluated by superadiabaticity with . The final value of subadiabaticity () near the base achieved in a statisticallystationary convection would be determined by a balance between a continuous supply of low entropy by downflows and suppressions of downflows which in turn limits the amount of low entropy supply. Note that if is relatively large, a relaxation effect by the vertical thermal conduction also enters this balance. The enhanced and extended subadiabatic layer in the lower convection zone may have a considerable impact on suppressing global convective amplitudes: it not only suppresses the vertical motions locally through the buoyant deceleration but it can also limit the global convective mode to shallower scales (Käpylä et al. 2017b) which then leads to a reduction of the global convective velocity.
In order to distinguish several thermal effects operating in the low regime, it is instructive to define the horizontal and vertical thermal diffusivities, and and discuss the thermal diffusion in each of the directions. We consider that the decrease in is essential for downflows to retain their large amount of entropy deficits and to achieve a mean entropy stratification which is more subadiabatic near the base. However, also has several important effects that can influence the proposed velocity suppression mechanism. Near the top and bottom boundaries, a decrease in prohibits the vertical relaxation of the mean entropy stratification, leading to a more subaidiabatic base and more superadiabatic surface. Moreover, near the surface layer where convection is continuously driven, the decrease in has another significant effect in increasing the Rayleigh number () and thus increasing the degree of turbulence there.
In the previous numerical studies of compressible thermal convection with larger than unity, the thermal diffusion has been treated isotropically (Warnecke et al. 2014; O’Mara et al. 2016; Käpylä et al. 2017a), and thus, all of the thermal effects described above are simultaneously included. However, recent MHD simulations revealed that the smallscale magnetism tends to make the smallscale motions highly anisotropic, suggesting that the SGS turbulent transport coefficients should also be anisotropic: More precisely, a main decrease in is inferred (Hotta et al. 2015a). In order to distinguish the different thermal effects and to investigate the influence of on the proposed velocity suppression mechanism, the anisotropy in the thermal diffusion is necessary. We thus conduct a set of high convection simulations to numerically confirm that the enhanced subadiabatic layer can suppress the global convective amplitude with anisotropic thermal diffusion included.
3 Numerical Model
3.1 Basic Equations
In this investigation, we solve the threedimensional hydrodynamic equations in a simplified system with Cartesian coordinate . In our definition, axis is directed vertically upward, antiparallel to the gravity. The basic equations consist of the equation of continuity, the equation of motion, the equation of entropy, and the equation of state,
(1)  
(2)  
(3)  
(4) 
Here, , , and denote timeindependent reference state values of density, pressure, and temperature, respectively. The reference state is assumed to be in an adiabaticallystratified hydrostatic equilibrium. The thermodynamic variables with subscript 1, , , and represent perturbations from the reference state that are small compared with background values so that the equation of continuity and the equation of state are linearized. Note that the entropy is normalized by the specific heat capacity at constant volume and the ideal gas is assumed for the ratio of specific heats, . denotes the gravitational acceleration and is assumed to be constant in space.
The reference state quantities are given in the following forms that are identical to Fan et al. (2003),
(5)  
(6)  
(7)  
(8) 
where , , , and denote the values of , , , and the pressure scale height evaluated at the bottom , respectively. Polytropic index takes an adiabatic value .
3.2 SubGridScale Diffusivities
denotes the viscous stress tensor and is the amount of dissipated energy which is converted from kinetic energy into internal energy. They are given by,
(9)  
(10) 
where and are the coefficient of kinetic viscosity and the Kronecker delta, respectively. Note that, for the sake of simplicity, we try to model the Maxwell stress on scales smaller than the grid resolution (subgridscale) by the viscous stress on the resolved scale . Here, represents an ensemble average on subgridscale so that , . Therefore, the coefficient should be regarded as a diffusive coefficient of an effective viscosity which mimics the smallscale magnetic tension force.
On the other hand, the thermal diffusivity is regarded as a subgridscale eddy diffusivity reflecting the heat transport by unresolved turbulent motions that are subject to a strong Lorentzforce feedback of superequipartition smallscale magnetic fields. Smallscale dynamo simulations of Hotta et al. (2015a) showed that the smallscale magnetism has a significant effect in suppressing the “horizontal” rms velocity, and as a result, the horizontal turbulent heat transport is greatly reduced. This means that the effective thermal diffusivity should be suppressed in a horizontal direction when we take into account the unresolved turbulent motions that are magnetized. We, therefore, introduce the anisotropic subgridscale thermal diffusion by expressing the thermal diffusivity tensor as,
(11) 
Here, and are the horizontal and vertical thermal diffusivities, respectively. Both and are assumed to have the same height dependence,
(12) 
The same functional form has been commonly employed in the Anelastic Spherical Harmonics simulations (Miesch et al. 2000; Brun et al. 2004). Although the subgridscale physics in our model is different from that of these previous studies, height dependence of the diffusivities is fixed just for simplicity.
The subgrid scale that we are considering in this study represents a scale smaller than that at which the magnetic energy excited by smallscale dynamos exceeds the kinetic energy of the turbulent flows (let us call this scale “SSD scale” hereafter). We compute the convective motions on scales larger than SSD scale, assuming that the smallscale magnetic fields have little direct influence on these scales. Typically, SSD scale is expected to lie between the inertial range of the turbulent flows (Hotta et al. 2015a, 2016) so that the moderate Reynolds number is enough for this study. It should be emphasized that we do not aim at any realistic solar magnetoconvection simulations with this model but rather study the physical processes of convective velocity suppression under a simplifying assumption that the effects of smallscale magnetism can be modeled as subgridscale diffusivities.
3.3 Radiative Energy Fluxes
In our model, the convectively unstable stratification is formed by injecting the artificial radiative energy flux from the bottom and extracting the same amount of energy flux through the upper boundary. For specifying the radiative heating term, the functional form presented in Featherstone & Hindman (2016) is adopted which can mimic the normalized radiative energy flux tabulated in a standard solar model. The radiative heating is assumed to be proportional to the background pressure,
(13) 
where the normalization factor is specified by the input energy flux as,
(14) 
The radiative energy flux is then defined as,
(15) 
The radiative cooling term at the surface is given in a similar form to Hotta et al. (2014) as,
(16)  
(17) 
Therefore, the cooling effect is localized near the top boundary and the thickness of the upper thermal boundary layer is mostly set by the local pressure scale height there.
3.4 Numerical Method
Case  

H1V1………  :  2.49  
H2V2………  :  3.48  
H6V6………  :  5.83  
H1V6………  :  4.21  
H6V1………  :  3.55 
Note. – Horizontal and vertical Prandtl number and are the free parameters in our model. is the volume and timeaveraged (nearly ) effective Reynolds number defined by equation (23). measures the entropy difference across the numerical domain averaged over nearly simulation time. denotes the minimum superadiabaticity that is usually achieved near the base. represents the depth where the superadiabaticity changes from negative to positive value. here denotes the volumeaveraged rms velocity in each case. is the kinetic energy in a stationary state integrated over the entire volume and averaged over nearly simulation time.
We solve the equations (1)(3) numerically with periodic boundary conditions for horizontal directions (at and ) and impenetrable, stressfree boundary conditions for lower and upper boundaries (). The free boundary condition is applied for density and entropy, i.e., , at the top and the bottom. The size of the numerical domain is set by and . The vertical domain spans density scale heights with density contrasts.
The magnitude of subgridscale diffusivities are measured by the Reynolds number at the base and the Prandtl number which, in our simulations, is defined distinctly for both horizontally and vertically . Here, the velocity scale is estimated based on the input energy flux as . The energy flux is determined by specifying the modified Mach number which is defined as . The values of the above dimensionless numbers for our simulations are given by , . With this numerical setup, the subsequent amplitudes of the dimensionless entropy or the superadiabatic gradient would be on the order of . The assumption of linearized equation of state is verified considering that . Prandtl numbers are treated as free parameters in this study and given as listed in Table 1. We should emphasize here that the kinetic viscosity is kept fixed in all simulations and we only change and . In fact, the SGS should also increase due to the Lorentzforce of the smallscale magnetic fields. However, the dynamical effect originating from the change in is not considered and only the thermal effect is investigated in this study, which may underestimate the degree of velocity suppression.
We use the fourthorder centereddifferencing method for space and the fourthorder RungeKutta scheme for the timeintegration (Vögler et al. 2005). The same artificial viscosity used in Rempel (2014) are implemented and added to velocity fields to stabilize numerical computations. The code is parallelized using message passing interface (MPI). For all our calculations, an uniform resolution of is used. Each simulation is initiated by giving a small random perturbation for vertical velocity with the other variables , , , and set to zero.
4 Results
4.1 Overview
First, let us review the time evolution of the system. In the following, we use the linear approximation for thermodynamic variables and also omit terms proportional to that are small for low Mach number flows. From equations (2) and (3), the energy equation is derived as,
(18) 
Under the boundary conditions described earlier, a global energy conservation law can be expressed as,
(19) 
Figure 2 shows the temporal evolution of total kinetic energy , thermal energy , and gravitational potential energy in Case H1V1 for a reference. The volume integrated energies are defined as,
(20)  
(21)  
(22) 
Note that, in our definition, the thermal energy represents a perturbation with respect to the energy contained by the adiabatic background. The internal energy of the system is expressed as . As prescribed by the equation (19), a conservation of can be well confirmed from Figure 2.
Since the surface cooling can quickly produce highly superadiabatic stratification near the top boundary, the thermal convection typically sets in at an early stage and a statistically stationary state is achieved after as shown in Figure 2. We then further evolve the system for at least one thermal diffusion time after the statistically stationary state is reached. From now on, we are going to discuss the results of this stage. The statistical properties of the convection are investigated by temporallyaveraging the data with a cadence of .
4.2 Isotropic Thermal Diffusion
In this section, the results of Cases H1V1, H2V2, and H6V6 are discussed where thermal diffusion is isotropic (). We try to make our argument that the change in the mean entropy stratification is mostly responsible for the convective velocity suppression in high regime by relating the suppressed convective velocities with the enhanced subadiabatic layers. Figure 3 clearly shows that the rms velocity systematically declines as is decreased ( is increased), as found by O’Mara et al. (2016). The dashed lines in Figure 3 represent vertical rms velocities that are directly subject to buoyancy accelerations. It is noteworthy that the vertical rms velocities are decreased in the bulk of the domain except for the top layer and correspondingly that the total rms velocities are suppressed mainly in the deeper layer. The effective Reynolds number is estimated in each case using the dependent rms velocities as,
(23) 
and presented in the fourth column of Table 1. In our model, the convection typically shows a moderate degree of turbulence with the typical value of which is of the same order as used in previous studies (O’Mara et al. 2016; Käpylä et al. 2017b).
Figure 4 shows the mean entropy normalized by for each case. Here, denotes the horizontal averaging. In this paper, we normalize the entropy (and superadiabaticity ) by to make a comparison with results with different input energy fluxes easier. It is obvious from Figure 4 that the highlysuperadiabatic thermal boundary layer is formed near the top boundary, and the bulk of the domain is more close to adiabatic compared with the top layer. Although the same form of surface cooling flux is imposed in all cases, the thickness of the upper thermal boundary layer varies according to due to the vertical thermal conduction near the surface. As a result, the entropy difference between the top boundary and the bulk of the domain increases as increases. However, dependence of the upper boundary layer is small compared with the previous studies (O’Mara et al. 2016), having no scaling relations between the thermal diffusivity , the thickness of the thermal boundary layer , and the entropy difference , such as or which can only hold for simulations imposing a diffusiontype upper boundary condition where the energy flux is released through the thermal conduction term (Featherstone & Hindman 2016).
The zoomin inset of Figure 4 manifests firstly that the mean entropy tends to be slightly subadiabatic near the bottom boundary and secondly that the subadiabaticity there increases as increases. These results can be interpreted that the subadiabatic layer is formed by continuous depositions of low entropy transported by downflows. When the thermal diffusivity is reduced, the amount of low entropy fluids that can be retained by downflows during the descent becomes large, which leads to an enhancement of the subadiabaticity near the base. The enhanced subadiabatic layer then suppresses downflow motions by alleviating the thermal contents of downflows and finally sets the net accumulation rate of the low entropy to the bottom. In a statisticallystationary state, the continuous accumulation of low entropy is compensated by the vertical thermal conduction near the bottom boundary.
Some caution should be taken here regarding the impenetrable lower boundary condition in our idealized model. It is true that this lower boundary condition indeed favorably affects the formation of this weaklysubadiabatic layer, helping the reflection and accumulation of low entropy around the base. However, the existence of the slightlysubadiabatic convection zone has been repeatedly reported even in the numerical simulations of convective overshoots where a stablystratified layer is included (Brummell et al. 2002; Käpylä et al. 2017b) and also predicted by nonlocal semianalytical model of solar overshoot layer (Xiong & Deng 2001; Rempel 2004). Moreover, recent numerical simulations of solar overshoot region, which aim to achieve more realistic parameter regimes by imposing a much lower energy flux, estimated the depth of solar overshoot layer to be less than a few percent of the local pressure scale height there (Hotta 2017), suggesting that the bottom of the solar convection zone may act as an impenetrable wall to a good approximation.
Figure 5 shows horizontallyaveraged energy fluxes across the numerical domain for Case H1V1 and H6V6. The definitions of the enthalpy flux , the kinetic energy flux , the thermal conductive flux , and the viscous dissipative flux are,
(24)  
(25)  
(26)  
(27) 
Due to the velocity suppression in Case H6V6, both the amplitudes of enthalpy flux and kinetic energy flux decrease. The reduction in the vertical thermal conductive flux near the surface in Case H6V6 is compensated by the verticallyupward peak shift of the enthalpy flux.
It should be noted that the enthalpy flux takes positive value throughout the numerical domain so that the thermal energy is transported vertically upward even near the base where the mean stratification is subadiabatic. Therefore, this subadiabatic layer formed near the base is not an overshooting layer where downflows are quickly decelerated by buoyancy and thus the enthalpy flux becomes negative. Rather it should be regarded as a result of nonlocal heat transport of downflow plumes, which cannot be described by the typical local mixing models assuming the enthalpy flux proportional to the local superadiabaticity. Recently, Brandenburg (2016) modified the mixinglength theory incorporating the effects of nonlocal heat transport by cold downflow plumes and showed that the enthalpy flux can be upward even in the subadiabatic region. We will discuss on this issue later in Section 5.1. Even though the subadiabaticity near the base is not strong enough to reverse the sign of enthalpy flux to be a overshooting layer, this layer has significant effects on the convective velocity amplitudes.
4.3 Anisotropic Thermal Diffusion
Next, we examine results of Case H1V6 and Case H6V1 where vertical and horizontal thermal diffusivities and are decreased independently so that several thermal effects can be separated. In Case H1V6, we drop only to the same value used in Case H6V6 while keeping identical to Case H1V1. In Case H6V1, on the other hand, only is dropped to the same level of Case H6V6 with unchanged from Case H1V1.
First, the overall convection patterns are explained. Figure 6 and Figure 7 show vertical velocities and entropy perturbations for Cases H1V1, H1V6, H6V1, and H6V6 from left to right panels with upper and lower panels presenting the horizontal cuts near the top and at the middle, respectively. In Figure 6, dependency of the vertical convection is more apparent in the horizontal cuts at the middle layer (panels (e)(h)): We can see that downflows become more pointwise when is decreased. In other words, convective structure is qualitatively changed from lanetype downflows to plumetype downflows as decreases. Moreover, it is also obvious that the overall amplitudes of upflows and downflows are reduced in Cases H6V1 and H6V6 where is reduced, which will be confirmed later in Figure 8.
As for entropy, it is clear from Figure 7(a)(d) that the amplitude of the entropy perturbation at the surface becomes large for Cases H1V6 and H6V6. This is because the suppression in makes the upper thermal boundary layer steeper, leading to a thinner highlysuperadiabatic surface layer and an enlargement of the absolute value of entropy difference. The change in , on the other hand, only plays a role in forming thinner downflow lanes near the surface. As we go into deeper convection zone, however, the entropy distributions are dominantly characterized by , as shown in Figure 7(e)(h): Downflow plumes are able to retain their low entropy at small scales in Cases H6V1 and H6V6 where horizontal diffusion is highly prohibited. In Cases H1V1 and H1V6, in contrast, efficient horizontal heat exchanges blur the entropy at deeper convection zone.
Figure 8 shows distributions of total, vertical, and horizontal rms velocities for Cases H1V1, H1V6, H6V1, and H6V6. We can confirm from Figure 8(a) that the significant reduction of rms velocity only occurs in Case H6V1 and H6V6 in which is decreased (see also the volumeaveraged rms velocity values tabulated in the eighth column of Table 1). The volumeintegrated kinetic energies are calculated and presented in the ninth column of Table 1. is highly reduced up to for Case H6V6 and for Case H6V1 with respect to that of Case H1V1, whereas in Case H1V6 the kinetic energy is suppressed only slightly (). Especially, it is obvious from Figure 8(b) that the values of vertical rms velocities are mostly determined by , except for the top surface where influences the convective properties by increasing the Rayleigh number . Therefore, it is concluded here that the prohibition of horizontal, not vertical, diffusive transport of entropy between warm upflows and cold downflows is responsible for the velocity suppression in high regime. We attribute this convective velocity suppression to the change of mean entropy profile as follows.
Figure 9 shows profiles of the mean entropy and corresponding superadiabaticity calculated by,
(28)  
for Case H1V6 and H6V1 at upper and lower panels, respectively. The amount of entropy difference between the top and the bulk is substantially influenced by as shown in Figure 9(a) and (c) (see also the fifth column of Table 1). For example, in Case H1V1 (black) and Case H6V1 (green) where and is relatively large, strong vertical thermal diffusion near the top efficiently alleviates the highlysuperadiabatic entropy gradient, resulting in the small entropy difference .
Figure 9(b) shows that the superadiabaticity profile is almost unaffected by a decrease in in the bulk of the domain. On the other hand, we can clearly observe from Figure 9(d) that the subadiabatic layer is enhanced and extended vertically upward in the lower portion of the domain and also that the mean stratification in the bulk becomes less superadiabatic from Case H1V1 (black) to Cases H6V1 (green) and H6V6 (red) as decreases. Thus, is only sensitive to and this is why huge reductions of convective velocity occur in Cases H6V1 and H6V6. Both the enhancement of subadiabaticity and weakening of superadiabaticity lead to the reduction of the net buoyancy acceleration and suppression of the convective amplitudes.
The different sensitivity of on and on comes from the fact that downflows can retain their cold entropy only when the horizontal thermal diffusion is inhibited, as clearly shown in Figure 10 where entropy fluctuations of upflows and downflows are plotted. Firstly, it is obvious that only the entropy contents of downflows are modified and the thermal properties of upflows are almost unchanged by . In Case H1V6 (blue), thermal fluctuation of downflows is mostly identical to that of Case H1V1 (black) in the bulk, which means that downflows tend to quickly lose their cold entropy via the horizontal thermal diffusion even though strong entropy deficits are generated at the top. In Case H6V1 (green), on the other hand, downflows can retain their entropy deficits against horizontal thermal diffusion so that they can greatly contribute to the overall accumulations of low entropy around the base. It is considered that the amount of negative entropy perturbation that can be transported by downflows and accumulated near the base should largely set the subadiabaticity there.
In fact, these results should be interpreted with some caution because the thermal diffusivities have a vertical dependence expressed by the equation (12). In principle, the profile of in a stationary state should be determined by both and ; the former one adjusts the amount of low entropy that is supplied to the base and the latter one sets the final entropy profile after the vertical thermal relaxation. Although the dependence of near the base may be emphasized in our model by originally prohibiting the vertical thermal conduction in the lower part in all cases, our conclusion that the decrease in is crucial for the enhancement of subadiabatic layer and for the convective velocity suppression is indeed supported by results of another set of simulations where spatiallyuniform is used (not shown here): The suppression in occurs only when are decreased. When the amplitude of the entropy perturbation is increased, low entropy material is transported more to the base. This process leads to an extension of the subadiabatic layer to the upper part of the convection zone and thus makes the bulk stratification more subadiabatic and less superadiabatic. On the other hand, effects of the vertical thermal diffusion are restricted to the lower and upper thermal boundary layers and does not largely affect the mean entropy stratification in the bulk which is more close to adiabatic. As a result, is almost unaffected by .
In order to quantitatively examine the effects of the enhanced subadiabatic layer, work density done by buoyancy force is calculated as,
(29) 
and presented in Figure 11. The buoyancy work becomes positive in the upper convection zone and negative in the lower convection zone. This general tendency is consistent with the fact that superadiabatic (subadiabatic) stratification accelerates (decelerates) thermal convection. The existence of the region near the bottom where the buoyancy work is negative can be thus regarded as a clear evidence of convection suppression by the subadiabatic layer. In Case H6V1 (green) and Case H6V6 (red), the positive buoyancy work is reduced due to the less superadiabatic stratification in the upper convection zone except for the surface region. Near the bottom, on the other hand, the absolute values of negative buoyancy work is decreased in these two cases where subadiabaticity is enhanced, which at first glance seems contradictory to our argument. We consider that the buoyancy force can do less negative work (to decelerate convection) just because the vertical velocity amplitudes are reduced in these two cases. Note here that the buoyancy works calculated by the above equation just represent the energy conversion rate between thermal and vertical kinetic energies in a statistically stationary state. A large amount of positive buoyancy work does not necessarily result in a large amount of kinetic energy or large convective amplitudes.
Vertical dotted lines shown in Figure 11 denote the heights from which mean stratifications change from subadiabatic to superadiabatic (see the seventh column of Table 1). Although the subadiabatic stratification eventually results in the negative buoyancy work in the lower convection zone, the deceleration actually sets in below the critical height so that there is a gap zone inbetween where the mean stratification is subadiabatic but buoyancy work is still positive . The origin of this nonlocalness comes from the fact that the thermal fluctuations possessed by downflows (larger density than surroundings for instance) must be gradually modified as downflows travel inside the subadiabatic zone and become less heavy (Hotta 2017). Figure 11 shows that the downflows are quickly decelerated after penetrating into the subadiabatic zone for Cases H1V1 (black) and H1V6 (blue) where the thermal contents of downflows are relatively small as shown in Figure 10. However, the nonlocalness become more substantial as downflows contain colder entropy from Case H1V1 (black) to Case H6V1 (green) to Case H6V6 (red). This is because the amount of thermal fluctuations necessary for reversing the sign of buoyancy work gets larger in Case H6V1 and H6V6, so that it becomes more difficult and needs longer distance to achieve for each downflow. Therefore, the deceleration region is localized near the bottom in Cases H6V1 and H6V6, and as a result, it helps for downflows to reach the base and to reinforce the subadiabaticity by accumulating the low entropy there. It should be emphasized again that since the subadiabaticity is on the order of the local Mach number square (see also sixth column of Table 1) and not strong enough to be an overshoot layer, downflows only feel gradual decelerations and thus the buoyancy works are not simultaneously responded.
In summary, we have shown up to here by comparing the results of Cases H1V1, H1V6, H6V1, and H6V6 that the suppression of is essential in establishing the enhanced and extended subadiabatic layer near the base and in decreasing the whole convective amplitudes. All of these results presented in this section are in qualitative accordance with the velocity suppression mechanism discussed in Section 2.
5 Discussion
5.1 Implications for Solar Convection
In the preceding sections, we have explained how the convective flow speed could be suppressed when the horizontal thermal diffusivity is decreased. Since the main focus of our research is to investigate the physical mechanism of velocity suppression and we do not aim at any realistic solar convection modeling, some caution is needed when discussing the applicability of our results to the solar convective conundrum. Huge discrepancy of convective amplitudes (by more than two orders of magnitude) between the helioseismic finding of Hanasoge et al. (2012) and the numerical simulation of Miesch et al. (2008) was found mainly on the low wavenumber subsurface flow (with the spherical harmonic degree ), which in other words suggests that the giant cells may not exist or their amplitudes may be far smaller than previously considered. Considering that our numerical simulations employ a simplified Cartesian box with periodic boundaries in horizontal directions, discussions on the globalscale () flow amplitudes are beyond the scope of this paper.
Nonetheless, our results have profound implications on the solar convective structure. In general, it is considered that the excess in the low wavenumber power may reflect too large convective amplitudes in a deeper convection zone (Lord et al. 2014). A main advantageous feature of our proposed thermal effect is that convective amplitudes in a deeper layer can be selectively suppressed along the enhancement of the subadiabaticity there. This may potentially alleviate the convective conundrum by suppressing the deep convection and decreasing the low wavenumber power of the subsurface horizontal flow.
Then, how much convective velocities can be reduced at most via the thermal effect? Can the thermal effect be efficient enough to fully resolve the discrepancies between observations and simulations? To derive some hints for answering these questions, another set of numerical simulations is conducted where only the horizontal Prandtl number is systematically increased from up to with the vertical Prandtl number fixed (). The numerical domain is horizontally restricted to and an uniform resolution of is used except for Case where a resolution is increased up to in order to keep the influence of the numerical diffusivity sufficiently small.
Figure 12 shows the results of rms velocities obtained for different values of . Localmaximum (localminimum) values of the rms velocities near the surface (bottom) are plotted on a doublelogarithmic scale. The thermal effect does not saturate within the parameter regime investigated where Prandtl numbers are kept still moderate () so that the rms velocities are found to monotonically decrease as increases. Figure 12 confirms that the rms velocities are more suppressed near the base than in the surface region. The scaling relations of maximum rms velocity near the surface and minimum rms velocity near the bottom are calculated as,
(30) 
Although this result shows a promising feature, it infers that the dependence of our proposed thermal effect is relatively weak and that a significantly large effective Prandtl number on the order of is required to lower the velocity amplitude by .
Several careful considerations are needed especially when trying to apply this result to the Sun or to the other numerical systems. First of all, we remind the reader that the SGS viscous diffusivity (and thus the Reynolds number ) is fixed in all calculations for simplicity. We must note that the scaling relation above is derived for the parameter regime employed in our numerical setup and that the scaling index may vary according to . For example, since we impose large viscous and thermal SGS diffusivities with the typical Reynolds number calculated as , the parameter regime studied in our simulations is laminar. It is expected that, if SGS is decreased and the Reynolds number (and thus Rayleigh number ) increases, the thermal effect would become ineffective and would cease to diminish being independent of the values of diffusivities (Featherstone & Hindman 2016).
Secondly, determining whether or not the thermal effect saturates for higher , or if do, determining the saturated value of or the critical would be another important issue that needs to be investigated when discussing the applicability of the scaling relation (equation (30)). Note that in the previous numerical study of O’Mara et al. (2016) the saturation and limitation of velocity suppression in high regime was mainly discussed in relation to the thickness of the upper thermal boundary layer which scales as owing to their conductivetype upper boundary condition: O’Mara et al. (2016) argued that their suppression mechanism would become unphysical when is decreased up to the point where approaches to the actual depth of the photospheric boundary of the Sun. In our model, on the other hand, the depth of the upper boundary layer is set by the artificial surface cooling function and is almost independent of . Therefore, the thermal effect discussed in this paper should not saturate in the same way as the model of O’Mara et al. (2016). Instead, it is expected that in our model the saturation occurs when the thermal diffusion becomes negligible and accumulations of low entropy by downflows becomes unable to further change the mean entropy stratification: If becomes small enough to make the vertical thermal conduction essentially ineffective, the subadiabaticity near the base is determined so that the strong buoyant decelerations can stop downflows before they reach the bottom boundary and limit their supply of low entropy to the base.
From another side, the mean superadiabatic stratification of the solar interior has been estimated based on the mixinglength theory. Although a great agreement between mixinglength models and solar surface convection simulations makes this model a helpful tool to describe the solar convection (e.g., Trampedach & Stein 2011), its reliability in the deep convection zone is still elusive: In fact, mixinglength models typically predict the existence of giant convective cells in the deep convection zone whose signals can hardly be captured by helioseismology (Hanasoge et al. 2012). As already pointed out in Section 4.2, the subadiabatic layer formed in the lower part of our numerical domain reflects the nonlocal effect of heat transport. The importance of the nonlocal treatment of mixinglength model has been repeatedly recognized mainly in the context of solar overshoot modeling. The nonlocal convective overshoot models naturally predict an extended weaklysubadiabatic layer above the radiation zone (Xiong & Deng 2001; Rempel 2004), which can typically extend up to depending the nonlocality of convection (Skaley & Stix 1991). Recently, Brandenburg (2016) modified the stellar mixinglength expression of the enthalpy flux by considering an additional term to incorporate the effects of nonlocal heat transport by strong downflow plumes, which enables even a weaklysubadiabatic region to transport the enthalpy upward. The importance of our study lies in that the effect of the enhanced subadiabatic layer on convection is connected for the first time to the change in the effective Prandtl number via the strongly nonlocal energy transport by convective downflows which has been widely recognized since Spruit (1997). Our results offer a possibility that, if downflows can retain their low entropy in effectively high regime and the nonlocality accordingly increases, the subadiabatic lower convection zone would be enhanced and extended upward in the lower part of the convection zone.
5.2 Possible Effects of Rotation
While we discussed so far the thermal effects on the convective velocity suppression in a nonrotating system for simplicity, it would be instructive to address some possible rotational effects on the proposed suppression mechanism. In general, Coriolis force tends to bend downflows into a longitudinal direction. It is thus expected that the subadiabatic layer near the bottom become more difficult to be formed due to the inefficient transport of low entropy by downflows. If the rotational effects are too strong (), the dependency of convective amplitudes is expected to diminish, because downflows are quickly distorted and the convectional structure should be dominated by the TaylorProudman state being independent of the thermal properties. Nonetheless, the proposed thermal effects may provide some important implications on the solar convection zone dynamics if we focus on a regime where rotational effects are moderate.
If the rotational effects are relatively weak (), the subadiabaticity at the bottom is expected to exhibit substantial latitudinal dependence since Coriolis forces can work effectively on downflows near the equatorial region whereas downflows near the polar region barely feel the rotational effects. As a result, a highlysubadiabatic layer should be formed near the pole and less subadiabatic (or even superadiabatic) layer should be formed near the equator, leading to a positive latitudinal gradient of superadiabaticity . This has a considerable influence on the thermal wind balance of the solar differential rotation. It has been argued that the negative latitudinal entropy gradient is needed for explaining the noncylindrical rotational profile of the observed differential rotation via the thermal wind balance (e.g., Miesch 2005). Although there are several theoretical explanations on the origin of this latitudinal entropy gradient (e.g., Kitchatinov & Rüdiger 1995; Rempel 2005; Masada 2011), the global convection simulations can hardly reproduce it in a selfconsistent manner and thus an adhoc latitudinal entropy variation has been commonly imposed at the bottom boundary in order to artificially break the TaylorProudman constraint (Miesch et al. 2006; Fan & Fang 2014). The latitudinal dependence of the superadiabaticity in the lower convection zone as a natural consequence of the proposed thermal effects in high regime may offer a promising mechanism to create the latitudinal entropy gradient at the base in a selfconsistent manner via the interaction with the clockwise (anti clockwise) meridional circulation in the southern (northern) hemisphere (Rempel 2005).
Another concern related to the rotational effects would be the angular momentum transport by the turbulent Reynolds stress (Kitchatinov & Rüdiger 1993). This issue is not only essential for our understanding of differential rotation but also for determining the meridional circulation structure which, despite its importance for fluxtransport dynamo model, varies significantly depending on theoretical models of the Reynolds stresses (Bekki & Yokoyama 2017), global parameters employed in 3D simulations (Passos et al. 2015; Featherstone & Miesch 2015), and inversion techniques of the helioseismic observations (Zhao et al. 2013; Rajaguru & Antia 2015). As already seen in Section 4.2, the convectional structure in high regime is qualitatively different from those in which has been commonly investigated in great detail: High thermal convection is characterized by strong downflowing plumes descending deeper into convection zone across several pressure scale heights, transporting highlyconcentrated cold entropy in a nonlocal way. It would be nontrivial, therefore, whether we can apply the knowledge on the turbulent Reynolds stress derived from calculations to the solar convection zone whose effective Prandtl number might be larger than unity. Further numerical work is needed in order to investigate how the properties of the turbulent momentum transport depend upon the Prandtl number as well as on the Rossby number .
6 Summary and Conclusions
In this paper, we have investigated one possible velocity suppression mechanism, motivated by the recentlyrecognized problem that the current solar convection simulations may be overestimating the amplitudes of deep convection (e.g., Hanasoge et al. 2016). We have especially focused on one interesting feature of the smallscale magnetism that it may decrease the effective thermal diffusivity by inhibiting the smallscale turbulent mixing of entropy between warm upflows and cold downflows (Hotta et al. 2015a). By conducting a set of thermal convection simulations where the effects of smallscale magnetism are incorporated as SGS diffusivities, we have shown in Section 4.2 that the convective velocities are suppressed as we decrease the effective thermal diffusivity through the enhancement of the subadiabatic layer which is formed near the base of the convection zone.
We further introduced the anisotropy of thermal diffusion in Section 4.3 to finally conclude that the decrease in the horizontal thermal diffusivity is critical for the convective velocity suppression, which is consistent with the results of smallscale dynamo simulations (Hotta et al. 2015a). This can be interpreted that the strong inhibition of horizontal entropy diffusion can promote the efficient transport of low entropy, and hence, promote the enrichment of the subadiabatic layer formation.
These results are synthesized to yield a picture of how the deep solar convection may be suppressed, which we outline as follows.
This is because the Lorentzforce of smallscale magnetism may act as an effective viscosity, and as a result, the effective Reynolds number for the flows on scales larger than the energy containing scale of the magnetic field is decreased. The presence of smallscale magnetic field, on the other hand, leads to smaller structure of entropy field because turbulent mixing of entropy is highly prohibited on scales smaller than the energy containing scale of the magnetic fields. This scale separation due to the smallscale magnetism between velocity spectrum (shift towards low wavenumber with effectively large ) and entropy spectrum (shift towards high wavenumber for effectively small ) can be roughly modeled as an increase in the effective Prandtl number .

Owing to a decrease in the effective thermal diffusivity , low entropy is conveyed almost adiabatically to the bottom, which results in the enhancement and extension of the subadiabatic layer near the base.
Note that the subadiabaticity here is not strong enough to be an overshoot region where downflows are quickly decelerated by buoyancy and therefore the enthalpy flux becomes downward. The enthalpy flux is in fact transported upward even inside the subadiabatic layer formed in the lower portion of the numerical domain, which means that this layer results from the nonlocal heat transport by strong downflow plumes (Brandenburg 2016; Käpylä et al. 2017b).

Downflows are subject to weak and gradual buoyant decelerations during the descent within the subadiabatic zone before reaching to the bottom beneath which a very thin and stiff overshoot layer is supposed to exist (Hotta 2017). Therefore, convective amplitude in the deeper subadiabatic layer is selectively suppressed, which may be observed as a reduction in the low wavenumber power of the horizontal velocity near the surface.
Although the velocity suppression mechanism described above is possible, the horizontal Prandtl number dependence of the convective rms velocities appears to be weak and still insufficient to explain the huge discrepancies between observations and simulations.
Further work is required to investigate whether this effect saturates or not at a finite effective Prandtl number, which is quite important when trying to apply this effect to the solar convective conundrum.
Future work will also focus on the effects of rotation on the degree of velocity suppression as discussed in Section 5.2 or on the applicability to the solar differential rotation problem especially to clarify whether the reduction in via the thermal effect can modify the Rossby number greatly enough to shift the differential rotation regime from antisolar to solarlike (Käpylä et al. 2014; Karak et al. 2015).
The authors thank an anonymous referee for helpful comments on the manuscript.
Y.B. further wish to thank Mark Miesch, Ben Brown, Mark Rast, and Sacha Brun for their thoughtful comments.
This work was supported by JSPS KAKENHI Grant Number 15H03640.
Y.B. also acknowledge financial support from the Leading Graduate Course for Frontiers of Mathematical Sciences and Physics (FMSP) of the University of Tokyo.
H.H. is also supported by MEXT/JSPS KAKENHI Grant Number JP16K17655 and JP16H01169.
Numerical computations were mostly carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan (CfCA), and in part on JAXA Supercomputer System generation 2 (JSS2).
Appendix
Comparison With Mhd Simulations
In this paper, we focus on the thermal effects of the smallscale magnetism that it can prohibit the effective thermal diffusivity. In fact, the highresolution simulations of smallscale dynamo should include both the dynamical effect (direct suppression by Lorentz force) and thermal effects. However, only few attention has been paid so far on the analysis of latter mechanism. We therefore present as an appendix the evidence that the thermal effect indeed exists in the smallscale dynamo simulations conducted by Hotta et al. (2015a). Their numerical domain consists of a simplified horizontallyperiodic box extending , but a realistic solar stratification, equation of state, and radiative diffusivity or energy flux values are used in a vertical direction extending .
Figure 13 shows the height dependence of the superadiabaticity for their runs H2048 and M2048 calculated by,
(31) 
where , , and are the dimensional entropy, background pressure scale height, and a specific heat at constant pressure, respectively. For making a comparison with our results easier, a typical superadiabaticity value which we use for normalizing is estimated as follows and overplotted in Figure 13 by a dotted line. The typical superadiabaticity should be on the order of Mach number square, according to the mixinglength argument. The modified Mach number at the base is derived in the same way as , where pressure and density at the base () are given by and , respectively. The typical convective velocity is calculated using the solar energy flux as with , which finally leads to a value . Although care must be taken in that the sound speed is artificially reduced in their calculations by a factor of , the comparison of the normalized superadiabaticity would be helpful to examine the influence of the change in .
First of all, it is obvious from Figure 13 that the subadiabatic layer formed near the base is enhanced and vertically extended with the inclusion of magnetic field. The degree of the subadiabaticity enhancement is similar to our results with Case: In both cases, the subadiabaticity near the base is increased by a factor of . Moreover, the vertically upward extension of the subadiabatic zone shifts the overall profile, which also leads to the decrease in the superadiabaticity in the bulk of the convection zone similarly to our results shown in Figure 9. The rms velocity for their MHD case is then suppressed by about with respect to HD case (see Figure 13 of Hotta et al. (2015a)), which is the same suppression degree seen in our model as shown in Figure 8. We therefore consider that in the real magnetized thermal convection, the thermal effect (velocity suppression via the enhanced subadiabatic layer due to the prohibition of the effective thermal diffusivity) may play a critical role in determining the convective amplitudes in addition to the dynamical effect.
It is also noteworthy that the superadiabaticity is not affected by the inclusion of magnetic fields near the highlysuperadiabatic surface, just like our result of Case H6V1. This, according to our discussion made in Section 4.2, indicates that the smallscale magnetic field does not inhibit the effective thermal conduction in a vertical direction: it only suppresses the “horizontal” effective thermal diffusivity, as already inferred by Hotta et al. (2015a). We thus make an argument that if we try to conduct LargeEddy Simulations with a SGS model on the effects of smallscale magnetic field, the anisotropy in the thermal diffusion term should be introduced and only the “horizontal” thermal diffusivity should be decreased.
References
 Bekki & Yokoyama (2017) Bekki, Y., & Yokoyama, T. 2017, ApJ, 835, 9
 Brandenburg (2016) Brandenburg, A. 2016, ApJ, 832, 6
 Brandenburg & Subramanian (2005) Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1
 Brummell et al. (2002) Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825
 Brun et al. (2004) Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073
 Cattaneo (1999) Cattaneo, F. 1999, ApJ, 515, L39
 Fan et al. (2003) Fan, Y., Abbett, W. P., & Fisher, G. H. 2003, ApJ, 582, 1206
 Fan & Fang (2014) Fan, Y., & Fang, F. 2014, ApJ, 789, 35
 Featherstone & Hindman (2016) Featherstone, N. A., & Hindman, B. W. 2016, ApJ, 818, 32
 Featherstone & Miesch (2015) Featherstone, N. A., & Miesch, M. S. 2015, ApJ, 804, 67
 Gastine et al. (2013) Gastine, T., Wicht, J., & Aurnou, J. M. 2013, Icarus, 225, 156
 Glatzmaier (1984) Glatzmaier, G. A. 1984, Journal of Computational Physics, 55, 461
 Greer et al. (2015) Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17
 Hanasoge et al. (2016) Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annual Review of Fluid Mechanics, 48, 191
 Hanasoge et al. (2012) Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proceedings of the National Academy of Science, 109, 11928
 Hotta (2017) Hotta, H. 2017, ApJ, 843, 52
 Hotta et al. (2014) Hotta, H., Rempel, M., & Yokoyama, T. 2014, ApJ, 786, 24
 Hotta et al. (2015a) —. 2015a, ApJ, 803, 42
 Hotta et al. (2015b) —. 2015b, ApJ, 798, 51
 Hotta et al. (2016) —. 2016, Science, 351, 1427
 Käpylä et al. (2014) Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43
 Käpylä et al. (2017a) Käpylä, P. J., Käpylä, M. J., Olspert, N., Warnecke, J., & Brandenburg, A. 2017a, A&A, 599, A4
 Käpylä et al. (2017b) Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017b, ArXiv eprints, arXiv:1703.06845
 Karak et al. (2015) Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&A, 576, A26
 Kitchatinov & Rüdiger (1993) Kitchatinov, L. L., & Rüdiger, G. 1993, A&A, 276, 96
 Kitchatinov & Rüdiger (1995) —. 1995, A&A, 299, 446
 Lord et al. (2014) Lord, J. W., Cameron, R. H., Rast, M. P., Rempel, M., & Roudier, T. 2014, ApJ, 793, 24
 Masada (2011) Masada, Y. 2011, MNRAS, 411, L26
 Miesch (2005) Miesch, M. S. 2005, Living Reviews in Solar Physics, 2, 1
 Miesch et al. (2008) Miesch, M. S., Brun, A. S., DeRosa, M. L., & Toomre, J. 2008, ApJ, 673, 557
 Miesch et al. (2006) Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618
 Miesch et al. (2000) Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593
 O’Mara et al. (2016) O’Mara, B., Miesch, M. S., Featherstone, N. A., & Augustson, K. C. 2016, Advances in Space Research, 58, 1475
 Ossendrijver (2003) Ossendrijver, M. 2003, A&A Rev., 11, 287
 Passos et al. (2015) Passos, D., Charbonneau, P., & Miesch, M. 2015, ApJ, 800, L18
 Rajaguru & Antia (2015) Rajaguru, S. P., & Antia, H. M. 2015, ApJ, 813, 114
 Rempel (2004) Rempel, M. 2004, ApJ, 607, 1046
 Rempel (2005) —. 2005, ApJ, 622, 1320
 Rempel (2014) —. 2014, ApJ, 789, 132
 Skaley & Stix (1991) Skaley, D., & Stix, M. 1991, A&A, 241, 227
 Spruit (1997) Spruit, H. 1997, Mem. Soc. Astron. Italiana, 68, 397
 Trampedach & Stein (2011) Trampedach, R., & Stein, R. F. 2011, ApJ, 731, 78
 Vögler & Schüssler (2007) Vögler, A., & Schüssler, M. 2007, A&A, 465, L43
 Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
 Warnecke et al. (2014) Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, ApJ, 796, L12
 Xiong & Deng (2001) Xiong, D. R., & Deng, L. 2001, MNRAS, 327, 1137
 Zhao et al. (2013) Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, Jr., T. L., & Hartlep, T. 2013, ApJ, 774, L29