HIGH-RESOLUTION CALCULATIONS OF THE SOLAR GLOBAL CONVECTION WITH THE REDUCED SPEED OF SOUND TECHNIQUE. I. THE STRUCTURE OF THE CONVECTION AND THE MAGNETIC FIELD WITHOUT THE ROTATION

The solar convection zone is filled with turbulent thermal convection due to its superadiabatic stratification, which is maintained by the radiative energy flux from the core of the Sun. Due to the influence of rotation, the anisotropic convection transports angular momentum as well as energy. This maintains mean flows such as the meridional flow and the differential rotation. These flows are thought to be important for the dynamo and related 11 year cyclic variations of the sunspot numbers (Parker 1955; Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Hotta & Yokoyama 2010a, 2010b). In addition, the turbulent thermal convection itself is also important for the generation of the magnetic field. Idealized simulations (e.g., Brandenburg et al. 1996; Cattaneo 1999) and the realistic calculations for the photosphere (Vögler & Schüssler 2007; Pietarila Graham et al. 2010) have suggested that chaotic motions have the ability to maintain magnetic energy through a small-scale dynamo process. Such dynamos have not been studied so far on the global scale, i.e., the entire convection zone.

The purpose of this study is to achieve multi-scale convection, including scales from 10 Mm (smaller than supergranulation) up to ∼200 Mm (global scale). We accomplish this by approaching the solar surface up to 0.99 R_☉ using unprecedented resolution, studying in particular the influence of the location of the upper boundary on the convective structure in the deeper parts of the convection zone. In addition, we use our approach to study the transport and the generation of the magnetic field by turbulent convection in the absence of rotation, i.e., we study the operation of a small-scale dynamo in the highly stratified convection zone.

The paper is organized as follows. We provide an overview of the reduced speed of sound technique (RSST) in Section 2, introduce our numerical setting in Section 3, and show the results in Section 4, present the properties based on the obtained spatial distributions in Section 4.1, discuss the energy balance using the integrated flux in Section 4.2, analyze the properties of the convection in cases without the magnetic fields using the spherical harmonic expansion in Section 4.3, analyze the cases with magnetic fields using spherical harmonic expansion and the probability density function in Section 4.4, investigate the transportation and the generation of the magnetic field in the convection zone in Section 4.5, and summarize our investigation in Section 5.

Developing numerical simulations of solar global convection is challenging due to the high resolution and the broad range of dynamical timescales that are required. One of the most paramount difficulties arises from the high speed of sound and the related low Mach-number flows throughout most of the convection zone. At the base of the convection zone, the speed of sound is about 200 km s⁻¹, while the speed of convection is thought to be 50 m s⁻¹ (e.g., Stix 2004). The time step must therefore be shorter due to the Courant–Friedrichs–Lewy (CFL) condition in an explicit fully compressible method, even when we are interested in the phenomena related to convection. In order to avoid this situation, the anelastic approximation is frequently adopted in which the mass conservation equation is replaced with ∇ · (ρ₀v) = 0, where ρ₀ is the reference density and v is the fluid velocity. In this approximation, the speed of sound is assumed to be infinite, and we need to solve the elliptic equation for the pressure, which filters out the propagation of the sound wave. Since the anelastic approximation is well applicable deep in the convection zone and the time step is no longer limited by the large speed of sound, the solar global convection has been investigated with this method in numerous studies about the differential rotation, the meridional flow, the global dynamo, and the dynamical coupling of the radiative zone (Miesch et al. 2000, 2006, 2008; Brun & Toomre 2002; Brun et al. 2004, 2011; Browning et al. 2006; Ghizaru et al. 2010).

There are, however, two drawbacks in the anelastic approximation. The first is the breakdown of the approximation near the solar surface. Since the convection velocity becomes higher and the speed of sound becomes lower in the near-surface layer (>0.98 R_☉), these have similar values and the anelastic approximation cannot be applied. Since the interaction between the small-scale convection generated in the near-surface layer (granulation and supergranulation) and large-scale convection (giant cell and global one) is an important topic, the connection between the near-surface layer and the global convection has been investigated (e.g., Augustson et al. 2011). The global calculation, however, which includes all such multiple scales, has not yet been achieved. The second drawback is the difficulty arising from the increase of the required resolution. The pseudo-spectral method based on the spherical harmonic expansion is frequently adopted, especially for solving the elliptic equation of the pressure. In this method, non-linear terms require the transformation of the physical variable from the real space to the spectral space and vice versa at every time step. The calculation cost of the transformation is estimated as $\mathcal {O}(N_\theta ^2N_\phi \log N_\phi)$ due to the absence of a fast algorithm for the Legendre transformation, which is as powerful as the fast Fourier transformation, where N_θ and N_ϕ are the maximum mode numbers in the latitude and the longitude, respectively. Thus, the computational cost of this method is significant, limiting the achievable resolution. Due to this, some numerical calculations of the geodynamo adopt the finite difference method in order to achieve high resolution (Kageyama et al. 2008; Miyagoshi et al. 2010). As explained above, when the near-surface layer is included in the calculation, the typical convection scale becomes small and a large number of grid points is required. The resolution is one of the keys to access the solar surface. We note that some studies using the finite difference have been done in the stellar or solar context using a moderate ratio of the speed of sound and the convection velocity with adjustments of the radiative flux and the stratification (Käpylä et al. 2011, 2012). Although this type of approach provides insights for the maintenance of the mean flow and the magnetic field, proper reproduction using the solar parameters, such as the stratification, the luminosity, the partial ionization effect and the rotation, and direct comparison with actual observations cannot be achieved.

The RSST (Rempel 2005, 2006; Hotta et al. 2012b) can overcome these drawbacks while avoiding the severe time step caused by the speed of sound. In the RSST, the equation of continuity is replaced with

$$\begin{eqnarray} &&\frac{\partial \rho }{\partial t}=-\frac{1}{\xi ^2}\nabla \cdot (\rho {\bf v}). \end{eqnarray} \tag{ 1 }$$

The speed of sound is then reduced ξ times, but the dispersion relationship for sound waves remains (where wave speed decreases equally for all wavelengths). This technique does not change the hyperbolic character of the equations that can be integrated explicitly. Due to this hyperbolicity, only the local communication is required, which decreases the communication overhead in parallel computing. A simple algorithm and low cost of communication are significant factors that enable high-resolution calculations. Hotta et al. (2012b) investigated the validity of the RSST in a thermal convection problem. It was concluded that the RSST is valid when the Mach number defined with the root mean square (rms) velocity and the reduced speed of sound $\hat{c_s}$ is smaller than 0.7. Another advantage of this method is its accessibility to the real solar surface with an inhomogeneous ξ. The Mach number varies substantially in the solar convection zone. When a moderate or zero reduction of the speed of sound is used in the near-surface layer while a large ξ around the bottom part of the convection zone is used, the properties of the thermal convection, even including the surface, can be properly investigated without undermining the physics. Hotta et al. (2012b) confirmed that the inhomogeneous ξ is valid again when the Mach number is less than 0.7.

We solve the three-dimensional magnetohydrodynamic equations with the RSST in the spherical geometry (r, θ, ϕ) as

$$\begin{eqnarray} &&\frac{\partial \rho _1}{\partial t} = {-}\frac{1}{\xi ^2}\nabla \cdot (\rho _0 {\bf v}), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \rho _0\frac{\partial {\bf v}}{\partial t} &=& {-}\rho _0({\bf v}\cdot \nabla){\bf v} - \nabla \left(p_1 + \frac{B^2}{8\pi }\right)\nonumber\\ &&+\nabla \cdot \left(\frac{{\bf BB}}{4\pi }\right) - \rho _1g{\bf e_r}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} && \frac{\partial {\bf B}}{\partial t}=\nabla \times \left({\bf v}\times {\bf B}\right), \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \rho _0 T_0\frac{\partial s_1}{\partial t} &=& {-}\rho _0 T_0({\bf v}\cdot \nabla)s_1\nonumber\\ &&+ \frac{1}{r^2} \frac{d}{dr} \left(r^2 \kappa _\mathrm{r}\rho _0 c_\mathrm{p}\frac{dT_0}{dr} \right) + \Gamma, \end{eqnarray} \tag{ 5 }$$

where ρ, p, s, T, v, and B are the density, the gas pressure, the specific entropy, the temperature, the fluid velocity, and the magnetic field, respectively. The subscript 1 denotes the fluctuation from the time-independent spherically symmetric reference state, which has subscript 0. g, κ_r, and Γ are the gravitational acceleration, the coefficient of the radiative diffusivity, and the surface cooling term, respectively. The setting of these values is explained in the following paragraph. In order to close the MHD system, we adopt a non-ideal equation of state as

$$\begin{eqnarray} &&p_1 = \left(\frac{\partial p}{\partial \rho }\right)_s \rho _1 + \left(\frac{\partial p}{\partial s }\right)_\rho s_1. \end{eqnarray} \tag{ 6 }$$

Since our upper boundary is at r = 0.99 R_☉ at maximum, in the near-surface region partial ionization is important (see Figure 1(f)) and is included in our treatment (see Appendix A for details) by using the OPAL repository with the solar abundance of X = 0.75, Y = 0.23, and Z = 0.02, where X, Y, and Z specify the mass fraction of hydrogen, helium, and other metals, respectively.

Zoom In Zoom Out Reset image size

Note that we do not have any explicit turbulent thermal diffusivity and viscosity (Miesch et al. 2000) in order to maximize the fluid and magnetic Reynolds number, but use the artificial viscosity introduced in Rempel et al. (2009). Our treatment ensures that the dissipated kinetic and magnetic energies with the artificial viscosity are converted to the internal energy. Note also that the rotation is not included in this study, since we focus on the connection of thermal convection between the small and large scales and its local dynamo action on the global scale. The combination of rotation and stratified convection leads to differential rotation and additional turbulent induction effects that enable large-scale dynamo action. It is difficult to distinguish the local dynamo action from global dynamo action if rotation is present. Since the rotation effect in the sun is significant, we will carry out studies of rotating convection in the future.

Figure 1 shows our reference state in comparison with Model S (Christensen-Dalsgaard et al. 1996). The reference stratification is determined by solving the one-dimensional hydrostatic equation and the realistic equation of state as

$$\begin{eqnarray} && \frac{dp_0}{dr}=-\rho _0 g, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} && \rho _0 = \rho _0(p_0,s_0), \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} && \frac{ds_0}{dr} = 0. \end{eqnarray} \tag{ 9 }$$

Equation (8) is calculated with the OPAL repository including partial ionization. We set an adiabatic stratification as a reference and initial state. The stratification becomes superadiabatic after convection develops as a consequence of radiative heating near the bottom and radiative cooling at the top as described below. The gravitational acceleration and the radiative diffusion are adopted from Model S. The boundary is set at r = 0.998 R_☉ with the values from Model S, and the equations are integrated inward.

On the actual Sun, the surface is continuously cooled by the radiation. Since our boundary is not located on the real solar surface (even though it is unprecedentedly closer to the real surface), we add an artificial cooling (Γ) as

$$\begin{eqnarray} &&\Gamma (r) = -\frac{1}{r^2}\frac{\partial }{\partial r}(r^2 F_\mathrm{s}), \end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray} && r^2F_\mathrm{s}(r) = r^2_\mathrm{min}F_\mathrm{r}(r_\mathrm{min}) \exp \left[-\left(\frac{r-r_\mathrm{max}}{d_\mathrm{c}}\right)^2\right], \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&F_\mathrm{r}(r) = -\kappa _\mathrm{r}\rho _0 c_\mathrm{p}\frac{dT_0}{dr}, \end{eqnarray} \tag{ 12 }$$

where r_min and r_max denote the location of the bottom and the top boundary, respectively. This procedure ensures that the radiative luminosity imposed at the bottom is released through the top boundary. Since the realistic simulation for the near-surface layer shows that the thickness of the cooling layer by the radiation is similar to the local pressure scale height (e.g., Stein et al. 2009), we adopt two pressure scale heights for the thickness of the cooling layer, i.e., d_c = 2H_p0(r_max) except for cases H1 and M1 (see Table 1), where H_p0 = p₀/(ρ₀g) is the pressure scale height.

We carry out three calculations named H0, H1, and H2, with different settings that are primarily hydrodynamic. In case H0, the top boundary is located at r = 0.99 R_☉, and the density contrast ρ₀(r_min)/ρ₀(r_max) exceeds 600. To our knowledge, this is the largest value achieved so far in numerical simulations of solar global convection. In cases H1 and H2, the top boundary is at r = 0.96 R_☉, and the density contrast is around 40. In case H1, the thickness of the cooling layer is the same as case H0, i.e., d_c = 3740 km, which is the two pressure scale heights at r = 0.99 R_☉, while case H2 adopts two scale heights at its top boundary (r_max = 0.96 R_☉) for the thickness (d_c). Our magnetic runs M0, M1, and M2 use the same settings and differ only through the initial magnetic field added. We note that the thickness of the cooling layer (d_c) has almost the same role as the value of the turbulent thermal diffusivity on the entropy that is adopted in the anelastic spherical harmonic simulation (Miesch et al. 2000, 2008).

In this study, the factor of the RSST is set to make the adiabatic reduced speed of sound uniform in space. The adiabatic speed of sound is defined as

$$\begin{eqnarray} &&c_\mathrm{s}(r) = \sqrt{\left(\frac{\partial p}{\partial \rho }\right)_s}. \end{eqnarray} \tag{ 13 }$$

Then the factor of the RSST is set as

$$\begin{eqnarray} &&\xi (r) = \xi _0 \frac{c_\mathrm{s}(r)}{c_\mathrm{s}(r_\mathrm{min})}. \end{eqnarray} \tag{ 14 }$$

In this study, we adopt ξ₀ = 120 for all the calculations, thus the reduced speed of sound ${\hat{c}_\mathrm{s}} \equiv c_\mathrm{s}/\xi = 1.88\ \mathrm{km\ s^{-1}}$ at all depths. Hotta et al. (2012b) suggest that the RSST is valid for the thermal convection under the criterion of $v_\mathrm{rms}/\hat{c}_\mathrm{s} < 0.7$, where v_rms is the rms convection velocity. Thus, we can properly treat the convection with v_rms < 1.3 km s⁻¹ in this model. It is also suggested in Hotta et al. (2012b) that the total mass is not conserved with the inhomogeneous ξ, and we cannot avoid long-term drift from the reference state, i.e., the mass continuously decreases or increases. In this study, however, we adopt a different way to avoid this type of long-term drift. When the equation of continuity is treated as

$$\begin{eqnarray} &&\frac{\partial }{\partial t}(\xi ^2 \rho _1) = -\nabla \cdot \left(\rho _0 {\bf v}\right), \end{eqnarray} \tag{ 15 }$$

the value $\hat{M}$ is conserved in a rounding error with appropriate boundary conditions, where

$$\begin{eqnarray} &&\hat{M} = \int _V \xi ^2 \rho _1 dV. \end{eqnarray} \tag{ 16 }$$

Although the radial distribution of the density is different from the original one, the fluctuation part remains small (e.g., ρ₁/ρ₀ ∼ 10⁻⁶), and this does not affect the character of the thermal convection. It is confirmed in Hotta et al. (2012b) that statistical features are not influenced by the inhomogeneous ξ.

We use the stress-free and the impenetrable boundary conditions for the fluid velocities v_r, v_θ, and v_ϕ. The free boundary condition is adopted for the density and the entropy, i.e., ∂ρ₁/∂r = ∂s₁/∂r = 0. The magnetic field is vertical at the top boundary, and the perfect conductor boundary condition is used at the bottom boundary.

ρ₁, B_r, B_θ, B_ϕ, and s₁ are zero initially. The fluid velocities v_r, v_θ, and v_ϕ have small random values. After the convection reaches a statistically steady state, a uniform magnetic field (B_ϕ = 100 G) is added. Although a net toroidal flux exists initially using this condition, it disappears through the upper boundary after around 75 days since horizontal magnetic field is zero there. Thus, we do not need to consider the influence of the net flux for our local dynamo study.

Under the setting explained above, we solve Equations (2)–(5) by using the fourth-order Runge–Kutta scheme for the time-integration and the fourth-order space-centered derivative (Vögler et al. 2005). The divergence-free condition, i.e., ∇ · B = 0, is maintained with the diffusion scheme for each Runge–Kutta loop (see the Rempel et al. 2009 appendix for details). In order to include all the spherical shell, we adopt the Yin–Yang grid (Kageyama & Sato 2004). The Yin–Yang grid is a set of two congruent spherical geometries. They are combined in a complemental way to cover a whole spherical shell. The boundary condition for each grid is calculated using the interpolation of the other grid. In all of our calculations, each grid covers 0.715 R_☉ < r < r_max, π/4 − δ_θ < θ < 3π/4 + δ_θ, and −3π/4  −  δ_ϕ < ϕ < 3π/4 + δ_ϕ, where δ_θ and δ_ϕ are the margins for the interpolation. Here δ_θ = 3Δθ/2 and δ_ϕ = 3Δϕ/2 are adopted for the interpolation with the third-order function, where Δθ and Δϕ are the grid spacing in the latitudinal and longitudinal directions. Here r = 0.715 R_☉ = r_min is the location of the base of the convection zone, and r_max is taken as a parameter (r_max = 0.99 R_☉ for the typical case). Although the boundary condition for the Yin grid is applied on the edge of the Yin grid, the boundary condition for the Yang grid is applied on the edge of the Yin grid in order to avoid the double solution in the overlapping area of the Yin–Yang grid. Figure 2 shows this configuration. The thick red lines show the location of the horizontal boundary for both Yin and Yang grids.

Zoom In Zoom Out Reset image size

The horizontal grid spacing is 1100 km at the top boundary and 375 km radially. With the Yin–Yang geometry, the number of grid points is 1024(N_θ) × 3072(N_ϕ) × 2. The last factor 2 indicates a Yin and Yang pair. The number of grid points in the radial direction is shown in Table 1. Since this resolution in Yin–Yang geometry has almost the same quality as that of 512(r) × 2048(θ) × 4096(ϕ) in ordinary spherical geometry in cases H0 and M0, we succeed in doubling the resolution in each direction from the previous study (Miesch et al. 2008). After the uniform magnetic field is added, the calculations are newly named M0, M1, and M2, which use the results of H0, H1, and H2, respectively. Using a hybrid MPI and automatic intra-node parallelization approach, the code scales efficiently up to 10⁵ cores with almost linear weak scaling, and achieves a 14% performance at maximum on the RIKEN K-computer in Japan. Since our code includes almost no global communication among cores, this linear scaling is expected to hold further with larger core counts.

4.1. Structure of Convection and the Magnetic Field

Figure 3 shows the rms velocities in cases H0, H1, and H2. The maximum rms velocity is 5 × 10² m s⁻¹ at the top boundary in case H0. Since the Mach number determined with the reduced speed of sound $\hat{c}_\mathrm{s}$ is always under 0.3 throughout the convection zone, the requirement for the validity of the RSST in this study is well satisfied (Hotta et al. 2012b). The difference between cases H1 and H2 is mostly a local effect. The thinner cooling layer in case H1 leads to stronger acceleration of plasma and causes a higher rms velocity close to the top boundary. In case H0, the rms velocity is significantly larger in the near-surface layer (>0.96 R_☉) since the energy has to be transported with a small background density. This has some influence on the velocity throughout the convection zone.

Zoom In Zoom Out Reset image size

Figure 4 shows the radial velocity (v_r) around the top boundary in case H0 (r_max = 0.99 R_☉, which is 7 Mm below the photosphere) in the orthographic projection. (The corresponding animation is provided in the online journal.) Note that since the shown location is close to the impenetrable top boundary, the value of the radial velocity is rather small. Since in the area near the upper boundary the pressure scale height is less than 2 Mm, the convection pattern shows significantly small cells of about ∼7 Mm. The typical cell size is slightly smaller than the supergranulation that is observed in the photosphere. This study is the first that well resolves the 10 Mm scale convection pattern in a calculation of the solar global convection zone. Figures 5(a)–(c) show the radial velocity at r = 0.99 R_☉, r = 0.95 R_☉, and r = 0.85 R_☉ in case H0 by using the orthographic projection. In deeper layers, the pressure scale height becomes larger (see Figure 1(e)), and the convection cell increases in size. A detailed analysis using the spherical harmonic expansion of the convective structure is shown in Section 4.3.

Download figure:

Video Standard image High-resolution image

Zoom In Zoom Out Reset image size

Figures 6(a), (d), and (g) show the zoomed-in contour of the radial velocity in case M0, i.e., after the inclusion of the magnetic field in H0. The region is indicated by the white rectangle in Figures 5(a)–(c). Figures 6(c), (f), and (i) show the vorticity (ω_r = (∇ × v)_r) at r = 0.99 R_☉, 0.95 R_☉, and 0.85 R_☉, respectively. As already seen in case H0 (Figure 5), it is clear that the scale of the thermal convection pattern significantly depends on the depth. In addition, the large-scale downflow is associated with the small-scale and strong vorticity in the deeper layer (especially in r = 0.85 R_☉). Figure 7(a) shows ρ₀[s₁(r, θ, ϕ = 0) − 〈s₁〉] in the meridional plane, where 〈s₁〉 is the horizontal average of the entropy. The low and high entropy materials correspond to the downflow and upflow, respectively. In the near-surface region, the convection structure shows a combination of broad upflows and narrow downflows ∼7 Mm forming at the top boundary. These small-scale downflows merge in the middle of the convection zone and build large-scale downflow. Although the overall structure of such convection is large, there is a superimposed turbulent pattern especially in the downflow region, which is shown in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 5(d), (e), and (f) show the results of case H1 with a different location of the top boundary (r_max = 0.96 R_☉) at r = r_max, r = 0.95 R_☉, and r = 0.85 R_☉, respectively. Since the location and the pressure scale height of the shown images is different between case H0 in Figure 5(a) (r = r_max = 0.99 R_☉) and H1 in Figure 5(d) (r = r_max = 0.96 R_☉), it is natural that the convective structures are much different, i.e., those in H0 have smaller-scale convection than in H1. It is more important that the structures at the same depth, 0.95 R_☉, are significantly different from each other between these cases (H0 in Figure 5(b) and H1 in Figure 5(e)). The small-scale downflow plumes penetrate near the surface layer and influence its structure (see also Figure 7(a)). When the downflow goes deeper, the influence becomes smaller. This is seen in the comparison of Figure 5(c) (r_max = 0.99 R_☉ and r = 0.85 R_☉) and Figure 5(f) (r_max = 0.96 R_☉ and r = 0.85 R_☉), where the difference of convection structure seems insignificant. Figures 5(g)–(i) show the results in case H2, in which the location of the top boundary is the same as in H1 but the thickness of the cooling layer is larger. The convective structure around the top boundary shows the largest scale (Figure 5(g)), while again the difference becomes smaller in the deeper layer (Figures 5(c), (f), and (i)). This is also shown in Figure 8 by the areas occupied by the upflow and the downflow, i.e., positive and negative radial velocities (v_r). Up to the middle of the convection zone (<0.9 R_☉), cases H0, H1, and H2 all show similar behavior, i.e., the fractional area by the upflow is larger (∼65%) but decreases below ∼0.85 R_☉ to equal to that of the downflow. Interestingly, this behavior is quantitatively the same in spite of the significant difference between H0 and H1 in the density contrast (see Table 1).

Zoom In Zoom Out Reset image size

4.2. Integrated Energy Flux

Figure 9 shows the integrated fluxes. The integrated enthalpy flux (L_e), the integrated kinetic flux (L_k), the radiative luminosity (L_r), and the luminosity form of the surface cooling (L_s) are defined as

$$\begin{eqnarray} L_\mathrm{e} &=& \int _s \left[\rho _0e_1 + p_1 - \frac{p_0\rho _1}{\rho _0}\right]v_r dS, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} L_\mathrm{k} &=& \int _s \frac{1}{2}\rho _0v^2 v_r dS, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} L_\mathrm{r} &=& \int _s F_\mathrm{r}(r) dS, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} L_\mathrm{s} &=& \int _s F_\mathrm{s}(r) dS. \end{eqnarray} \tag{ 20 }$$

The radiative flux (F_r) and surface cooling flux (F_s) are defined in Equations (11) and (12), respectively. The derivation of the enthalpy flux is given in Appendix B. Figure 9(a) shows the integrated fluxes in case H0 (r_max = 0.99 R_☉ and d_c = 3740 km). The total integrated flux L_t = L_e + L_k + L_r + L_s is almost constant along the depth. This indicates that the convection zone in our calculation is in energy equilibrium. This takes about 50 days. Note that since we do not use a conservative form for the total energy nor estimate the energy flux contributions caused by the artificial diffusivities, the total flux is not completely constant. The enthalpy flux transports twice the solar luminosity upward at maximum, and the kinetic flux transports almost the same amount of the energy as the solar luminosity downward at maximum. Although the kinetic energy flux is frequently ignored in the one-dimensional mixing length model (e.g., Stix 2004), our result shows the importance of the kinetic flux. This has been already suggested by Miesch et al. (2008). Figures 9(b) and (c) show the results in cases H1 and H2, respectively. The integrated fluxes show almost the similar behavior as those in case H0, but the maximum absolute values of the enthalpy and kinetic flux are smaller. Since these absolute values gradually decrease from H0 to H2, we conclude that both the thickness of the cooling layer and the location of the upper boundary contribute to this result. It is possible that when the upper boundary becomes closer to the real solar surface and the cooling layer becomes thinner, the absolute values of the enthalpy flux and the kinetic flux become even larger than our case H0.

Zoom In Zoom Out Reset image size

Figure 10 shows the integrated enthalpy flux and the kinetic flux transported by upflow (L_eu, L_ku) and downflow (L_ed, L_kd), respectively. Note that the enthalpy flux by upflow and downflow is estimated with the perturbation from the reference state. Regarding the enthalpy flux, both upflow and downflow transport energy upward. The downflow transports most of the energy (>70%). We note that the enthalpy flux of the upflow shows the negative value near the bottom boundary since the cool fluid bounces at the boundary and moves upward. Regarding the kinetic energy flux, upflows (downflows) transport energy upward (downward). The larger kinetic energy flux of the downflows make the kinetic energy flux negative. These results show that downflows play a key role in the transports of energy in the solar convection zone. While we find a significant difference in the overall amplitude among cases H0 to H2, the ratio of contributions from up- and downflows does not change much despite the significant difference in the density contrast.

Zoom In Zoom Out Reset image size

4.3. Analysis Using Spherical Harmonics for the Hydrodynamic Cases

In this section, the results of the analysis using the spherical harmonics expansion are shown. We focus on the following question. How do the location of the upper boundary and the thickness of the surface cooling layer influence the convective structure throughout the convection zone?

A real function f(θ, ϕ) can be expressed in spherical harmonics as

$$\begin{eqnarray} &&f(\theta,\phi) = \sum ^{l_\mathrm{max}}_{l=0}\sum ^l_{m=-l} f_{lm}Y_{lm}(\theta,\phi), \end{eqnarray} \tag{ 21 }$$

where Y_lm(θ, ϕ) is the spherical harmonics for degree l and order m. In the analysis the absolute total value along m without m = 0,

$$\begin{eqnarray} &&\bar{f}_l = \sqrt{\sum _{m=-l}^l |f_{lm}|^2}, \end{eqnarray} \tag{ 22 }$$

is shown. The value is normalized in order to satisfy the relation

$$\begin{eqnarray} &&\frac{\int _\Omega (f(\theta,\phi)-\langle f(\theta,\phi)\rangle)^2\sin \theta d\theta d\phi }{4\pi } = \sum _{l=1}^{l_\mathrm{max}} \bar{f}_l^2. \end{eqnarray} \tag{ 23 }$$

Our spherical harmonic analysis is performed using the freely available software archive SHTOOLS (shtools.ipgp.fr).

Figures 11 and 12 show the spectra of the radial velocity (v_r) and the latitudinal velocity (v_θ), respectively, as a function of the horizontal wavelength (L_h = 2πr/l), where l is the spherical harmonic degree, i.e., the horizontal wavenumber. The black, blue, and red lines show the results in case H0 (r_max = 0.99 R_☉ and d_c = 3740 km), H1 (r_max = 0.96 R_☉ and d_c = 3740 km), and H2 (r_max = 0.96 R_☉ and d_c = 18780 km), respectively. The black line in Figure 11(a) shows a peak around L_h ∼ 7–8 Mm (see also Figure 4). This peak moves to the larger scale L_h with increasing depth. At r = 0.80 R_☉, the peak is around L_h ∼ 300 Mm (Figure 11(f), black line). This reflects the variation of the pressure scale height in the solar model, H_p = 1.9 Mm at r = 0.99 R_☉ and H_p = 44 Mm at r = 0.8 R_☉ (see Figure 1(e)). Again the peak is in the smaller scale (∼7 Mm) at the bottom boundary r = r_min = 0.715 R_☉. This is caused by the collision between the downflow and the impenetrable bottom boundary. In this process, a thin boundary layer whose thickness is determined by the numerical resolution is formed. The scale of turbulence near the boundary (Figure 11(f)) is a consequence of this boundary layer. A comparison of black lines in Figures 11(a) and 12(a) indicates that the peak of the latitudinal velocity is at a much larger scale (∼400 Mm) than the radial velocity, even close to the upper boundary. We note that almost all the spectra have this feature around L_h ∼ 400 Mm, i.e., the peak or the bended feature in which the power-law index varies. The length of the scale 400 Mm is twice the thickness of the convection zone, i.e., the radial extent of our computational domain. This result is consistent with our previous study (Hotta et al. 2012a), which argues that the typical scale of the convection is determined by the pressure (or density) scale height or the height of the computational domain. Figure 12(a) shows that the horizontal velocity is more likely to be influenced by the large-scale structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The influences from the location of the boundary and the thickness of the cooling layer are discussed by comparing the black, blue, and red lines in Figures 11 and 12. The common feature is that the spectra do not depend on these two factors in the deep layer (r ⩽ 0.85 R_☉: panels (d)–(f)). This suggests that the small-scale convection caused by the short pressure scale height or the thin cooling layer cannot influence the convection in the deeper region. On the top boundary (Figures 11(a) and 12(a)), the small-scale convection is suppressed gradually from case H0 to H2. This is confirmation of our understanding from the appearance of the convection pattern in Section 4.1. In the near-surface layer (r = 0.95 R_☉), the difference still remains unchanged. Moving the location of the top boundary upward increases the amplitude of the fluid velocity at all scales (black line). Reducing the thickness of the cooling layer while keeping the location of the boundary unchanged leads to an increase of velocity on small scales (<50 Mm), while larger scales are unaffected (blue line). We point out that when the surface layer (0.96 R_☉ < r < 0.99 R_☉) is included, the spectrum of the latitudinal velocity (v_θ) is flat from the middle to the small scale (10 Mm  <  L_h  <  40 Mm in Figure 12(b)). This is caused by the penetration of the corresponding scale plume from the near-surface layer.

4.4. Analysis Using Spherical Harmonics and the Probability Density Function for the Magnetohydrodynamic Cases

In this subsection, we analyze the results of the magnetohydrodynamic calculation in cases M0, M1, and M2 using spherical harmonics and probability density function. We note that cases M0, M1, and M2 use the same parameters as cases H0, H1, and H2, respectively, with the magnetic field. We focus on the influence of the location of the boundary layer and the thickness of the cooling layer on the structure of the magnetic field.

Figure 13 shows the time evolution of the magnetic energy (B²/(8π)) averaged over the simulation domain at each time step. The initial linear growth stops around 10 days after the input of the seed field in every case. Even after that, the magnetic energy continues to increase gradually with a rather small growth rate until it saturates after around 150 days in all cases. A comparison among the cases is performed using the data from t = 115 days to t = 162 days in which the generation of the magnetic field is almost saturated. Basically, the differences among the cases are insignificant, although case M0 saturates with a slightly smaller average magnetic energy. Since the equipartition magnetic field strength in the near-surface region is smaller due to the small density (ρ₀), the increase of the volume in case M0 causes the slight decrease in the average magnetic energy. This conclusion can be supported by the fact that the average value over r_min < r < 0.96 R_☉ in M0 is larger, as shown by the dashed line in Figure 13. The growth rate of this energy (dashed line) is higher than those in M1 and M2 at t < 50 days. The timescale of the convection in the near-surface layer is short. The generated magnetic field is transported downward (see also the discussion about the pumping in Section 4.5).

Zoom In Zoom Out Reset image size

Figure 14 shows the spectra of the magnetic energy B²/(8π) at selected depths. Similar to the results introduced in the previous sections, the differences can be seen in the near-surface layer, and these differences become insignificant as we go to the deeper layers.

Zoom In Zoom Out Reset image size

Figure 15 shows the probability density functions (PDFs) for the three components of velocity (v_r, v_θ, and v_ϕ) and the three components of the magnetic field (B_r, B_θ, and B_ϕ), the radial vorticity ω_r, the horizontal divergence ζ, and the temperature perturbation T₁ in case M0 at t = 115 days. Although the rotation is not included in this study, some features of the PDFs are similar to findings from previous studies including rotation (Brun et al. 2004; Miesch et al. 2008). In this study, the PDF is the normalized histogram on a horizontal surface, corrected for the grid convergence at the poles. Figure 15(a) shows the significant asymmetry in the radial velocity. This reflects the asymmetry between up- and downflows due to stratification (see also Figures 8 and 10). The horizontal velocities (v_θ and v_ϕ) almost show a Gaussian distribution (see also Figure 16). The magnetic fields have high intermittency compared with the velocities (Brandenburg et al. 1996). Despite our asymmetric initial condition for the longitudinal magnetic field B_ϕ = 100 G, the PDFs show a close to symmetric distribution peaked at zero after a sufficiently long temporal evolution. The maximum value of the strength of the magnetic field is around 10⁴ G. The PDF for the radial vorticity is similar to that of the magnetic field, i.e., with high intermittency. This can be expected based on the similarity between the induction and vorticity equations. The horizontal divergence has similarity to the radial vorticity ω_r in the convection zone with high intermittency, while ζ shows the asymmetry near the top boundary is similar to the radial velocity v_r. The temperature perturbation has significant asymmetry, which also reflects the asymmetry between the upflow and downflow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to evaluate the influence of the location of the boundary condition and the thickness of the cooling layer, we investigate the moments of the PDF, in particular the kurtosis $\mathcal {K}$ and the skewness $\mathcal {S}$ as

$$\begin{eqnarray} && \mathcal {K}=\frac{1}{\sigma ^4}\int (x-\langle x\rangle)^4f(x)dx, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} && \mathcal {S}=\frac{1}{\sigma ^3}\int (x-\langle x\rangle)^3f(x)dx, \end{eqnarray} \tag{ 25 }$$

where f(x) is the PDF, x is each variable, and σ is the standard deviation as

$$\begin{eqnarray} && \sigma = \sqrt{\int (x-\langle x\rangle)^2f(x)dx}. \end{eqnarray} \tag{ 26 }$$

The kurtosis $\mathcal {K}$ and the skewness $\mathcal {S}$ denote the intermittency and the asymmetry of the distribution, respectively. We note that the PDF is normalized as ∫f(x)dx = 1. For example, the Gaussian PDF is characterized by $\mathcal {K}=3$ and $\mathcal {S}=0$. Figure 16 shows the distribution of the kurtosis and the skewness for the velocity and the magnetic field. As noted above, the horizontal velocities (v_r and v_ϕ) have almost a Gaussian distribution, i.e., $\mathcal {K}\sim 3$ and $\mathcal {S}\sim 0$ throughout the convection zone. The radial velocity has high intermittency and asymmetry in the convection zone. The magnetic field has intermittency and almost symmetric distribution. These features are common among cases M0, M1, and M2. Quantitatively, the kurtosis and skewness agree with each other in all cases. While the values in the near-surface region (r > 0.85 R_☉) are influenced by the two factors considered, the location of the upper boundary and the thickness of the cooling layer, in the deeper region the values converge.

4.5. Generation and Transportation of the Magnetic Field

In this subsection, we investigate the generation and transport of the magnetic field by turbulent thermal convection. The global structure of the mutual interaction between the plasma and the magnetic field is our paramount interest.

Figure 6 shows the radial velocity v_r, the radial magnetic field B_r, and the radial vorticity ω_r at different depths (r = 0.99 R_☉: panels (a)–(c), r = 0.95 R_☉: panels (d)–(f), r = 0.85 R_☉: panels (g)–(i)). As introduced in Section 4.1, there is good agreement between the downflow and regions with a large amplitude of the radial vorticity. This means that downflows, especially in the deeper region, include most of the turbulent small-scale horizontal motions. We can also see the preferential association of a strong magnetic field with downflows and regions with strong vorticity in the constant-depth plane (middle and right columns of Figure 6) and also in the meridional plane (Figure 7). In order to investigate this aspect quantitatively, we estimate the correlation between the radial velocity v_r and the absolute value of the magnetic field B, i.e., 〈v_r, B〉 in Figure 17(a), where our definition of the correlation between quantity A and B is defined as

$$\begin{eqnarray} &&\langle A,B\rangle = \frac{\int AB dS}{\sqrt{\int A^2dS}\sqrt{\int B^2 dS}}. \end{eqnarray} \tag{ 27 }$$

Figure 17(a) shows the negative value in most of the convection zone. This means that the magnetic field is preferentially found in downflows. This is also seen in the joint PDFs of v_r and B in Figure 18. The figure shows the asymmetric distribution about the v_r = 0 axis in the convection zone. A strong magnetic field is more likely to be located in downflow regions. We note that a symmetric distribution is found at r = r_min, which corresponds to 〈v_r, B〉 ∼ 0 in Figure 17(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are two possibilities for the preference of the magnetic field to downflow regions. One is that the magnetic field is generated uniformly in space and transported to the downflow region by converging motion. The other is that the magnetic field is generated in the downflow region. In order to answer this question, we define and evaluate the generation rate of magnetic energy by the stretching (W_str) and the compression (W_cmp) as

$$\begin{eqnarray} && W_\mathrm{str} = \frac{\bf B}{4\pi }\cdot \left[ (\bf B\cdot \nabla) \bf v\right], \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} && W_\mathrm{cmp} = -\frac{B^2}{4\pi }(\nabla \cdot {\bf v}). \end{eqnarray} \tag{ 29 }$$

The estimations for the horizontal averages, 〈W_str〉 and 〈W_cmp〉, are shown in Figure 17(b). It is clear that the value of the stretching 〈W_str〉 is much larger than that of the compression throughout the convection zone and that the generation of the magnetic field is basically done by the stretching in the turbulent motion. This is also seen by the realistic calculations in the photosphere (Pietarila Graham et al. 2010), which suggest that 95% of the gain of the magnetic energy is done by the stretching. We also estimate the correlation between the radial velocity v_r and the energy generation rate by the stretching W_str (Figure 17(c)). The distribution of this correlation is similar to that of 〈v_r, B〉, i.e., the effective stretching prefers the downflow region. Thus we conclude that the reason the strong magnetic field prefers the downflow region is that the magnetic field is more likely to be generated there. The correlation between B and W_str is also estimated. This shows a larger value of (>0.4) throughout the convection zone (Figure 17(d)). The features discussed above are basically common among the studied cases.

In order to investigate magnetic field transport in the convection zone, we evaluate the Poynting flux given by

$$\begin{eqnarray} &&F_\mathrm{m} = \frac{c}{4\pi } (E_\theta B_\phi - E_\phi B_\theta ),\end{eqnarray} \tag{ 30 }$$

where the electric field is defined as E = −(v × B)/c and c is the speed of light (see Brun et al. 2004). Figure 19 shows the horizontally integrated Poynting flux L_m as a function of the depth. Since the absolute value (∼10³¹ erg s⁻¹) is much smaller than the solar luminosity (L_☉ = 3.84 × 10³³ erg s⁻¹), the Poynting flux does not contribute significantly to the total energy flux balance. The Poynting flux is negative in most of the convection zone, since a strong magnetic field is concentrated in downflow regions. This is suggested by previous studies as the turbulent pumping effect (see, e.g., Tobias et al. 1998, 2001). Within the convection zone, the flux has a positive value only in the thin layer (∼0.01 R_☉) close to the bottom. This is caused by the bounced motion from the bottom boundary. The magnetic energy is transported upward, causing the magnetic flux in the upflow region. This can be one of the reasons why the absolute values of the correlation 〈v_r, B〉 is small in the deeper region (Figure 17(a)) and why the distribution of the joint PDF in the bottom (Figure 18(f)) is symmetric.

Zoom In Zoom Out Reset image size

In the following discussion, we investigate the scale of the magnetic field generated. Figure 20 shows the spectra of the kinetic energy (black lines: E_kin = ρ₀v²/2) and the magnetic energy (red lines: E_mag = B²/(8π)) in case M0 averaged from t = 173 days to t = 237 days in which the generation of the magnetic field is close to being saturated. The dashed black and red lines show the kinetic energy without the magnetic field and the magnetic energy at t = 5.8 days, respectively. In the upper convection zone (⩾0.85 R_☉), the spectra of the magnetic energy peaks at the smallest scale, which is typical for the kinematic phase of a local dynamo. Finding this feature in the saturated phase indicates that the local dynamo is likely not very efficient for our resolved scales (>7 Mm). We see also no indication that the kinetic energy spectrum changed due to the presence of the dynamo near the top of the domain. This situation is different in the lower half of the domain. Figures 20(e) and (f) show some peak shift of the magnetic energy to larger scales and some feedback on thermal convection, i.e., the kinetic energy is suppressed on the smallest scales. Here, the magnetic energy slightly exceeds the kinetic energy near the smallest scales. This shows that a more efficient local dynamo can be achieved to some extent for our resolution in the lower part of the convection zone (<0.85 R_☉).

Zoom In Zoom Out Reset image size

Figure 21 shows that the spectra of the horizontal divergence (ζ) and the radial vorticity (ω_r) show a similar distribution to that of the magnetic field with the peak at small scales. Moffatt (1961) suggested that the power spectrum varies as D_b(k)∝k²D_v(k) when the magnetic field is proportional to (∇ · v) or (∇ × v), where D_b(k) and D_v(k) are magnetic and kinetic spectra. This was recently confirmed through solar observations in the photosphere using the Hinode satellite (Katsukawa & Orozco Suárez 2012). Since in our calculation there is similarity between the magnetic field and the vorticity or the divergence, the magnetic energy peaks at smaller scales than does the kinetic energy. We note that differences can be seen between the shear of the velocity (ω_r and ζ) and the magnetic energy (E_mag) at the base of the convection zone (Figures 20(f) and 21(f)), where the feedback from the magnetic field is stronger.

Zoom In Zoom Out Reset image size

Even though the numerical resolution does not vary much with depth, the effectiveness of the local dynamo is not expected to be depth independent. There are two reasons for this. The intrinsic convective scale varies with depth, and a downward Poynting flux exists almost everywhere in the convection zone. In the upper half of the convection zone the divergence of the Poynting flux provides an energy sink, while at the same time the dynamo is not very efficient on the resolved scale (a similar discussion can be found in Stein et al. 2003). In the lower convection zone (<0.85 R_☉), the radial gradient of the Poynting flux is negative, thus the magnetic energy is accumulated. At the same time the intrinsic scale of convection is larger and better resolved, leading to a more efficient dynamo. It should be noted that Vögler & Schüssler (2007) found that even in near-surface areas the small-scale convection reproduced with high resolution has a timescale short enough to amplify the magnetic field against the pumping effect.

In order to confirm our idea about the generation and transportation of the magnetic field in our calculation, we evaluate the effective shear for the magnetic field,

$$\begin{eqnarray} &&f_\mathrm{eff} = \frac{\langle W_\mathrm{str}\rangle }{\langle B^2/(8\pi)\rangle }. \end{eqnarray} \tag{ 31 }$$

f_eff has the unit of s⁻¹ and indicates the timescale of the amplification. Figure 22(a) shows the time evolution of log f_eff and indicates that regardless of the phase of the generation of the magnetic field, the larger value of f_eff is seen in the upper convection, which reflects the short timescale of the thermal convection there. At later times with a stronger magnetic field, the effective shear f_eff is reduced, and after t ∼ 150, it remains constant. Figure 22(b) shows the ratio to the value at t = 5.8 days, i.e., f_eff(t)/f_eff(t = 5.8 days). The suppression of the effective shear depends on the depth. It is more suppressed in the deeper region. This might be caused by the feedback of the magnetic field on the velocity seen in Figure 20 (the solid and dashed black lines) as well as a misalignment between the shear and magnetic field around the base of the convection zone. The ineffective suppression of f_eff in the upper convection zone confirms our idea presented in the previous paragraph. The local dynamo saturates there with little non-linear feedback (suppression of f_eff) since the pumping effect works well in the upper region and the dynamo is not very efficient on the resolved scales to begin with.

Zoom In Zoom Out Reset image size

Figure 23 also shows the variation of dynamo efficiency in the convection zone. Figure 23(a) shows the equipartition magnetic field $B_\mathrm{eq} = \sqrt{4\pi v_\mathrm{rms}^2}$ and the rms magnetic field B_rms as functions of the depth. The solid and dotted lines show the values at the downflow and upflow regions, respectively. Basically, both B_eq and B_rms increase with the depth and, Figure 23(b) shows an increase of B_rms/B_eq with depth, which is caused by the ineffectiveness of the local dynamo in the upper region due to the pumping effect and our insufficient resolution.

Zoom In Zoom Out Reset image size

We carry out the high-resolution calculations of the solar global convection which resolve 10 Mm scale convection smaller than the supergranulation using the RSST (described in Section 2). The RSST leads to a simple algorithm and requires only local communication in the parallel computing. In addition, this method has the capability to access the real solar surface without undermining the important physics. This enables us to capture near-surface small-scale convection while keeping a global domain. Our main conclusions are given as follows. (1) Small-scale convection is excited close to the surface (>0.9 R_☉), when we expand our domain upward to 0.99 R_☉ to capture the near-surface layers with small pressure scale heights. (2) In deeper convection zones (<0.9 R_☉), the convection flow is not influenced by the location of the top boundary and the assumed thickness of the thermal boundary layer. Changing the position of the top boundary or thickness of the cooling layer does not lead to significant differences in the convective structure and properties of the local dynamo, in terms of the power-spectrum, the probability density function, and the local dynamo. (3) Using a large eddy simulation approach, we can achieve small-scale dynamo action and maintain a field of about 0.15 B_eq–0.25 B_eq throughout the convection zone. (4) The overall dynamo efficiency varies significantly in the convection zone as a consequence of the downward directed Poynting flux and the depth variation of the intrinsic convective scales. For a fixed numerical resolution, the dynamo-relevant scales are better resolved in the deeper convection zone and are therefore less affected by numerical diffusivity, i.e., the effective Reynolds numbers are larger.

Recently, Hanasoge et al. (2012) suggested that there is a significant difference between the calculation by Miesch et al. (2008) and flow velocities inferred from local helioseismology. The observed convection amplitude is much smaller than that in the hydrodynamic calculation. We expected that the modification to the location of the upper boundary may help to resolve this issue. Conclusion 1 above, however, suggests that the difference between the calculation and the observation becomes larger when we move the top boundary further upward. We see, nevertheless, some indication that the kinetic energy is reduced on the largest scales when a local dynamo is present. At the same time, we cannot fully rule out the possibility that unresolved small-scale turbulence near the top boundary contributes to this issue, since the intrinsic scale of the convection near our top boundary is still not well resolved.

Conclusion 2 is one of the most important results in this study, since it suggests that previous calculations such as by Miesch et al. (2008) are physically reasonable in the deeper convection zone even if the top boundary condition is placed significantly below the solar surface.

We summarize our simulation results in a schematic picture shown in Figure 24. The conclusion that the dynamo is ineffective near the top comes from the comparison of the kinematic and saturated spectra and the only moderate reduction of the effective shear in the upper convection zone (see Section 4.5). Since several aspects, in particular with regard to the local dynamo, require higher resolution before being directly applicable to the solar convection zone, Figure 24 does not necessarily apply to the real Sun. High-resolution simulations of a local dynamo in the solar photosphere (Vögler & Schüssler 2007) suggest that efficient dynamo action is possible even in the presence of the pumping effect. These simulations, however, use a grid resolution of about a factor of 100 larger than our setup, which is currently inapplicable for global-scale convection simulations. Therefore, the dynamo rms field strength of 0.15 B_eq–0.25 B_eq can likely be considered a lower limit.

Zoom In Zoom Out Reset image size

One issue we cannot address in this study is the problem regarding the small magnetic Prandtl number (Pm ∼ 10⁻³) in the solar convection zone since we adopt numerical diffusivities that assume that the magnetic Prandtl number is around unity. Several authors argued that smaller magnetic Prandtl numbers make local dynamos less efficient (e.g., Schekochihin et al. 2004; Boldyrev & Cattaneo 2004). In other words, a small magnetic Prandtl number requires a larger magnetic Reynolds number in order to achieve a super-critical dynamo. While this can be a significant problem for numerical simulations with rather moderate Reynolds numbers, this is less likely an issue in the solar convection zone with Reynolds numbers as large as Rm ∼ 10¹¹ (Brandenburg 2011). Thus, we believe that our approach relying only on numerical diffusivities can capture the physics of the local dynamo in the solar convection zone.

Rotation is not included in this study, although it plays a key role in organizing large-scale flows and enabling large-scale dynamo action. Studying the near-surface small-scale convection in a setup with rotation is of particular interest to us for investigation of the origin of the near-surface shear layer observed on the Sun. Our investigation of this is a work in progress and will be presented in a forthcoming paper. Further progress can be made by using non-uniform grids, such as nested grids or adaptive mesh refinement, which are useful in order to overcome the significant differences in spatial and temporal dynamic scales. These methods are already implemented in our numerical code (Hotta & Yokoyama 2012).

The authors thank M. Miesch for helpful comments on the draft. H.H. is supported by Grant-in-Aid for JSPS Fellows. The National Center for Atmospheric Research is sponsored by the National Science Foundation. Part of the results was obtained by using the K computer at the RIKEN Advanced Institute for Computational Science (Proposal number hp120287). This research was partly conducted using the Fujitsu PRIMEHPC FX10 System (Oakleaf-FX) in the Information Technology Center, The University of Tokyo. The authors are grateful to Ryoji Matsumoto for managing our computational resources. We have greatly benefited from the proofreading/editing assistance from the GCOE program. This work is partially supported by "Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures" and "High Performance Computing Infrastructure" in Japan. The authors thank the referee for helpful suggestions.

In the numerical simulations of the convection in the near-surface layer, the ordinary tabular equation of state is widely used (Vögler et al. 2005; Rempel et al. 2009). However, this is not a viable approach for our current simulations since the deviations from the reference state are small, e.g., ρ₁/ρ₀ ∼ 10⁻⁶ around the base of the convection zone. Thus, we adopt another way to treat the ionization effect in the near-surface layer. The fluctuations from the reference state are calculated as

$$\begin{eqnarray} && p_1 = \left(\frac{\partial p}{\partial \rho }\right)_s \rho _1 + \left(\frac{\partial p}{\partial s }\right)_\rho s_1, \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} && T_1 = \left(\frac{\partial T}{\partial \rho }\right)_p \rho _1 + \left(\frac{\partial T}{\partial p }\right)_\rho p_1, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} && e_1 = \left(\frac{\partial e}{\partial \rho }\right)_T \rho _1 + \left(\frac{\partial e}{\partial T }\right)_\rho T_1. \end{eqnarray} \tag{ A3 }$$

The first derivatives, such as (∂p/∂ρ)_s, are described by the background variables, ρ₀(r), p₀(r)... and are regarded as functions of the depth r. In the OPAL routine (Rogers et al. 1996), the values (∂e/∂ρ)_T, (∂e/∂T)_ρ, (∂p/log ρ)_T, and (∂p/log T)_ρ are provided for given ρ₀, T₀, and the mass fraction of hydrogen (X), helium (Y), and other metals (Z). We show the relations between the OPAL-provided variable and the required variable derived from the first law of thermodynamics (Mihalas & Mihalas 1984)

$$\begin{eqnarray} && \left(\frac{\partial p}{\partial \rho }\right)_s = \frac{c_\mathrm{p}}{\kappa _\mathrm{t}\rho _0 c_\mathrm{v}}, \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} && \left(\frac{\partial p}{\partial s}\right)_\rho = \frac{\beta _\mathrm{p}T_0}{\kappa _\mathrm{t}c_\mathrm{v}}, \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} && \left(\frac{\partial T}{\partial \rho }\right)_p = -\frac{1}{\rho _0 \beta _\mathrm{p}}, \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} && \left(\frac{\partial T}{\partial p}\right)_\rho = T_0{/}\left(\frac{\partial p}{\partial \log T}\right)_\rho, \end{eqnarray} \tag{ A7 }$$

where β_p, c_v, c_p, and κ_t are the coefficient of thermal expansion, the specific heat at constant volume and pressure, and the coefficient of isothermal compressibility, respectively, defined as

$$\begin{eqnarray} && \beta _\mathrm{p} \equiv -\left(\frac{\partial \log \rho }{\partial T}\right)_p = \frac{1}{T_0}\left(\frac{\partial p}{\partial \log T}\right)_\rho {/}\left(\frac{\partial p}{\partial \log \rho }\right)_T,\nonumber\\&&\quad \end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray} && c_\mathrm{v} = \left(\frac{\partial e}{\partial T}\right)_\rho, \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} && c_\mathrm{p} = c_\mathrm{v} -T_0\beta _\mathrm{p}\left[\left(\frac{\partial e}{\partial \rho }\right)_T -\left(\frac{p}{\rho ^2}\right)\right], \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} &&\kappa _\mathrm{t} = \left(\frac{\partial \log \rho }{\partial p}\right)_T. \end{eqnarray} \tag{ A11 }$$

When the equation of the entropy is given as Equation (5), the total energy is not conserved even mathematically. In this appendix, the deviation in the conservation of the total energy caused by the RSST is introduced. Here, we ignore the magnetic field, the gravity, and the radiation for the sake of simplicity. The equation of continuity for RSST is defined as Equation (2). Thus, there are two relations as

$$\begin{eqnarray} \frac{D\rho }{Dt} &=&\frac{\partial \rho }{\partial t} + {\bf v}\cdot \nabla \rho \nonumber \\ &=&{-}\frac{1}{\xi ^2}\nabla \cdot (\rho {\bf v}) + {\bf v}\cdot \nabla \rho \nonumber \\ &=& \left(1-\frac{1}{\xi ^2}\right)\nabla \cdot (\rho {\bf v})-\rho \nabla \cdot {\bf v} \end{eqnarray} \tag{ B1 }$$

and

$$\begin{eqnarray} \rho \frac{D A}{Dt} &=& \rho \frac{\partial A}{\partial t} + \rho {\bf v}\cdot \nabla A\nonumber \\ &=&\frac{\partial }{\partial t}(\rho A) + \nabla \cdot \left(\rho A {\bf v}\right) -A\frac{\partial \rho }{\partial t} - A\nabla \cdot (\rho {\bf v})\nonumber \\ &=&\frac{\partial }{\partial t}(\rho A) + \nabla \cdot \left(\rho A {\bf v}\right) -A\left(1-\frac{1}{\xi ^2}\right)\nabla \cdot (\rho {\bf v}). \nonumber\\ \end{eqnarray} \tag{ B2 }$$

There is deviation from the relation with the original equation of continuity. The deviations are proportional to (1 − 1/ξ²)∇ · (ρv). Therefore, the deviation between the conservative form and the primitive form of the equation of motion is

$$\begin{eqnarray} &&\left(1-\frac{1}{\xi ^2}\right){\bf v}\nabla \cdot \left(\rho {\bf v}\right). \end{eqnarray} \tag{ B3 }$$

Then, the deviation between the equation of entropy and the conservative form for the total energy is introduced as

$$\begin{eqnarray} \rho T\frac{Ds}{Dt}&=& \rho \frac{De}{Dt}+\rho p\frac{D}{Dt}\left(\frac{1}{\rho }\right) \nonumber \\ &=& \frac{\partial }{\partial t}\left(\rho e\right) + \nabla \cdot \left(\rho e {\bf v}\right) + p\nabla \cdot {\bf v}\nonumber \\ &&- \left(e+\frac{p}{\rho }\right)\left(1-\frac{1}{\xi ^2}\right)\nabla \cdot (\rho {\bf v}) \nonumber \\ &=& Q, \end{eqnarray} \tag{ B4 }$$

where Q includes thermal diffusivity and surface cooling. When we start from the primitive form of the equation of motion, the equation of the kinetic energy is expressed as

$$\begin{eqnarray} \frac{\partial }{\partial t}\left(\frac{1}{2}\rho v^2\right) &+& \nabla \cdot \left(\frac{1}{2}\rho v^2{\bf v}\right) = -{\bf v}\cdot \nabla p\nonumber\\ &+&\frac{1}{2} v^2\left(1-\frac{1}{\xi ^2}\right)\nabla \cdot (\rho {\bf v}). \end{eqnarray} \tag{ B5 }$$

Thus, the conservative form for the total energy is expressed as

$$\begin{eqnarray} \frac{\partial }{\partial t}\left(\rho e +\frac{1}{2}\rho v^2\right) &=& {-}\nabla \cdot \left[\left(\rho e + p+\frac{1}{2}\rho v^2\right){\bf v}\right]\nonumber \\ &&+ \left(e + \frac{p}{\rho } + \frac{1}{2} v^2\right)\nonumber\\ &&\times\left(1-\frac{1}{\xi ^2}\right)\nabla \cdot (\rho {\bf v}) +Q. \end{eqnarray} \tag{ B6 }$$

If the speed of sound is much faster than the fluid velocity, the anelastic approximation ∇ · (ρv) = 0 is well achieved. In our calculations over 200 days, we could not see any coherent trend in the energy.

We derive the first-order radial enthalpy flux. According to Equation (B6), the original form of the radial enthalpy flux is expressed as

$$\begin{eqnarray} &&F_\mathrm{e} = \left(\rho e + p\right)v_r. \end{eqnarray} \tag{ C1 }$$

In the statistical steady state, horizontal integration and the time average ensure no net transport of mass along the depth, i.e.,

$$\begin{eqnarray} &&\int _s \rho v_rdS = 0. \end{eqnarray} \tag{ C2 }$$

In this study, the equation of state for the perfect gas is not adopted, and the form of the enthalpy flux is slightly different from previous studies (e.g., Miesch et al. 2008; Käpylä et al. 2011) in which the enthalpy flux is simply expressed as F_e = c_pρ₀T₁v_r. The integrated radial enthalpy flux is calculated as

$$\begin{eqnarray} L_\mathrm{e} &=& \int _s \left(e + \frac{p}{\rho }\right)\rho v_r dS\nonumber \\ &=& \int _s \left(e + \frac{p_0+p_1}{\rho _0 + \rho _1}\right)\rho v_r\nonumber dS\\ &\sim & \int _s \left[e + \frac{p_0+p_1}{\rho _0} \left(1-\frac{\rho _1}{\rho _0}\right)\right]\rho v_r\nonumber dS \\ &\sim & \int _s \left(e_1 + \frac{p_1}{\rho _0}-\frac{p_0\rho _1}{\rho _0^2}\right)\rho v_r dS\nonumber \\ &\sim & \int _s \left(\rho _0 e_1 + p_1 - \frac{p_0\rho _1}{\rho _0}\right)v_r dS. \end{eqnarray} \tag{ C3 }$$

The final line of the equation is then adopted for the enthalpy flux in this study.