STELLAR DYNAMO MODELS WITH PROMINENT SURFACE TOROIDAL FIELDS

One of the most compelling problems in modern dynamo theory is the formulation of a realistic coupling between the internal magnetic field and the external field in the atmosphere. In fact, the boundary conditions for the electric and the magnetic field at the stellar surface put severe constraints on the allowed coronal field configurations, and it is often necessary to resort to very crude approximations for the latter.

The standard textbook boundary condition employed in mean-field dynamo theory considers a current-free field in the region r ≥ R, where R is the stellar radius, so that ∇ × B = 0 in this domain. Although in the solar case this assumption is motivated by the possibility of describing the almost rigid rotation of the coronal holes in the lower corona (Nash et al. 1988), it be might incorrect to extend its validity in more active stars.

In fact, as current-free fields represent the states of minimum energy under the constraint that the normal component of the field at the photosphere is fixed, they cannot provide the additional energy required to sustain a significant activity level. Recent measurements of Faraday rotation in the solar corona support the evidence for large-scale coronal currents (Spangler 2007), an essential ingredient to explain coronal heating in terms of Joule dissipation. On the other hand, on much smaller scales, the presence of currents is unavoidable in order to explain the twisted field structure of filaments and prominences.

Force-free magnetic fields, defined by ∇ × B = α_ff(x)B, where α_ff(x) is a scalar function, can be more appealing from the physical point of view, at least for very low plasma-β values. However, recent investigations based on direct numerical simulations have shown that the free magnetic energy and the efficiency of coronal heating via current dissipation are still very limited in these models (Peter et al. 2015).

An important consequence of current-free boundary conditions is that the toroidal field must identically vanish in the all domain r ≥ R. In order to illustrate this point in detail, it is convenient to introduce the scalar potential ${\rm{\Phi }}(r,\vartheta ,\varphi )$ so that the toroidal field can be written as

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{{\rm{T}}}=-\displaystyle \frac{1}{\sin \theta }\displaystyle \frac{\partial {\rm{\Phi }}}{\partial \varphi }{{\boldsymbol{e}}}_{\vartheta }+\displaystyle \frac{\partial {\rm{\Phi }}}{\partial \vartheta }{{\boldsymbol{e}}}_{\varphi }.\end{eqnarray} \tag{ 1 }$$

It is not difficult to show that if Φ is a single valued function, ${\rm{\nabla }}\times {{\boldsymbol{B}}}_{{\rm{T}}}=0$ in a volume always necessarily implies ${{\boldsymbol{B}}}_{{\rm{T}}}\equiv 0$ everywhere (Krause & Raedler 1980). In particular, if the field configuration is spherically symmetric, the azimuthal component of the magnetic field must vanish at the surface, so that $\tfrac{\partial {\rm{\Phi }}}{\partial \vartheta }{| }_{r=R}=0$, as imposed in most of the dynamo models (see Moss & Sokoloff 2009 for an interesting discussion on this issue).

Spectropolarimetric observations of photospheric magnetic fields in solar-like stars have revealed surface toroidal field which are mostly axisymmetric and have a predominant toroidal component (Petit et al. 2005, 2008; Fares et al. 2010, 2013; See et al. 2015) implying that most of the magnetic energy resides in the toroidal field. This is the case of the solar-like stars like HD72905, with 82% of the magnetic energy stored in the toroidal field which is nearly completely axisymmetric (97%), or of the G8 dwarf ξ Boo A with 81% of toroidal energy of which 97% is due to the axisymmetric component, or the case of HD56124 with 90% of the energy in axisymmetric field configurations, and roughly the same strength of poloidal and toroidal component. The situation is even more dramatic if one considers M dwarfs like WX UMa or AD Leo where nearly all the energy is stored in an axisymmetric field with prominent photospheric non-zero toroidal component (see See et al. 2015 for a detailed discussion).

While in the case of the Sun a similar observational strategy has confirmed that the magnetic energy is mostly (>90%) poloidal (Vidotto 2016), it is clear that the boundary conditions based on current-free coronal field might not be correct for other, more active stars. In fact, it has been further noticed that fast-rotating, low-mass stars have on average stronger surface toroidal fields than solar-mass slow rotators (See et al. 2016). Can this fact be explained as an enhanced dynamo action due a αΩ mechanism? Indeed, as strong toroidal fields can alter the average atmospheric structure, there is no reason to assume that a current-free field is still a reasonable approximation for discussing the magnetic energy budget in these objects.

The idea proposed in this Letter is that on large scales and on timescales of the order of the stellar cycle, the external field can be considered harmonic so that

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\boldsymbol{B}}+{k}^{2}{\boldsymbol{B}}=0,\end{eqnarray} \tag{ 2 }$$

where the wavenumber k, assumed to be real, determines the characteristic spatial length of the field.

As is well known, in the Sun, various MHD instabilities trigger multiple modes harmonics on different length scales (Nakariakov & Verwichte 2005; Arregui et al. 2013). In more active stars, one can thus argue that the dynamo waves from the interior trigger global oscillations of the coronal field that then propagate according to (2).

It is important to note that linear force-free fields with ${\alpha }_{{ff}}^{2}={k}^{2}$ are solutions of (2) (note that strictly speaking α_ff is a pseudoscalar, while k is a scalar), although the converse is not true in general (Chandrasekhar & Kendall 1957). In this investigation, it turns out that the coronal fields obtained by coupling (2) with dynamo solutions in the interior are approximately force-free in the sense of Warnecke & Brandenburg (2010) as $\langle {({\boldsymbol{J}}\times {\boldsymbol{B}})}^{2}\rangle \ll \langle {{\boldsymbol{B}}}^{2}\rangle \langle {{\boldsymbol{J}}}^{2}\rangle$.¹

It will then be shown that standard MHD continuity conditions at the surface naturally allow for non-zero surface toroidal fields that are directly extrapolated from the interior field produced by the dynamo. In this Letter, several solutions obtained for solar-like stars will therefore be discussed and analyzed. It turns out that depending on the strength of the differential rotation and the surface meridional circulation, it is possible to obtain surface fields whose energies distribution between toroidal and poloidal components is consistent with the observations.

Let us assume that the field periodically evolves with a characteristic cycle frequency ω so that ${\boldsymbol{B}}={e}^{-i\omega t}{\boldsymbol{ \mathcal B }}$. Equation (2) thus reads

$$\begin{eqnarray}&&({{\rm{\nabla }}}^{2}+{k}^{2}){\boldsymbol{ \mathcal B }}=0,\end{eqnarray} \tag{ 3 }$$

where the wavenumber k is considered to be real. Clearly we must assume that Rk ≪ 1, otherwise the typical spatial structure of the field would be too short to be consistent with our quasi-homogeneous approximation. As we shall see, as long as Rk ≪ 1, our results are not quantitatively dependent on the value of k.

As usual, the boundary conditions are the continuity of the normal component of the magnetic field and the tangential component of the electric field across the stellar surface:

$$\begin{eqnarray}&&[[{\boldsymbol{n}}\cdot {\boldsymbol{B}}]]=0,\ \ [[{\boldsymbol{n}}\times {\boldsymbol{E}}]]=0,\end{eqnarray} \tag{ 4 }$$

where n is the normal to the surface. In spherical symmetry, the following decomposition for the magnetic field ${\boldsymbol{ \mathcal B }}$ can be used

$$\begin{eqnarray}&&{\boldsymbol{ \mathcal B }}=-{\boldsymbol{r}}\times {\rm{\nabla }}{\rm{\Psi }}-{\rm{\nabla }}\times ({\boldsymbol{r}}\times {\rm{\nabla }}{\rm{\Phi }})\equiv {{\boldsymbol{ \mathcal B }}}_{{\rm{T}}}+{{\boldsymbol{ \mathcal B }}}_{{\rm{P}}},\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Psi }}={\rm{\Psi }}(r,\vartheta ,\varphi )$ and ${\rm{\Phi }}={\rm{\Phi }}(r,\vartheta ,\varphi )$ are scalar functions, ${{\boldsymbol{ \mathcal B }}}_{{\rm{T}}}$ is the toroidal component, and ${{\boldsymbol{ \mathcal B }}}_{{\rm{P}}}$ is the poloidal one (see Krause & Raedler 1980 for details). The vector Helmholtz (Equation (3)) decouples in the two scalar equations

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}+{k}^{2}{\rm{\Phi }}=0\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Psi }}+{k}^{2}{\rm{\Psi }}=0\end{eqnarray} \tag{ 6b }$$

and variable separation in (6) gives

$$\begin{eqnarray}&&{\rm{\Phi }}=R\displaystyle \sum _{n}^{\infty }[{A}_{n}{j}_{n}(\xi x)+{B}_{n}{y}_{n}(\xi x)]{P}_{n}(\cos \theta )\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \sum _{n}^{\infty }[{C}_{n}{j}_{n}(\xi x)+{D}_{n}{y}_{n}(\xi x)]{P}_{n}(\cos \theta ),\end{eqnarray} \tag{ 7b }$$

where x is the normalized stellar radius x = r/R, ξ = kR, ${P}_{m}(\cos \theta )$ are the Legendre polynomials, and j_n(x) and y_n(x) are the spherical Bessel functions. It should be noted that k = 0 in (6) does not necessarily imply Ψ ≡ 0 in the all r ≥ R domain, as it must instead hold in the case of a current-free field (vacuum boundary condition).

The A_n–D_n constants are complex numbers whose value must be determined by the boundary conditions (4) imposed at the stellar surface and at some finite outer radius $r={R}_{\mathrm{out}}$ where a transition to a wind-dominated field topology occurs (Parker 1958). For our purposes it will be sufficient to assume that at $r={R}_{\mathrm{out}}$ the solution is radially dominated so that²

$$\begin{eqnarray}&&{B}_{\theta }(r={R}_{\mathrm{out}})={B}_{\phi }(r={R}_{\mathrm{out}})=0.\end{eqnarray} \tag{ 8 }$$

In the stellar interior, the field is described by the mean-field dynamo equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{U}}\times {\boldsymbol{B}}+\alpha {\boldsymbol{B}})-{\rm{\nabla }}\times (\eta {\rm{\nabla }}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 9 }$$

where, as usual, α is a pseudoscalar function representing the turbulent α-effect, U is the mean flow, and η is the turbulent (eddy) diffusivity. As before ${\boldsymbol{B}}={e}^{-i\omega t}{\boldsymbol{ \mathcal B }}$, and we apply the fundamental decomposition (5) for ${\boldsymbol{ \mathcal B }}$. Finally, we write

$$\begin{eqnarray}&&{\rm{\Phi }}=R\displaystyle \sum _{n}^{\infty }{\phi }_{n}(x){P}_{n}(\cos \theta )\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \sum _{n}^{\infty }{\psi }_{n}(x){P}_{n}(\cos \theta ),\end{eqnarray} \tag{ 10b }$$

where ψ_n and ϕ_n are the complex eigenfunctions of antisymmetric parity of the linear operator (9) (see Rädler 1973 for details).

The field at the inner boundary is assumed to be a perfect conductor, but the boundary conditions (4) at the surface imply the continuity of ϕ_n, ψ_n, and their derivatives across the boundary. It is not difficult to show that in order for the interior solution to be consistent with the external field in (7) the following relation must hold at x = 1:

$$\begin{eqnarray}&&\displaystyle \frac{d{\phi }_{n}}{{dx}}+\left[\displaystyle \frac{{\gamma }_{n}{y}_{n+3/2}(\xi )+{j}_{n+3/2}(\xi )}{{\gamma }_{n}{y}_{n+1/2}(\xi )+{j}_{n+1/2}(\xi )}-n\right]{\phi }_{n}=0,\end{eqnarray} \tag{ 11 }$$

where γ_n = A_n/B_n are determined by imposing the outer boundary condition (8) on $r={R}_{\mathrm{out}}$. A similar equation can be obtained for the ψ_n components.

Let us specialize our formalism to the case of a solar-like star. We assume the following form for the velocity field ${\boldsymbol{U}}={\boldsymbol{u}}(r,\theta )+r\sin \theta {\rm{\Omega }}(r,\theta ){{\boldsymbol{e}}}_{\phi }$, and we model the solar-like differential rotation in the following way:

$$\begin{eqnarray}&&{\rm{\Omega }}(r,\theta )={{\rm{\Omega }}}_{c}+\displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(\displaystyle \frac{x-{x}_{c}}{d}\right)\right]({{\rm{\Omega }}}_{s}(\theta )-{{\rm{\Omega }}}_{c}),\end{eqnarray} \tag{ 12 }$$

where x_c is the location of the convection zone in units of the stellar radius, d = 0.02, Ω_c is the uniform angular velocity of the radiative core, ${{\rm{\Omega }}}_{s}(\theta )={{\rm{\Omega }}}_{0}-d{\rm{\Omega }}{\cos }^{2}\theta$. Here, $d{\rm{\Omega }}={{\rm{\Omega }}}_{{eq}}-{{\rm{\Omega }}}_{c}$ is the surface differential rotation (in principle, obtained from observations). For actual calculations, we fixed x_c = 0.7 and Ω_c/Ω_eq = 0.9, but $d{\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{eq}}$ is allowed to vary. The radial profile of the turbulent diffusivity is assumed to be the following:

$$\begin{eqnarray}&&\eta ={\eta }_{c}+\displaystyle \frac{1}{2}({\eta }_{t}-{\eta }_{c})\left[1+\mathrm{erf}\left(\displaystyle \frac{x-{x}_{x}}{d}\right)\right],\end{eqnarray} \tag{ 13 }$$

with η_c/η_t = 10⁻¹. The α effect is proportional to $\cos \vartheta$, and its radial profile is assumed to be uniformly distributed in all the convection zone (see Bonanno et al. 2002 for details). The meridional circulation u is obtained from the stream function S(r), which as explained in Bonanno (2013), S(r) can be obtained from an underlying stellar model. The flow as usual is equatorward at the equator and poleward at the surface.

As usual, let us introduce the following dynamo numbers, ${C}_{{\rm{\Omega }}}={R}^{2}{{\rm{\Omega }}}_{0}/{\eta }_{t}$, ${C}_{\alpha }=R{\alpha }_{0}/{\eta }_{t}$, ${C}_{u}={{RU}}_{0}/{\eta }_{t}$, ${C}_{\omega }={R}^{2}\omega /{\eta }_{t}$, where Ω₀ is the rotation rate at the equator and U₀ is the maximum strength of u_θ at the bottom of the convection zone. The resulting eigenvalue problem can be conveniently solved by inverting a block-diagonal complex matrix (Rädler 1973; Bonanno 2013) that must be truncated to the desired numerical accuracy. For calculations we used the kinematic dynamo code CTDYN, developed by the author, and extensively tested in Jouve et al. (2010). Our results are summarized in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Summary of the Numerical Simulations

${ \mathcal M }$	${C}_{{\rm{\Omega }}}$	${C}_{u}$	$d{\rm{\Omega }}/{{\rm{\Omega }}}_{\mathrm{eq}}$	kR	${C}_{\alpha }$	${C}_{\omega }$	${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$
A	1000	0	0.1	0.1	7.87	$\infty$	0.67
B	1000	0	0.1	0.4	7.88	$\infty$	0.67
C	1000	0	0.3	0.1	11.24	$\infty$	0.67
D	1000	100	0.3	0.1	7.78	$\infty$	0.37
E	2000	0	0.1	0.1	9.49	$\infty$	0.63
F	2000	0	0.3	0.1	10.97	79.02	0.91
G	2000	0	0.3	vac	11.94	69.98	0
J	4000	0	0.1	0.1	10.11	62.93	0.85
K	4000	0	0.1	0.5	10.04	63.03	0.86
L	4000	100	0.1	0.1	5.57	$\infty$	0.68
M	4000	200	0.1	0.1	5.13	$\infty$	0.37

Note. In particular, it shows how ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$ depends on ${C}_{{\rm{\Omega }}}$ and various other input parameters ("vac" stands for current-free boundary condition). ${ \mathcal M }$ is the model name, and in particular for model J the toroidal and poloidal field are displayed in Figure 1 at various values of the cycle phase.

Download table as:  ASCIITypeset image

Let us first stress that the ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$ at the surface is non-zero, and it is now an increasing function of ${C}_{{\rm{\Omega }}}$, an effect clearly expected for an αΩ dynamo. The surface differential rotation also plays an important role because it strongly influences the value of ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$, as it can be deduced by looking at models E and F, for instance. On the contrary, it should also be noted that ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$ decreased as the flow is increased (see models J, L, and M, for instance). This fact has a clear physical interpretation: as the flow is poleward below the surface, the toroidal belts become more and more confined below the surface at high ${C}_{u}$. It is also reassuring to notice that as long as kR ≪ 1, our results are not sensitive to the choice of ξ. On the other hand, with the new boundary conditions, the value of the critical ${C}_{\alpha }$ is in general smaller than with the standard current-free boundary conditions, as can be observed by looking at models F and G. We also verified that our results are not qualitatively dependent on the location of the external boundary ${R}_{\mathrm{out}}$. Moreover, low-mass stars with a more extended convection zone have in general higher values of ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$, if ${C}_{{\rm{\Omega }}}$ is large enough, as discussed in the model of Figure 2.

Zoom In Zoom Out Reset image size

If the dynamo action is instead driven by the meridional circulation (as in mean-field models of the solar dynamo), the new boundary condition does not alter the internal dynamo action because the toroidal field is localized at the bottom of the convection zone.

In the case of a generic αΩ advection-dominated dynamo action, we observe a change of sign of the dimensionless ratio ${\boldsymbol{J}}\cdot {\boldsymbol{B}}/\langle {{\boldsymbol{B}}}^{2}\rangle$ from positive (in the northern hemisphere) in the turbulent zone in the interior to negative in the exterior. The opposite happens in the southern hemisphere as one can see in Figure 3. A similar change of sign between the turbulent zone and the exterior has also been observed in Warnecke et al. (2011), although in our mean-field models the current helicity is mostly negative in the outer layers in the northern hemisphere.

Zoom In Zoom Out Reset image size

The description of the external field in terms of solution of the Helmholtz equation allowed us to extrapolate the internal toroidal field generated by the dynamo on the photosphere, and finally to make contact with the observations. The assumption beyond this idea is the possibility of treating the corona as an external passive medium with effective macroscopic dielectric properties, if averaged over long enough timescales. Although this approach oversimplifies the complex physics of the corona, in our opinion, it represents a significant improvement of the current-free boundary conditions for which ${{\boldsymbol{B}}}_{{\rm{T}}}\equiv 0$ on the surface, at least in some classes of very active stars. The resulting dynamo numbers are in general smaller than the standard critical dynamo numbers; the ratio ${E}_{{\rm{T}}}/{E}_{\mathrm{tot}}$ increases with ${C}_{{\rm{\Omega }}}$ and with a more extended convection zone and decreases with ${C}_{u}$. Indeed, fast-rotating stars with a larger convection zone should approach a cylindrical rotation law in the interior with a smaller surface meridional circulation.

One can therefore argue that the general increase of the surface toroidal energy in low-mass fast-rotating stars finds its natural explanation in an underlying αΩ dynamo mechanism.

A detailed study of all the parameter space and a comparison with the global topologies inferred from observations will be discussed in a longer paper. We also plan to extend this investigation including non-axisymmetric solutions of higher azimuthal modes.

I would like to thank the colleagues of the MHD group of the Leibniz Institute for Astrophysics in Potsdam for important comments and hospitality. I am also indebted to Rim Fares for clarifications on Zeeman–Doppler imaging and on magnetic field topology reconstruction and to the anonymous referee for the constructive criticism.