THE SUBSURFACE-SHEAR-SHAPED SOLAR αΩ DYNAMO

It is widely believed that the 11 year sunspot activity is produced and organized by large-scale magnetic fields generated somewhere in the deep convection zone. Most of the solar dynamo models suggest that the toroidal magnetic field that emerges on the surface and forms sunspots is generated near the bottom of the convection zone, in the tachocline or just beneath it in a convection overshoot layer (see, e.g., Choudhuri et al. 1995; Rüdiger & Brandenburg 1995; Dikpati & Charbonneau 1999; Bonanno et al. 2002; Tobias & Weiss 2007). The belief in a deep-seated solar dynamo comes from the fact that this region is sufficiently stable and can store magnetic flux despite the magnetic-flux–tube buoyancy effect (Parker 1975; Spiegel & Weiss 1980; van Ballegooijen 1982; Spruit & Roberts 1983; van Ballegooijen & Choudhuri 1988; Choudhuri 1990). The tachocline represents a strong radial shear of the angular velocity. Yet, turbulent diamagnetism (see, e.g., Zeldovich 1957 or Kitchatinov & Rüdiger 1992) pumps the magnetic fields from the intensively mixed interior of convection zone to its boundaries. This effect can substantially amplify the toroidal magnetic fields near the convection zone boundaries (see, e.g., Krivodubskij 1987; Guerrero & de Gouveia Dal Pino 2008).

However, an attention was drawn to a number of theoretical and observational problems concerning the deep-seated dynamo models (Brandenburg 2005, 2006). A renewed discussion of the place of the solar dynamo can be found, e.g., in papers by Brandenburg (2005) and Tobias & Weiss (2007). In particular, there are some arguments that the subsurface angular velocity shear could play an important role in the dynamo distributed in the convection zone. This shear layer becomes an important ingredient of the flux-transport dynamo models as well (see, e.g., Guerrero & de Gouveia Dal Pino 2008).

In this Letter, we discuss the importance of the surface boundary conditions for the dynamo models, which include the subsurface shear layer. The boundary conditions commonly used in the dynamo models correspond to a perfect conductor at the bottom of the convection zone and vacuum boundary conditions at the top. Both the vacuum and perfect conductor boundary conditions can be regarded as a mathematically convenient idealization. The top boundary conditions play a particularly important role because they control the escape of the dynamo-generated magnetic fields from the Sun.

The perfect conductor boundary condition is usually identified as "closed" (e.g., Choudhuri 1984), because in this case there is no penetration of the generated magnetic flux to the outside. For the axisymmetric magnetic fields all magnetic-field flux is closed inside the dynamo region. The vacuum boundary condition is identified as "open." In this case the poloidal field lines are open to the outside, and the corresponding poloidal magnetic flux "freely" escapes. Also, the strength of the toroidal magnetic field goes smoothly to zero at the boundary. This means that the vacuum boundary condition does not allow the toroidal field to penetrate into the surface. With such boundary condition it is hardly possible to form sunspots from the near-surface large-scale toroidal magnetic fields.

Bearing in mind the dynamical nature of magnetic fields on the solar surface one can model the near-surface behavior by using a combination of the "open" and "closed" types of the boundary conditions. Various consequences of this idea were explored (see, e.g., Choudhuri 1984; Tavakol et al. 1995, 2002; Covas et al. 1998; Kitchatinov & Mazur 1999; Kitchatinov et al. 2000; Käpylä et al. 2010). Here, we apply this approach to a solar dynamo model that extends from the bottom of the convection zone to the top, including the region of the strong subsurface rotational shear. We show that allowing the toroidal magnetic flux to penetrate into the surface brings the butterfly diagram of the toroidal large-scale magnetic field in the upper convection zone and also the phase relations between the different components of the dynamo-generated magnetic field in agreement with solar-cycle observations.

The evolution of the axisymmetric magnetic field (B being the azimuthal component of the magnetic field, A is proportional to the azimuthal component of the vector potential) is governed by the following equations:

$$\begin{equation} \frac{\partial A}{\partial t} = r\sin \theta \mathcal {E}_{\phi }, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial B}{\partial t} = {-}\!\sin \theta \left(\frac{\partial \Omega }{\partial r}\frac{\partial A}{\partial \mu }-\frac{\partial \Omega }{\partial \mu }\frac{\partial A}{\partial r}\right)+\frac{1}{r}\frac{\partial r\mathcal {E}_{\theta }}{\partial r}+\frac{\sin \theta }{r}\frac{\partial \mathcal {E}_{r}}{\partial \mu }, \end{equation} \tag{ 2 }$$

where

$$\begin{eqnarray} r\sin \theta \mathcal {E}_{\phi } & = & \psi _{\eta }\eta _{T}\big\lbrace f_{2}^{(d)}+2f_{1}^{(a)}\big\rbrace \frac{\partial ^{2}A}{\partial r^{2}}\nonumber\\ &&+\psi _{\eta }\eta _{T}\big\lbrace f_{2}^{(d)}+2f_{1}^{(a)}\big\rbrace \frac{(1-\mu ^{2})}{r^{2}}\frac{\partial ^{2}A}{\partial \mu ^{2}}\nonumber\\ && + \psi _{\alpha }C_{\alpha }\eta _{T}Gf_{10}^{(a)}r\mu \sin \theta f(\theta)B\nonumber \\ && + G\!\left(\!2f_{1}^{(a)}\mu \sin ^{2}\theta \frac{\partial A}{\partial \mu }{-}\big(f_{3}^{(a)}-f_{1}^{(a)}\sin ^{2}\theta \big)\frac{\partial A}{\partial r}\right)\nonumber\\ \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \mathcal {E}_{r} & = & \frac{\eta _{T}\psi _{\eta }}{r}\left\lbrace \big(f_{2}^{(d)}+2f_{1}^{(a)}(1-\mu ^{2})\big)\frac{\partial \sin \theta B}{\partial \mu }\right.\nonumber \\ && + \left.2f_{1}^{(a)}\mu \sin \theta \frac{\partial }{\partial r}rB+2rGf_{1}^{(a)}\mu \sin \theta B\right\rbrace \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} r\mathcal {E}_{\theta } & = & \eta _{T}\psi _{\eta }\bigg\lbrace \big(f_{2}^{(d)}+2f_{1}^{(a)}\mu ^{2}\big)\frac{\partial rB}{\partial r}\nonumber\\ && + 2f_{1}^{(a)}\mu \sin \theta \frac{\partial }{\partial \mu }(\sin \theta B)\nonumber \\ && - rG\big(f_{1}^{(a)}\cos 2\theta +f_{3}^{(a)}\big)B\bigg\rbrace. \end{eqnarray} \tag{ 5 }$$

These equations are similar to those used by Pipin & Seehafer (2009) and Seehafer & Pipin (2009). We use the same notations for the functions and parameters as in the paper of Pipin (2008, hereafter P08). Here, G = ∂_rlog ρ is the density stratification scale. Functions f^(a,d)_1,2,3,10 depend on the Coriolis number Ω* = 2τ_cΩ₀; functions ψ_η,α describe magnetic quenching and depend on $\beta =B/\sqrt{\mu _{0}\rho \bar{u^{2}}}$. For reference, these functions are given in the Appendix. The parameter C_α controls the strength of the α-effect. In the presented model the α-effect is distributed in the bulk of the convection zone. For a more clear demonstration of the boundary condition impact, we confine the α-affect in a low-latitude region where the radial gradient of the angular velocity is positive in the most part of the solar convection zone. Similarly to Dikpati et al. (2004), we specify the confinement function:

$$\begin{equation} f(\theta)=(1+e^{30(|\theta -\pi /2|-\pi /6)})^{-1}. \end{equation} \tag{ 6 }$$

In the radial direction the α-effect depends on the density stratification, G, and function of the Coriolis number f^(a)₁₀(Ω*). We introduce parameter C_η to control the turbulent diffusion coefficient, η_T = C_ηη⁽⁰⁾_T, where $\eta _{T}^{(0)}=\tau _{c}\bar{u^{2}}/3$. The internal parameters of the solar convection zone are given by Stix (2002). At the top of the solar convection zone the stratification is strongly deviates from adiabatic, and also the turbulence parameters vary sharply. For this reason we confine the integration domain between 0.71 R_☉ and 0.972 R_☉ in radius, and it extends from pole to pole in latitude. The differential rotation profile, Ω = Ω₀f_Ω(x, μ) (shown in Figure 1(a)) is a slightly modified version of the analytical approximation proposed by Antia et al. (1998):

$$\begin{eqnarray} f_{\Omega }(x,\mu) & = & \frac{1}{\Omega _{0}} [\Omega _{0}+55\left(x-0.7)\phi (x,x_{0})\phi (-x,-0.96)\right.\nonumber\\ && - 200(x-0.95)\phi (x,0.96)\nonumber \\ &&+ \left(21P_{3}\left(\mu \right)+3P_{5}\left(\mu \right)\right]\left(\frac{\mu ^{2}}{j_{p}\left(x\right)} +\frac{1-\mu ^{2}}{j_{e}\left(x\right)}\right)\nonumber \\ j_{p} & = & \frac{1}{1+\exp \left(\frac{0.709-x}{0.02}\right)}, \quad j_{e}= \frac{1}{1+\exp \left(\frac{0.692-x}{0.01}\right)}, \end{eqnarray} \tag{ 7 }$$

where Ω₀ = 2.87 ×  10⁻⁶{rms⁻¹ is the equatorial angular velocity of the Sun at the surface, x = r/R_☉, ϕ(x, x₀) = 0.5[1 + tanh [100(x − x₀)]], x₀ = 0.71. The distribution of the Coriolis number, turbulent diffusivity, and the rms convection velocity are shown in Figure 1(b). The radial profile of the α-effect is shown in Figure 1(c).

Zoom In Zoom Out Reset image size

At the bottom of the integration domain we apply the perfect conductor ("closed") boundary conditions: $\mathcal {E}_{\theta }=0,\, A=0$. The boundary conditions at the top are defined as the following. Bearing in mind the idea of the partial escape of the toroidal flux from the Sun discussed in Introduction, we explore a combination of the "open" and "closed" boundary conditions at the top, controlled by a parameter δ. For the toroidal field we use condition

$$\begin{equation} \delta \frac{\eta _{T}}{r_{e}}B+\left(1-\delta \right)\mathcal {E}_{\theta } = 0. \end{equation} \tag{ 8 }$$

This is similar to the boundary condition discussed by Kitchatinov et al. (2000). For the poloidal field, we apply a combination of the local condition A = 0 and condition of smooth transition from the internal poloidal field to the external potential (vacuum) field:

$$\begin{equation} \delta \left(\left.\frac{\partial A}{\partial r}\right|_{r=r_{e}} -\left.\frac{\partial A^{(\rm vac)}}{\partial r}\right|_{r=r_{e}}\right)+\left(1-\delta \right) A= 0, \end{equation} \tag{ 9 }$$

where the external potential field is

$$\begin{equation} A^{(\rm vac)}\left(r,\mu \right)= \sum a_{n}\left(\frac{r_{e}}{r}\right)^{n}\sqrt{1-\mu ^{2}}P_{n}^{1}\left(\mu \right), \end{equation} \tag{ 10 }$$

P¹_n(μ) is the associated Legendre polynomial of degree n. For the numerical implementation of Equation (9), we take a one-side finite difference approximation for the radial derivative at an angular mesh point μ_j:

$$\displaystyle \left.\frac{\partial A_{j}}{\partial r}\right|_{r=r_{e}} =\frac{3A_{N\, j}-4A_{N-1\, j}+A_{N-2\, j}}{2h_{r}},$$

where h_r is the radial discretization interval, and consider the expansion (10) at the top boundary r = r_e: $A_{N\, j}=\sum a_{n}\sqrt{\vphantom{A^A}\smash{\hbox{${1-\mu _{j}^{2}}$}}}P_{n}^{1}\left(\mu _{j}\right)$. Then, we define matrices $M_{nj}^{\left(a\right)}=\sqrt{1-\mu _{j}^{2}}P_{n}^{1}\left(\mu _{j}\right)$ and $\tilde{M}_{jn}=n\sqrt{\vphantom{A^A}\smash{\hbox{${1-\mu _{j}^{2}}$}}}P_{n}^{1}\left(\mu _{j}\right)$, with μ_j being collocation points of P¹_n (see Boyd 2001). This procedure allows us to express coefficients a_n in Equation (10) via the values of potential A at the grid points: a_n = M^(a)−1_njA_N j. Substituting this in Equation (9) and solving it we get

$$A_{N\, j}=\delta \big(3E_{jk}+2h_{r}\tilde{M}_{jk}M_{nk}^{(a)-1}\big)^{-1}(4A_{N-1\, jk}-A_{N-2\, k}),$$

where E is a unit diagonal matrix.

Parameter δ in the boundary conditions describes a transition between the "closed" (δ = 0) to "open" (δ = 1) boundaries. Physically, it controls penetration of the dynamo-generated fields into the outer atmosphere.

Decreasing δ in Equations (8) and (9) results in stronger tangential and weaker radial large-scale magnetic fields at the surface. While the strong toroidal magnetic field is a desired feature of the model, the weak radial magnetic field decreases the efficiency of the radial subsurface shear to produce large-scale toroidal magnetic fields. In fact, the simulations reveal that the critical dynamo number, C_α, is greater when the penetration parameter δ is smaller. For this reason, we consider the case of a small deviation from the vacuum ("open") boundary conditions. Moreover, in order to match the dynamo period to the solar cycle we choose the magnetic diffusivity parameter C_η = 0.05, which is significantly lower than the value predicted by the mixing length theory. To demonstrate the effect of the new boundary conditions with the field penetration we show for comparison in Figures 2 and 3 the results of two runs for δ = 0.95 (corresponding to a partial penetration of toroidal field) and δ = 1 (the vacuum boundary conditions).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results show that allowing the large-scale toroidal magnetic field to penetrate into the surface layers of the Sun changes the direction of the latitudinal migration of the toroidal field activity and produces the magnetic butterfly diagram in a good qualitative agreement with solar-cycle observations. The dynamo wave penetrates into the surface and propagates along the iso-surface of angular velocity in the subsurface shear layer. This is in agreement with the Yoshimura rule (Yoshimura 1975).

Dikpati et al. (2002) explored generation of toroidal magnetic fields by the Ω-effect in the subsurface shear layer in the Babcock–Leighton-type dynamo models. They found that the phase relation between the subsurface toroidal magnetic field and the surface radial magnetic field is inconsistent with observations. We believe that this inconsistency was due to the fact that in their model the source of the surface poloidal magnetic fields was related with the bottom of convection zone. Therefore, the subsurface toroidal field that is generated in the subsurface shear layer does not contribute directly to the generation of the poloidal magnetic field.

Both of our simulation runs shown in Figures 2 and 3 were started with initial conditions representing equally mixed symmetrical and antisymmetrical (relative to the equator) magnetic field components. The evolution retains only the dipole-like parity configurations in both cases though the relaxation time in the penetration case of δ = 0.95 is much longer than in the case of the pure vacuum boundary conditions. If we relax the confinement of the α-effect in latitude, i.e., f(θ) = 1 instead of Equation (6), the general patterns of Figure 2 are hold except that the maximum of the toroidal magnetic field is shifted to higher latitude ≈40°. Therefore, we can conclude that the αΩ-dynamo model with the boundary conditions that allow a small partial penetration of the toroidal field into the outer layers of the Sun can robustly reproduce the solar-cycle butterfly diagram for the near-surface large-scale magnetic-field evolution. These results demonstrate the importance of the subsurface rotational shear layer in the solar dynamo mechanism.

The authors thank F. Busse, N. Mansour, and M. Stix for useful comments. This work was supported by the NASA LWS NNX09AJ85G grant and partially by the RFBR grant 10-02-00148-a.

Here, we give the definitions of the functions which were used in the model. The given functions describe the efficiency of the Coriolis force and the mean magnetic field to act on the stratified turbulence and to produce the dynamo α-effect, anisotropy of magnetic diffusion, turbulent magnetic pumping, magnetic quenching of the turbulent effects, etc. These effects are discussed in Pipin (2008). Functions f^(a,d)_1,2,3,10 depend on the Coriolis number Ω* = 2τ_cΩ₀; functions ψ_η,α describe magnetic quenching and depend on $\beta =B/\ssty\sqrt{\mu _{0}\rho \bar{u^{2}}}$:

$$\begin{eqnarray*} f_{1}^{(a)} & = & \frac{3}{4\Omega ^{*\,2}}\left((\Omega ^{*2}+3)\frac{\arctan \Omega ^{*}}{\Omega ^{*}}-3\right),\\ f_{2}^{(d)} & = & \frac{3}{4\Omega ^{*\,2}}\left(((\varepsilon -1)\Omega ^{*\,2}+3\varepsilon +1)\frac{\arctan (\Omega ^{*})}{\Omega ^{*}}-(3\varepsilon +1)\right),\\ f_{3}^{(a)} & = & \frac{3}{4\Omega ^{*\,2}}\left(\left(\left(\varepsilon -1\right)\Omega ^{*\,2}+\varepsilon -3\right)\frac{\arctan \Omega ^{*}}{\Omega ^{*}}+3-\varepsilon \right),\\ f_{10}^{(a)} & = & {-}\frac{1}{3\Omega ^{*\,3}}\bigg(3(\Omega ^{*2}+1)(\Omega ^{*2}+\varepsilon -1)\frac{\arctan \Omega ^{*}}{\Omega ^{*}}\nonumber\\ &&-((2\varepsilon +1)\Omega ^{*2}+3\varepsilon -3)\bigg),\\ \psi _{\eta } & = & \frac{3}{16\beta ^{2}}\left(((4\beta ^{2}+3)\varepsilon +4\beta ^{2}-1)\frac{\arctan \left(2\beta \right)}{2\beta }+1-3\varepsilon \right)\\ \psi _{\alpha } & = & \frac{5}{128\beta ^{4}}\left(16\beta ^{2}-3-3(4\beta ^{2}-1)\frac{\arctan \left(2\beta \right)}{2\beta }\right). \end{eqnarray*}$$

The parameter ε measures the ratio between the turbulent energies of the kinetic and magnetic fluctuations $\varepsilon =\bar{b^{2}}/ \mu _{0}\rho \bar{u^{2}}$, in the background turbulence (in the absence of the mean fields). Note, in notation of Pipin (2008), the turbulent diffusivity and α-effect quenching functions are defined as follows: ψ_η = ϕ₃ + ϕ₂ − 2ϕ₁ and ψ_α = −3/4ϕ^(a)₆, respectively. Expressions for ϕ_1,2,3 and ϕ^(a)₆ are given in Pipin (2008).