A SIMULATION OF CONVECTIVE DYNAMO IN THE SOLAR CONVECTIVE ENVELOPE: MAINTENANCE OF THE SOLAR-LIKE DIFFERENTIAL ROTATION AND EMERGING FLUX

Figure 1 shows the magnetic and kinetic energies in the statistically steady convective flows in the simulation domain over a time span of about 74 yr. The total magnetic energy E_m maintained by the dynamo is about 10% of the total kinetic energy E_k of the convective envelope. The energy of the azimuthally averaged (mean) magnetic field E_{m, mean} only constitutes a small fraction of E_m, oscillating from about 1% to 10% of E_m.

Zoom In Zoom Out Reset image size

Figure 2 shows the depth variation of the mean entropy gradient established in the CZ domain in the statistical steady state. The entropy gradient reaches a value of about 4.26 × 10⁻⁶ erg⁻¹ K⁻¹ cm⁻¹ at the top boundary at about 0.97 R_s, which is of a similar order of magnitude as the entropy gradient (∼10⁻⁵ erg⁻¹ K⁻¹ cm⁻¹) at this depth in the solar model of JCD. Figure 3 shows the various horizontally integrated energy fluxes (normalized to the solar luminosity L_s) through the domain as a function of radius established in the statistical steady state. These are, respectively, the integrated radiative diffusive heat flux (red curve): L_rad = 4πr²F_rad, the convective enthalpy flux (black curve):

$$\begin{equation} L_{\rm conv} = 4 \pi r^2 \rho _0 c_p \langle v_r T_1 \rangle, \end{equation} \tag{ 17 }$$

the conductive energy flux by thermal diffusion (yellow curve):

$$\begin{equation} L_{\rm cond} = - 4 \pi r^2 K \rho _0 T_0 \frac{d \langle s_1 \rangle }{d r}, \end{equation} \tag{ 18 }$$

the kinetic energy flux (blue curve):

$$\begin{equation} L_{\rm kin} = 4 \pi r^2 \frac{\rho _0}{2} \langle v^2 v_r \rangle, \end{equation} \tag{ 19 }$$

the viscous energy flux (black dashed curve):

$$\begin{equation} L_{\rm vis} = - 4 \pi r^2 \langle v_i {\cal D}_{ir} \rangle, \end{equation} \tag{ 20 }$$

the Poynting flux (green curve):

$$\begin{equation} L_{\rm poyn}= - 4 \pi r^2 \left\langle \left(\frac{1}{4 \pi } ({\bf v} \times {\bf B}) \times {\bf B} \right)_r \right\rangle, \end{equation} \tag{ 21 }$$

the resistive energy flux (cyan curve):

$$\begin{equation} L_{\rm res}= 4 \pi r^2 \left\langle \left(\frac{1}{4 \pi }(\eta \nabla \times {\bf B}) \times {\bf B} \right)_r \right\rangle, \end{equation} \tag{ 22 }$$

and the sum of all the energy fluxes, L_tot, is shown as the dash–dotted curve in Figure 3. In the above, 〈〉 denotes averaging over the spherical shell surface and time. The gradual decline of L_tot reflects a numerical deviation from exact energy conservation, which is mainly caused by the numerical diffusion of the magnetic field due to the Alfvén wave-upwind scheme used for advancing the induction equation (Fan 2008; Stone & Norman 1992). This numerical dissipation of magnetic energy is not being put back into the thermal energy in the entropy equation and results in a loss of the total energy and hence a loss of about 13% of the total energy flux exiting the domain at the top compared to the total energy flux (L_s) entering the domain from the bottom. We note that the explicit resistive dissipation of the magnetic field (due to η), and both the explicit viscous dissipation (due to ν) and the numerical diffusion of momentum are put into the thermal energy in the entropy equation as resistive and viscous heating to maintain energy conservation. From Figure 3, it can be seen that the enthalpy flux of the resolved convection transports about 66% of the solar luminosity in the middle of the CZ, and due to the high thermal diffusivity K, thermal conduction also transports a substantial fraction of the solar luminosity (about 36% at the middle of the CZ). The kinetic energy flux of the convective flows is downward and peaks at about 16% of the solar luminosity. The energy fluxes due to the Poynting flux (mostly downward), resistive, and viscous transport are all much smaller.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the depth variation of the peak downflow (solid black curve), and the rms speed v_rms (dash–dotted black curve), of the statistical steady convective flows in the domain. Note that in computing v_rms, we take out the azimuthally averaged velocity components and only sum up the azimuthally fluctuating parts of the velocity components. Also shown are the corresponding magnetic field strength in equipartition with the peak downflow speed (solid red curve) and the rms speed (dash–dotted red curve). It can be seen that the equipartition field strength B_eq corresponding to the peak down flow speed reaches ≈63 kG, while B_eq corresponding to the rms speed is ∼10 kG for the deep and mid-convection zone, and decrease to about 5000 G near the top boundary at about r = 0.971 R_s. Following Käpylä et al. (2012), we compute the following non-dimensional numbers characterizing the convective flows. The Reynolds number R_e = u_rms/νk_f ranges from about 130 at the bottom to about 50 at the top, and with a mid-convection zone value of about 128, where k_f = 2π/(r_o − r_i) and u_rms = ((3/2)〈v_r² + v_θ²〉) is the rms velocity averaged over each depth, omitting the contribution from the azimuthal velocity. The Coriolis number C_O = (2Ω/u_{rms, all} k_f) = 1.3, where u_{rms, all} = ((3/2)〈v_r² + v_θ²〉) with the averaging 〈〉 done for the entire domain. We can compare the values of these non-dimensional numbers with the corresponding ones in Käpylä et al. (2012): R_e = 36 and C_O = 7.6. It appears the convective flow in our dynamo simulation is moderately more turbulent as characterized by the larger R_e, especially in the deeper layers of the CZ. Our Coriolis number C_O is significantly lower, indicating that our convective dynamo is operating in a significantly less rotationally dominant regime. If we were to scale their typical rms velocity to be similar to ours u_{rms, all} ≈ 100 m s⁻¹, then their C_O would imply a significantly more rapidly rotating stellar envelope (with the solar CZ depth) at about five times the solar rotation rate.

Zoom In Zoom Out Reset image size

Figure 5(a) shows the latitude-time variation of the mean (azimuthally averaged) toroidal magnetic field at a depth near the bottom of the CZ. The mean toroidal magnetic field tends to be of opposite signs for the two hemispheres, and exhibits an irregular cyclic behavior with oscillations of the field strength on timescales ranging from about 5 to about 15 yr and undergoes irregular sign/polarity reversals. The strongest mean toroidal field is concentrated near the bottom of the CZ (see Figure 5(b)), peaking at about 7 kG. Figure 5(c) shows a shell-slice of B_ϕ at a depth near the bottom of the CZ, at a cycle maximum phase indicated by the green line in Figure 5(a). It shows that strong toroidal fields B_ϕ of a preferred sign (opposite for the two hemispheres) are concentrated in individual channels or filaments in each hemisphere, reaching a peak field strength of about 30 kG, which exceeds the field strength in equipartition with the local rms convective speed (B_eq ≈ 13 kG), but is below the equipartition field strength corresponding to the peak down flow speed (B_eq ≈ 63 kG). Thus, these strong field filaments are not passively advected by convective flows, but would be pinned down by the strong down flows if in their paths.

Zoom In Zoom Out Reset image size

Figure 6(a) shows the time and azimuthally averaged rotation rate in the convective envelope self-consistently maintained in the convective dynamo simulation. It shows a solar-like differential rotation profile (e.g., Thompson et al. 2003) with a faster rotation rate at the equator than at the polar region by about 30% of the mean rotation rate, and more conical shaped iso-rotation contours in the mid-latitude zone. The time and azimuthally averaged mean meridional flow pattern is shown in Figure 6(b) in terms of the mass flux function f where $\rho _0 \langle {\bf v} \rangle = \nabla \times [(f / r \sin \theta) {\hat{\bf \phi }} ]$. The meridional circulation has a complex multi-cell structure with a counter-clockwise (clockwise) cell pattern in the low latitude region of the northern (southern) hemisphere, i.e., a poleward near-surface flow in the low latitude region.

Zoom In Zoom Out Reset image size

Interestingly, we find that the presence of the magnetic field is necessary for the self-consistent maintenance of the solar-like differential rotation profile in the current parameter regime. We have carried out a hydrodynamic simulation (hereafter referred to as the HD case) which is identical to the present convective dynamo simulation except that the magnetic field is set to zero, and found that a very different differential rotation profile (Figure 6(f)) is established in the statistical steady state of the HD simulation. It shows a significantly larger differential rotation with a faster rotation rate in the polar region than at mid-latitudes and the equator. The iso-rotation contours are also more cylindrical. The meridional flow (Figure 6(g)) shows a more prominent counter-clockwise (clockwise) cell pattern in the northern (southern) hemisphere in the middle depths in the CZ, with much weaker reversed cells in the near surface layer.

To compare the transport of angular momentum in the dynamo and the HD cases, we show in Figures 6(c)–(e) the meridional profile of the angular momentum flux density in the r_⊥ direction (perpendicular to and away from the rotational axis) due, respectively, to the Reynolds stress of the rotationally influenced convection (panel (c)):

$$\begin{equation} {\rm RS} = \rho _0 r_{\perp } \left\langle v^{\prime }_{r_{\perp }} v^{\prime }_{\phi } \right\rangle, \end{equation} \tag{ 23 }$$

the viscous stress (panel (d)):

$$\begin{equation} {\rm VS} = \rho _0 \nu r_{\perp }(\langle S_{\phi r } \rangle \sin \theta + \langle S_{\phi \theta } \rangle \cos \theta), \end{equation} \tag{ 24 }$$

and the Maxwell stress (panel (e)):

$$\begin{equation} {\rm MS}=- r_{\perp } \frac{1}{4 \pi } \langle B_{\phi } B_{r_{\perp }} \rangle, \end{equation} \tag{ 25 }$$

where 〈〉 denotes time and azimuthal averages and ' denotes the azimuthally varying component. The meridional profiles of the angular momentum flux density RS and VS for the corresponding HD case are shown in Figures 6(h) and (i). We find that there is a significant difference in the angular momentum flux density by the Reynolds stress between the dynamo and HD cases. The RS for the dynamo case shows an overall more outward transport in its meridional distribution compared to that for the corresponding HD case. Near the lower boundary in the mid-latitude range, the presence of the concentrated magnetic field (which helps to damp the convective downflows) results in a more enhanced outward RS flux at the lower boundary layer. The angular momentum flux density MS due to the Maxwell stress in the dynamo case (Figure 6(e)) is found to oppose RS, and its strength is greatest at the lower boundary layer where the magnetic field concentrates. The angular momentum flux density VS simply acts to reduce the differential rotation as expected for both the dynamo and the HD cases. The difference between the dynamo and the HD cases is more clearly seen by evaluating the net angular momentum fluxes in the r_⊥ direction integrated over individual concentric cylinders of radii r_⊥ centered on the rotation axis, as shown in Figure 7(a) for the dynamo case and Figure 7(b) for the corresponding HD case. It can be seen that the net angular momentum flux due to RS is outward throughout (except near the top boundary) for the dynamo simulation (black curve in Figure 7(a)), which drives a faster rotation in the outer equatorial region. The net angular momentum flux due to RS is mainly counteracted by the net flux due to MS by the Maxwell stress, with the remaining difference balanced by the significantly smaller net fluxes due to VS and the meridional flow. In contrast, the HD simulation shows a significant inward angular momentum transport due to the Reynolds stress (black curve in Figure 7(b)) across the inner cylinders in the high to mid-latitude region. This drives a faster rotation in the polar region. Thus, it appears that the presence of the magnetic field alters the convective flows such that the resulting Reynolds stress from the convective motions produces a more outward (away from the rotation axis) net transport of the angular momentum needed to drive a solar-like differential rotation. The Rossby number R_O = v_{rms, all}/(ΩH_p) is about 0.74 for the dynamo case and 0.96 for the HD case, where v_{rms, all} = 125 m s⁻¹ for the dynamo case and 157 m s⁻¹ for the HD case, is the rms velocity (with the azimuthally averaged mean flow velocity taken out) averaged over the entire volume, and H_p is the pressure scale height at the bottom of the CZ. The Rossby number measures the importance of the Coriolis force in the force balance. The lower Rossby number in the dynamo simulation shows that the magnetic fields suppress the convective motions so that they are more rotationally constrained.

Zoom In Zoom Out Reset image size

A recent systematic study by Gastine et al. (2014) of rotating stellar convection considering a wide range of models shows that the differential rotation profile transitions from being solar-like, with a faster rotating equator, to being anti-solar, with a faster polar rotation rate, at a value of about 1 for the Rossby number. This result is found to be quite general, independent of the detailed model setup (presence of a magnetic field, thickness of the convective layer, density stratification). Our HD case with R_O = 0.96 appears to be very close to the transition, and a reduction of ∼23% of the overall rms velocity and R_O by the presence of the magnetic field in the dynamo case is able to significantly alter the angular momentum transport by the rotationally constrained convection, leading to a transition into the solar-like differential rotation. We note that even though the HD case is quite close to the transition, its anti-solar differential rotation appears to be a stable solution not dependent on the history, i.e., not one of two bistable states (Gastine et al. 2014). We have arrived at the statistically steady HD solution with the anti-solar differential rotation by either starting from the initial setup with a seed velocity field as described in Section 2, or by starting from the statistically steady state dynamo solution (for which the differential rotation is solar-like) and zero out the magnetic field. The solar-like differential rotation obtained in the convective dynamo simulation also appears to be a stable solution that is not dependent on the history. The differential rotation at the pole and equator remain statistically steady without exhibiting any systematic drift for the ≳ 74 yr period (comparable to the maximum viscous timescale near the bottom of CZ) we have run after the dynamo solution has reached a statistical steady state. Also, as we start the dynamo simulation from the initial setup described in Section 2, we find that the convective dynamo goes through an earlier phase of anti-solar differential rotation (for ∼6 yr) before it evolves toward the solar-like differential rotation profile as the mean entropy gradient and the convective energy flux settle down to their statistical steady state. Thus, it appears that the solar-like differential rotation is the preferred stable solution in the dynamo case.

The transition to a solar-like differential rotation can alternatively be achieved in the non-magnetic hydro simulations by simply increasing the viscosity to reduce R_O. We have run another hydrodynamic simulation (hereafter referred to as the HVHD case, meaning "high viscosity hydro") where we increase the viscosity ν by 5 times (with the same ρ₀^−1/2 depth dependence) compared to the HD case (or the dynamo case). The resulting rms velocity of the statistical steady convection reached is v_{rms, all} = 120 m s⁻¹ and the Rossby number R_O = 0.71, much closer to those of the dynamo case. The bottom row panels of Figure 6 show, respectively, the resulting differential rotation profile (Figure 6(j)), meridional circulation (Figure 6(k)), the angular momentum flux density in the r_⊥ direction, RS (Figure 6(l)) due to the Reynolds stress, and VS (Figure 6(m)) due to the viscous stress. The integrated net (outward) angular momentum fluxes across concentric cylinders of radius r_⊥ centered on the rotation axis is shown in Figure 7(c). It is found that a solar-like differential rotation profile (Figure 6(j)) with faster rotating equator and with a more conical iso-rotation contours in mid-latitude zones is established, although the contrast of rotation rate between the equator and the polar region is larger, about 44% of the mean rotation rate (compared to about 32% in the dynamo case). The mean meridional circulation (Figure 6(k)) shows a counter-clockwise (clockwise) cell pattern in the low-latitude region of the northern (southern) hemisphere, i.e., a poleward near-surface flow in the low latitude region, similar to the dynamo case (Figure 6(b)). We find that the angular momentum transport in the r_⊥ direction due to the Reynolds stress for the HVHD case is very similar to that for the dynamo case, both in the meridional profile of the flux density (Figure 6(l) compared to Figure 6(c)) as well as in the integrated net flux across the constant r_⊥ concentric cylinders (Figure 7(c) compared to Figure7(a)). However, this similar outward net angular momentum flux by the Reynolds stress is now balanced almost entirely by the transport due to the viscous stress, with the absence of the Maxwell stress which is the major component that balances the angular momentum flux by the Reynolds stress in the dynamo case. The comparison between the dynamo and the HVHD cases suggests an effective role of enhanced viscosity played by the magnetic fields, which (1) suppresses the large-scale convective motions such that they are more rotationally constrained (lower R_O) to produce an outward transport of the angular momentum by the Reynolds stress, necessary to drive a solar-like differential rotation, and (2) takes up the main role to balance the Reynolds stress transport with the Maxwell stress instead of the viscous stress.

Further in Figure 8, we show the various horizontally integrated energy fluxes through the domain for the HD case (upper panel) and the HVHD case (lower panel), in comparison with the energy fluxes shown in Figure 3 for the dynamo case. It can be seen that the dynamo case and the HVHD case show a similar convective energy flux L_conv (reaching about 66% L_s in the dynamo case and about 60% L_s in the HVHD case). In both the dynamo and the HVHD cases, the downward kinetic energy flux L_kin (reaching about 16% L_s in the dynamo and 10% L_s in the HVHD case) is significantly reduced compared to the HD case (reaching about 45% L_s). In fact, in the HD case (upper panel in Figure 8), the downward kinetic energy flux is so large that the outward convective energy flux L_conv exceeds the solar luminosity in the middle of the CZ to counter it. The result here indicates again the similar role played by the magnetic fields and the enhanced viscosity in suppressing the downward convective flows.

Zoom In Zoom Out Reset image size

In the convective dynamo, the large-scale mean toroidal field as shown in Figure 5(b) is produced by the latitudinal differential rotation shearing a dipolar poloidal mean field. The reason that the mean toroidal field is concentrated towards the bottom of the CZ is mainly due to a downward advective transport of the magnetic energy in the bulk of the CZ, as represented by 〈v_r(B_θ² + B_ϕ²)/8π〉 shown in Figure 9(a). This causes the distribution of the magnetic energy (for both the mean field and the small-scale field) to be strongly concentrated toward the bottom (see Figure 9(b)). The decrease with depth of the magnetic diffusivity η would also promote stronger fields toward the bottom but is less important here because of the small magnitude of η and the long diffusive timescale: (Δr)²/η ∼ 10 yr, using the peak η value near the surface and the depth of the CZ domain Δr. The advective timescale for the downward magnetic energy transport across the CZ is Δr/u_m ∼ 0.5 yr is significantly shorter, where we have used u_m = 〈v_r(B_θ² + B_ϕ²)〉/〈B_θ² + B_ϕ²〉 evaluated at the middle of the CZ as a measure of the transport speed. Thus, the advective transport acts more quickly.

Zoom In Zoom Out Reset image size

In the midst of magneto-convection, we find occasional active regions like flux emergence events in the top layer of the simulation domain. Such an example is shown in Figure 10, where panels (a), (b), (c), and (d) show, respectively, snapshots of B_r, B_ϕ, v_r, and v_ϕ at a constant r slice at the depth of 30 Mm below the photosphere. The location of the emerging bipolar region is indicated by an arrow in the panels. It is characterized by a diverging bipolar pattern in B_r (panel (a)) and the emergence of a strong toroidal field patch reaching a peak field strength of 9800 G (panel (b)) (also see the animation in the online journal). The emerging region corresponds to an up flow region in v_r (panel (c)), but the upward velocity is not significantly different from that of other up flow convective cells. The zonal velocity v_ϕ of the emerging region shows a diverging pattern, and when averaged over the emerging region, is ∼100 m s⁻¹ faster than the mean zonal velocity of that latitude. Figure 10(e) shows the subsurface 3D magnetic field configuration in the convective envelope by showing field lines traced from randomly seeded points throughout the volume. The field lines are colored based on their azimuthal field B_ϕ as indicated by the color table. It can be seen that relatively more coherent bundles of strong toroidal flux are embedded in the turbulent magnetic fields. In Figure 10(f), regions of strong field strength where the Alfvén speed (v_a) exceeds the rms convective velocity (v_rms) for the corresponding depth is outlined with the equipartition iso-surfaces (with v_a/v_rms = 1), which are again colored based on the B_ϕ value on the iso-surfaces. There is a systematic preference for these strong flux regions to be green or of negative B_ϕ (red or of positive B_ϕ) in the northern (southern) hemisphere. The arrows in Figures 10(e) and (f) mark the toroidal flux bundle with super-equipartition field strength that gives rise to the emerging region.

Figure 11 shows a more zoomed in view of the thermodynamic properties of the emerging region at the same depth as that shown in the upper four panels of Figure 10. We see that there is a systematic reduction of density (i.e., buoyant; see Figure 11(b)) and pressure (Figure 11(d)) in the emerging region compared to the surrounding, although the reduction magnitude is rather moderate compared to the fluctuations seen in strong downflow lanes and in strong vertical flux tubes in the downflow lanes. The temperature change in the emerging region compared to the surrounding is smaller, partly due to the large thermal conduction. Averaged over the emerging region (area enclosed in the yellow contour in Figure 11(a), 〈ρ₁/ρ₀〉 ≈ −1.6 × 10⁻⁵, 〈p₁/p₀〉 ≈ −1.1 × 10⁻⁵, and 〈T₁/T₀〉 ≈ 0.5 × 10⁻⁵. It can be seen that the temperature change is relatively small compared to the density and pressure change in the emerging region. This suggests that the buoyancy or reduction in density is mainly due to the reduction in gas pressure provided largely by the presence of the magnetic pressure, instead of mainly due to an increase in temperature. In other words, the buoyancy contribution is more from the magnetic buoyancy than the thermal buoyancy.

Zoom In Zoom Out Reset image size

As can be seen from Figure 10(e), the emerging flux region is the apex of a roughly east-west oriented (toroidal) super-equipartition flux bundle which remain relatively coherent for some distance, before the two ends connect in complex ways to other flux systems. The following end of the coherent flux bundle extends into the middle of the CZ. The fact that the emerging flux has a prograde zonal speed of ∼100 m s⁻¹ relative to the mean zonal speed of the latitude indicates that it is not a toroidal flux tube rising in isolation from the bottom of the CZ. If it were, it would have a retrograde flow due to angular momentum conservation as is found in many previous studies of isolated rising flux tubes in the rotating solar CZ (e.g., Caligari et al. 1995; Fan 2008, 2009). The emerging flux bundle must have well mixed with the local plasma through reconnections, and is continually sheared and amplified by the differential rotation and the local flows against resistive dissipation. The sequence of images in Figure 12 show the sub-surface development of the super-equipartition emerging flux bundle (marked by the arrow) over a nine-day period prior to the time of the flux emergence event shown in Figure 10. It shows that local shear in the upper CZ contributes significantly to the development of the emerging flux bundle. The left column images show iso-volumes of super-equipartition fields with the surface of the volume colored by B_ϕ. The right column images show representative field lines traced from the iso-volume corresponding to the emerging flux bundle. It can be seen that a segment of a super-equipartition flux bundle in the middle of the CZ is sheared and stretched in the prograde direction into a hairpin turn with the upper side of the hairpin forming the emerging flux bundle reaching the top boundary.

Zoom In Zoom Out Reset image size

We have done a statistical study of the super-equipartition emerging fields. For a time period of about one year centered at the cycle maximum phase (green line in Figure 5(a)) and at an interval of 12 hr, we find in the shell slice at 30 Mm depth all the area where the emerging horizontal field exceeds $\sqrt{2}$ times the field strength in equipartition with the rms convective velocity of that depth. For each pixel (or grid point) of the selected emerging field area, we compute a tilt angle of the horizontal field vector based on the local B_ϕ and B_θ. The resulting tilt angle distribution of all the pixels is shown in Figure 13. The quadrant of the tilt angle is such that, if the sign of the azimuthal field B_ϕ is consistent with Hale's polarity rule of the cycle, i.e., negative (positive) in the northern (southern) hemisphere, then the tilt angle falls in quadrants I and IV. If the horizontal field vector is tilted clockwise (anti-clockwise) from the cycle preferred azimuthal field direction in the northern (southern) hemisphere by an acute angle, consistent with the mean tilt of solar active regions, then the tilt angle falls in quadrant I. We find from Figure 13 that there is a preference for Hale's polarity rule for the emerging azimuthal field by a ratio of 2.4 to 1 in the area. For those pixels satisfying Hale's rule, the mean tilt angle is 75, with an estimated uncertainty of 16. Thus, the super-equipartition emerging fields have a statistically significant mean tilt similar to the active region mean tilt (e.g., Stenflo & Kosovichev 2012).

Zoom In Zoom Out Reset image size