A study of global magnetic helicity in self-consistent spherical dynamos

ABSTRACT

Magnetic helicity is a fundamental constraint in both ideal and resistive magnetohydrodynamics. Measurements of magnetic helicity density on the Sun and other stars are used to interpret the internal behaviour of the dynamo generating the global magnetic field. In this note, we study the behaviour of the global relative magnetic helicity in three self-consistent spherical dynamo solutions of increasing complexity. Magnetic helicity describes the global linkage of the poloidal and toroidal magnetic fields (weighted by magnetic flux), and our results indicate that there are preferred states of this linkage. This leads us to propose that global magnetic reversals are, perhaps, a means of preserving this linkage, since, when only one of the poloidal or toroidal fields reverses, the preferred state of linkage is lost. It is shown that magnetic helicity indicates the onset of reversals and that this signature may be observed at the outer surface.

Keywords:

• Magnetohydrodynamics
• magnetic helicity
• self-consistent convective dynamo

1. Introduction

Magnetic helicity is an invariant of ideal magnetohydrodynamics (MHD) that describes field line topology weighted by magnetic flux (Moffatt 1969, Berger and Field 1984). This quantity is important for both laboratory (Taylor 1986) and astrophysical plasmas (Pevtsov et al. 2014, MacTaggart and Hillier 2020) since it remains approximately conserved in the weakly resistive limit (Berger 1984, Faraco and Lindberg 2019, Faraco et al. 2022). Magnetic helicity has been studied, in detail, in the context of mean field dynamo models (see Brandenburg and Subramanian 2005, Cameron et al. 2017, Brandenburg 2018, and references within). In self-consistent convection-driven dynamos in spherical shells, however, magnetic helicity has received little attention. On the scales typically considered in such global simulations, it is not possible to reach the high magnetic Reynolds number range needed to maintain an approximately constant value of the global magnetic helicity. This does not mean, however, that magnetic helicity is not a useful measure for understanding the nature of self-consistent spherical dynamos. In this work, rather than considering the role that magnetic helicity plays in turbulence (such as its influence on the α-effect in mean-field models), we will focus on what the global magnetic helicity can tell us about the global topology of the magnetic field, namely how the toroidal and poloidal fields link with each other and how variations in this linkage manifest themselves, particularly at the outer surface (which can be observed in stars).

Some initial inroads have already been made in this direction in observational studies of global magnetic helicity. Recent examples of work in this area, which focus on the Sun, include Pipin and Pevtsov (2014), Hawkes and Berger (2018), Pipin et al. (2019). In particular, Hawkes and Berger (2018) show that helicity flux during solar minima is a good predictor of the subsequent solar maxima, and can improve upon predictions based on the polar magnetic field alone. This result seems to be stronger when the helicity flux measured in only one hemisphere is considered, a calculation that can be performed due to the equatorial symmetry of helicity density at the solar surface. Magnetic helicity has also begun to be measured in stars other than the Sun (e.g. Lund et al. 2020, 2021). Again, only the surface density of magnetic helicity can be measured, and at limited resolution. Despite this, scaling laws relating the strength of the surface helicity density to that of the toroidal field have been found (Lund et al. 2021). Further, the results of the latter study suggest that the field topology is different for stars on different parts of their evolutionary track. Thus, it may be possible to use magnetic helicity density, in conjunction with other observables, to characterise different stages of stellar, and hence dynamo, evolution.

The purpose of this note is to investigate global magnetic helicity in various types of self-consistent convective dynamos in spherical shells. Our approach for this pilot investigation is not to try to model the Sun or a particular star or planet, but rather to analyse from a general standpoint some typical solutions to a well-studied model for spherical dynamos (Simitev and Busse 2005, Busse and Simitev 2011) that are representative of various known dynamo regimes. The three particular solutions that we consider increase in complexity and allow us to assess the usefulness of magnetic helicity as a prediction and analysis tool when confronted with increasingly chaotic spatial and temporal flow and field structures.

The paper is set up as follows. First, the main model equations are introduced, together with brief remarks on the numerical method used. This is followed by a description of magnetic helicity in spherical shells. We then consider each particular dynamo solution in turn, with a focus on what information is provided by magnetic helicity. The report ends with a summary of the results and a discussion of their interpretation.

2. Main equations

2.1. Mathematical model formulation

Following an established model of convection-driven spherical dynamo action (e.g. Busse and Simitev 2006, Simitev and Busse 2009, 2012) we consider a spherical shell that rotates about a fixed axis with a constant angular speed Ω. A static state is assumed to exist with the temperature distribution (1a) $T_S=T_0-\text{β}{d}^2{r}^2/2+\mathrm{\Delta }T\eta r^{-1}(1-\eta {)}^{-2},$(1a) (1b) $\text{β}=q/\left(3\kappa c_p\right),$(1b) (1c) $T_0=T_1-\mathrm{\Delta }T/(1-\eta ),T_S\left(r_i\right)=T_1,T_S\left(r_o\right)=T_2,\mathrm{\Delta }T=T_2-T_1,$(1c) where $T_1$ and $T_2$ are constant temperatures at the inner and outer spherical boundaries, $\eta =r_i/r_o$ is the ratio of the inner radius $r_i$ to the outer radius $r_o$, q is a uniform heat source density, κ is the thermal diffusivity, d is the shell thickness and r is the magnitude of the position vector to a point in the spherical shell.

Within the spherical shell, we solve the equations of magnetohydrodynamics (MHD) under the Boussinesq approximation (2a) $\begin{array}{rl}\mathbf{\nabla }\mathbf{\cdot }\mathit{u}& =0,\mathbf{\nabla }\mathbf{\cdot }\mathit{B}=0,\end{array}$(2a) (2b) $\begin{array}{rl}(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\mathit{u}\mathbf{\cdot }\mathbf{\nabla })\mathit{u}& =-\mathbf{\nabla }\pi -\tau \mathit{k}\times \mathit{u}+\vartheta \mathit{r}+{\mathrm{\nabla }}^2\mathit{u}+\mathit{B}\mathbf{\cdot }\mathbf{\nabla }\mathit{B},\end{array}$(2b) (2c) $\begin{array}{rl}P(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\mathit{u}\mathbf{\cdot }\mathbf{\nabla })\vartheta & =[R_i+R_e\eta r^{-3}(1-\eta {)}^{-2}]\mathit{r}\mathbf{\cdot }\mathit{u}+{\mathrm{\nabla }}^2\vartheta ,\end{array}$(2c) (2d) $\begin{array}{rl}{P}_m(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\mathit{u}\mathbf{\cdot }\mathbf{\nabla })\mathit{B}& =P_m\mathit{B}\mathbf{\cdot }\mathbf{\nabla }\mathit{u}+{\mathrm{\nabla }}^2\mathit{B}.\end{array}$(2d) In the above equations, $\mathit{u}$ is the velocity field, $\mathit{B}$ is the magnetic field, ϑ is the deviation from the static temperature distribution, π is an effective pressure including all terms that can be represented as a gradient, $\mathit{r}$ is the position vector and $\mathit{k}$ is the unit vector in the positive vertical direction about which the shell rotates.

Above, the Boussinesq approximation is used, where the density ρ is constant everywhere except in the buoyancy term. There the density varies linearly with respect to a constant background density ${\rho }_0$ and has the form (3) $\rho ={\rho }_0(1-\alpha \vartheta ),$(3) where α is a constant and represents the specific thermal expansion coefficient. The non-dimensional parameters in equations (2a(1b) $\text{β}=q/\left(3\kappa c_p\right),$(1b) ) are (4) $\begin{array}{r}{R}_i=\frac{\alpha \gamma \text{β}{d}^6}{\nu \kappa },R_e=\frac{\alpha \gamma \mathrm{\Delta }T{d}^4}{\nu \kappa },\tau =\frac{2\Omega d^2}{\nu },P=\frac{\nu }{\kappa },P_m=\frac{\nu }{\lambda },\end{array}$(4) which are the internal and external thermal Rayleigh numbers $R_i$ and $R_e$, the Coriolis number τ, the Prandtl number P and the magnetic Prandtl number $P_m$, respectively. Of the quantities not yet defined, λ is the magnetic diffusivity, ν is the viscosity and γ is a constant related to the gravitational acceleration $\mathit{g}=-d\gamma \mathit{r}$.

To solve the above equations, the velocity and magnetic field are first split into poloidal and toroidal parts, (5a) $\begin{array}{rl}\mathit{u}& ={\mathit{u}}_P+{\mathit{u}}_T=\mathbf{\nabla }\times (\mathbf{\nabla }\times \mathit{r}v)+\mathbf{\nabla }\times \mathit{r}w,\end{array}$(5a) (5b) $\begin{array}{rl}\mathit{B}& ={\mathit{B}}_P+{\mathit{B}}_T=\mathbf{\nabla }\times (\mathbf{\nabla }\times \mathit{r}h)+\mathbf{\nabla }\times \mathit{r}g,\end{array}$(5b) where $\mathit{r}=r{\mathit{e}}_r$. The scalar quantities are decomposed into spherical harmonics, (6a) $\begin{array}{rl}v& =\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l{V}_l^m(r,t){\mathrm{P}}_l^m(\cos \theta ){\mathrm{e}}^{\mathrm{i}m\varphi },w=\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l{W}_l^m(r,t){\mathrm{P}}_l^m(\cos \theta ){\mathrm{e}}^{\mathrm{i}m\varphi },\end{array}$(6a) (6b) $\begin{array}{rl}g& =\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l{G}_l^m(r,t){\mathrm{P}}_l^m(\cos \theta ){\mathrm{e}}^{\mathrm{i}m\varphi },h=\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l{H}_l^m(r,t){\mathrm{P}}_l^m(\cos \theta ){\mathrm{e}}^{\mathrm{i}m\varphi },\end{array}$(6b) where ${\mathrm{P}}_l^m\left(\mathbf{\cdot }\right)$ denotes the associated Legendre polynomials of first kind and $(r,\theta ,\varphi )$ are spherical coordinates. Once all the scalars are decomposed as above, the radial functions are expanded in terms of Chebychev polynomials and the series are truncated to a chosen resolution $(n_r,n_{\theta },n_{\varphi })$. The MHD equations are solved using a pseudospectral method and a combination of Crank-Nicolson and Adams-Bashforth time integration schemes. Further details of the numerical method can be found in Tilgner (1999) and an open source version of the code is available at Silva and Simitev (2018). For brevity, we do not repeat this description and instead guide the reader to these other works. The calculations performed in this work have been run with the resolutions $(n_r,n_{\theta },n_{\varphi })=(33,64,129)$ and $(n_r,n_{\theta },n_{\varphi })=(41,96,193)$. Azimuthally averaged components of the fields v, w, h and g will be indicated by an overbar.

To complete the specification of the above mathematical model, we require boundary conditions. We assume fixed temperatures at $r=r_i$ and $r=r_o$ and consider both cases with no-slip conditions (7a) $v=\frac{\mathrm{\partial }v}{\mathrm{\partial }r}=w=0,$(7a) and cases with stress-free conditions, (7b) $v=\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{r}^2}=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{w}{r}\right)=0.$(7b) For the magnetic field, we assume electrically insulating conditions at $r=r_i$ and $r=r_o$. At these locations, the poloidal function h matches the function $h^{\left(e\right)}$ which describes the potential fields outside the spherical shell (7c) $g=h-h^{\left(e\right)}=\frac{\mathrm{\partial }}{\mathrm{\partial }r}(h-h^{\left(e\right)})=0.$(7c)

2.2. Magnetic helicity

The magnetic field is generated and sustained by thermal convective motions inside the spherical shell (e.g. Busse and Simitev 2005). The field emanates outside of the spherical fluid shell where, in the absence of sustaining fluid motion, it takes the form of a freely decaying potential field. Thus, the magnetic field is not closed and we must consider relative magnetic helicity (Berger and Field 1984) instead of classical helicity (Woltjer 1958, Moffatt 1969). Assuming h is continuous at the inner and outer boundaries, Berger (1985) showed that the relative helicity in the spherical shell volume V has the appealing form (8) $H=2{\int }_V\mathit{L}h\mathbf{\cdot }\mathit{L}g\mathrm{d}V=2{\int }_V\left({\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{l}}^{-1}{\mathit{B}}_P\right)\mathbf{\cdot }{\mathit{B}}_T\mathrm{d}V,$(8) where $\mathit{L}=-\mathit{r}\times \mathbf{\nabla }$. Note that the definition of $\mathit{L}$ is chosen to match the particular magnetic field decomposition used in equation (5b(4) $\begin{array}{r}{R}_i=\frac{\alpha \gamma \text{β}{d}^6}{\nu \kappa },R_e=\frac{\alpha \gamma \mathrm{\Delta }T{d}^4}{\nu \kappa },\tau =\frac{2\Omega d^2}{\nu },P=\frac{\nu }{\kappa },P_m=\frac{\nu }{\lambda },\end{array}$(4) ). The representation in (8(7b) $v=\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{r}^2}=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{w}{r}\right)=0.$(7b) ) can be interpreted in terms of the global mutual linkage of the poloidal and toroidal magnetic fields (see also Berger and Hornig 2018). This is a generalised form of linkage as the poloidal and toroidal fields occupy the same region of space. However, as will be made clear in subsequent analysis, the linkage interpretation of magnetic helicity will prove very useful.

By simple vector algebra, we have that ${\mathit{B}}_T=\mathit{L}g$. There is no simple relation between ${\mathit{B}}_P$ and $\mathit{L}h$, but, in spherical coordinates, the latter can be expressed as (9) $\mathit{L}h=(0,\frac{1}{\sin \theta }\frac{\mathrm{\partial }h}{\mathrm{\partial }\varphi },-\frac{\mathrm{\partial }h}{\mathrm{\partial }\theta }).$(9) A test to confirm that the equation for the helicity H has been coded correctly is described in the Appendix.

3. Dynamo solutions

We now investigate the behaviour of global magnetic helicity in various dynamo solutions. Three cases of increasing complexity are considered as summarised in table 1.

Table 1. The three dynamo solutions considered in the study.

	Case 1: Steady dynamo	Case 2: Quasi-periodically reversing dynamo	Case 3: Aperiodically-reversing dynamo
η	0.35	0.4	0.4
$R_i$	0	3.5 $\times {10}^6$	$8.5\times {10}^5$
$R_e$	${10}^5$	0	0
τ	2000	$3\times {10}^4$	$3\times {10}^4$
P	1	0.75	0.1
$P_m$	5	0.65	1
Boundary conditions	no-slip, fixed temperatures, insulating outer space	stress-free, fixed temperatures, insulating outer space	stress-free, fixed temperatures, insulating outer space

3.1. Case 1: Steady dynamo

To begin with, we will consider a laminar dynamo solution that has been used as a benchmark case by the community for validating and testing the accuracy of numerical codes (Christensen et al. 2001, Matsui et al. 2016). The non-dimensional parameters and boundary conditions are listed in the first column of table 1.

The main feature of this dynamo is that it develops a steady-state for both the magnetic and velocity fields. The steady velocity profile exhibits columnar convection, as shown in figure 1. Focusing now on the magnetic field, there is a global dipole with inversions at the equator. This situation, together with other details, is displayed in figure 2.

Figure 1. Plots of the velocity field for Case 1 at a particular time (this is a steady flow so other times just represent rigid rotations of this solution) (a) shows lines of constant ${\overline{u}}_{\varphi }$ (left half) and $r\sin \theta {\mathrm{\partial }}_{\theta }\overline{v}$ (right half) in a meridional plane, visualising the toroidal and poloidal fields respectively; (b) shows lines of constant $u_r$ at $r=r_i+0.5$. There is a distinct pattern of columnar convection at both the equator and near the poles; and (c) shows lines of constant $r{\mathrm{\partial }}_{\varphi }v$ in the equatorial plane. Blue is negative, red is positive and green is zero in this figure and in all subsequent contour plots. (Colour online).

Figure 2. Plots of the magnetic field (top row) and magnetic helicity (bottom row) for Case 1 at a particular time. The top panel (a) shows meridional lines of constant ${\overline{B}}_{\varphi }$ (left half) and $r\sin \theta {\mathrm{\partial }}_{\theta }\overline{h}$ (right half) visualising the toroidal and poloidal fields respectively; (b) shows lines of constant $B_r$ at $r=r_i+0.96$; and (c) displays equatorial streamlines, $r{\mathrm{\partial }}_{\varphi }h$ = const. The bottom panel (d) shows contour lines of azimuthally averaged helicity density (left half) and the helicity density in a particular slice $\varphi =\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$ (right half); (e) shows the helicity density at $r=r_i+0.96$; and (f) shows the helicity density at the equatorial plane.

In figure 2, the top row reveals the behaviour of the magnetic field in (a) the meridional plane, (b) a near-outer spherical surface, and (c) the equatorial plane. The left-most plot shows the toroidal and poloidal components of the magnetic field. Both quantities shown in the hemispheres are meridional averages. In the left-hand hemisphere, displaying the large-scale behaviour of the toroidal field, the magnetic field pattern is perfectly antisymmetric about the equator. In the right-hand hemisphere, revealing the large-scale behaviour of the poloidal field, a dominant dipolar field is present, except near the surface at the equator. Here, smaller and weaker (compared to the polar fields) positive and negative polarities exist, as can be seen from the surface radial magnetic field plot in (b). The radial component of the magnetic field shows a dipolar symmetry, which means $B_r=0$ on the equatorial plane. This magnetic profile, which rotates rigidly with the sphere, is simple enough to provide a complete interpretation of the behaviour of magnetic helicity.

The bottom row of figure 2 displays the contour lines of (d) azimuthally averaged helicity density (left half) and unaveraged helicity density in a particular slice (right half), (e) the helicity density at the outer surface, and (f) the helicity density in the equatorial plane. Figure 2(d) reveals how the sign of the helicity density in different parts of the spherical shell depends on the linkage of the toroidal and poloidal fields. In order to clarify the signs of the helicity density displayed in this figure, we sketch an idealised linkage of the toroidal and poloidal fields in figure 3.

Figure 3. A representation of the linkage of toroidal and poloidal fields. Colours correspond to those in figure 2(a). The (Gauss) linkage in (a) is $-1$ and in (b) is $+1$. These signs correspond with those of the helicity density in figure 2(d). (Colour online).

The representation of linkage in figure 3 is idealised because the toroidal and poloidal fields are space-filling in the spherical shell and are not the union of two disjoint regions. However, due to the relative simplicity of the field structure (clearly separated regions of positive and negative toroidal field), this idealisation serves to indicate the correct sign of linkage. Figure 3(a) shows the linkage of a negative (blue) part of the toroidal field with the poloidal field and figure 3(b) shows the linkage of a positive (red) part of the toroidal field with the poloidal field. The Gauss linking numbers for these links are $-1$ and $+1$ respectively, corresponding to the negative and positive patches of helicity density in figure 2(d). Just to reiterate, this is only a simple representation of the linkage, but it is effective in showing how the magnetic helicity reveals information about the connectivity of poloidal and toroidal fields on large scales.

Figure 2(d) displays the distribution of helicity density near the surface. This shows a more detailed morphology compared to the azimuthal average plots, but there is still a clear dominance of positive helicity density in the north and negative helicity density in the south. As in the other plots, there is perfect antisymmetry about the equator. Thanks to this property, together with this being a steady solution, H = 0.

This simple dynamo case is a clear example showing that although the total magnetic helicity is zero, this value is due to the cancellation of positive and negative large-scale structures. In particular, the signs of these structures depend on the linkage of poloidal and toroidal fields. Magnetic helicity, however, is not only a measure of linkage but is weighted by magnetic flux. This consideration will become important for understanding cases exhibiting more irregular temporal behaviour and spatial morphology as discussed below.

3.2. Case 2: Quasi-periodically reversing dynamo

We now consider a time-dependent quasi-periodic dynamo solution that regularly reverses the signs of both the poloidal and toroidal parts of the magnetic field. The non-dimensional parameters for this case, together with the boundary conditions, are listed in the second column of table 1.

Figure 4 displays the typical structure of the velocity field. Although there is columnar convection like Case 1, the velocity is now evolving chaotically in time and exhibits much smaller-scale structures.

Figure 4. Typical structures of the velocity field for Case 2. (a) shows lines of constant ${\overline{u}}_{\varphi }$ in the left half and streamlines $r\sin \theta {\mathrm{\partial }}_{\theta }\overline{v}$ = const in the right half, all in the meridional plane; (b) shows lines of constant $u_r$ at $r=r_i+0.5$; and (c) shows streamlines, $r{\mathrm{\partial }}_{\varphi }v$ = const. in the equatorial plane. These images correspond to the time t = 5.538 in the simulation. (Colour online).

Unlike in Case 1, where the field structure was fixed in time and interpreting the helicity density in terms of linkage was straightforward, here the situation is more involved. Further, the magnetic Prandtl number $P_m$ is an order of magnitude smaller for this case compared to Case 1, so the total magnetic helicity is not expected to be strongly conserved. This is indeed the case upon examination of figure 5(a).

Figure 5. Time series of magnetic helicity for Case 2. (a) shows the total helicity as a function of time. (b) shows the helicity density integrated in the northern hemisphere (solid) and the southern hemisphere (dashed).

Since magnetic helicity is not conserved in time for this dynamo, we cannot make any clear causal conclusions based on helicity conservation. However, this does not mean that magnetic helicity does not reveal useful information. To reveal the structure of magnetic helicity in this dynamo solution, we restrict its calculation to northern and southern hemispheres separately. The results are shown in figure 5(b). There is a clear wave solution in both hemispheres, with a complete cycle taking about 0.02 time units. This result indicates that although there can be local changes in the magnetic helicity density in a hemisphere, there is no complete change in the linkage of hemispheric toroidal and poloidal fields. This global preservation of linkage arises from the nearly co-temporal reversal of the global toroidal and poloidal fields, leaving little time for a different large-scale linkage to develop. The reversals of the poloidal and toroidal fields are indicated in figure 6(a). Since the global toroidal and poloidal fields reverse together, the linkage, and thus the overall sign, of magnetic helicity in each hemisphere remains the same.

Figure 6. Time series related to reversals. (a) displays reversals of the poloidal and toroidal fields. The poloidal and toroidal fields are represented by their dominant spectral expansion components $H_1^0$ (solid line) and $G_1^0$ (dashed line), respectively. Both of these quantities are calculated at $r=r_i+0.5$. (b) displays the total magnetic energy $E_B$, which oscillates at a frequency similar to the hemispheric helicities displayed in figure 5.

Although the integrated magnetic helicity density does not change sign in each hemisphere, it does oscillate. Magnetic helicity can change due to variations in field linkage and magnetic field strength, the latter indicated by the magnetic energy shown in figure 6(b), and both these effects cause the oscillation shown in figure 5(b). The global magnetic helicity, in each hemisphere or in the entire spherical shell, is not an observable quantity. The magnetic helicity density at the surface is observable, however, and, for this dynamo solution, provides an indication of when a reversal develops. Figure 7 displays a time-latitude plot of the magnetic helicity density near the surface.

Figure 7. A latitude vs time plot of the magnetic helicity density at the surface for Case 2 with a nonlinear colourmap: indigo (−400, −30); maroon (−30, 30) and seismic (blue and red) (30, 400).

At t = 5.54, there is the clear positive/negative split in the north/south hemispheres. This time corresponds to the peaks of the magnetic helicity magnitudes in the hemispheres (see figure 5(b)). As the hemisphere magnitudes decrease to their minima, this is represented in figure 7 by, first, a more mixed pattern of linkage (i.e. increased negative magnetic helicity density in the north and vice versa in the south), and then by a decrease in the strength of the magnetic helicity density. These two behaviours can be seen just before and after t = 5.55 in figure 7, and repeat for all the cycles shown. The weakening of the density corresponds to the time when the poloidal and toroidal fields reverse and this is followed by the return of these fields to their peak values, resulting in a new phase of dominating positive/negative magnetic helicity density in the north/south hemispheres.

4. Case 3: Aperiodically-reversing dynamo

We now consider a dynamo solution in which aperiodic reversals of the global magnetic field occur. The general behaviour of this dynamo is rather more chaotic compared to the previous two cases considered, and there is no clear precursor pattern for global reversals. The non-dimensional parameters and boundary conditions used for this case are listed in the third column of table 1. For completeness, we show, in figure 8, typical profiles of velocity components for this dynamo solution. Like Case 2, the velocity has a columnar but chaotic morphology.

Figure 8. Typical structures of the velocity field in the case. (a) shows lines of constant ${\overline{u}}_{\varphi }$ in the left half and streamlines $r\sin \theta {\mathrm{\partial }}_{\theta }\overline{v}$ = const. in the right half, all in the meridional plane; (b) shows lines of constant $u_r$ at $r=r_i+0.5$ and (c) shows streamlines, $r{\mathrm{\partial }}_{\varphi }v$ = const. in the equatorial plane. These images correspond to the time t = 13.52 in the simulation.

This particular dynamo solution has been studied in Busse and Simitev (2008). Therefore, rather than repeating the description given in that work, we now focus on the behaviour of magnetic helicity at a reversal. Unlike Case 2, the reversals for this dynamo solution do not occur at regular intervals. That being said, despite the differences between these two cases, the magnetic helicity can again be used to interpret the behaviour of the global magnetic field during a reversal. Figure 9(a) displays dominant components of the poloidal and toroidal scalars $H_1^0$ and $G_1^0$ at a reversal and figure 9(b) displays the integrals of the magnetic helicity density in both the north and south hemispheres.

Figure 9. Time series during a reversal. $H_1^0$ (solid) and $G_1^0$ (dashed) are displayed in (a). Both quantities are evaluated at $r=r_i+0.5$. The blue vertical lines indicate the reversal times of these quantities. The northern (solid) and southern (dashed) helicities are shown in (b). Reversal times are indicated by green vertical lines. Just before the final green line, there are very small reversals about the H = 0 axis, and these are not shown. (Colour online).

Apart from the lack of periodicity, the other striking difference compared to Case 2 is that the hemispheric helicities change sign. The blue lines in Figures 9(a) and (b) indicate the reversals of $H_1^0$ and $G_1^0$ (both evaluated at $r=r_i+0.5$). The first two green lines in figure 9(b) indicate the reversals of $H\left(V_N\right)$ and $H\left(V_S\right)$. The helicity reversals occur long before those of $H_1^0$ and $G_1^0$. The third green line in figure 9(b) marks the time at which the signs of both the hemispheric helicities return to their original values (from the start of the time series). This occurs after the reversal of $H_1^0$ and before that of $G_1^0$. Some care is required in interpreting the reversal times as the helicities are dependent on the magnetic field throughout an entire hemisphere, whereas the values $H_1^0$ and $G_1^0$ presented here are evaluated at a single radius. Despite this, however, the results indicate that an oscillation in the hemispheric helicities, resulting in sign changes, precedes the global reversal. In other words, the linkage of the poloidal and toroidal fields changes, flipping and returning to its original state, leading to a global reversal. This suggests that the topology of the magnetic field changes and reaches an unstable state, in the sense of this state not lasting a long time. The way the magnetic field returns to a stable (longer lasting) configuration, and thus a stable topological state, is through a reversal. This pattern of helicity reversing before a global reversal is also true for the other reversals (labelled 2 and 3) that occur in the solution (for the time span we have simulated it). These results are displayed in table 2.

Table 2. The helicity and magnetic field reversal times for the other observed reversals.

Reversals	$H_1^0$	$G_1^0$	$H\left(V_N\right)$	$H\left(V_S\right)$	$H\left(V_N\right)$ and $H\left(V_S\right)$
2	26.3232	26.374	26.199	26.3313	26.3504
3	31.2784	31.3534	31.169	31.2532	31.3047

The quantities discussed so far in this section are available in simulations but not to observers, who only have access to data at the outer surface. Figure 10 displays the magnetic helicity density near the outer surface in a latitude-time plot.

Figure 10. A latitude vs time plot of the magnetic helicity density at the surface for Case 3 with a nonlinear colourmap: indigo (−500, −30); maroon (−30, 30; shown) and seismic (blue and red) (30, 500).

In figure 10 there is also an indication in the behaviour of the magnetic helicity density, at the surface, that reversals are taking place. The reversal of $H_1^0$ in figure 9, indicating the reversal of the poloidal field, occurs at $t\approx 13.45$. From about $t\approx 13.28$, there is a signature of stronger patches of positive magnetic helicity density near the equator. At around $t\approx 13.45$, there is an abrupt change with much weaker magnetic helicity density at all latitudes. This then switches to a phase of strong positive/negative magnetic helicity density, just south/north of the equator. This signature (a change in sign) is indicative of the reversal of the poloidal field but not the toroidal. When the toroidal field reverses at $t\approx 13.552$, the magnetic helicity density distribution returns to its original state. Therefore, even with just information about the magnetic helicity density at the surface, it is possible to identify the reversals of both the global poloidal and toroidal fields.

5. Discussion

In this note, we have analysed three dynamo solutions in rotating spherical shells and have investigated the behaviour of magnetic helicity in each. Our approach has been to adopt a well-studied model (Boussinesq MHD) for convection-driven dynamos in rotating spherical shells and consider dynamo solutions of increasing complexity – steady (Case 1), periodically-reversing (Case 2), and aperiodically-reversing (Case 3). Although magnetic helicity is not conserved in the final two cases, it does provide important information about the magnetic field in all cases. It has been shown, as indicated by its definition, that magnetic helicity can provide clear information about the linkage of toroidal and poloidal magnetic fields. For the cases with reversals, we have described how magnetic helicity relates to the reversal of both poloidal and toroidal magnetic fields and that changing states of the magnetic helicity density on the surface (an observable quantity in stars) can indicate the onset of reversals.

Since magnetic helicity is not strictly conserved in the simulations we have considered, we cannot claim that reversals are linked causally to magnetic helicity. However, what our results do suggest is that there are perhaps preferred states of global magnetic linkage for particular dynamo solutions. When the dynamo moves away from such states, it can return via a global reversal. For example, if there is a reversal of the large-scale toroidal field, this changes the global field linkage. To return to the original linkage, the poloidal field reverses, thus completing the global magnetic field reversal and returning the magnetic helicity distribution to what it was before any reversal took place. Therefore, even if magnetic helicity is not the cause of reversals in these simulations, it is a clear signature of the global poloidal-toroidal linkage, which is intimately linked to reversals. Furthermore, maps of the magnetic helicity density at the surface, which are measurable quantities in stellar observations, indicate the onset of reversals. This result is robust for dynamo solutions with different reversal mechanisms.

As mentioned previously, the main focus on magnetic helicity in dynamos has been the role it plays in the α-effect, in mean-field models. Thus, there is an assumed scale separation of magnetic helicity. Our simulations focus predominantly on the large-scale magnetic helicity, which alone (without the small-scale contribution) is not conserved. Nevertheless, our results do tie in qualitatively with some results from mean-field studies. For example, in a dynamo model based on the Babcock-Leighton α-effect, Choudhuri et al. (2004) found that at the start of cycles, helicities tend to be opposite of the preferred hemispheric trends. Temporary changes in the helicity hemisphere pattern have also been recorded in observations (Zhang et al. 2010).

The interpretation of our results in terms of “preferred” states of linkage may also provide a new interpretation of the origin of the hemisphere rule of magnetic helicity. Our results also suggest that imbalances in this rule can be caused by changes in the linkage of large-scale field (i.e. larger than the small-scale field considered in mean-field models). More work is needed to develop this area, but the interpretation of helicity in terms of toroidal/poloidal linkage may help to develop existing mean-field models which attempt to explain hemisphere imbalances (e.g. Yang et al. 2020).

In future work, in order to perform a closer comparison with solar models and observations, we will move beyond Boussinesq MHD and study the global magnetic helicity in anelastic models. This will enable us not only to mimic solar parameters more closely, but also the parameters of other specific stars.

Acknowledgements

Numerical computations were performed using the DiRAC Extreme Scaling service at the University of Edinburgh, operated by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This equipment was funded by BEIS capital funding via STFC capital grant ST/R00238X/1 and STFC DiRAC Operations grant ST/R001006/1. DiRAC is part of the National e-Infrastructure. PG acknowledges PhD sponsorship TRAW6PAA8 from the School of Mathematics and Statistics, the University of Glasgow.

Appendix

In order to verify the correct implementation of the helicity formula in the code, the following test was performed. In a magnetically closed spherical shell, linear force-free solutions exist where (A1) $\mathbf{\nabla }\times \mathit{B}=\lambda \mathit{B},$(A1) with constant λ. The values of λ that satisfy equation (A1(8) $H=2{\int }_V\mathit{L}h\mathbf{\cdot }\mathit{L}g\mathrm{d}V=2{\int }_V\left({\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{l}}^{-1}{\mathit{B}}_P\right)\mathbf{\cdot }{\mathit{B}}_T\mathrm{d}V,$(8) ) in the given domain form a discrete spectrum of eigenvalues for the curl operator. The corresponding eigenfunction of each eigenvalue has the form (in spherical coordinates) (A2a) $\begin{array}{rl}{B}_r& =\frac{1}{r^{3/2}}\left[c_1{\mathrm{J}}_{3/2}\right(\lambda r)+c_2{\mathrm{Y}}_{3/2}(\lambda r\left)\right]\cos \theta ,\end{array}$(A2a) (A2b) $\begin{array}{rl}{B}_{\theta }& =-\frac{1}{2{\lambda }^{3/2}{r}^3}\{\sin (\pi -\lambda r)[c_1(\sqrt{\frac{2}{\pi }}{\lambda }^2{r}^2-\sqrt{\frac{2}{\pi }})+\sqrt{\frac{2}{\pi }}{c}_2\lambda r]\\ & +\cos (\pi -\lambda r)[c_2(\sqrt{\frac{2}{\pi }}{\lambda }^2{r}^2-\sqrt{\frac{2}{\pi }})-\sqrt{\frac{2}{\pi }}{c}_1\lambda r]\}\sin \theta ,\end{array}$(A2b) (A2c) $\begin{array}{rl}{B}_{\varphi }& =\frac{\lambda }{2\sqrt{r}}\left[c_1{\mathrm{J}}_{3/2}\right(\lambda r)+c_2{\mathrm{Y}}_{3/2}(\lambda r\left)\right]\sin \theta ,\end{array}$(A2c) where ${\mathrm{J}}_{3/2}\left(x\right)$ and ${\mathrm{Y}}_{3/2}\left(x\right)$ are Bessel functions of the first and second kind, respectively, and $c_1$ and $c_2$ are constants determined from the boundary conditions (Cantarella et al. 2000).

For a closed magnetic field, we have $B_r=0$ on the inner and outer boundaries of the spherical shell. The two conditions can be written as the matrix-vector system (A3) $\left[\begin{array}{cc}{\displaystyle {\mathrm{J}}_{3/2}\left(\lambda r_i\right)}& {\mathrm{Y}}_{3/2}\left(\lambda r_i\right)\\ {\mathrm{J}}_{3/2}\left(\lambda r_o\right)& {\mathrm{Y}}_{3/2}\left(\lambda r_o\right)\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$(A3) For a non-trivial solution, we require that (A4) ${\mathrm{J}}_{3/2}\left(\lambda r_i\right){\mathrm{Y}}_{3/2}\left(\lambda r_o\right)-{\mathrm{Y}}_{3/2}\left(\lambda r_i\right){\mathrm{J}}_{3/2}\left(\lambda r_o\right)=0,$(A4) and we consider the smallest values of λ satisfying equation (A4(A2c) $\begin{array}{rl}{B}_{\varphi }& =\frac{\lambda }{2\sqrt{r}}\left[c_1{\mathrm{J}}_{3/2}\right(\lambda r)+c_2{\mathrm{Y}}_{3/2}(\lambda r\left)\right]\sin \theta ,\end{array}$(A2c) ). With λ found, $c_1$ and $c_2$ are readily determined. For $r_i=2/3$ and $r_o=5/3$, $\lambda =3.383384$ and we take $c_1=1$ and $c_2=1.974996$.

Since this particular magnetic field is force-free, it satisfies $g=\lambda h$. Using this property, the magnetic helicity can be written as (A5) $H=2\lambda {\int }_V|\mathit{L}h{|}^2\mathrm{d}V=\frac{2}{\lambda }{\int }_V|\mathit{L}g{|}^2\mathrm{d}V.$(A5) For this force-free field, the above relations represent an alternative way to calculate H compared to the general equation (8(7b) $v=\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{r}^2}=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{w}{r}\right)=0.$(7b) ), and provide a useful test that the general formula has been coded correctly. For the values used in this work, we find, for the above force-free solution, that H = 6.91 using either of equations (8(7b) $v=\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{r}^2}=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{w}{r}\right)=0.$(7b) ) and (A5(A3) $\left[\begin{array}{cc}{\displaystyle {\mathrm{J}}_{3/2}\left(\lambda r_i\right)}& {\mathrm{Y}}_{3/2}\left(\lambda r_i\right)\\ {\mathrm{J}}_{3/2}\left(\lambda r_o\right)& {\mathrm{Y}}_{3/2}\left(\lambda r_o\right)\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$(A3) ).

In figure A1(b), the meridional plot of the azimuthal average of the helicity density is equal to that of a specific slice (shown on the right-had side). This is because the magnetic field is symmetric and, also, not dependent on ϕ in this example. However, this result is connected to a more general one regarding the magnetic helicity of any magnetic field in a spherical shell. That is, (A6) ${\int }_0^{2\pi }\mathit{L}h\mathbf{\cdot }{\mathit{e}}_{\theta }\mathrm{d}\varphi ={\int }_0^{2\pi }\mathit{L}g\mathbf{\cdot }{\mathit{e}}_{\theta }\mathrm{d}\varphi =0,$(A6) for given values of r and θ. The proof of this result can be most easily seen by expanding the integrands in terms of their spectral decompositions, e.g. $\mathit{L}h\mathbf{\cdot }{\mathit{e}}_{\theta }=\sum _{l=0}^{\mathrm{\infty }}\sum _{m=-l}^l\frac{\mathrm{i}m}{\sin \theta }{H}_l^m(r,t){\mathrm{P}}_l^m(\cos \theta ){\mathrm{e}}^{\mathrm{i}m\varphi }.$ The azimuthal average is found by setting m = 0, hence confirming (A6(A4) ${\mathrm{J}}_{3/2}\left(\lambda r_i\right){\mathrm{Y}}_{3/2}\left(\lambda r_o\right)-{\mathrm{Y}}_{3/2}\left(\lambda r_i\right){\mathrm{J}}_{3/2}\left(\lambda r_o\right)=0,$(A4) ). Thus, plots of the azimuthal average of the helicity density depend only on $\left(\mathit{L}h\mathbf{\cdot }{\mathit{e}}_{\varphi }\right)\left(\mathit{L}g\mathbf{\cdot }{\mathit{e}}_{\varphi }\right)$. This constraint is another way to check that magnetic helicity has been calculated correctly in the code.

Figure A1. (a) shows meridional lines of constant ${\overline{B}}_{\varphi }$ (left half) and $r\sin \theta {\mathrm{\partial }}_{\theta }\overline{h}$ (right half); (b) shows averaged helicity density (left half) and unaveraged helicity density (right half); (c) shows the helicity density in the equatorial plane. (Colour online).