## Abstract

The solar convection zone is filled with turbulent convection in highly stratified plasma. Several theoretical and observational studies suggest that the numerical calculations overestimate the convection velocity. Since all deep convection zone calculations exclude the solar surface due to substantial temporal and spatial scale separations, the solar surface, which drives the thermal convection with efficient radiative cooling, has been thought to be the key to solve this discrepancy. Thanks to the recent development in massive supercomputers, we are successful in performing the comprehensive calculation covering the whole solar convection zone. We compare the results with and without the solar surface in the local domain and without the surface in the full sphere. The calculations do not include the rotation and the magnetic field. The surface region has an unexpectedly weak influence on the deep convection zone. We find that just including the solar surface cannot solve the problem.

## INTRODUCTION

The solar convection zone (CZ), the outer 30% of the solar interior, is filled with turbulent convection. Because of substantial contrast in density and temperature, the temporal and spatial scales of the convection vary drastically. For example, the temporal and spatial scales at the base of the CZ are 1 month and 100,000 km, respectively, whereas at the surface, they are minutes and 1000 km. It is difficult to cover these scale separations in numerical simulations, and there have been no comprehensive calculations covering the whole CZ.

An important requirement for such a comprehensive calculation is an understanding of the thermal convection itself. Several investigations indicated recently that convection in numerical calculations is much faster than in reality. For example, the polar region tends to rotate faster, called “antisolar” differential rotation, in high-resolution deep convection calculations (*1*, *2*), whereas the equator region rotates faster, called “solar-like” differential rotation, in the real Sun. This is certainly caused by the high convection velocity in the numerical calculations. To achieve the solar-like differential rotation with the solar rotation rate in numerical calculations, there are currently two options for numerical calculations. One is reducing the luminosity (*2*), and the other is adopting large viscosity and/or large thermal conductivity (*1*). The purpose of these is to suppress the convection velocity. While it is known that the magnetic field has a role in suppressing the convection velocity (*1*), even with the magnetic field, we still need large thermal conductivity to achieve the solar-like differential rotation. The magnetic field itself cannot solve the problem. We see discussions about antisolar and solar-like differential rotations in a lot of literatures (*3*–*5*). Helioseismology and supergranulation studies give similar arguments, i.e., the observational constraint on the convection velocity is smaller than that reproduced in the numerical calculations (*6*, *7*). One possible solution to this issue is the inclusion of the solar surface. For the deep convection, an artificial top boundary is set around 0.98*R*_{⊙} (*8*–*10*), where *R*_{⊙} = 696 Mm is the solar radius. The near-surface layer is calculated separately. The bottom boundary of the surface calculations is typically located at a depth of 6 to 20 Mm, i.e., 0.99 − 0.97*R*_{⊙ }(*11*, *12*). The granulation properties in the calculations match the observations well. As the calculation with the deepest bottom boundary is located at a depth of 49 Mm, i.e., 0.93*R*_{⊙} (*7*), our understanding of the influence of the surface layer on the deep CZ is limited. The radiative cooling at the surface occurs on a short time scale (~10 s), and the cooling layer is thin (~100 km). Considering this fact, Spruit (*13*) argued that it is probable that the deep convection simulations with the artificial top boundary, which is settled at a deep layer (~0.98*R*_{⊙}), overestimate the driving in the deep layer. Because the low-entropy plumes generated with the efficient radiative cooling descend to the deep layer, the deep layer is not necessarily superadiabatic. In the deep convection calculation, artificial cooling forces the stratification to be superadiabatic and drives a large-scale flow into the deep layer (*14*). There is no calculation that includes the photosphere, but there are several numerical challenges to mimicking the existence of the photosphere and its radiative cooling using an open top boundary condition with a cool plume (*15*) and enforcing a thin and strong cooling layer (*16*). These references also report that the existence of strong near-surface cooling leads to a less superadiabatic stratification in the deep layer and suppression of the convection velocity at large scales. To address the influence of the photospheric cooling in a self-consistent manner, a calculation that covers the whole CZ with realistic radiative cooling is indispensable.

## RESULTS

Thanks to the recent development of supercomputers, we can carry out a comprehensive calculation of the whole CZ with a new code R2D2 [Radiation and RSST (reduced speed of sound technique) for Deep Dynamics]. We show three cases, which we call Local/Photosphere (LP), Local/Deep (LD), and Spherical/Deep (SD). The case LP is a comprehensive calculation including the photosphere in a local Cartesian box. The vertical domain of the calculation box is from the base of the CZ, *x*/*R*_{⊙} = 0.71, to about 700 km above the τ = 1 surface (*x*/*R*_{⊙} = 1.001), where *x* denotes the vertical direction. Around the surface, we solve the radiative energy transfer equation. We use a nonuniform vertical grid spacing, which is 48 km around the top boundary and about 1000 km around the base of the CZ. We adopt a combination of the linear and full equations of state using the OPAL repository (*17*). Around the base of the CZ, we adopt the diffusion approximation for the radiative energy transfer because of the large optical depth there. For the case LD, the top boundary is at *x*/*R*_{⊙} = 0.992, and the artificial cooling layer is added around the top boundary. The other conditions, grid spacing, and horizontal domain size are the same as in the case LP. We prepare 512 and 384 grid points for the vertical direction for the cases LP and LD, respectively. The horizontal extent of the cases LP and LD is 196.608 Mm, and the grid spacing is 192 km. Similar resolution for the photosphere is adopted in previous studies of the solar surface (*7*). In the comparison between the cases LP and LD, we can investigate the influence of the photosphere with realistic radiative cooling on the deep layer. Our method for including the solar surface is validated with the result of the numerical code RAMENS (Radiation Magnetohydrodynamics Extensive Numerical Solver) (*18*), which is consistent with the solar observations. The calculation box is limited to enable the simulation to run for a sufficiently long time. The effect of this limitation is investigated with the case SD. The case SD covers the whole sphere of the CZ with the top boundary at *x*/*R*_{⊙} = 0.99 and the artificial cooling layer. For all cases, we adopt the stratification from Model S (*19*). The top and bottom boundaries are closed for the vertical velocity and stress free for the horizontal velocities. The symmetric boundary condition is adopted for the entropy perturbation at both boundaries. Only for the case LP, the positive vertical velocity is symmetric, and the negative vertical velocity is closed at the top boundary. This top boundary condition is adopted in the previous photospheric calculations (*12*). To carry out this comprehensive calculation, one of the important breakthroughs is the RSST. In this study, we limit the maximum speed of sound to 20 km s^{−1}. This allows us to take a large time step and reduce the number of temporal integrations (*20*). The magnetic field and the rotation are not included in this study. For the cases LP and LD, we calculate the convection for 75 days (1800 hours). For the analysis, we average the values between 800 and 1800 hours. For the case SD, we calculate 456 days and average the values between 216 and 456 days. The detail of the calculation setting is shown in Materials and Methods.

Figure 1 shows the three-dimensional (3D) structure of the thermal convection for the cases LP (Fig. 1A) and SD (Fig. 1B). In Fig. 1A, the perturbation of the entropy (*s*′) normalized with the temporal average of the root mean square (RMS) entropy (*s*_{rms}) is shown, where *s*′ is the deviation from the horizontal and temporal average. A hierarchical pattern is seen in the entropy perturbation. The small-scale plumes around the surface merge repeatedly and construct a large-scale convection in the deep interior. Figure 1B shows the vertical velocity at 0.98*R*_{⊙}. The black rectangle in Fig. 1B indicates the numerical domain for the cases LP and LD.

Figure 2 shows the emergent intensity and horizontal slices of the vertical velocity (*v*_{x}) for the case LP. The emergent intensity is shown in Fig. 2A. The 1 Mm scale granulation is well reproduced in the intensity image. A similar granulation pattern is observed in the vertical velocity at the τ = 1 surface (Fig. 2B). In the deep region, convection patterns at the larger scales of 10 Mm (0.99*R*_{⊙} (Fig. 2C) and 30 Mm (0.98*R*_{⊙} (Fig. 2D) are seen. By comparing the black rectangle in Figs. 1B and 2D, it is found that the overall convection structures are similar between the cases despite the different settings. In deeper regions (Fig. 2, E and F), the cool plumes merge further and, last, a box-scale (~200 Mm) convection cell is observed.

Figure 3 shows comparisons of the cases LP and LD. The influence of the surface radiative cooling on the deep CZ is investigated. Figure 3 (A and B) shows the energy fluxes in the cases LP and LD, respectively. The overall structure in the deep layer (<0.98*R*_{⊙}) is the same. Compared with the thick artificial cooling layer in the case LD, the case LP has a very thin (~100 km) cooling layer (green lines). Both cases show an almost constant total energy flux. This indicates that the energy transport is statistically steady over the period for averaging. The convection velocity increases because of a substantial reduction in the near-surface density. The behavior of the energy flux shows an additional effect. The realistic radiative cooling causes a cool plume with a large amplitude [], which descends to the deep layer and leads to the asymmetry between the upflow and downflow and the large downward kinetic flux in the near-surface layer (>0.98*R*_{⊙}; blue lines). To compensate for the negative kinetic flux, the enthalpy flux (red lines) needs to be increased. When the entropy perturbation is increased with a fixed energy flux, the convection velocity tends to decrease. In this case, however, the increase in the amplitude of the negative kinetic flux increases both the enthalpy flux and the convection velocity in the near-surface layer. The RMS velocities for the cases LP (red) and LD (black) are shown in Fig. 3C. Although the convection velocity in the near-surface layer (>0.98*R*_{⊙}) is increased as explained, in the deep layer, the convection velocities are almost the same between the cases. We note that the convection velocity is higher than the other solar global calculations (*21*). This is caused by the absence of the rotation. When the rotation is included, the convection is bent, and the convection RMS velocity is suppressed (*3*). Figure 3D shows the superadiabaticity in the cases LP and LD. In the near-surface layer, the superadiabaticity in the LP case is larger than that in the case LD. In the deep layer, the superadiabaticity values also coincide for both cases. This is not expected in (*13*–*16*), where the superadiabaticity is thought to decrease when surface radiative cooling is included.

Figure 4 shows the kinetic energy spectra in the cases LP (red lines), LD (black lines), and SD (blue lines). Figure 4 (A and B) shows the kinetic energy spectra for the vertical flow at *x*/*R*_{⊙} = 0.96 and 0.72, respectively. Figure 4 (C and D) shows the horizontal energy spectra in the same manner. We adopt the definition of the energy spectra shown in (*22*). In Fig. 4 (B and D), the corrected energy spectra considering the geometry are also plotted with a dashed line for the case LP. Despite the different geometry and domain in the calculations, the kinetic energy spectra in the overlapping range match well (*ℓ* ~ 20 − 1000). This comparison indicates that the previous global calculations without the photosphere would not change substantially even with the photosphere included. In addition, the result of the case LP would not change with increasing horizontal domain size because the larger domain simulation (SD) shows a consistent result. The green line in the Fig. 4C shows the upper limit of the convection velocity constrained by the time-distance helioseismology (*6*). While the helioseismic constraint itself is under debate (*23*), this is an indication that all our calculations still overestimate the convection velocity.

To evaluate the influence of the photosphere on the deep CZ, we estimate the mixing length parameter. We adopt the method suggested in (*24*), where the mixing length *L*_{ml} is calculated as(1)where *F*_{down} is the horizontally averaged mass flux with negative vertical velocity. We normalize the mixing length *L*_{ml} with the local pressure scale height *H*_{p} = (*d* log *p*/*dx*)^{−1}. This is the mixing length parameter, i.e., α_{ml} = *L*_{ml}/*H*_{p}. This distribution of the mixing length parameter is shown in Fig. 5. We note that the mixing length parameter increases toward the bottom and diverges around *x* = 0.75*R*_{⊙}. In addition, negative superadiabaticity, i.e., the subadiabatic layer, is observed around the base of the CZ, while the enthalpy flux is positive (Fig. 3D). This indicates that the mixing length description is not appropriate around the base of the CZ. We note that the negative superadiabaticity is a common feature in convection calculation even with the radiation zone and the overshoot layer (*25*, *26*).

## DISCUSSION

The reason that the solar surface has a very weak influence on the deep convection structure, i.e., including the photosphere does not change the deep convection, is that efficient mixing around the near surface is achieved. This leads to a situation where the convection properties are mostly determined by the temperature gradient in nearby layers. This argument is supported by the estimation of the mixing length parameter. The parameter in Fig. 5 shows an order of unity (~1.8) around the surface, which indicates efficient mixing. In particular, in the near-surface layer, the pressure scale height is very short (e.g., 2 Mm = 0.003*R*_{⊙} at *x*/*R*_{⊙} = 0.99), and the influence of the photosphere disappears within a short distance.

The case SD is almost same as the result of our calculation without the rotation (*27*). We also carried out the calculation with the solar rotation in the previous study (*2*). It is noted that just adding the solar rotation in the calculation in (*27*) shows the antisolar differential rotation, and the luminosity is decreased by a factor of 18 to suppress convection velocity and achieve the solar-like differential rotation. The energy spectra in cases LP and LD are consistent with that in case SD. Considering this fact, we expect that even when we include the photosphere for the global calculation, the thermal convection constructs the antisolar differential rotation. Our state-of-the-art calculation reveals that the inclusion of the photosphere does not resolve the problem of the high convection velocity in the deep solar interior at the currently affordable numerical resolution.

While the case LP has a relatively low resolution, i.e., 192-km horizontal grid spacing, for resolving the convection in the photosphere, the validity of the calculation is checked with the high-resolution calculation, i.e., 24-km horizontal grid spacing, using RAMENS as shown in fig. S1. The convection velocity and overall structure are mainly determined by the imposed luminosity and the stratification, and the resolution has little influence on our investigated range. These results match well with other calculations that are consistent with the solar photospheric observation. While the smallest eddy in the real solar photosphere is estimated to be 0.01 cm with a Reynolds number of 10^{13} (*28*) and the calculation covering all the scales is unfeasible, we have not found any substantial influence of the low resolution on the overall convection structure at the photosphere. We note that this statement does not deny the possibility of unknown influence of the unresolved small-scale turbulence.

In this study, the rotation and the magnetic field are not included. The possible influence of these is discussed here. The rotation tends to bend the convective flow and makes the typical convection scale small. This indicates that the rotation itself would decrease the influence of the solar surface. There is a possibility that the magnetic field has a key role in solving the problem. The magnetic field tends to suppress the small-scale convection via the Lorentz force. This decreases the mixing between the upflow and downflow and effectively increases the mixing length (*29*), and this possibly increases the influence of the solar surface. Calculation for achieving this situation is challenging since increasing the mixing length by the Lorentz force requires an efficient small-scale dynamo. The efficient small-scale dynamo at the solar surface requires a much higher resolution, i.e., grid spacing of 8 km (*12*). In addition, saturated small-scale dynamo in the deep CZ requires a longer calculation time of several hundred days (*29*).

## MATERIALS AND METHODS

### Equations

In this section, we explain the model settings for cases LP and LD. The setting for the case SD is identical to that of the case H0 in (*27*), but the number of the grid is 512(*r*) × 768(θ) × 2304(φ) × 2 in the Yin-Yang grid. Especially for the case LP, we developed the code R2D2. We solved the hydrodynamic equations in 3D Cartesian geometry (*x*, *y*, *z*), where *x* denotes the vertical direction and the others are the horizontal directions. The equations are expressed as(2)(3)(4)where ρ, **v**, *s*, *T*, *p*, and ξ denote the density, velocity, specific entropy, temperature, pressure, and reduction factor for the speed of sound, respectively. We solved the entropy equation instead of the total energy equation because we needed to deal with a substantially small entropy gradient in the deep CZ. We did not use the linear assumption (e.g., ρ_{1} << ρ_{0}) for the perturbation because it is not valid around the photosphere. The equations were solved with the fourth-order space-centered derivative and the four-step Runge-Kutta method (*30*). We also adopted a slope-limited artificial diffusion to all variables to stabilize the calculation (*12*). In this calculation, we did not have any explicit viscosity and thermal conductivity, and the calculations were solely stabilized by the artificial viscosity. Since we adopted the same scheme for the velocity and the entropy, the effective Prandtl number was expected to be around unity. We used a uniform horizontal grid spacing of 192 km with 1024 grid points for each horizontal direction (*y* and *z*). For the vertical direction (*x*), we used a nonuniform grid spacing with 512 grid points. The distribution of the vertical grid points is shown in Fig. 6A. To calculate the radiative cooling and heating (*Q*_{rad}), we adopted the diffusion approximation around the base of the CZ, where the optical depth is large(5)where is the radiative diffusivity calculated from the OPAL opacity (*17*). We made a modification considering the geometric effect, whereby our radiative diffusivity can be expressed as a function of the original radiative diffusivity *κ*_{r} as . We switched off the diffusion approximation at *x*/*R*_{⊙} = 0.995, where the radiative energy transport is not important (see the green line in Fig. 3A). Around the solar surface, in which the optical depth becomes small and the diffusion approximation is no longer valid, we solved the gray radiation transfer(6)where *I*, τ, and *S* are the intensity, optical depth, and source function, respectively. After calculating the intensity *I*, the radiative cooling rate was calculated in two ways, depending on the optical depth. In the optically thick region(7)(8)was used. In the optically thin region(9)(10)was used, where κ, μ, and ω are the opacity, the unit vector of the direction of the radiative intensity, and the surface element of the unit sphere, respectively. Then, the cooling rate around the surface is(11)where τ_{0} = 0.1. We adopted the assumption of a local thermodynamic equilibrium (LTE), where the source function is expressed with the Planck function. The LTE assumption is widely accepted in the CZ and the photosphere because of its high density, while non-LTE should be adopted in the chromosphere. The validity of the LTE assumption has been confirmed with the comparison between the numerical calculations and the observations (*31*). To reduce the numerical cost, we solved the radiative transfer only along the vertical direction in each horizontal grid under the assumption of a plane parallel atmosphere. Thus, the energy transport is only vertically upward and downward. Typical radiation codes solve 9 to 17 rays (*32*). We started solving the radiative transfer from the layer 3.072 Mm below to the top boundary, where the optical depth is already large and the intensity is solely determined by the local temperature. Near the photosphere, the horizontal energy exchange by the radiation causes diffusion of the subgranular-scale structure. This process is complemented by the diffusion approximation for the horizontal direction when τ > 5. The numerical diffusion of the hydrodynamic solver, referred to as slope-limited diffusion above, partially contributes to this process. When the resolution is high, i.e., grid spacing is smaller than 24 km in our numerical scheme, our method of the radiative transfer without horizontal direction enhances small-scale thermal features. The influence of the assumption on the radiative transfer, however, is not substantial with the coarse horizontal grid (192 km) in this paper. The above approximations for the treatment near the solar surface were validated by a comparison with the standard photospheric convection codes in the Supplementary Materials.

The perturbation, ρ_{1}/ρ_{0}, *p*_{1}/*p*_{0}, and *s*_{1}/*s*_{0}, changes by five orders of magnitude between the surface () and the deep CZ (). To calculate the pressure using the density and entropy, we switched between the linear and full equations of state to deal with this substantial difference. In the deep CZ, the linear equation of state is(12)where (∂*p*/∂ρ)_{s} and (∂*p*/∂*s*)_{ρ} are regarded as functions of *x*. In the near-surface region, we used the table equation of state and calculated the perturbation of the pressure.(13)

Then, the perturbation of the pressure was determined in response to the amplitude of the perturbation of the density and entropy as(14)(15)(16)where *c*_{t0} = −2 and *dc*_{t} = 0.1.

We used the RSST to deal with the fast speed of sound and relax the CFL (Courant-Friedrichs-Lewy) condition (*20*). We adopted an inhomogeneous reduction rate of the speed of sound ξ. The distribution of ξ is shown in Fig. 6B. The resulting reduced speed of sound is shown in Fig. 6C. ξ does not depend on time. In the case LP (red line), ξ is defined as(17)where(18)(19)(20)

The original speed of sound *c*_{s} was calculated from the time-independent background stratification. We set ξ_{0} = 120, *x*_{RSST} = *x*_{max} −5.7 × 10^{8} cm, and *d*_{RSST} = 5 × 10^{7} cm, where *x*_{min} and *x*_{max} are the locations of the bottom and top boundaries, respectively. This expression allows the speed of sound around the photosphere to be unchanged, i.e., ξ = 1. In the case LD (black line), the reduced speed of sound was set to be uniform . ξ is defined as(21)

The choice of ξ was determined by the effective Mach number, i.e., the RMS velocity (Fig. 3C) divided by the reduced speed of sound. Figure 6D shows the effective Mach number. When the RSST was used, the effective Mach number was kept smaller than 0.2, which can maintain the low Mach number situation. In addition, a too small Mach number causes notably diffusive features in the calculations, since the artificial viscosity is proportional to the characteristic velocity in our numerical code. Thus, the choice is a compromise of the validity and the numerical diffusivity. We note that the original speed of sound can be calculated as the product of Fig. 6 (B and C), which is 225 km s^{−1} around the base of the CZ.

In the case LD, we added an artificial cooling layer around the top boundary. This is included in the radiative cooling as(22)(23)where *F*_{0} = 6.32 × 10^{10} erg s^{− 1} cm^{− 2} is the solar radiative flux in the photosphere. We adopted the thickness of the cooling layer *d*_{art} = 3.74 × 10^{8} cm.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/1/eaau2307/DC1

Validation of the near-surface calculation

Energy flux

Fig. S1. Comparison of the results in this study with the results using RAMENS (black, case LP; red, using RAMENS).

Movie S1. 3D structure of the entropy perturbation.

Movie S2. Emergent intensity of thermal convection.

Movie S3. Scan movie of the vertical velocity.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We are grateful to M. Rempel, P. Käpylä, and R. Matsumoto for the insightful comments on the manuscript. We thank S. Hanasoge for providing the data on the helioseismic constraint. The results were obtained using the K computer at the RIKEN Advanced Institute for Computational Science (proposal nos. hp180042, hp170239, hp170012, hp160026, hp160252, and ra000008).

**Funding:**This work was supported by MEXT/JSPS KAKENHI grant nos. JP18H04436, JP16K17655, JP16H01169, and JP15H05816. This research was supported by MEXT as “Exploratory Challenge on Post-K Computer” (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System) and the Project for Solar-Terrestrial Environment Prediction.

**Author contributions:**H.H. carried out the numerical calculations and analyses. H.H. and H.I. developed the numerical scheme and code for the calculations. K.K. supervised the whole project. All the authors participated in the project discussion and manuscript preparation.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in this paper are present in the paper and/or Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).