TOROIDAL FIELD REVERSALS AND THE AXISYMMETRIC TAYLER INSTABILITY

We solve the MHD equations in axisymmetric spherical coordinates in the radial range from 0.15 R_☉ to 0.95 R_☉. The radiative region extends from 0.15 R_☉ to 0.71 R_☉ and the convection zone extends from 0.71 R_☉ to 0.90 R_☉ and we impose an additional stably stratified region from 0.90 R_☉ to 0.95 R_☉. The reference state thermodynamic variables are taken from a polynomial fit to the standard solar model from Christensen-Dalsgaard et al. (1996). Therefore, the density ranges from approximately 60 gm cm⁻³ at the bottom of the radiation zone to 10⁻² gm cm⁻³ at the top of the domain. The stable stratification is that given by the standard solar model and the superadiabaticity of the convection zone is set to 10⁻⁷. The grid resolution is 512 latitudinal grid points by 1500 radial grid points, with 500 radial zones dedicated to the radiative interior and 1000 zones dedicated to the tachocline and convection zone.

Because this model is not fully three dimensional (3D), we cannot reproduce the Reynolds stresses required to set up the observed differential rotation (Thompson et al. 2003) in the convection zone and tachocline. We therefore impose this differential rotation through a forcing term on the azimuthal component of the momentum equation:

$$\begin{eqnarray} {\frac{\partial {}}{\partial {t}}}\left(\frac{U_{\phi }}{r \sin \theta }\right)&=&\,{-}\frac{(\nabla \cdot {\bf u}{U_{\phi }})_{\phi }}{r\sin \theta }+\frac{((\nabla \times {\bf B})\times {\bf B})_{\phi }}{\bar{\rho }\mu r\sin \theta }\nonumber \\ & &+\frac{(2{\bf u}\times \Omega)_{\phi }}{r\sin \theta }+\frac{\nu }{r\sin \theta }\nabla ^{2}u_{\phi }+\frac{\nu u_{\phi }}{r^{3}\sin ^{3}\theta }\nonumber\\ &&+\frac{1}{\tau }\left(\Omega ^{\prime }(r,\theta)-\frac{u_{\phi }(r,\theta)}{r\sin \theta }\right), \end{eqnarray} \tag{ 1 }$$

where the last term on the right-hand side represents the forcing term. Ω'(r, θ) is the prescribed (observed) differential rotation relative to the constant Ω and τ is the e-folding time of the forcing which we set to 10⁴.¹ The rotation profile is given by

$$\begin{eqnarray} \Omega ^{\prime }(r,\theta)&=&\Omega \qquad\qquad\qquad\qquad\quad\ \ \ \ \mbox{for $r < 0.67R_{\odot }$}\nonumber\\ \Omega ^{\prime }(r,\theta)&=&(A+B\cos ^{2}\theta +C\cos ^{4}\theta -\Omega)\nonumber\\ &&\times {\rm erf}\left(\frac{2(r-0.67\,R_{\odot })}{0.05\,R_{\odot }}\right) \quad\ \ \ \ \ \mbox{for $0.67<r<0.72$}\nonumber\\ \Omega ^{\prime }(r,\theta)&=&(A+B\cos ^{2}\theta\nonumber\\ && +C\cos ^{4}\theta -\Omega) \qquad\quad\quad\quad\ \ \mbox{for $r>0.72$},\nonumber\\ \end{eqnarray} \tag{ 2 }$$

where Ω is the constant value set to 441 nHz, A is 456 nHz, B is −42 nHz, and C is −72 nHz (Thompson et al. 2003).

Initially, we run a purely hydrodynamic model until convection is sufficiently established at which point we imposed a dipolar field with the form

$$\begin{equation} B_{1}=B_{o}r^{2}\left(1.0-\frac{r}{0.71\, R_{\odot }}\right)^{2}, \end{equation} \tag{ 3 }$$

so that initially all field lines close within the radiative interior, but overlap the bottom of the tachocline. We ran five models, which we label Models A, B, C, D, and E with initial field strengths, B_o, of 10 G, 40 G, 80 G, 200 G, and 4000 G, respectively, with all other parameters the same. Figure 1 shows a schematic of the initial setup and model, while Figure 2 shows the time-averaged azimuthal velocity as forced in Equation (2), along with the time-averaged meridional flow that results from this differential rotation profile (Figure 2(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following we discuss the primary results from Model B (B_o = 40 G) and leave the discussion of the variation in field strength for Section 3.3. Because of the slight overlap of the imposed dipolar poloidal field with the tachocline and because of magnetic diffusion, it does not take long for the dipolar poloidal field to be stretched into toroidal field by the differential rotation in the tachocline and convection zone. Initially, this induces a toroidal field which has two signs each in the northern and southern hemispheres, see Figure 3. This structure can be understood by looking at the azimuthal component of the magnetic induction equation:

$$\begin{eqnarray} {\frac{\partial {B_{\phi }}}{\partial {t}}}&=&rB_{r}{\frac{\partial {}}{\partial {r}}}\left(\frac{u_{\phi }}{r}\right)-ru_{r}{\frac{\partial {}}{\partial {r}}}\left(\frac{B_{\phi }}{r}\right)+\frac{\sin \theta B_{\theta }}{r}{\frac{\partial {}}{\partial {\theta }}}\left(\frac{u_{\phi }}{\sin \theta }\right) \nonumber \\ &-& \frac{\sin \theta u_{\theta }}{r}{\frac{\partial {}}{\partial {\theta }}}\left(\frac{B_{\phi }}{\sin \theta }\right)+B_{\phi }u_{r}h_{\rho }+\eta \nabla ^{2}B_{\phi }-\eta \frac{B_{\phi }}{r^{2}\sin ^{2}\theta }. \nonumber\\ \end{eqnarray} \tag{ 4 }$$

The structure of the differential rotation gives ∂u_ϕ/∂r negative at high latitudes and positive at low latitudes and the first term on the right-hand side of Equation (4) induces oppositely signed toroidal field at high and low latitudes, with asymmetry about the equator. This is seen in Figure 3(a).

Zoom In Zoom Out Reset image size

However, the first term on the right-hand side of Equation (4), induction due to radial gradient in the rotation rate, is dominant for only a short time. The latitudinal toroidal field gradient becomes unstable (see Figure 6, Section 3.2) and the dominant induction term becomes the fourth term on the right-hand side of Equation (4), advection of toroidal field by latitudinal flow. This leads to a single dominant signed toroidal field in each hemisphere, asymmetric about the equator. In the following I will refer to positive toroidal field (white in Figure 3) as the dominant sign in the Northern hemisphere, and negative toroidal field (black in Figure 3) as the minority field sign in the Northern hemisphere and conversely, in the Southern hemisphere. These particular signs arise naturally from the equatorward directed u_θ coupled with a toroidal field which is stronger at the poles than the equator. Subsequently, and for most of the simulation, the advection of toroidal field by meridional flow (fourth term on the right-hand side of Equation (4)) is dominant within the tachocline (Figure 3(b)). As can be seen from Figures 3(b) and (c), a "reversal" takes place in which the dominant field is replaced by the minority field in the bulk of the convection zone and tachocline. This reversal proceeds by the advection of minority signed field toward the equator by meridional circulation at the base of the convection zone.

The evolution in time of the reversal process is better demonstrated in a time–latitude plot, as seen in Figure 4(a). Shown there is the toroidal field amplitude as a function of time and latitude. One can see the initial phase of two signs in each hemisphere giving way to a dominant sign in each hemisphere. Reversals subsequent to that shown in Figure 3 are also seen here. During the first reversal (around 7 × 10⁷s, or around 30 rotation periods), the northern and southern hemispheres are relatively in phase, however later reversals are stronger in the southern hemisphere and the coherence between north and south is not strictly maintained. As is seen in Figure 4, these reversals are seen both in the tachocline and at the top of the convection zone. We note that no reversal occurs in the radiative interior, nor does any reversal occur in the poloidal field.

Zoom In Zoom Out Reset image size

Except for the very initial phases, the dominant term in the tachocline responsible for toroidal field evolution is the advection of toroidal field by the meridional circulation. The time-averaged equatorward latitudinal velocity at the base of the convection zone is ≈3 × 10³ cm s⁻¹, which would imply an advection time from pole to equator of ≈3 × 10⁷ s (≈30 rotation periods). Looking at Figure 4, one can determine an approximate time for toroidal field to propagate from the pole to the equator to be ≈1.5 × 10⁷ to 2 × 10⁷ s during a reversal. This slight discrepancy is likely due to the difference between the time-averaged latitudinal velocity compared to the typical latitudinal velocity of the convective cells, which at the base of the convection zone is ≈10⁴–10⁵ cm s⁻¹. This indicates that the bulk of the latitudinal transport is due to time-varying convection rather than the time-averaged meridional circulation.² In the animation³ of the images shown in Figures 3(a)–(d) one can see that the toroidal field is mostly advected by the convective eddies, therefore leading to a slightly more rapid equatorward propagation than would be implied by the mean latitudinal flow.

Over the duration of the simulation, the poloidal field does little else than diffuse and be advected by the overlying convection zone, as can be seen in Figures 3(e)–(h). As this model is axisymmetric, there can clearly be no regeneration of the poloidal field and hence, no dynamo. Therefore, one would expect the reversals of field sign to be due to fluctuations of the azimuthal velocity superimposed on the imposed differential rotation. In Figure 5, we show the azimuthal velocity as a function of time and latitude. One can see the general behavior of the (imposed) observed differential rotation with fast equator (prograde flow) and slow poles (retrograde flow), this is maintained throughout the simulation. However, one can also see weak signs of prograde flow at both poles, intermittent in the tachocline but persistent in the convection zone. Likewise one can see bands of prograde azimuthal flows starting at mid-latitudes with slight propagation toward the equator, although on much shorter timescales than the magnetic reversals, contrary to observations. This short timescale structure with equatorward propagation is also seen in the Figures 6–8, in the Tayler instability panel and, at some level, in the ratio of toroidal to poloidal field indicating that there is some interplay between magnetism and azimuthal flow variations, we are currently pursuing this interaction.

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

Because of the lack of obvious variation in the azimuthal flow on the same timescales of the magnetic field it appears that the reversals must be due to a local change in sign and amplification in toroidal field followed by advection. This can happen if the persistent minority sign seen at high latitudes is amplified and becomes dominant in the region.⁴

There has been significant discussion about the possible role of magnetic instabilities in stellar interiors (Wright 1973; Tayler 1975; Acheson 1978; Spruit 1999; Parfrey & Menou 2007). Although it is generally believed that convective or turbulent motions provide the required α effect to regenerate poloidal field from toroidal field, it was suggested by Spruit (2002), that convection is not needed because this effect could be provided by an instability in the field itself. He argued that the most relevant instability in stellar interiors was likely to be the Tayler instability. Parfrey & Menou (2007) argued that the diffusive magneto-rotational instability (MRI) was likely at high latitudes in the tachocline and that while this region overlaps with the region of Tayler instability the MRI is likely to be more dynamically important because of its exponential growth. They argued that this instability could disrupt magnetic fields, thus explaining the lack of sunspots at high latitude.

Here we find that the toroidal field reversals observed in Figure 4 are consistent with times when the toroidal field is unstable to the axisymmetric Tayler instability as described by the simple instability criterion provided in Goossens et al. (1981) and Spruit (1999), namely,

$$\begin{equation} \cos \theta {\frac{\partial {}}{\partial {\theta }}}\left(\ln \left(\frac{B^{2}}{\sin ^{2}\theta }\right)\right)>0. \end{equation} \tag{ 5 }$$

In Figure 6(b), we show the left-hand side of the inequality of (5). One can clearly see that reversals of the field are coincident with a Tayler instability that proceeds even at low latitudes. The instability does not appear in the bulk of the radiative interior where the toroidal field strength is substantially lower due to weaker differential rotation and the (stabilizing) poloidal field is correspondingly stronger. The instability is strongest in the tachocline and is present, but much less coherent in the convection zone. The bulk of the time the instability is confined to higher latitudes as expected (Parfrey & Menou 2007). In Figure 6(d), which shows the ratio of the local toroidal to poloidal field energy at high latitudes, one can see an almost oscillatory behavior of the instability. At these high latitudes the dominant toroidal field grows due to differential rotation, once the local ratio reaches a value around 10–20, the instability disrupts the field growth and magnetic flux destruction ensues. Normally, the field decay is weak and growth due to differential rotation is re-established shortly thereafter. However, occasionally, field destruction caused by the instability is more severe and the toroidal field energy drops precipitously. During these times the minority signed field can become dominant and be advected equatorward. Therefore, times of field reversals as seen in Figures 6–8 are associated with minima of the toroidal field energy relative to the poloidal energy and these minima are induced by the Tayler instability.

Indeed, one can see in Figure 6(d) that the slight north–south asymmetry in reversal coincides with a slight discrepancy in timing between flux destruction in the Northern and Southern hemispheres. This can be seen in all of the composite figures (Figures 6–8), when flux destruction occurs there is a reversal, when there is weak or no destruction, no reversal ensues.

The stability condition expressed in (5) is for a purely toroidal field. It has been shown (Wright 1973; Braithwaite & Spruit 2004) that mixed field configurations are more stable. Therefore, in addition to the condition expressed in (5) the toroidal field must also overcome the stabilizing effect of the poloidal field. This is what we see in the bottom panels of Figures 6–8. Substantial flux destruction only occurs when the criterion in (5) is met, along with the ratio of local toroidal to poloidal magnetic energy reaching some critical value. The critical value appears to be around 10–20, but clearly is not strict and other factors, such as local turbulent dissipation could become relevant. It should be noted that the non-axisymmetric Tayler instability and MRI are likely more relevant in real stars, this could lead to more regular flux destruction and possibly more regular reversals, although whether they would proceed like the ones here is unclear.

In the preceding section, we concentrated on Model B and found that it showed equatorward propagating toroidal field reversals. These reversals were precipitated by a Tayler instability and were sometimes coherent between Northern and Southern hemispheres, but not always. We ran four additional models with varying initial field strengths. All of the models, except Model A, show some evidence of reversal, however, only Model B shows reversals which are coherent in Northern and Southern hemispheres. Furthermore, only Model B and E show reversals of any substantial duration. Model A, which shows no reversal does show dominant toroidal field which weakens and then increases in time, but this is not accompanied by amplification and advection of minority field.

Figures 7 and 8 show the toroidal magnetic field, Tayler instability criterion, and ratio of toroidal to poloidal magnetic field energy for Model E and C, respectively. One can see the initial two signed field structures giving way to a single dominant field in each hemisphere, similar to Model B. Other equatorial propagating field reversals are also seen, particularly in the Southern hemisphere for Model E and the Northern hemisphere for Model C. One should note that the local ratio of the toroidal to poloidal field strengths in both of these models is similar to that seen in Model B, despite the vast difference in initial poloidal field strength; again implying the ratio of the field strengths plays some role in determining the growth of the Tayler instability. It is unclear why Model B shows symmetric reversals, while the others do not, it is also not clear why the various reversals are so different both in latitudinal excursion and in duration. These simulations indicate that weakening of the dominant toroidal field by the Tayler instability is robust, however, accompanying solar-like reversals may be possible only in limited parameter space of magnetic field strength. The dynamics that dictate that parameter space are presently unclear.