Issue 
A&A
Volume 539, March 2012



Article Number  A30  
Number of page(s)  7  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201118268  
Published online  22 February 2012 
Research Note
Numerical calculation of convection with reduced speed of sound technique
^{1} Department of Earth and Planetary Science, University of Tokyo, 731 Hongo, Bunkyoku, 1130033 Tokyo, Japan
email: hotta.h@eps.s.utokyo.ac.jp
^{2} High Altitude Observatory, National Center for Atmospheric Research, Boulder, CO, USA
Received: 13 October 2011
Accepted: 31 December 2011
Context. The anelastic approximation is often adopted in numerical calculations with low Mach numbers, such as those including stellar internal convection. This approximation requires socalled frequent global communication, because of an elliptic partial differential equation. Frequent global communication is, however, negative factor for the parallel computing performed with a large number of CPUs.
Aims. We test the validity of a method that artificially reduces the speed of sound for the compressible fluid equations in the context of stellar internal convection. This reduction in the speed of sound leads to longer time steps despite the low Mach number, while the numerical scheme remains fully explicit and the mathematical system is hyperbolic, thus does not require frequent global communication.
Methods. Two and threedimensional compressible hydrodynamic equations are solved numerically. Some statistical quantities of solutions computed with different effective Mach numbers (owing to the reduction in the speed of sound) are compared to test the validity of our approach.
Results. Numerical simulations with artificially reduced speed of sound are a valid approach as long as the effective Mach number (based on the lower speed of sound) remains less than 0.7.
Key words: methods: numerical / stars: interiors / convection / Sun: interior
© ESO, 2012
1. Introduction
Turbulent thermal convection in the solar convection zone plays a key role for the maintenance of large scale flows (differential rotation, meridional flow) and solar magnetic activity. The angular momentum transport of convection maintains the global mean flows. Global flows relate to the generation of the global magnetic field, i.e. the solar dynamo. Differential rotation bends the preexisting poloidal field and generates the strong toroidal field (the Ω effect), and mean meridional flow transports the magnetic flux equatorward at the base of the convection zone (Choudhuri et al. 1995; Dikpati & Charbonneau 1999). The internal structures of the solar differential rotation and the meridional flow are revealed by the helioseismology (see review by Thompson et al. 2003). Some mean field studies have reproduced these global flows (Kichatinov & Rüdiger 1993; Küker & Stix 2001; Rempel 2005; Hotta & Yokoyama 2011). These studies, however, used some kinds of parameterization of the turbulent convection, i.e. turbulent viscosity and turbulent angular momentum transport. Thus, a selfconsistent thorough understanding of global structure requires the detailed investigation of the turbulent thermal convection. Turbulent convection is also important when describing the magnetic field itself. The strength of the next solar maximum in the predictions of the mean field model depends significantly on the turbulent diffusivity (Dikpati & Gilman 2006; Choudhuri et al. 2007; Yeates et al. 2008). In addition, the turbulent diffusion has an important influence on the parity of solar global field, the strength of polar field and so on (Hotta & Yokoyama 2010a,b).
There have been numerous LES numerical simulations of the solar and stellar convection (Gilman 1977; Gilman & Miller 1981; Glatzmaier 1984; Miesch et al. 2000, 2006; Brown et al. 2008) and their magnetic fields (Gilman & Miller 1981; Brun et al. 2004; Brown et al. 2010, 2011). In these studies, the anelastic approximation is adopted to avoid any of the difficulties caused by the high speed of sound. At the base of the convection zone, the speed of sound is about 200 km s^{1}. In contrast, the speed of convection is thought to be 50 m s^{1} (Stix 2004). The time step must therefore be shorter owing to the CFL condition in an explicit fully compressible method, even when we are interested in phenomena related to convection. In the anelastic approximation, the equation of continuity is treated as (1)where ρ_{0} is the stratified background density and v denotes the velocity. The anelastic approximation assumes that the speed of sound is essentially infinite, resulting in an instantaneous adjustment of pressure to flow changes. This is achieved by solving an elliptic equation for the pressure, which filters out the propagation of sound waves. As a result, the time step is only limited by the much lower flow velocity. However, owing to the existence of an elliptic the anelastic approximation has a weak point. When performing the numerical calculation using parallel computing, we must employ the frequent global communication. At the present time, the efficiency of scaling in parallel computing is saturated with about 2000–3000 CPUs in solar global simulations with the pseudospectral method (Miesch, priv. comm.). Higher resolution, however, is needed to understand the precise mechanism of the angular momentum and energy transport by the turbulent convection and the behavior of magnetic field especially in thin magnetic flux tube.
In this paper, we test the validity of a different approach to circumvent the severe numerical time step constraints in low Mach number flows. We use a method in which the speed of sound is reduced artificially by transforming the equation of continuity to (see, e.g., Rempel 2005) (2)where t denotes the time. Using this equation, the effective speed of sound becomes ξ times smaller, but the dispersion relationship for sound waves remains otherwise unchanged (wave speed decreases equally for all wavelengths). Since this technique does not change the hyperbolic character of the underlying equations, the numerical treatment can remain fully explicit, thus does not require global communication in parallel computing. We refer to the technique as the Reduced Speed of Sound Technique (RSST). This technique was used previously by Rempel (2005, 2006) in mean field models of solar differential rotation and nonkinematic dynamos, which essentially solve the full set of timedependent axisymmetric MHD equations. Those solutions were however restricted to the relaxation toward a stationary state or very slowly varying problems on the timescale of the solar cycle. Here we apply this approach to thermal convection, where the intrinsic timescales are substantially shorter. To this end, we study two and threedimensional (3D) convection, in particular the latter will be nonstationary and turbulent.
The detailed setting of test calculation is given in Sect. 2. The results of our calculations are given in Sect. 3. We summarize our paper and give discussion of the RSST in Sect. 4.
2. Model
In this section, we describe the detailed settting of this study.
2.1. Equations
The two or threedimensional equation of continuity, equation of motion, equation of energy, equation of state, are solved in Cartesian coordinate (x,z) or (x,y,z), where x and y denote the horizontal direction and z denotes the vertical direction. The basic assumptions underlying this study are as:

1.
the timeindependent hydrostatic reference state;

2.
that the perturbations caused by thermal convection are small, i.e. ρ_{1} ≪ ρ_{0} and p_{1} ≪ p_{0}, here ρ_{0} and p_{0} denote the reference state values, whereas ρ_{1} and p_{1} are the fluctuations of density and pressure, respectively. Thus a linearized equation of state is used as Eq. (6);

3.
that the profile of the reference entropy s_{0}(z) is a steady state solution of the thermal diffusion equation ∇·(Kρ_{0}T_{0}∇s_{0}) = 0 with constant K.
The formulations are almost identical to those of Fan et al. (2003), and given by where T_{0}(z), and s_{0}(z) denote the reference temperature and entropy, respectively, and e_{z} denotes the unit vector along the zdirection, γ is the ratio of specific heats, with the value for an ideal gas being γ = 5/3, and s_{1} denotes the fluctuation of entropy from reference atmosphere. We note that the entropy is normalized by a value of the specific heat capacity at constant volume c_{v}. The quantity g is the gravitational acceleration, which is assumed to be constant. The quantity Π denotes the viscous stress tensor (7)and ν and K denote the viscosity and thermal diffusivity, respectively, and ν and K are assumed to be constant throughout the simulation domain. We assume for the reference atmosphere a weakly superadiabatically stratified polytrope where ρ_{r}, p_{r}, T_{r}, H_{r}, and δ_{r} denote the values of ρ_{0}, p_{0}, T_{0}, H_{0} (the pressure scale height), and δ (the nondimensional superadiabaticity) at the bottom boundary z = 0. Since δ ≪ 1, the value m is nearly equal to the adiabatic value, meaning that m = 1/(γ − 1). The strength of the diffusive parameter ν and K is expressed in terms of the nondimensional parameters: the Reynolds number Re ≡ v_{c}H_{r}/ν and the Prandtl number Pr ≡ ν/K, where the velocity unit v_{c} ≡ (8δ_{r}gH_{r})^{1/2}. We note that in this paper the unit of time is H_{r}/v_{c}. In all calculations, we set Pr = 1.
Parameters for numerical simulations.
2.2. Boundary conditions and numerical method
We solve Eqs. (3)–(6) numerically. At the horizontal boundaries (x = 0, L_{x} and y = 0, L_{y}), periodic boundary conditions are adopted for all variables. At the top and the bottom boundaries, impenetrative and stressfree boundary conditions are adopted for the velocities and the entropy is fixed to At both the top and bottom boundaries (z = 0 and z = L_{z}), we set p_{1} in the ghost cells such that the right hand side of the zcomponent of Eq. (4) is zero at the boundary (which is between ghost cells and domain cells), where the ghost cells are the cells beyond the physical boundary.
We adopt the fourthorder spacecentered difference for each derivative. The first spatial derivatives of quantity q is given by (18)where i denotes the index of the grid position along a particular spatial direction. The numerical solution of the system is advanced in time with an explicit fourthorder RungeKutta scheme. The system of partial equations can be written as (19)for U_{n + 1}, which is the value at t_{n + 1} = (n + 1)Δt is calculated in four steps The longest allowed time step, Δt_{max}, is determined by the CFL criterion. When both advection and diffusion terms are included in calculation, the time step reads (24)where (25)and c_{tot} is the total wave speed (26)where the effective speed of sound is expressed as (27)The time step determined by the diffusion term is (28)and c_{ad} and c_{v} are safety factors of order unity.
Using δ_{r} = 1 × 10^{4}, the original speed of sound is about 35v_{c} at the bottom and 12v_{c} at the top boundary, respectively. In all calculations, Δx ~ 2.3 × 10^{2}H_{r}. Thus, if we use c_{ad} = 1 and c_{v} = 1, Δt_{ad} = 6.4 × 10^{4} H_{r}/v_{c} and Δt_{v} = 1.5 × 10^{1} H_{r}/v_{c}. For all the ξ values considered in this paper, the time step remains restricted by the (reduced) speed of sound, thus the calculation is about ξ times more efficient with RSST.
3. Results
3.1. Twodimensional study
We carried out four twodimensional (2D) calculations with the RSST and one calculation without approximation (case 1–5). The values of the free parameters are given in Table 1. The superadiabaticity (δ) is 1 × 10^{6} at the bottom and about 2 × 10^{5} at the top boundary. It is almost the same superadiabaticity value as at the base of the solar convection zone. Figure 1 shows the timedevelopment of entropy. In the beginning, nonlinear timedependent convection occurred (top panel), which changed to a steady state at later times (bottom panel). For the steady state, we compared the root mean squared (rms) velocity with different ξ. The rms velocity was defined as (29)Figure 2 shows our results where the value of ξ is from 1 to 80. The effective Mach number is defined as . If we have larger ξ value than 80, we cannot reach a stationary state, since there are some shocks generated by supersonic convection. The discussion of unsteady convection is given in the next session of 3D calculations. Even though the Mach number reaches 0.6 using ξ = 80 (panel a), the horizontal and vertical rms velocities are almost identical to those calculated with ξ = 1 (without RSST). The ratio of the rms velocities for each ξ and ξ = 1 are shown in Fig. 2e. The difference is always a few percent. This result is unsurprising since in the stationary state the equation of continuity becomes (30)whose solution must not depend on the value of ξ. We note that the cell size is to some degree affected by the aspect ratio of the domain. This does not affect our conclusions since this influence is the same for all considered values of ξ. To confirm the robustness of our conclusions, we repeated this experiment with a wider domain (26.16H_{r} × 2.18H_{r} instead of 8.72H_{r} × 2.18H_{r}) in which we find four steady convection cells with about 10% different rms velocities. In addition we find a dependence on ξ that is similar to that shown in Fig. 2, i.e. the solutions show differences of only a few percent as long as ξ < 80. In this section, we confirm that the RSST is valid for the 2D stationary convection when the effective Mach number is less than 1, i.e. ξ < 80 with δ_{r} = 1 × 10^{6}. If we use δ_{r} = 1 × 10^{4} (result is not shown), the criterion becomes ξ < 8 for the 2D calculation.
Fig. 1 Timedevelopment of entropy in a 2D calculation with parameters of case 1. (Top panel) t = 20. (Bottom panel) t = 400. 
Fig. 2 Some quantities in 2D calculation. a) Maximum Mach number of each time step. b) Distribution of Mach number estimated with the rms velocity. c) Distribution of horizontal rms velocity v_{h}. d) Distribution of vertical rms velocity. e) Ratio of the rms velocities for each ξ and ξ = 1. 
As a next step, we investigate the dependence of the linear growth rate on ξ during the initial (timedependent) relaxation phase toward the final stationary state. Figure 3 shows the linear growth of maximum perturbation density ρ_{1} with different ξ. Black and red lines show the results with δ_{r} = 1 × 10^{6} and δ_{r} = 1 × 10^{4}, respectively. These calculation parameters are not given in Table 1. In the calculation with δ = 1 × 10^{6}, the growth rate decreases for values of ξ > 300, whereas it occurs at ξ > 30 in the calculation with δ_{r} = 1 × 10^{4}. The reason can be explained as follows. In the convective instability, upflow (downflow) generates a positive (negative) entropy perturbation and then a negative (positive) density perturbation is generated by the sound wave. If the speed of sound is fairly slow, the generation mechanism of density perturbation is ineffective. In the calculation with δ_{r} = 1 × 10^{6} (δ_{r} = 1 × 10^{4}), the growth rate with ξ = 200 (ξ = 20), however, is almost same as that with ξ = 1, even though flow with ξ = 200 (ξ = 20) is expected to be the Mach number Ma = 1.3, i.e. supersonic convection flow in the saturated state.
Fig. 3 Behavior in linear phase. Dependence of growth rate on ξ is shown. Black and red lines show the results with δ_{r} = 1 × 10^{6} and 1 × 10^{4}, respectively. 
3.2. Threedimensional study
We now investigate the suitability of the RSST for 3D unsteady thermal convections (case 6–13). The value of superadiabaticity at the bottom boundary is 1 × 10^{4} and at the top 2 × 10^{3}. Although this value is relatively large compared to solar value, the expected speed of convection is much smaller than the speed of sound, hence is small enough to allow us to investigate the validity of the RSST. The entropy of 3D convections with ξ = 1, 20, and 80 at t = 100H_{r}/v_{c} are shown in Fig. 4. The convection is completely unsteady and turbulent (animation is provided). The appearance of convection with ξ = 20 is almost the same as that with ξ = 1. We verify this using the Fourier transformation and autodetection technique. However, the appearance of convection with ξ = 80 (bottom panel) differs completely from the others. This difference is most clearly visible in the animation of Fig. 4 that is provided with the online version.
Fig. 4 Snapshots of entropy of the 3D convection. Top, middle, and bottom panels correspond to ξ = 1, 10, and 80, respectively. Left (right) panels show entropy at top (bottom) boundary. (Animation is provided, and the difference between ξ = 1 and 80 is most clearly visible in the animation of Fig. 4, which is provided with the online version.) 
Fig. 5 Some quantities averaged in time between t = 100 and 200H_{r}/v_{c} in 3D calculations. a) Distribution of Mach number estimated with the rms velocity. b) Distribution of vertical rms velocity v_{h}. c) Distribution of horizontal rms velocity. d) Distribution of rms power density of pressure and buoyancy. e) Distribution of rms power density of inertia. 
We estimate the rms velocities for different ξ by calculating the average of values between t = 100 to 200H_{r}/v_{c} (Fig. 5: panel b and c). Without the RSST, i.e. ξ = 1 and using δ_{r} = 1 × 10^{4}, the Mach number is 1 × 10^{2} at the bottom and 4 × 10^{2} at the top boundary. The rms velocities at ξ = 40 and 80 differ from those at the ξ = 1 by less than 15% and 30%, respectively. When we adopt ξ = 40 and 80, the Mach number estimated by rms velocity exceeds unity (Fig. 5: panel a), i.e. supersonic convection. This supersonic downflow frequently generates shocks and positive entropy perturbation, thus downflow is decelerated. This is the reason why the rms velocities with large ξ(=40, 80) are small. When ξ = 5, 10, 15, and 20, however, the rms velocities are in good agreement with that for ξ = 1. We estimate the rms power density of pressure, buoyancy, and inertia (Fig. 5). This results display almost the same trend as for the rms velocities. The power profiles ξ = 5, 10, 15, 20 agree with each other and those with ξ = 40, 80 disagree with those with ξ = 1. These results show that the RSST is a valid technique at least for ξ = 20, at which the Mach number is around 0.7. We note that the convection pattern among our 2D and 3D cases varies substantially. As a consequence, the ξ values for which the validity of RSST breaks down are also different in the 2D and 3D setups.
Fig. 6 Comparison of horizontal velocity spectra with different ξ. a) z = 2.18H_{r} b) z = 1.09H_{r} c) z = 0. Velocities are averaged in time between t = 100 and 200H_{r}/v_{c}. 
Fig. 7 Detection of convective cell. a) Original contour of density perturbation at z = 2.18H_{r}. b) Distribution of detected convection cell. Color and label (#n) show each convective cell. 
In Fig. 6, we compare the average spectral amplitudes for different values of ξ. There is no significant difference between the spectral amplitudes for different values of ξ in the range from 1 to 20.
We investigate the distribution of cell sizes at the top boundary. This value is strongly related to the turbulent diffusivity and transport of either angular momentum or energy. The method for detecting the cell is as follows. At the boundary of the cell i.e., the region of downflow, the perturbation of density is positive and has a high value. When the density in a region exceeds a threshold, the region is regarded as the boundary of a convective cell. When a region is surrounded by one continuous boundary region, the region is defined as one convective cell. Figure 7 shows the detected cells. Each color and each label (#n) correspond to each detected cell. We estimate the size of all cells and compare the distribution of cell sizes for different ξ. The results are shown in Fig. 8. The cell size distribution follows a power law from about 0.01 to 10 and there is no dependence on ξ in the range from 1 to 20. Although when using this type of technique the size of cells tends to be large when neglecting smaller cells, our conclusion is not wrong, since all autodetections are equally affected.
Fig. 8 Distribution of convective cell size with different ξ is shown. The cell size distribution is averaged in time between t = 100 and 200H_{r}/v_{c}. 
For the above two investigations, i.e. the Fourier analysis and study of cell detection, we conclude that the statistical features are unaffected by the RSST method as long as the effective Mach number (computed with the reduced speed of sound) does not exceed values of about 0.7, which corresponds to ξ = 20 in our setup.
To confirm our criterion that the RSST is valid if the Mach number is smaller than 0.7, we conduct calculations with larger superadiabaticity, i.e. δ_{r} = 1 × 10^{3} (cases 14–18). If our criterion is valid, required ξ with larger superadiabaticity must decrease. The results are shown in Fig. 9. Using δ_{r} = 1 × 10^{3}, the calculations with ξ = 10, 15, and 20 generate a supersonic convection flow near the surface (Fig. 9a). It is clear that the results for ξ = 10, 15, and 20 differ from those with ξ = 1 and 5. The calculation with ξ = 5 corresponds to a flow whose Mach number is around 0.6. Thus, our criterion is not violated even for different values of the superadiabaticity.
Fig. 10 Dependence of [∇·(ρ_{0}v)] _{rms} on ξ. In case 13, the value of ξ is 20 at the bottom and 4.5 at the top boundary. 
Owing to our altered equation of continuity the primitive and conservative formulation of Eqs. (3) to (5) are no longer equivalent. For example, the equation of motion in conservative form is expressed as (31)where F denotes the pressure gradient, gravity, and Lorentz force. With some transformations, we can obtain (32)If the equation of continuity is satisfied, i.e. (∂ρ/∂t = −∇·(ρv)), the primitive form is obtained as (33)However, using the RSST these two forms are no longer equivalent. We used the primitive formulation at the expense that energy and momentum are not strictly conserved; however, the consistency of our results for different values of ξ strongly indicates that this is not a serious problem for the setup we considered. Alternatively, we could also implement our modified equation of continuity in a conservative formulation. This would ensure that density, momentum, and energy are strictly conserved at the expense of a modified set of primitive equations. Figure 10 shows the dependence of [∇·(ρ_{0}v)] _{rms} (=[ξ^{2}∂ρ_{1}/∂t] _{rms}) on ξ. Using ξ = 5, this term is almost the same as that with ξ = 1. Although the deviation becomes large as ξ increases, it is not proportional to ξ^{2}.
In the previous discussion, we kept ξ constant in the entire computational domain. In the solar convection zone, the Mach number varies however substantially with depth, from ~1 in the photosphere to <10^{7} at the base of the convection zone, the speed of sound itself varies from about 7 km s^{1} in the photosphere to 200 km s^{1} at the base of the convection zone. A reduction in the speed of sound is therefore most significant in the deep convection zone, but not in the near surface layers. This could be achieved with a depth dependent ξ. Even if we use a conservative form of the equation of continuity, (34)the result must not be same as the result without any approximation of the statistical steady state. When the value averaged during a statistical steady state is expressed as ⟨ a ⟩ , where a is a physical value, the equation of continuity becomes (35)with inhomogeneous ξ. The solution of Eq. (35) is different from the solution of the original equation of continuity in the statistical steady state, i.e. 0 = ∇·(ρ_{0} ⟨ v ⟩ ). Thus, the statistical features such as rms velocity are not reproduced with an inhomogeneous ξ. If we use the nonconservative form of the equation of motion given in Eq. (2), this problem does not occur. Thus, we investigate the nonuniformity of ξ in case 13 using the nonconservative form, i.e., Eq. (3). To keep the Mach number uniform for all heights, we use ξ = 20/(δ(z)/δ_{r})^{1/2}, i.e. ξ = 20 at the bottom and ξ = 4.5 at the top boundary. We note that the ratio of the speed of convection to the speed of sound is roughly estimated as . In this setting, the value of is 1 × 10^{2} at the bottom and 4 × 10^{2} at the top boundary. The same analysis as that for uniform ξ is done (see Figs. 5, 6, and 10). There is no significant difference between ξ = 20/(δ(z)/δ_{r})^{1/2} and ξ = 1. Although mass is not conserved locally for a nonhomogeneous ξ using a nonconservative form of the equation of continuity, the horizontally averaged vertical mass flux is approximately zero in the statistically steady convection and the conservation total mass is not broken significantly. We conclude that an inhomogeneous ξ is valid under the previously obtained condition, i.e. the Mach number is less than 0.7.
4. Summary and discussion
We have applied RSST (described in Sect. 1) to 2D and 3D simulations of low Mach number thermal convection and confirmed the validity of this approach as long as the effective Mach number (computed with the reduced speed of sound) stays below 0.7 everywhere in the domain. The overall gain in computing efficiency that can be achieved therefore depends on the maximum Mach number that was present in the setup of the problem.
Since the Mach number is estimated to be 10^{4} in the base of the solar convection zone (Stix 2004), a several thousand times longer time step can be taken using the RSST. Therefore, the RSST and parallel computing with a large number of CPUs will make it possible to calculate largescale solar convection with high resolution in the near future.
Compared to the anelastic approximation, there are three major advantages in RSST:

1.
It can be easily implemented into any fully compressible code(regardless of the numerical scheme or grid structure) since itonly requires a minor change in the equation of continuity andleaves the hyperbolic structure of the equations unchanged.

2.
Owing to the explicit treatment, it does not require any additional communication overhead, which makes it suitable for massive parallel computations.

3.
The base of the convection zone and the surface of the sun where the anelastic approximation is broken can be connected, using a spacedependent ξ.
Overall we find that RSST is a very useful technique for studying low Mach number flows in stellar convection zones as it substantially alleviates any stringent timestep constraints without adding any computational overhead.
Acknowledgments
The authors want to thank Dr. M. Miesch for his helpful comments. Numerical computations were in part carried out on Cray XT4 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan. We would like to acknowledge highperformance computing support provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. The National Center for Astmospheric Reseach (NCAR) is sponsored by the National Science Foundation. This work was supported by the JSPS Institutional Program for Young Researcher Overseas Visits and the Research Fellowship from the JSPS for Young Scientists. We have greatly benefited from the proofreading/editing assistance from the GCOE program. We thank the referee for helpful suggestions that helped to improve our paper.
References
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2008, ApJ, 689, 1354 [CrossRef] [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
 Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri, A. R., Schussler, M., & Dikpati, M. 1995, A&A, 303, L29 [NASA ADS] [Google Scholar]
 Choudhuri, A. R., Chatterjee, P., & Jiang, J. 2007, Phys. Rev. Lett., 98, 131103 [Google Scholar]
 Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 [NASA ADS] [CrossRef] [Google Scholar]
 Dikpati, M., & Gilman, P. A. 2006, ApJ, 649, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., Abbett, W. P., & Fisher, G. H. 2003, ApJ, 582, 1206 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1977, Geophys. Astrophys. Fluid Dynam., 8, 93 [Google Scholar]
 Gilman, P. A., & Miller, J. 1981, ApJS, 46, 211 [Google Scholar]
 Glatzmaier, G. A. 1984, J. Comp. Phys., 55, 461 [Google Scholar]
 Hotta, H., & Yokoyama, T. 2010a, ApJ, 709, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., & Yokoyama, T. 2010b, ApJ, 714, L308 [NASA ADS] [CrossRef] [Google Scholar]
 Hotta, H., & Yokoyama, T. 2011, ApJ, 740, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Kichatinov, L. L., & Rüdiger, G. 1993, A&A, 276, 96 [NASA ADS] [Google Scholar]
 Küker, M., & Stix, M. 2001, A&A, 366, 668 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miesch, M. S., Elliott, J. R., Toomre, J., et al. 2000, ApJ, 532, 593 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Brun, A. S., & Toomre, J. 2006, ApJ, 641, 618 [Google Scholar]
 Rempel, M. 2005, ApJ, 622, 1320 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2006, ApJ, 647, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Stix, M. 2004, The Sun: an introduction, 2nd edn., Astronomy and astrophysics library (Berlin: Springer), ISBN: 3540207414 [Google Scholar]
 Thompson, M. J., ChristensenDalsgaard, J., Miesch, M. S., & Toomre, J. 2003, ARA&A, 41, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., Nandy, D., & Mackay, D. H. 2008, ApJ, 673, 544 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Timedevelopment of entropy in a 2D calculation with parameters of case 1. (Top panel) t = 20. (Bottom panel) t = 400. 

In the text 
Fig. 2 Some quantities in 2D calculation. a) Maximum Mach number of each time step. b) Distribution of Mach number estimated with the rms velocity. c) Distribution of horizontal rms velocity v_{h}. d) Distribution of vertical rms velocity. e) Ratio of the rms velocities for each ξ and ξ = 1. 

In the text 
Fig. 3 Behavior in linear phase. Dependence of growth rate on ξ is shown. Black and red lines show the results with δ_{r} = 1 × 10^{6} and 1 × 10^{4}, respectively. 

In the text 
Fig. 4 Snapshots of entropy of the 3D convection. Top, middle, and bottom panels correspond to ξ = 1, 10, and 80, respectively. Left (right) panels show entropy at top (bottom) boundary. (Animation is provided, and the difference between ξ = 1 and 80 is most clearly visible in the animation of Fig. 4, which is provided with the online version.) 

In the text 
Fig. 5 Some quantities averaged in time between t = 100 and 200H_{r}/v_{c} in 3D calculations. a) Distribution of Mach number estimated with the rms velocity. b) Distribution of vertical rms velocity v_{h}. c) Distribution of horizontal rms velocity. d) Distribution of rms power density of pressure and buoyancy. e) Distribution of rms power density of inertia. 

In the text 
Fig. 6 Comparison of horizontal velocity spectra with different ξ. a) z = 2.18H_{r} b) z = 1.09H_{r} c) z = 0. Velocities are averaged in time between t = 100 and 200H_{r}/v_{c}. 

In the text 
Fig. 7 Detection of convective cell. a) Original contour of density perturbation at z = 2.18H_{r}. b) Distribution of detected convection cell. Color and label (#n) show each convective cell. 

In the text 
Fig. 8 Distribution of convective cell size with different ξ is shown. The cell size distribution is averaged in time between t = 100 and 200H_{r}/v_{c}. 

In the text 
Fig. 9 The results with δ_{r} = 1 × 10^{3}. The format is the same as Fig. 5. 

In the text 
Fig. 10 Dependence of [∇·(ρ_{0}v)] _{rms} on ξ. In case 13, the value of ξ is 20 at the bottom and 4.5 at the top boundary. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.