Generation of Solar-like Differential Rotation

In this subsection, the overall convection and magnetic field are discussed. Figures 2–7 show the overall structure of the radial velocity and the radial magnetic field. Figures 2, 4, and 6 show the radial velocity v_r at r = 0.95, 0.9, and 0.85R_⊙, respectively. Figures 3, 5, and 7 show the radial magnetic field B_r at r = 0.95, 0.9, and 0.85R_⊙, respectively. The results from Low (panels (a), (d)), Middle (panels (b), (e)), and High (panels (c), (f)) cases are shown in these figures. The panels (d), (e), and (f) show zoomed views indicated by a white dashed box in panel (a). The radial velocity at r = 0.95R_⊙ (Figure 2) shows a typical convection pattern, i.e., thin concentrated downflows surrounded by broad upflows. Two effects cause this pattern. The first effect is stratification. Because the solar convection zone is gravitationally stratified, the upper layer has a lower gas pressure. A rising fluid parcel expands because of the stratification, while the descending parcel contracts. This asymmetry of the upflows and downflows causes the typical convection pattern. In addition, we should see a boundary effect at this depth. A wall exists at r = 0.96R_⊙ where the radial motion stops. This process leads to diverging and converging motions in the upflows and downflows, respectively. The convection patterns in the Low case (Figure 2(a) and (d)) are similar to previous calculations (Miesch et al. 2008), i.e., the smallest scale is the downflow lane. In the High case, we can see smaller-scale structures even in the downflow lanes (Figure 2(f)). The banana-cell, the north–south aligned convection cell, cannot be seen in all the cases at this depth.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the radial magnetic field B_r at r = 0.95R_⊙. The magnetic field strength increases from the Low case to the High case. This tendency is also seen in the other depth (see Figures 5 and 7). In all the cases, the radial magnetic field is swept up to the downflow region. This concentration is also seen in the previous calculation in the deep interior (Brun et al. 2004) and the photosphere (Vögler et al. 2005). While the previous simulations and the Low case in this study typically show a sheet-like magnetic flux aligned to the downflow lane, we can occasionally observe a blob-shaped magnetic flux (a notable one is indicated by the dashed orange circle in Figure 3(f)). This structure shows significantly superequipartition magnetic field strength and low gas pressure. The convection is suppressed in this region. This structure is important for magnetic field generation (see discussion at Section 3.6). At r = 0.9R_⊙, the convection shows a larger-scale pattern. The small-scale convection around the top boundary is merged to construct the larger scale while increasing the pressure/density scale height in the deep region (see Stein & Nordlund 1998; Lord et al. 2014). The banana-cell-like feature begins to appear in this depth. In the deeper layer (r = 0.85 R_⊙, middle of the convection zone), the flow pattern shows the banana-cell-like features more clearly than the upper layers (Figure 6). In the mixing length theory, the convection velocity v_c scales as ${\rho }_{0}{v}_{{\rm{c}}}^{3}\sim {L}_{\odot }/4\pi {r}^{2}$ (Biermann 1948), where L_⊙ is the solar luminosity. This dependence indicates that the convection velocity decreases in the deeper layers with the larger density ρ₀. The convection timescale is τ ∼ H_p /v_c. These relations mean that the convection timescale increases in the deeper layers by increasing the pressure scale height and decreasing the convection velocity. As a result, the convection tends to obey the rotation influence (Coriolis force) and shows the banana-cell in the deep layers. The magnetic field distribution is chaotic at this depth. The strong magnetic field tends to be located at the downflow plume, but the coincidence between the downflow and the strong magnetic field is worse than the upper layers.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this subsection, we discuss statistical properties of the convection and magnetic fields. Here, we define the statistical values of a quantity Q, the longitudinal average 〈Q〉, the longitudinal rms ${Q}_{(\mathrm{RMS})}^{{\prime} }$, and the latitudinally averaged longitudinal rms Q_(RMS) as follows.

$$\begin{eqnarray}&&\langle Q\rangle (r,\theta )=\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{Qd}\phi \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{Q}_{(\mathrm{RMS})}^{{\prime} }(r,\theta )=\sqrt{\displaystyle \frac{1}{2\pi }{\int }_{0}^{2\pi }{\left(Q-\langle Q\rangle \right)}^{2}d\phi }\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{Q}_{(\mathrm{RMS})}(r)=\sqrt{\displaystyle \frac{1}{2}{\int }_{0}^{\pi }Q{{\prime} }_{(\mathrm{RMS})}^{2}\sin \theta d\theta }\end{eqnarray} \tag{ 11 }$$

We note that we define spherical average $\widetilde{Q}$ and spherical rms Q_(rms) in Section 3.7 differently from the current definition.

Figure 8 shows the longitudinal rms velocity (panel (a)) and the ratio of the magnetic energy E_mag(RMS) to the kinetic energy E_kin(RMS), where the energies are defined as follows:

$$\begin{eqnarray}&&{E}_{\mathrm{kin}(\mathrm{RMS})}=\displaystyle \frac{1}{2}{\rho }_{0}{v}_{(\mathrm{RMS})}^{2}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{E}_{\mathrm{mag}(\mathrm{RMS})}=\displaystyle \frac{{B}_{(\mathrm{RMS})}^{2}}{8\pi }.\end{eqnarray} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

Figure 8(a) shows that the convection velocity is suppressed in the higher resolution. This is a general tendency of the high-resolution simulations (e.g., Hotta et al. 2015b). The horizontal velocity (dashed lines) is more suppressed than the radial velocity (solid lines). A key to understanding convection suppression is the magnetic field. Figure 8(b) shows the ratio of the magnetic energy to the kinetic energy. While in the Low case (blue line), the magnetic energy is smaller than the kinetic energy throughout the convection zone, and we observe a superequipartition magnetic field in the High case (green line). The ratio exceeds 2.5 at maximum. This result indicates the strong influence of the magnetic field on the convective flow. The reason why the High case has such a strong magnetic field is discussed in Section 3.6, and the suppression mechanism of the convection velocity by the magnetic field is shown in Section 3.7.

Figure 9 shows the kinetic and magnetic energy spectra. We show the results at the layers at r = 0.73 (panel (a)), 0.85 (panel (b)), and 0.9R_⊙ (panel (c)). The definition of the spherical harmonic expansion for an arbitrary quantity Q(θ, ϕ) is

$$\begin{eqnarray}&&{\widehat{Q}}_{m}^{{\ell }}={\int }_{0}^{\pi }{\int }_{0}^{2\pi }Q(\theta ,\phi ){Y}_{m}^{{\ell }}(\theta ,\phi )d\phi d\theta ,\end{eqnarray} \tag{ 14 }$$

where ${Y}_{m}^{{\ell }}$, ℓ, m are the spherical harmonics, the spherical harmonic degree, and the spherical harmonic order, respectively. We use a normalization that satisfies

$$\begin{eqnarray}&&{Q}_{(\mathrm{RMS})}^{2}=\sum _{{\ell }}\sum _{m,m\ne 0}\displaystyle \frac{{\widehat{Q}}_{m}^{{\ell }}{\widehat{Q}}_{m}^{{\ell }* }}{r},\end{eqnarray} \tag{ 15 }$$

where * denotes the complex conjugate. The kinetic ${\widehat{E}}_{\mathrm{kin}}({\ell })$ and magnetic ${\widehat{E}}_{\mathrm{mag}}({\ell })$ energy spectra are calculated as

$$\begin{eqnarray}&&{\widehat{E}}_{\mathrm{kin}}({\ell })=\displaystyle \frac{1}{2}{\rho }_{0}\sum _{m}\widehat{{\boldsymbol{v}}}\cdot {\widehat{{\boldsymbol{v}}}}^{* }\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\widehat{E}}_{\mathrm{mag}}({\ell })=\displaystyle \frac{1}{8\pi }\sum _{m}\widehat{{\boldsymbol{B}}}\cdot {\widehat{{\boldsymbol{B}}}}^{* }.\end{eqnarray} \tag{ 17 }$$

The tendency of the kinetic energy spectra is almost the same among the different layers. While the large-scale (ℓ < 10) energy does not change from the Low to Middle cases (blue and green lines, respectively), the energy is significantly reduced in the High case (green line). This reduction is one of the main topics in this paper. The relation between the kinetic and magnetic energies depends on the resolution. When the magnetic energy surpasses the kinetic energy, we expect an efficient small-scale dynamo (e.g., Rempel 2014). In the Low case, while the magnetic energy exceeds the kinetic energy in the deep layer (r = 0.73R_⊙, panel (a)) in the small scale (ℓ ∼ 40), this clear excess cannot be seen at the shallower layer (r = 0.9R_⊙, panel (c)). The inefficient small-scale dynamo in a shallower layer is a common feature in the global dynamo calculation (e.g., Hotta et al. 2014). Because the shallower layer has a smaller energy injection scale of the convection because of a small pressure/density scale height and a short timescale for downward magnetic energy transport, we need a high resolution to resolve the small-scale dynamo (see discussion by Stein & Nordlund 2002; Vögler & Schüssler 2007). This difficulty of the small-scale dynamo is solved in the Middle case (orange line). While the turnover scale depends on the layer depth, the excess of the magnetic field in the small scales is achieved in all the layers in the Middle case. The situation drastically changes in the High case (green line). The kinetic energy is reduced in all the scales but especially in the large scale (ℓ < 30). This significant suppression is seen at all depths. As a result, the magnetic energy exceeds or is comparable to the kinetic energy in all scales.

Zoom In Zoom Out Reset image size

Figure 10 shows the differential rotation 〈Ω〉/2π and the meridional flow 〈 v _m〉 = 〈v_r 〉 e _r + 〈v_θ 〉 e _θ . The angular velocity is defined as Ω = Ω₀ + Ω₁ and ${{\rm{\Omega }}}_{1}={v}_{\phi }/(r\sin \theta )$. While the Low case shows the fast pole (panel (a)), we reproduce the fast equator in the High case as shown in HK21. The reason why we have the fast equator in the high-resolution calculation is discussed in Section 3.9. Also, the differential rotation succeeds in avoiding the Taylor–Proudman constraint, i.e., ∂Ω/∂z ≠ 0, where z is the direction of the rotational axis. This topic is discussed in Section 3.8. The meridional flow structure also depends on the resolution (Figure 10(d), (e), (f)). In the Low case, an anticlockwise flow is dominant, and we can see a clear poleward flow around the surface. We can observe a tiny clockwise cell around the base of the convection zone. In the Middle case, the meridional flow is separated around the tangential cylinder of the base of the convection zone. Anticlockwise and clockwise flow cells are seen in low and high latitudes, respectively. In the High case, a clockwise meridional flow is dominant throughout the convection zone. The poleward flow around the base of the convection zone is an essential feature for the fast equator (see Section 3.9). The poleward meridional flow around the surface becomes weak in Middle and High cases. Note that we can recover clear poleward meridional flow when the top boundary is closer to the real solar surface (Hotta et al. 2015a). In a high-resolution calculation, we have already checked this tendency and will introduce it in a future publication (H. Hotta, K. Kusano & T. Sekii, 2022, in preparation).

Zoom In Zoom Out Reset image size

In this subsection, we evaluate several nondimensional parameters for comparisons with the previous studies. Since we do not use any explicit diffusivities (viscosity ν, magnetic diffusivity η, and thermal conductivity κ), the effective diffusivities, ν_eff, η_eff, and κ_eff need to be evaluated. The evaluation procedure is shown in Appendix E. We evaluate the effective viscosity from the kinetic energy spectra. Since we use the same numerical scheme for the magnetic field and the entropy as the velocity, the Prandtl number Pr = ν_eff/κ_eff, and the magnetic Prandtl number Pm = ν_eff/η_eff are assumed to be unity. The mean rms velocity ${\overline{v}}_{\mathrm{RMS}}$ is defined as

$$\begin{eqnarray}&&{\overline{v}}_{\mathrm{RMS}}={\int }_{{r}_{\min }}^{{r}_{\max }}{v}_{(\mathrm{RMS})}{r}^{2}{dr}/{\int }_{{r}_{\min }}^{{r}_{\max }}{r}^{2}{dr}.\end{eqnarray} \tag{ 18 }$$

The nondimensional numbers are defined as follows:

• 1.  
  Reynolds number $(\mathrm{Re})$ (Featherstone & Miesch 2015).

  $$\begin{eqnarray}&&\mathrm{Re}\,=\,\displaystyle \frac{{\overline{v}}_{\mathrm{RMS}}d}{{\nu }_{\mathrm{eff}}}\end{eqnarray} \tag{ 19 }$$

  where $d={r}_{\max }-{r}_{\min }$ is the radial extent of the computational domain.
• 2.  
  Rayleigh number (Ra) (Gastine et al. 2014).

  $$\begin{eqnarray}&&\mathrm{Ra}=\displaystyle \frac{{g}_{0}{d}^{3}{\rm{\Delta }}s}{{c}_{{\rm{p}}0}{\nu }_{\mathrm{eff}}{\kappa }_{\mathrm{eff}}}\end{eqnarray} \tag{ 20 }$$

  where g₀ and c_p0 are the gravitational acceleration and the heat capacity at constant pressure at the outer boundary $r={r}_{\max }$. ${\rm{\Delta }}s=\max ({\tilde{s}}_{1})-\min ({\tilde{s}}_{1})$, where ${\tilde{s}}_{1}$ is the horizontally averaged entropy perturbation.
• 3.  
  Flux Rayleigh number (Ra_F ) (Featherstone & Hindman2016b).

  $$\begin{eqnarray}&&{\mathrm{Ra}}_{F}=\displaystyle \frac{{g}_{0}{F}_{0}{d}^{4}}{{c}_{{\rm{p}}0}{\rho }_{{\rm{t}}}{T}_{{\rm{t}}}{\nu }_{\mathrm{eff}}{\kappa }_{\mathrm{eff}}}\end{eqnarray} \tag{ 21 }$$

  where F₀, ρ_t , and T_t are the energy flux density, the background density ρ₀, and the background temperature T₀ at $r={r}_{\max }$, respectively.
• 4.  
  Ekman number (Ek) (Gastine et al. 2014).

  $$\begin{eqnarray}&&\mathrm{Ek}=\displaystyle \frac{{\nu }_{\mathrm{eff}}}{{{\rm{\Omega }}}_{0}{d}^{2}}\end{eqnarray} \tag{ 22 }$$
• 5.  
  Convective Rossby number (Ro_c ) (Gastine et al. 2014).

  $$\begin{eqnarray}&&{\mathrm{Ro}}_{c}=\sqrt{\displaystyle \frac{{\mathrm{RaEk}}^{2}}{\Pr }}\end{eqnarray} \tag{ 23 }$$
• 6.  
  Rossby number (Ro) (Featherstone & Miesch 2015).

  $$\begin{eqnarray}&&\mathrm{Ro}=\displaystyle \frac{{\overline{v}}_{\mathrm{RMS}}}{2{{\rm{\Omega }}}_{0}d}\end{eqnarray} \tag{ 24 }$$
• 7.  
  Local Rossby number (Ro_ℓ ) (Christensen & Aubert 2006).

  $$\begin{eqnarray}&&{\mathrm{Ro}}_{{\ell }}=\displaystyle \frac{{\overline{{\ell }}}_{{\rm{u}}}{\overline{v}}_{\mathrm{RMS}}}{\pi {{\rm{\Omega }}}_{0}d},\end{eqnarray} \tag{ 25 }$$

  where ${\overline{{\ell }}}_{{\rm{u}}}$ is the mean spherical harmonic degree defined as

  $$\begin{eqnarray}&&{\overline{{\ell }}}_{{\rm{u}}}={\int }_{0}^{{{\ell }}_{\max }}{\ell }{\widehat{E}}_{{\rm{h}}}({\ell })d{\ell }/{\int }_{0}^{{{\ell }}_{\max }}{\widehat{E}}_{{\rm{h}}}({\ell })d{\ell }.\end{eqnarray} \tag{ 26 }$$

  ${\widehat{E}}_{{\rm{h}}}$ is the kinetic energy spectra of the horizontal velocity. We evaluate it at r = 0.83R_⊙.

All these values are summarized in Table 1. Thanks to a large number of grid points and the slope-limited artificial viscosity, the effective diffusivities are significantly reduced, and the Reynolds and Rayleigh numbers reach huge values. Convective and ordinary Rossby numbers Ro_c and Ro are typical values for the solar simulations (Featherstone & Miesch 2015; Mabuchi et al. 2015). These are values for the transition between the fast equator and the fast pole. In these values, the effect of small-scale turbulence realized by the high resolution is not considered. The local Rossby number using the kinetic energy spectra is a way to consider the small-scale turbulence (Christensen & Aubert 2006). Gastine et al. (2014) suggest that the transition between the fast equator and the fast pole occurs around Ro_ℓ ∼ 1. Our local Rossby numbers are much larger than the critical value due to the small-scale turbulence. This result indicates that only with the angular momentum transport by the Reynolds stress, i.e., the turbulence, we cannot maintain the solar-like differential rotation. This issue is again discussed in Section 4.

In this subsection, we discuss the generation mechanism of the magnetic field, especially the superequiparition magnetic field achieved in the High case (Figure 8), i.e., the magnetic energy is larger than the kinetic energy. We analyze the magnetic energy equation to investigate the mechanism. The equation is as follows.

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{{B}^{2}}{8\pi }\right) & = & \mathop{\underbrace{-\displaystyle \frac{{\boldsymbol{B}}}{4\pi }\cdot \left[\left({\boldsymbol{v}}\cdot {\rm{\nabla }}\right){\boldsymbol{B}}\right]}}\limits_{{T}_{{\rm{m}}(\mathrm{ADV})}}\\ & & \,\mathop{\underbrace{+\displaystyle \frac{{\boldsymbol{B}}}{4\pi }\cdot \left[\left({\boldsymbol{B}}\cdot {\rm{\nabla }}\right){\boldsymbol{v}}\right]}}\limits_{{T}_{{\rm{m}}(\mathrm{STR})}}\\ & & \mathop{\underbrace{-\displaystyle \frac{{B}^{2}}{8\pi }\left({\rm{\nabla }}\cdot {\boldsymbol{v}}\right)}}\limits_{{T}_{{\rm{m}}(\mathrm{CMP})}}\end{array}\end{eqnarray} \tag{ 27 }$$

There are three contributions to change the magnetic energy, which are advection T_m(ADV), stretching ${T}_{{\rm{m}}(\mathrm{STR})}$, and compression T_m(CMP). Figure 11(a) shows the spherically averaged terms in Equation (27). The solid lines indicate the contribution from the advection T_m(ADV). Around the top boundary, the strong magnetic field is concentrated in the downflow region, and the magnetic energy is transported downward. As a result, the advection contribution is negative and positive in the near-surface layer and the deep convection zone, respectively. Because the higher resolution shows a stronger magnetic field, this effect increases with increased resolution. Next, the dashed lines are the contribution by the stretching term ${T}_{{\rm{m}}(\mathrm{STR})}$. In most of the convection zone, the amplitude decreases in the higher resolutions. Because the magnetic field increases, the Lorentz feedback is amplified, and the production rate of the magnetic field decreases. Figure 11(b) shows the probability density function (PDF) between the radial velocity v_r and the stretching term ${T}_{{\rm{m}}(\mathrm{STR})}$. The result indicates that the main contribution of the stretching occurs at the downflows (v_r < 0). While the net contribution of the stretching is positive, we can see a significant negative contribution (energy transfer from magnetic to kinetic energies) in the downflow region. The dependence of the stretching ${T}_{{\rm{m}}(\mathrm{STR})}$ on the resolution (Figure 11(a)) indicates that the stretching is not responsible for the superequipartition magnetic fields in the High case. Finally, we discuss the compression term T_m(CMP) shown with the dotted lines in Figure 11(a). The amplitude of the compression term monotonically increases with increased resolution. Also, Figure 11(c) shows the PDF between the radial velocity v_r and the compression term T_m(CMP). Similar to the stretching, the important compression occurs at the downflow region, but the negative contribution of T_m(CMP) is not significant. The reason why the fluid can overcome and compress the strong magnetic field to amplify the field strength is shown in Figure 12. PDFs between the magnetic pressure (energy) and perturbation gas pressure are shown. The perturbation gas pressure is defined as

$$\begin{eqnarray}&&p{{\prime} }_{1}={p}_{1}-\langle {p}_{1}\rangle ,\end{eqnarray} \tag{ 28 }$$

and is the deviation from the longitudinal average. We define the gas pressure inside and outside the strong magnetic field as _pi and p_e, respectively. The perturbation gas pressure defined here is approximated as $p{{\prime} }_{1}\sim {p}_{{\rm{i}}}-{p}_{{\rm{e}}}$. The black dashed line in Figure 12 indicates $p{{\prime} }_{1}=-{B}^{2}/8\pi$, i.e., the magnetic pressure is balanced with the gas pressure. In the Low case (Figure 12(a)), the PDF distributes rather uniformly. Even the small magnetic energy ( ∼10⁸ dyn cm⁻²) has large perturbation gas pressure ( < −4 × 10⁷ dyn cm⁻²). In the High case (Figure 12), the magnetic field strength is amplified, and most of the strong magnetic field distributes on the $p{{\prime} }_{1}=-{B}^{2}/8\pi$ line. Also, the region with a weak magnetic field and the low gas pressure disappears. These results indicate that the gas pressure maintains the superequipartition magnetic field realized in the High case. Because the solar interior is in a low Mach number situation, the internal energy is huge compared with the kinetic and magnetic energies. If the magnetic field is balanced with the gas pressure (internal energy), the magnetic field strength can be a superequipartition to the kinetic energy. The region with the weak magnetic field and the low gas pressure disappears in the High case, indicating that the dynamo is efficient enough to amplify all the small-scale fields once the magnetic field enters the low gas pressure region. This process where the internal energy amplifies the magnetic field is similar to the explosion process (Moreno-Insertis et al. 1995; Rempel & Schüssler 2001; Hotta et al. 2012a). In the explosion process, the rising motion in the superadiabatic stratification leads to an entropy difference between the inside and outside of the flux tube. Figure 13 shows the PDFs between (a) the perturbation density and magnetic pressure and (b) the perturbation entropy and the magnetic pressure. While the density well correlates with the magnetic pressure, the entropy does not. This result indicates that the entropy does not contribute to the amplification and that the process achieved in this study is different from the explosion process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Also, as shown in Hotta et al. (2015b), the absolute amplitude of entropy perturbation increases with increased resolution. This increase also occurs in this study (see Section 3.7), and this tends to lower the gas pressure in the downflow region (s₁ < 0) because the linearized equation of state is expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{1}}{{p}_{0}}=\gamma \displaystyle \frac{{\rho }_{1}}{{\rho }_{0}}+\displaystyle \frac{{s}_{1}}{{c}_{{\rm{v}}}}.\end{eqnarray} \tag{ 29 }$$

Again, Figure 13(b) shows that the correlation between the entropy perturbation and the magnetic pressure is not good, and this fact indicates that the increase of the entropy perturbation does not contribute to amplifying the magnetic field.

We also investigate the location of the strong magnetic field amplification. Figure 14 shows the PDF between the gas pressure perturbation p'₁ and the radial velocity v_r . When we compare the Low (panel (a)) and High (panel (c)) cases, the low gas pressure region appears, especially at the downflow region. Considering the result shown, we can draw an overall picture of the amplification process of the superequipartition magnetic field. A schematic picture is shown in Figure 15. In a hydrodynamic case without the magnetic field (left panel), when a fluid parcel in the upper layer descends, the parcel has a low pressure compared with the external fluid in the lower layer because of the stratification. This pressure imbalance is instantaneously relaxed by the sound wave. In a magnetic case, especially with an efficient small-scale dynamo like the High case, the situation changes. When a fluid parcel goes down to the lower layer, the small-scale magnetic field is involved. The low gas pressure inside the magnetic field p_i and the magnetic pressure B²/8π are balanced with the external gas pressure p_e, i.e., p_i + B²/8π = p_e. Thus, the magnetic energy is amplified by the compression, i.e., maintained by the internal energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also discuss the spatial scale of magnetic field amplification. The spectral magnetic energy is expressed by

$$\begin{eqnarray}&&{\widehat{E}}_{\mathrm{mag}}({\ell })=\displaystyle \frac{1}{8\pi }\widehat{{\boldsymbol{B}}}({\ell })\cdot {\widehat{{\boldsymbol{B}}}}^{* }({\ell }),\end{eqnarray} \tag{ 30 }$$

where $\widehat{}$ and * denote the spherical harmonic transform and the complex conjugate, respectively. Then, the time evolution of ${\widehat{E}}_{\mathrm{mag}}({\ell })$ can be written as (see details in Pietarila Graham et al. 2010; Rempel 2014)

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\widehat{E}}_{\mathrm{mag}}({\ell })={\widehat{T}}_{{\rm{m}}(\mathrm{STR})}+{\widehat{T}}_{{\rm{m}}(\mathrm{ADV})}+{\widehat{T}}_{{\rm{m}}(\mathrm{CMP})},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&{\widehat{T}}_{{\rm{m}}(\mathrm{STR})}=\displaystyle \frac{1}{8\pi }\widehat{{\boldsymbol{B}}}\cdot {\widehat{\left({\boldsymbol{B}}\cdot {\rm{\nabla }}\right){\boldsymbol{v}}}}^{* }+c.c.,\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\widehat{T}}_{{\rm{m}}(\mathrm{ADV})}=-\displaystyle \frac{1}{8\pi }\widehat{{\boldsymbol{B}}}\cdot {\widehat{\left({\boldsymbol{v}}\cdot {\rm{\nabla }}\right){\boldsymbol{B}}}}^{* }+c.c.,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\widehat{T}}_{{\rm{m}}(\mathrm{CMP})}=-\displaystyle \frac{1}{8\pi }\widehat{{\boldsymbol{B}}}\cdot {\widehat{\left({\boldsymbol{B}}{\rm{\nabla }}\cdot {\boldsymbol{v}}\right)}}^{* }+c.c.,\end{eqnarray} \tag{ 34 }$$

where c. c. indicates the complex conjugate expression. Each term in the spectral magnetic energy evolution at r = 0.9R_⊙ is shown in Figure 16. The magnetic energy transfer by the advection ${\widehat{T}}_{{\rm{m}}(\mathrm{ADV})}$ does not depend on the resolution in a middle (ℓ ∼ 10²) to large scale (ℓ ∼ 1). The advection term contribution ${\widehat{T}}_{{\rm{m}}(\mathrm{ADV})}$ is typically negative because the downward magnetic energy transport is dominant at this height. Around the smallest scale in each resolution, ${\widehat{T}}_{{\rm{m}}(\mathrm{ADV})}$ is positive. The dominant magnetic energy production source is the stretching ${\widehat{T}}_{{\rm{m}}(\mathrm{STR})}$, but the production rate decreases with increased resolution, especially at the small scale because the magnetic field strength and the resulting Lorentz feedback increase. Meanwhile, the compression contribution ${\widehat{T}}_{{\rm{m}}(\mathrm{CMP})}$ increases with the resolution at a middle scale (ℓ ∼ 100). The peak scale of the compression does not depend on the resolution. This result also supports our presented explanation of the amplification mechanism of the strong magnetic field. A complex small-scale magnetic field is concentrated at the downflow region. The field is strong enough to suppress the turbulent stretching, but the compression can still work.

Zoom In Zoom Out Reset image size

In this subsection, we discuss the driving mechanism of the thermal convection. In particular, the mechanism in which the large-scale convection is suppressed in the High case is discussed.

For the discussion in this subsection, we additionally define statistical values, spherical average $\tilde{Q}$, spherical rms Q_(rms), spherical correlation [Q₁ Q₂], and normalized spherical correlation $\overline{{Q}_{1}{Q}_{1}}$ as follows.

$$\begin{eqnarray}&&\widetilde{Q}(r)=\displaystyle \frac{1}{4\pi }{\int }_{S}{QdS}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{Q}_{(\mathrm{rms})}(r)=\sqrt{\displaystyle \frac{1}{4\pi }{\int }_{S}{\left(Q-\widetilde{Q}\right)}^{2}{dS}}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\left[{Q}_{1}{Q}_{2}\right](r)=\displaystyle \frac{1}{4\pi }{\int }_{S}{Q}_{1}{Q}_{2}{dS}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\overline{{Q}_{1}{Q}_{2}}(r)=\displaystyle \frac{\left[{Q}_{1}{Q}_{2}\right]}{{Q}_{1(\mathrm{rms})}{Q}_{2(\mathrm{rms})}}\end{eqnarray} \tag{ 38 }$$

We note that the spherical rms Q_(rms) defined in Equation (36) is different from the longitudinal rms Q_(RMS) defined in Equation (11). Figure 17 shows the superadiabaticity δ in different cases. The superadiabaticity is defined as follows:

$$\begin{eqnarray}&&\delta =-\displaystyle \frac{{H}_{p}}{{c}_{p}}\displaystyle \frac{d\tilde{s}}{{dr}}.\end{eqnarray} \tag{ 39 }$$

Zoom In Zoom Out Reset image size

We observe a thermal convectively stable region (δ < 0) in all cases. This layer is common in an effectively high Prandtl number convection (Bekki et al. 2017; Hotta 2017; Käpylä 2019). The effective high Prandtl number is achieved with the strong small-scale magnetic field. In a high Prandtl number regime, the thermal structure does not diffuse, and low entropy material is accumulated at the base of the convection zone. Brandenburg (2016) also shows that a nonlocal convection can cause this type of subadiabatic layer in his analytical model. This process results in the convectively stable region (δ < 0). Bekki et al. (2017) show that when the stable region is achieved around the base of the convection zone, the large-scale flow is suppressed, because the convection driving scale in the deeper layer is larger because of the large pressure/density scale height. Because the stable region expands and the absolute value of superadiabaticity ∣δ∣ increases with the resolution, this effect should contribute to suppressing the large-scale convection. The difference of the superadiabaticity, however, between the Low and Middle cases is larger than that between the Middle and High cases, while the large-scale flow is significantly suppressed only in the High case. This indicates that the main reason for the large-scale suppression is not the change of the superadiabaticity.

The basic value that determines the convection velocity is the energy flux. In the solar convection zone, the energy flux is fixed by the efficiency of the nuclear fusion. Figure 18 shows different types of fluxes. Definitions of the enthalpy F_e, kinetic F_k, Poynting F_m, radiative F_r, and total F_t flux densities are (see Hotta et al. 2014)

$$\begin{eqnarray}&&{F}_{{\rm{e}}}=\left({e}_{1}+\displaystyle \frac{{p}_{1}}{{\rho }_{0}}-\displaystyle \frac{{p}_{0}}{{\rho }_{0}^{2}}{\rho }_{1}\right)\rho {v}_{r},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{F}_{{\rm{k}}}=\displaystyle \frac{1}{2}\rho {v}^{2}{v}_{r},\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&{F}_{{\rm{m}}}=\ \displaystyle \frac{1}{4\pi }\left[\left({B}_{\theta }^{2}+{B}_{\phi }^{2}\right){v}_{r}-\left({v}_{\theta }{B}_{\theta }+{v}_{\phi }{B}_{\phi }\right){B}_{r}\right],\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{F}_{{\rm{r}}}={F}_{\mathrm{rad}}+{F}_{\mathrm{art}},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{F}_{{\rm{t}}}={F}_{{\rm{e}}}+{F}_{{\rm{k}}}+{F}_{{\rm{m}}}+{F}_{{\rm{r}}},\end{eqnarray} \tag{ 44 }$$

where e is the internal energy calculated with the OPAL repository. F_rad and F_art are defined at Equations (7) and (8), respectively. For convenience, the sum of the physics-based radiation flux density F_rad and an artificial energy flux density F_art is called the radiative flux density F_r in this study.

Zoom In Zoom Out Reset image size

In Figure 18, we integrate the flux densities over the full sphere and evaluate each corresponding luminosity (flux). The enthalpy flux L_e (magenta) slightly decreases with increased resolution. The decrease is more moderate than expected from the convection velocity suppression (Figure 8). In the mixing length theory, the enthalpy flux scales as ${L}_{{\rm{e}}}\propto {v}_{{\rm{c}}}^{3}$, and the suppression of the convection velocity v_c should have a strong influence on the enthalpy flux. This deviation from the mixing length theory is essential to investigate the suppression mechanism of the convection velocity. The slight decrease of the enthalpy flux can be compensated for by the kinetic flux. As is usual, the kinetic flux is negative because the downflow has larger kinetic energy. This is reduced because of the convection velocity suppression. The Poynting flux has minor contributions, but the flux has a negative value. This downward Poynting flux is also caused because the downflow region has larger magnetic energy. Figure 19 shows the 2D energy flux density distribution in the High case. The enthalpy flux density does not show a significant dependence on the latitude (Figure 19(a)). 〈F_k〉 is always negative in all latitudes. The inward kinetic energy flux is most effective at the equator and the poles (Figure 19(b)). Latitudinal variation is most prominent in the Poynting flux (Figure 19(c)). 〈F_m〉 is positive and negative at the low and high latitudes, respectively.

Zoom In Zoom Out Reset image size

In this paragraph, we discuss why the energy flux (especially the enthalpy flux) is maintained even with the suppressed convection velocity. With the equation of state for the perfect gas, the enthalpy flux density can be expressed as

$$\begin{eqnarray}&&{F}_{{\rm{e}}}\sim {\rho }_{0}{c}_{p}\left[{v}_{r}{T}_{1}\right].\end{eqnarray} \tag{ 45 }$$

Because the background density, ρ₀, and heat capacity at constant pressure, c_p , do not change in a low Mach number situation, the correlation between the radial velocity, v_r , and the temperature perturbation, T₁, determines the enthalpy flux. Figure 20 shows analysis to this end. Figure 20(a) shows the spherical rms for the radial velocity v_r(rms). As discussed, the convection velocity is suppressed. Figure 20(b) shows the spherical rms for the temperature perturbation T_1(rms). T_1(rms) increases with increased resolution. The magnetic field is amplified in the higher-resolution simulations that suppress the mixing between the upflow and downflow. This process increases the temperature perturbation (see also Hotta et al. 2015b). In addition, the latitudinal temperature difference increases because of the presented process (see Section 3.8). The increased latitudinal temperature difference also contributes to increasing the spherical rms for the temperature T_1(rms). Figure 20(c) shows the normalized spherical correlation between the radial velocity v_r and the temperature perturbation T₁. The correlation decreases with the increase in the resolution. This correlation should be good when the flow obeys the thermal convection. In the high-resolution simulations, small-scale turbulence, which does not behave as the thermal convection, increases, and the correlation decreases. As a result, the dimensional correlation [v_r T₁], which directly determines the energy flux, stays the same among different resolutions (Figure 20). As a summary, the suppressed convection velocity and the worse normalized correlation are compensated by the increase in temperature perturbation to maintain the energy flux.

Zoom In Zoom Out Reset image size

We also discuss the energy flux from the viewpoint of the spatial scale. Figure 21(a), (b), and (c) show the spectra of the radial velocity, the temperature perturbation, and these correlations, respectively. As discussed already, the radial velocity is suppressed in all the scales (Figure 21(a)). The increase of the temperature perturbation in the High case is mainly seen in the small scales (ℓ > 40, Figure 21(b)). These results support our interpretation of the increase of the temperature perturbation in the High case. We expect the suppression of the mixing by the magnetic field to increase the temperature perturbation effectively. This process is most effective on a small scale. The combination of the decrease of the radial velocity and the increase of the temperature perturbation in the small scales leads to a situation where the correlation $\widehat{{v}_{r}{T}_{1}}$ in the small scale (ℓ > 40) stays the same (Figure 21(c)). In addition, the higher-resolution calculation has a long tail of the correlation on a smaller scale. This result indicates that a significant fraction of the energy is transported by the small-scale turbulence in the High case. To evaluate the importance of the small scale in energy transport, we calculate a value S_ℓ defined as follows.

$$\begin{eqnarray}&&{S}_{{\ell }}=\displaystyle \frac{{\sum }_{{\ell }^{\prime} ={\ell }}^{{{\ell }}_{\max }}\widehat{{v}_{r}{T}_{1}}({\ell }^{\prime} )}{{\sum }_{{\ell }^{\prime} =0}^{{{\ell }}_{\max }}\widehat{{v}_{r}{T}_{1}}({\ell }^{\prime} )}\end{eqnarray} \tag{ 46 }$$

Zoom In Zoom Out Reset image size

S_ℓ shows the fraction of the correlation from ℓ to ${{\ell }}_{\max }$ to the total correlation. Figure 21(d) shows the dependence of S_ℓ on the resolution. S_ℓ reaches unity around ℓ ∼ 5 in the Low and Middle cases, while ℓ ∼ 20 is enough for S_ℓ to reach unity in the High case. This indicates that, in the High case, a significant fraction of the energy is transported by the middle to small scales (ℓ > 20), and the large scale cannot transport the energy. We conclude that this is the main reason why the large-scale convection is suppressed in the High case.

We also investigate the convection driving mechanism in the viewpoint of the scale. In the analyses, we assume the background density is constant in time. Similar to the spectral magnetic energy ${\widehat{E}}_{\mathrm{mag}}$ discussed in Section 3.6, the spectral kinetic energy ${\widehat{E}}_{\mathrm{kin}}$ evolution equation can be written as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\widehat{E}}_{\mathrm{kin}}={\widehat{T}}_{{\rm{k}}(\mathrm{ADV})}+{\widehat{T}}_{{\rm{k}}(\mathrm{BUO})}+{\widehat{T}}_{{\rm{k}}(\mathrm{LOR})},\end{eqnarray} \tag{ 47 }$$

where

$$\begin{eqnarray}&&{\widehat{T}}_{{\rm{k}}(\mathrm{ADV})}=-\displaystyle \frac{1}{2}{\rho }_{0}\widehat{{\boldsymbol{v}}}\cdot {\widehat{{\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{v}}}}^{* }+c.c.,\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\widehat{T}}_{{\rm{k}}(\mathrm{BUO})}=-\displaystyle \frac{1}{2}{\rho }_{0}\widehat{{\boldsymbol{v}}}\cdot {\widehat{-{\rho }_{1}{\boldsymbol{g}}+{\rm{\nabla }}p}}^{* }+c.c.,\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{\hat{T}}_{{\rm{k}}(\mathrm{LOR})}=\frac{1}{8\pi }\hat{{\boldsymbol{v}}}\cdot {\hat{\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)\times {\boldsymbol{B}}}}^{* }+c.c.\end{eqnarray} \tag{ 50 }$$

The result at r = 0.9R_⊙ is shown in Figure 22. We note that while the Coriolis force should affect the spectral analysis, the amplitude is 1–2 orders of magnitude smaller than the other values, and we do not include it in our discussion. The general tendency is that the buoyancy drives the thermal convection and the advection, and the Lorentz force reduces the kinetic energy in almost all the scales. Also, both the energy production and the suppression on the large scale are small in the High case. For all the contributions to the kinetic energy transfer, the velocity is multiplied. In the High case, the kinetic energy in the large scale is reduced, and the reduction of the energy transfer seems an obvious result. To investigate the effective importance of the large-scale suppression, we normalize the kinetic energy transfer by $\sqrt{\widehat{v}{\widehat{v}}^{* }/r}$ (Figure 22(b)). The normalized kinetic energy transfer by the buoyancy ${\widehat{T}}_{{\rm{k}}(\mathrm{BUO})}$ is reduced only in the High case. We also observe the suppression of the Lorentz force contribution. These results indicate that the suppression of the large-scale kinetic energy is caused by the suppression of the buoyancy. The magnetic field on a large scale does not directly contribute to the large-scale suppression.

Zoom In Zoom Out Reset image size

In this subsection, we discuss the force balance on the meridional plane, especially about the Taylor–Proudman constraint. The differential rotation in the High case does not obey the Taylor–Proudman constraints, i.e., ∂Ω/∂z ≠ 0, where z indicates the direction of the rotational axis. To address this aspect, we need to analyze the vorticity equation (e.g., Miesch & Hindman 2011). The longitudinal component of the vorticity equation in a steady state, ∂/∂t = 0, is written as (see also Balbus et al. 2009, for more detailed discussions about the balance)

$$\begin{eqnarray}\begin{array}{ccc}\mathop{\underbrace{-2r\sin \theta {{\rm{\Omega }}}_{0}\frac{\partial \langle {{\rm{\Omega }}}_{1} \rangle }{\partial z}}}\limits_{{P}_{\mathrm{COR}}} & = & \mathop{\underbrace{{{\rm{\nabla }}\times {\boldsymbol{v}}\times {\boldsymbol{\omega }}}_{\phi }}}\limits_{{P}_{\mathrm{ADV}}}\\ & & +\mathop{\underbrace{\frac{g}{{\rho }_{0}r}{\left(\frac{\partial \rho }{\partial s}\right)}_{p}\frac{\partial \langle {s}_{1} \rangle }{\partial \theta }}}\limits_{{P}_{\mathrm{BAR}}}+{\underbrace{{\unicode{x0232A}{\rm{\nabla }}\times \left[\frac{1}{4\pi {\rho }_{0}}\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)\times {\boldsymbol{B}}\right]\unicode{x02329}}_{\phi }}}_{{P}_{\mathrm{MAG}}}\end{array}.\end{eqnarray} \tag{ 51 }$$

We show each term in the equation in Figure 23. To suppress realization noise, especially in P_ADV ad P_MAG, we use a Gaussian filter with a width of 5 × 5 grid points. We also show the spherically averaged 1D profile of each term in Figure 24. We use a Gaussian filter with a width of five grid points also for the 1D profile. The raw data are shown with transparent lines. The results clearly show that the deviation from the Taylor–Proudman theorem (P_COR) is mainly balanced by the baroclinic term (P_BAR). We see a significant deviation from the Taylor–Proudman theorem around the top boundary. This is maintained both by the advection (P_ADV) and the magnetic field (P_MAG). While this tendency is important to discuss the near-surface shear layer, we leave this for our future publication for the near-surface layer (H. Hotta, K. Kusano & T. Sekii 2022, in preparation). In this paper, we focus on the discussion about the Taylor–Proudman theorem in the middle of the convection zone. The result shows that the Coriolis force is balanced with the baroclinic term, i.e., the latitudinal entropy gradient. Hotta (2018) argues that the efficient small-scale dynamo and generated magnetic field help construct the entropy gradient. As shown in Section 3.7, the temperature perturbation increases with increased resolution. In addition, the convection velocity is reduced in the higher resolutions (Figure 8). The Coriolis force bends a warm upflow (cold downflow) poleward (equatorward). Both the high-resolution effects (increasing the temperature perturbation and reducing the convection velocity) enhance this process. Figure 25 shows the entropy and the temperature distributions in the High case. We succeed in reproducing the negative entropy and temperature gradient in the whole convection zone. Miesch et al. (2006) enforce the entropy gradient at the bottom boundary to avoid the Taylor–Proudman constraint (see also Miesch et al. 2008; Fan & Fang 2014). Also, Brun et al. (2011) maintains the entropy gradient by a dynamical coupling of the convection and radiation zones (see also Rempel 2005). In their studies, maintaining the negative entropy gradient in the near-surface equator is difficult, and the differential rotation tends to be the Taylor–Proudman type topology in the near-surface equator region (for examples, see Figures 10 and 13 of Brun et al. 2011). In our simulations, the entropy gradient is generated by the turbulent process throughout the convection zone, and the differential rotation can avoid the Taylor–Proudman constraint, which is consistent with the observations (e.g., Schou et al. 1998). Figure 26 shows the resolution dependence of the entropy and the temperature gradient. It is clearly shown that the entropy and the temperature gradient increase with resolution. This result also indicates that the magnetic field maintains the entropy and temperature gradients because the magnetic strength increases with the resolution. The temperature difference between the equator and the pole at the base of the convection zone is 8 K in the High case. This value corresponds to the AB3 case in Miesch et al. (2006) with which they argue their most solar-like profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Recently Matilsky et al. (2020) show that the fixed flux boundary condition, which is similar to that in this study, is easier to generate the non-Taylor–Proudman differential rotation than the fixed entropy boundary condition. They show that the entropy gradient is generated by the anisotropic enthalpy flux caused by the Busse column in the low latitude. Since we use the fixed flux boundary condition, our calculation should also be benefited by the numerical setting.

This subsection discusses the angular momentum transport, and we explain why the equator is rotating faster than the polar region. To discuss the angular momentum transport, we should start from the angular momentum conservation law that is approximated as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left({\rho }_{0}\langle { \mathcal L }\rangle \right)={G}_{\mathrm{REY}}+{G}_{\mathrm{MER}}+{G}_{\mathrm{MAG}},\end{eqnarray} \tag{ 52 }$$

where

$$\begin{eqnarray}&&{G}_{\mathrm{REY}}=-{\rm{\nabla }}\cdot \left({\rho }_{0}\lambda \langle {{\boldsymbol{v}}}_{{\rm{m}}}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle \right),\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{G}_{\mathrm{MER}}=-{\rm{\nabla }}\cdot \left({\rho }_{0}\langle {{\boldsymbol{v}}}_{{\rm{m}}}\rangle \langle { \mathcal L }\rangle \right),\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{G}_{\mathrm{MAG}}=-{\rm{\nabla }}\cdot \left(-\lambda \displaystyle \frac{\langle {{\boldsymbol{B}}}_{{\rm{m}}}{B}_{\phi }\rangle }{4\pi }\right),\end{eqnarray} \tag{ 55 }$$

and $\lambda =r\sin \theta$. ${ \mathcal L }=\lambda {u}_{\phi }=\lambda {v}_{\phi }+{\lambda }^{2}{{\rm{\Omega }}}_{0}$ is the specific angular momentum. G_REY, G_MER, and G_MAG are the angular momentum by the turbulence, mean flow (meridional flow), and magnetic field. We define B _m = B_r e _r + B_θ e _θ . Because the large-scale magnetic field 〈 B 〉 is weak in this study (see Table 1), we do not divide the magnetic contribution G_MAG to turbulent component ${\boldsymbol{B}}^{\prime}$ and large-scale component 〈 B 〉.

At first, we discuss how to transport the angular momentum equatorward in the High (and Middle) cases. To this end, we evaluate the temporal evolution of the latitudinal angular momentum flux density at θ = π/4. The latitudinal angular momentum flux densities are as follows:

$$\begin{eqnarray}&&{F}_{\mathrm{tur}}={\rho }_{0}\lambda \langle {v}_{\theta }^{{\prime} }{v}_{\phi }^{{\prime} }\rangle ,\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{mer}}={\rho }_{0}\langle {v}_{\theta }\rangle \langle { \mathcal L }\rangle ,\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{mag}}=-\lambda \displaystyle \frac{\langle {B}_{\theta }{B}_{\phi }\rangle }{4\pi },\end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{tot}}={F}_{\mathrm{tur}}+{F}_{\mathrm{mer}}+{F}_{\mathrm{mag}}.\end{eqnarray} \tag{ 59 }$$

These fluxes are radially averaged at θ = π/4 in Figure 27. Although the turbulent angular momentum transport (orange line, F_tur) is positive (equatorward) in all cases, the resulting differential rotation is different in all cases. This indicates that the equatorward angular momentum transport cannot be the reason why we have the fast equator in the High case. A prominent difference can be seen in the angular momentum transport by the meridional flow (blue line, F_mer). F_mer is negative (poleward) in the initial phase ( <600 days) in all cases. This poleward angular momentum transport leads to the fast pole in the initial phase (see Figure 35 in Appendix C). While F_mer stays almost negative in the Low case, the other cases clearly show positive F_mer in the latter phase. This is the reason why we have a fast equator in the Middle and High cases. Because of the low Mach number situation, ${\rm{\nabla }}\cdot \left({\rho }_{0}{{\boldsymbol{v}}}_{{\rm{m}}}\right)=0$ is approximately satisfied. This leads to ∫ρ₀ v_θ rdr ∼ 0 at an arbitrary latitude with the closed boundary condition for the radial velocity. Because the specific angular momentum is

$$\begin{eqnarray}&&\langle { \mathcal L }\rangle ={r}^{2}{\sin }^{2}\theta \left(\langle {{\rm{\Omega }}}_{1}\rangle +{{\rm{\Omega }}}_{0}\right),\end{eqnarray} \tag{ 60 }$$

the deeper layers (small r) tend to have smaller angular momentum than near-surface layers (large r). The equatorward meridional flow in the middle of the convection zone is the direct reason for accelerating the equator. Due to the mass conservation, the fast meridional flow, which overcomes the poleward angular momentum transport near the surface, requires the poleward meridional flow around the base of the convection zone. The poleward meridional flow around the base of the convection zone is the primary key to why we have the fast equator in the High case.

Zoom In Zoom Out Reset image size

Gyroscopic pumping is useful in understanding the maintenance mechanism of the meridional flow. Gyroscopic pumping is the angular momentum conservation law in a steady state.

$$\begin{eqnarray}&&{\rho }_{0}\langle {{\boldsymbol{v}}}_{{\rm{m}}}\rangle \cdot {\rm{\nabla }}\langle { \mathcal L }\rangle \sim -{\rm{\nabla }}\cdot \left({\rho }_{0}\lambda \langle {{\boldsymbol{v}}}_{{\rm{m}}}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle -\lambda \displaystyle \frac{\langle {{\boldsymbol{B}}}_{{\rm{m}}}{B}_{\phi }\rangle }{4\pi }\right)\end{eqnarray} \tag{ 61 }$$

In this study and the solar case, the differential rotation is weak, i.e., Ω₁/Ω₀ ∼ 0.1, and the angular momentum $\langle { \mathcal L }\rangle$ does not change significantly even after the differential rotation is constructed. Thus, the gyroscopic pumping indicates that the angular momentum transport by the Reynolds stress and the magnetic field determines the topology of the meridional flow. Figure 28 shows each term in Equation (52). From this figure, we can discuss two topics: one is the generation mechanism of the poleward meridional flow around the base of the convection zone, and the other is the acceleration mechanism of the near-surface equator.

Zoom In Zoom Out Reset image size

At first, we discuss the generation mechanism of the meridional flow. Because of the poleward meridional flow around the base of the convection zone, the angular momentum increases (Figure 28(e) and (f)). This is compensated for by the magnetic angular momentum transport (G_MAG, Figure 28(h) and (i)). In other words, the poleward meridional flow around the base of the convection zone is maintained by the magnetic angular momentum transport. The magnetic angular momentum transport decreases the angular momentum around the base of the convection zone, and the poleward meridional flow increases it as compensation. While the turbulent angular momentum transport tends to increase the angular momentum around the base of the convection zone (Figure 28(a), (b), and (c)), this is not enough to compensate for the decrease by the magnetic angular momentum transport (see also Figure 34 for the sum of G_REY and G_MER).

As for the increase of the angular velocity in the near-surface equator, the magnetic angular momentum has the main contribution. The major difference in the differential rotation between the Middle and High cases is the angular velocity in the near-surface equator (Figure 10). The High case has a large angular velocity there, which is more consistent with the solar observation. It is apparent that this increase in the angular velocity in the High case is caused by the magnetic angular momentum transport (G_MAG, Figure 28(i)). In all cases, the near-surface equator is accelerated by G_MAG, but the amplitude of G_MAG increases in the High case because the magnetic field strength increases with the resolution (see Section 3.6).

In summary, the equatorward latitudinal angular momentum transport is done by the meridional flow constructed by the magnetic angular momentum transport. To have the large angular velocity in the near-surface equator, we need additional contributions by the magnetic field, which is stronger in the higher resolutions. For the angular momentum transport, both the strength and correlation are important. We analyze the result in this regard in the following paragraph.

Figure 29 shows the correlation between the velocities and the magnetic fields. For both the Reynolds stress $\langle {v}_{i}^{{\prime} }{v}_{j}^{{\prime} }\rangle$ and the Maxwell stress 〈B_i B_j 〉, the main contribution is the radial transport. The distributions of G_REY and G_MAG are roughly explained by the radial angular momentum transport. The Reynolds stress transports the angular momentum radially inward, while the Maxwell stress transports it in the opposite direction. The radially inward angular momentum transport is a usual result with a high Rossby number (weak rotational influence) situation (see Gastine et al. 2013; Hotta et al. 2015a; Featherstone & Miesch 2015; Karak et al. 2015). While the Rossby number (Ro) decreases from the Low to High cases (see Table 1), the radially inward angular momentum transport does not change much or even increase (see also normalized correlation in Figure 36). This result indicates the weak influence of rotation on the small-scale flow in all the cases. The local Rossby number (Ro_ℓ ), which increases from the Low to High cases, measuring the rotational influence, is not small enough to maintain the solar-like differential rotation by the Reynolds stress. As explained in the previous paragraph, the essential reason for the poleward meridional flow around the base of the convection zone and the large angular velocity around the near-surface equator is the magnetic angular momentum transport. Figure 29 shows that the negative correlation 〈B_r B_ϕ 〉, i.e., the radially outward, magnetic angular momentum transport, is responsible for both of these. The radially outward transport decreases and increases the angular momentum at the base and the top of the convection zone, respectively.

Zoom In Zoom Out Reset image size

The main possible reasons for the magnetic field correlation are the shear and the alignment to the flow. The shear term of the induction equation is written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{r}}{\partial t}=\displaystyle \frac{{B}_{\phi }}{r\sin \theta }\displaystyle \frac{\partial {v}_{r}}{\partial \phi }+[...],\end{eqnarray} \tag{ 62 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{\phi }}{\partial t}={B}_{r}\displaystyle \frac{\partial {v}_{\phi }}{\partial r}+[...].\end{eqnarray} \tag{ 63 }$$

Thus, the shear of the flow can correlate the magnetic field components. In addition, the magnetic induction equation in the high conductivity limit is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times \left({\boldsymbol{v}}\times {\boldsymbol{B}}\right).\end{eqnarray} \tag{ 64 }$$

This means that, when the magnetic field is parallel to the velocity, v × B = 0, the magnetic field does not evolve more. Conversely, the magnetic field tends to be parallel to the velocity. To understand the origin of the negative correlation of 〈B_r B_ϕ 〉, Figure 30 shows the PDF of (a) B_r B_ϕ versus $\partial {v}_{r}/\partial \phi /r\sin \theta$, (b) B_r B_ϕ versus ∂v_θ /∂r, and (c) B_r B_ϕ versus ${v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }$ at r = 0.9R_⊙ in the High case. While we do not see clear correlation between B_r B_ϕ and shears (Figure 30(a) and (b)), B_r B_ϕ and ${v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }$ correlate well. This indicates that the origin of the negative 〈B_r B_ϕ 〉 is not the flow shear but the negative correlation of velocities $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$. Figure 29(a), (b), and (c) shows that $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$ is negative at all latitudes. This is caused by the Coriolis force. The Coriolis force in the longitudinal equation of motion is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{\phi }}{\partial t}=[...]-2{{\rm{\Omega }}}_{0}\left({v}_{r}\sin \theta +{v}_{\theta }\cos \theta \right).\end{eqnarray} \tag{ 65 }$$

Thus, the radial velocity, which is the source of the thermal convection, is bent by the Coriolis force, and the negative $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$ is caused. Because the magnetic induction equation only suggests that the magnetic field tends to be parallel to the velocity, it is possible that 〈B_r B_ϕ 〉 is the origin of the $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$. To confirm the origin of $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$, we compare the hydro case (High-HD) and the magnetic case (High) in Figure 31 with PDFs. Figure 31 shows PDFs of (a) ${v}_{r}^{{\prime} }$ versus ${v}_{\phi }^{{\prime} }$; (b) B_r and B_ϕ from the High case are shown. Figure 31(c) shows the correlation of ${v}_{r}^{{\prime} }$ versus ${v}_{\phi }^{{\prime} }$ from the High-HD case. Even in the hydro case, we see a similar correlation between ${v}_{r}^{{\prime} }$ and ${v}_{\phi }^{{\prime} }$ (Figure 31(c)) to the magnetic case (Figure 31(a)). This result shows that the magnetic field is not the main origin of $\langle {v}_{r}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$, but the velocity is the origin of 〈B_r B_ϕ 〉.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We decompose the Reynolds stress to radial r and colatitudinal θ components. We note that the decomposition of parallel and perpendicular to the rotational axis directions, i.e., z and λ directions, are also useful. Since the constant surface of the specific angular momentum is parallel to the cylindrical surface (constant λ) in the leading order, it is difficult for the meridional flow to transport the angular momentum through the constant λ surface, i.e., $\int {\rho }_{0}\langle {v}_{\lambda }\rangle \langle { \mathcal L }\rangle {dz}\sim 0$ due to the anelastic approximation ∫ρ₀〈v_λ 〉dz ∼ 0. Thus the meridional flow mainly transports the angular momentum in the z direction. This tendency indicates that $\langle {v}_{z}^{{\prime} }{v}_{\phi }^{{\prime} }\rangle$ and $\langle {v}_{\lambda }^{{\prime} }{v}_{\phi }\rangle$ are responsible for the generation of the meridional flow and the differential rotation, respectively.

In this study, we analyze the Reynolds stress just as the velocity correlation $\langle {v}_{i}^{{\prime} }{v}_{j}^{{\prime} }\rangle$. The stress includes the diffusive part, i.e., so-called turbulent viscosity and nondiffusive part, so-called Λ effect (Ruediger 1980). If we can distinguish these two from the Reynolds stress, we can directly evaluate the anisotropy.

In this study, compression is an important process to amplify the magnetic field. In the process, we can use the internal energy, which is about 10⁶ times larger than the kinetic energy in the deep convection zone. We have not reached numerical convergence, i.e., the higher resolutions show a stronger magnetic field (see Figure 8(b)). Featherstone & Hindman (2016b) show that the kinetic energy can converge to a specific value by increasing the Rayleigh number in their hydrodynamic run around Ra ∼ 10⁵. Hindman et al. (2020) carry out a similar survey with the rotation and find the saturation around Ra ∼ 10⁷ with the solar rotation case. Our Rayleigh number is larger than these critical Rayleigh numbers. These results indicate that when we carry out a hydrodynamic run, the kinetic energy is expected to be converged. In addition, our result shows that the magnetic field generation requires further resolution to be numerically converged than the hydrodynamic models. Currently, we cannot conclude the magnetic field strength in the real Sun, but it is most probably stronger than our simulation result. The strength may be determined by a balance between the generation of compression and the destruction by the small-scale turbulence. At the same time, provided our suggested mechanism of the magnetic angular momentum transport is correct, we do not expect that the magnetic field strength in the real Sun is much stronger than our simulation because our differential rotation is similar to the observational results. In our simulation, the magnetic field directly determines the differential rotation topology. If our magnetic field strength is completely different from reality, the differential rotation is also away from reality. This is not the case in the simulation. Of course, there is a possibility that our suggested mechanism is incorrect. We need to perform higher-resolution simulations to reach numerical convergence and to understand magnetic field strength in reality and the validity of the mechanism.

We note that we use the RSST in our calculation. When the Alfvén velocity exceeds the reduced speed of sound, the amplification efficiency should be decreased. In this study, the maximum magnetic pressure is B²/(8π) ∼ 1 × 10⁸dyncm⁻² (see Figure 12), and the effective gas pressure evaluated with the reduced speed of sound at r = 0.9R_⊙ is ${\rho }_{0}{c}_{{\rm{s}}}^{2}/{\xi }^{2}\sim 2.4\times {10}^{9}\,\mathrm{dyn}\,{\mathrm{cm}}^{-2}$. These values indicate that we can ignore the influence of the RSST on the compression in this study. We also emphasize that, even if the RSST influenced the result, it would weaken the magnetic field strength. Our conclusion, i.e., a strong magnetic field constructs the differential rotation, should be robust.

Figure 33 shows a comparison of kinetic energy spectra of the longitudinal velocity E_ϕ between the simulations and an observation. We follow the definition of the spectra of Gizon & Birch (2012), where

$$\begin{eqnarray}&&\displaystyle \frac{{\int }_{V}{v}_{\phi }^{2}/2{dV}}{{\int }_{V}{dV}}=\sum _{{\ell }\gt 0}\displaystyle \frac{{E}_{\phi }}{r}.\end{eqnarray} \tag{ 66 }$$

To exclude the contribution of the differential rotation, we exclude m = 0 mode, where m is the azimuthal wavenumber. The integration is carried out in the whole computational domain. The magenta line shows the result from the High case in this study. The blue line shows the result from Hotta et al. (2019). In the calculation, the horizontal extent is restricted to 200 Mm, but it covers the whole convection zone vertically from the base to the photosphere. As suggested by Hotta et al. (2019), the existence of the photosphere does not change the energy spectra in the deeper layers, and the magenta and blue lines are consistent. The black line shows the result from another global calculation (Miesch et al. 2008) at r = 0.98R_⊙. The orange line indicates the upper limit suggested by the local helioseismology (Hanasoge et al. 2012). While we still have a large discrepancy between the simulation and the observation, the difference is relaxed. In our simulation, the large-scale convection is suppressed because the small-scale turbulence can efficiently transport the energy. Also, in this regard, the higher resolution possibly changes the result more. Meanwhile, recently, the helioseismology results have been revised (Proxauf 2021). Our simulation results in the High case are highly consistent with the revised result of Greer et al. (2015). Currently, we cannot conclude if our convective velocity is correct or not. A more detailed comparison between the simulations and observation is needed.

Zoom In Zoom Out Reset image size

Miesch et al. (2012) evaluate the lower limit of the convective velocity from the dynamical balance for the differential rotation. The evaluated value is not consistent with the local helioseismology (Hanasoge et al. 2012). In their study, they do not consider the magnetic contribution for the construction of the differential rotation. In this study, we find that the magnetic field is a dominant contribution. This means that the convection velocity has large freedom. The Rossby number does not solely determine the differential rotation. One remaining restriction on the convection velocity is the energy flux. The solar luminosity L_⊙ is determined; thus, there should be a lower limit on the convection velocity to transport the required energy. Our simulation also shows that the temperature perturbation increases with the resolution. If the temperature perturbation increases, the lower limit on the convective velocity should be relaxed. At the same time, a significantly large temperature perturbation should be detected by the local helioseismology with a mean travel time. Future observations for the convection velocity as well as the temperature perturbation will contribute to solving the problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Currently, the local helioseismology for the meridional flow is still controversial. Zhao et al. (2013) indicate the double cell flow with the poleward meridional flow around the base of the convection zone. On the other hand, Gizon et al. (2020) show the equatorward meridional flow around the base. In this regard, our result is more consistent with Zhao et al. (2013)'s result. This discrepancy is caused by the difference in observations. Zhao et al. (2013) adopt Solar Dynamics Observatory (SDO) data, and Gizon et al. (2020) use both Solar and Heliospheric Observatory (SoHO) and Global Oscillation Network Group (GONG) data. Gizon et al. (2020) find that SoHO and GONG data are consistent, but SDO data show a systematic difference from these two data. We should also note that the observations still have a tiny sensitivity in the deep convection zone because it requires long enough, separated two endpoints of Δ ∼ 45° for evaluating the travel time (Giles 2000). The observation has not accomplished enough precise observations for these separated two endpoints. For example, Gizon et al. (2020) show meridional flow results with and without the data with Δ > 30°, but the result does not change. This indicates that the data with Δ > 30° are not used for their inversion because of the large error, and the equatorward meridional flow is caused by the constraint of the mass conservation. This result indicates that we cannot conclude that our meridional flow is inconsistent with Gizon et al. (2020)'s result. In order to check whether our meridional flow model is compatible with the helioseismic observations, it is needed to compute the seismic travel times based on this solution and to compare them with the observations. Observations from different viewing angles, such as the Solar Orbiter (Müller et al. 2013), also enable us to understand the whole topology of the meridional flow, which should be of significant impact on the understanding of the convection and magnetic fields in the solar convection zone.