© ESO 2019

1. Introduction

Magnetic helicity is a measure of the topological complexity associated with the average linking of magnetic field lines (Moffatt 1969). In ideal magnetohydrodynamics, magnetic helicity is an ideal invariant in any magnetically closed region, and is known to be well conserved even in non-ideal conditions (Berger 1984). Indeed, the rate of change of magnetic helicity has been shown to give a lower bound to the rate of change of magnetic energy (Berger 1984).

The Sun’s corona is not magnetically closed, so the helicity in the corona can change due to boundary motions, even in a purely ideal approximation. If B = ∇ × A, then the magnetic helicity in the coronal volume V is given by

$$\begin{array}{c}\hfill H={{\displaystyle \int }}_V\mathit{A}·\mathit{B}{\mathrm{d}}^{3}x,\end{array}$$(1)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}H}{\mathrm{d}t}=-2{{\displaystyle \int }}_V\mathit{E}·\mathit{B}{\mathrm{d}}^{3}x-2{{\displaystyle \oint }}_S\mathit{A}\times \mathit{E}·\widehat{\mathit{r}}{\mathrm{d}}^{2}x,\end{array}$$(2)

where E is the electric field and S denotes the solar surface r = R_⊙. In the highly conducting corona, the volume term is usually negligible, but the surface term is not.

A difficulty with the magnetically open corona is that both H and dH/dt depend on the choice of the (non-unique) vector potential A. In Solar Physics, the standard resolution of this difficulty is to consider the relative helicity (Berger & Field 1984), with a potential reference field. For a spherical boundary S, this corresponds to the same formula H but with a specific gauge condition on A (DeVore 2000). Specifically, $\mathit{A}\times \widehat{\mathit{r}}={\mathit{A}}_{\mathrm{p}}\times \widehat{\mathit{r}}$ on S, where A_p satisfies ∇_h ⋅ A_p = 0 and $\widehat{\mathit{r}}·{\mathit{A}}_{\mathrm{p}}=0$ on S. With this definition, the helicity H for a given coronal B is uniquely defined, and we shall use H to denote this relative helicity throughout this paper.

Although H cannot be measured directly without knowing the three-dimensional magnetic field in the corona, the flux of helicity through the lower boundary S is determined entirely by an integral over this surface,

$$\begin{array}{c}\hfill -2{{\displaystyle \oint }}_S{\mathit{A}}_{\mathrm{p}}\times \mathit{E}·\widehat{\mathit{r}}{\mathrm{d}}^{2}x.\end{array}$$(3)

This integral depends purely on B_r and E. It follows that the flux of (relative) helicity out of the Sun is well-defined and independent of the coronal magnetic field B above the surface, and indeed of whether one defines an outer boundary of the corona or not. It is this flux that we estimate in this paper. Whilst H is conserved only for the global corona and not for each hemisphere separately, it is nevertheless more informative to consider the two hemispheres separately, since the helicity injection is roughly equal-and-opposite in each hemisphere (as we show in Sect. 3). In principle, coronal magnetic field lines may connect across the equator and thus cancel helicity between the hemispheres in the overall H integral. Estimation of such cross-equatorial cancellation requires extrapolation of the coronal magnetic field, which is beyond the scope of this paper. However, we expect that only a relatively small proportion of helicity is lost in this way in typical extrapolation models (cf. Yeates & Hornig 2016).

In this paper, we consider the global-scale injection of helicity due to surface flux transport, from large-scale solar rotation and meridional flow as well as small-scale supergranular diffusion. Previously, Berger & Ruzmaikin (2000) estimated that 4 × 10⁴⁶ Mx² of helicity were injected into each coronal hemisphere by solar rotation over a 22 year period from 1976 to 1998 (solar cycles 21 and 22), using low resolution magnetogram data from Wilcox Solar Observatory. This calculation was extended to 2018 by Hawkes & Berger (2018), who found, in particular, that the injected helicity in solar cycle 23 was lower than in the previous two cycles. Here we use surface flux transport modelling to extend this calculation to multiple solar cycles, considering magnetic flux distributions with higher spatial resolution and including also the helicity flux from meridional flow and supergranular diffusion.

DeVore (2000) took a different approach to estimating the helicity flux from differential rotation, based on an MHD calculation of the shearing of a single bipolar active region. By extrapolating to the number of active regions in a full solar cycle, the author gave a final estimate of 10⁴⁶ Mx² of helicity injected in this way during solar cycle 21. However, this approach neglects the contribution from weaker high-latitude magnetic fields. In this paper, the use of a surface flux transport model allows us to consistently include both the injection measured by DeVore (2000) as well as the contribution from large-scale fields that was captured by the approach of Berger & Ruzmaikin (2000).

In the approaches discussed so far, as well as in this paper, the electric field E on the solar surface is approximated by imposing an analytical flow and/or diffusion (see Sect. 2). By contrast, LaBonte et al. (2007) and Georgoulis et al. (2009) have used correlation tracking in high-cadence magnetogram sequences to estimate the electric field pattern in individual observed active regions, and hence the helicity flux. Using magnetograms from SoHO/MDI, the authors have estimated the helicity flux in the 393 largest active regions during solar cycle 23. Accounting for unobserved regions and uncertainties in the correlation tracking, they estimate a total helicity flux of 6.6 × 10⁴⁵ Mx² over solar cycle 23. Only about 20% of this flux came from differential rotation. However, owing to the need for high-resolution data, it has not been possible to apply this technique outside of active regions. This result is in agreement with that of Pevtsov (2008), who found a mean helicity of 1.7 × 10⁴³ Mx² per active region (based on a sample size of around 160), by fitting linear force free models to a series of magnetograms. Tziotziou et al. (2012) found a similar range of helicities in their analysis. Much of this helicity in young active regions arises from emergence and small-scale motions within the region, so would be present in addition to the large-scale generation studied in this article.

One can also try to put less direct constraints on solar helicity flux. From the viewpoint of the solar interior, Brandenburg & Sandin (2004) and Brandenburg (2009) suggest that the solar dynamo must shed 10⁴⁶ Mx² of magnetic helicity per cycle if it is to avoid catastrophic alpha quenching caused by a build up of small-scale magnetic helicity. From the viewpoint of the heliosphere, Bieber & Rust (1995) estimated a helicity ejection rate of 2 × 10⁴⁵ Mx² through coronal mass ejections by considering toroidal magnetic flux, while DeVore (2000) gave a higher estimate of 10⁴⁶ Mx² by modelling the magnetic structure of interplanetary magnetic clouds. More recently, Démoulin et al. (2016) extrapolated data from 107 observed magnetic clouds to estimate a total ejection rate of 2.5 × 10⁴⁶ Mx² over solar cycle 23, and a similar value was obtained independently by Lowder & Yeates (2017) through non-potential modelling of flux rope formation and ejection in the low corona.

In summary, all of these estimates are broadly consistent, but they are all approximations. It is a question of ongoing interest to study how much helicity is contributed by different physical processes, and particularly how this varies between solar cycles. In Sect. 2 we describe the two surface flux transport models used in this paper to estimate the helicity flux, before presenting the results in Sect. 3 and our conclusions in Sect. 5.

2. Methods

The helicity flux through the photosphere cannot be observed directly, so instead we need to estimate ${\mathit{A}}_{\mathrm{p}}\times \mathit{E}·\widehat{\mathit{n}}$ on this spherical surface indirectly. In this paper, we do so using the surface flux transport model. Two different implementations are used, as will be outlined in this section. The surface flux transport model (Jiang et al. 2014) evolves the photospheric radial magnetic field, B_r(R_⊙, θ, ϕ, t) through

$$\begin{array}{c}\hfill \frac{\partial B_{\mathrm{r}}}{\partial t}=-\widehat{\mathit{r}}·\mathrm{\nabla }\times \mathit{E}+S(\theta ,\varphi ,t)-D(\theta ,\varphi ,t)\end{array}$$(4)

where S(θ, ϕ) is a source term representing the emergence of new magnetic flux, and D(θ, ϕ, t) is an additional decay term that accounts for the effect of radial (inward) diffusion that cannot be treated self-consistently in this two-dimensional model (Baumann et al. 2006). The electric field takes the form

$$\begin{array}{cc}\hfill E_{\theta }(\theta ,\varphi ,t)& =-v_{\varphi }{B}_{\mathrm{r}}+\frac{\eta }{R_{\odot }sin\theta }\frac{\partial B_{\mathrm{r}}}{\partial \varphi }\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill E_{\varphi }(\theta ,\varphi ,t)& =v_{\theta }{B}_{\mathrm{r}}-\frac{\eta }{R_{\odot }}\frac{\partial B_{\mathrm{r}}}{\partial \theta }·\hfill \end{array}$$(6)

This represents the evolution of B_r with (a) advection by large-scale horizontal flows v(θ) and (b) horizontal (turbulent) diffusion due to supergranular convection, with coefficient η.

To estimate the helicity flux, we take the electric field (5), (6) as an approximation of the horizontal electric field on the solar photosphere. We thus neglect any electric field arising from flux emergence, as well as from localised flows within individual active regions. The systematic observations required to account accurately for these effects are simply not available for the time periods considered here. We also neglect the contribution from the decay term (we estimate at the end of this subsection that its contribution would be an order of magnitude smaller than that from supergranular diffusion).

In addition, we need to calculate A_p from B_r at each time. This is done by expressing ${\mathit{A}}_{\mathrm{p}}=-\mathrm{\nabla }\psi \times \widehat{\mathit{r}}$, whence ψ satisfies a two-dimensional Poisson equation

$$\begin{array}{c}\hfill B_{\mathrm{r}}=-\widehat{\mathit{r}}·\mathrm{\nabla }\times \mathrm{\nabla }\times (\psi \widehat{\mathit{r}})=-{\mathrm{\nabla }}_{\mathrm{h}}^2\psi .\end{array}$$(7)

This equation is solved at each time by a five-point finite-difference method that reduces to a Fast Fourier Transform in ϕ and a series of tridiagonal eigenvalue problems in θ (cf. van Ballegooijen et al. 2000). The two flux transport models differ both in the choice of imposed flow parameters v and η and in the form of the source term S, as we now describe.

2.1. Model driven by observed active regions (W18)

For solar cycles 21–23, we use the flux transport model of Yeates et al. (2015), in which the source term S(θ, ϕ, t) comprises individual active regions assimilated from synoptic magnetogram observations. These regions are selected from US National Solar Observatory radial-component magnetograms based on a flux threshold ( ≥ 39.8 G). Regions with sufficiently balanced flux are inserted instantaneously into the simulation by replacing the corresponding pixels of B_r, correcting to ensure that flux balance is maintained. The rest are discarded. This technique allows us to represent individual emerging regions, and without approximating them as idealised magnetic bipoles. The disadvantage is that systematic magnetogram observations are required, which are not available for earlier cycles.

The particular simulation used here, henceforth denoted W18, was published by Whitbread et al. (2018), who determined the parameters through systematic optimization against the observed magnetic butterfly diagram for cycles 21–23, using a genetic algorithm. The initial condition was set to

$$\begin{array}{c}\hfill B_{\mathrm{r}}(\theta ,\varphi )=B_0{|cos\theta |}^{7}cos\theta .\end{array}$$(8)

The azimuthal velocity was v_ϕ(θ) = R_⊙Ω(θ)sinθ, with angular velocity

$$\begin{array}{c}\hfill \mathrm{\Omega }(\theta )={\mathrm{\Omega }}_0+{\mathrm{\Omega }}_2{cos}^2\theta +{\mathrm{\Omega }}_4{cos}^4\theta ,\end{array}$$(9)

while the meridional velocity took the functional form

$$\begin{array}{c}\hfill v_{\theta }(\theta )=-\frac{v_0}{v_{\max }}{sin}^{\mathrm{p}}\theta cos\theta ,\end{array}$$(10)

where v_max = max_θ (sin^pθcosθ). The decay term in this model took the form of a simple exponential decay D(θ, ϕ, t) = B_r/τ.

The optimization by Whitbread et al. (2018) determined simultaneously the values of the parameters B₀, v₀, p, η and τ for the time period 1976 May 1–2008 October 23, selecting B₀ = 6.7 G, v₀ = 9.2 m s⁻¹, p = 2.33, η = 466.8 km² s and τ = 10.1 yr. The rotation profile Ω(θ) was fixed with Ω₂ = −2.396 deg day⁻¹ and Ω₄ = −1.787 deg day⁻¹ as determined by Snodgrass & Ulrich (1990). The original W18 simulation was performed in the Carrington frame, but for computing the helicity flux, we include not only the contribution from differential rotation in the Carrington frame but also that from uniform rotation of this frame with respect to the fixed stars. Thus we set Ω₀ = 14.713 deg day⁻¹ (the sidereal rotation rate at the equator). The v_θ and v_ϕ profiles are shown in Fig. 1.

Fig. 1. Velocity profiles as a function of latitude for meridional velocity vθ (left), and differential rotation vϕ (right), in W18 and J11. The black dashed curve for vϕ shows the profile used by Berger & Ruzmaikin (2000) and also by Hawkes & Berger (2018).

For a simple comparison of the magnitude of the various terms, note that the peak flow speed of rotation is |v_ϕ| ≈ 2 km s⁻¹, compared to the maximum meridional flow speed of |v_θ| ≈ 9.2 m s⁻¹. To calculate an effective diffusion speed, we must define a charecteristic length scale; given the signed nature of helicity, we choose this to be of the order of a unipolar structure, which corresponds to 10 degrees, or approximately 0.17 R_⊙. In this case, we calculate an effective diffusion speed of |v_η| ≈ η/0.17 R_⊙ = 3.97 m s⁻¹. We do not calculate the helicity flux solely attributable to the decay term D(θ, ϕ, t) – calculating its effective decay speed (at the same characteristic length scale) gives |v_D| ≈ 0.17 R_⊙/τ = 0.37 m s⁻¹, which is an order of magnitude below that of |v_η|. We see in Sect. 3 that this ordering determines the relative contribution of each of these terms to the helicity flux.

2.2. Model driven by statistical active regions (J11)

In order to estimate the helicity flux over multiple solar cycles, we employ a flux transport simulation by Jiang et al. (2011a), who simulated the solar magnetic field evolution continuously over the period 1700–2009. Henceforth denoted J11, this has several differences compared to the W18 model. Most significantly, the flux emergence term S(θ, ϕ, t) consists of bipolar active regions determined semi-synthetically from statistical distributions modulated by observed sunspot numbers (Jiang et al. 2011b). This allows the model to extend for many solar cycles prior to the magnetograph era. The flow profiles are also different: J11 uses a lower supergranular diffusivity η′ = 250 km² s⁻¹, compensated by a slightly faster meridional flow with a slightly sharper gradient at the equator, given by

$$\begin{array}{c}\hfill v{\prime }_{\theta }(\theta )=\{\begin{array}{cc}-v_0\prime sin(2.4\lambda )\hfill & \lambda <1.31\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array},\end{array}$$(11)

where λ  =  π/2 − θ and v′₀  =  11 m s⁻¹. The differential rotation profile was slightly different from W18, with Ω′₂  =   − 2.30 deg day⁻¹ and Ω′₄  =   − 1.62 deg day⁻¹. Again we set Ω₀′  =  14.713 deg day⁻¹. Finally, the decay term D(θ, ϕ, t) in J11 uses a more sophisticated form where different spherical harmonic components of the solution decay at different rates, commensurate with the diffusive decay of a uniform field in the solar convection zone. The strength of convection zone diffusivity used here gives an effective decay time of about 20 yr, a little longer than the W18 model.

To compute the helicity flux, we once again compute A_p from B_r′, using the same finite-difference method as for W18. The electric field, however, is determined using v_ϕ′, v_θ′ and η′.

Figure 2a shows that J11 produces a larger unsigned magnetic flux in each cycle than W18, commensurate with the lower diffusivity. On the other hand, Fig. 2b shows that the relative strength of their axial dipole moments (and hence polar fields) varies from cycle to cycle, commensurate with the differences both in transport and particularly in the source term.

Fig. 2. Time variation of unsigned magnetic flux (a), axial dipole moment (b), and hemispheric helicity fluxes (c)–(f) for the W18 and J11 models over solar cycles 21–23. The net helicity flux FN/S (c) is the sum of those due to solar rotation Fvϕ, N/S (d), meridional velocity Fvθ, N/S (e) and supergranular diffusion Fη, N/S (f) for both W18 (blue shades) and J11 (red shades).

3. Results for solar cycles 21–23

We first consider the period covered by both flux transport models: solar cycles 21–23.

3.1. Relative contributions of different terms

Figure 2c shows the net helicity flux out of the solar surface in each hemisphere as a function of time, while Figs. 2d–f show the separate contributions from differential (sidereal) rotation, meridional velocity and supergranular diffusion, for both W18 and J11. Throughout this paper, plots describing W18 are shaded blue, while those describing J11 are shaded red. The contributions from the Northern hemisphere N are defined as

$$\begin{array}{cc}\hfill F_{v_{\varphi },\mathrm{N}}(t)& =-2{{\displaystyle \int }}_N{A}_{\mathrm{p}\varphi }{v}_{\varphi }{B}_r{\mathrm{d}}^{2}x,\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill F_{v_{\theta },\mathrm{N}}(t)& =-2{{\displaystyle \int }}_N{A}_{\mathrm{p}\theta }{v}_{\theta }{B}_r{\mathrm{d}}^{2}x,\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill F_{\eta ,\mathrm{N}}(t)& =2\eta {{\displaystyle \int }}_N{\mathit{A}}_{\mathrm{p}}·\mathrm{\nabla }{B}_r{\mathrm{d}}^{2}x,\hfill \end{array}$$(14)

with their sum giving the net hemispheric flux

$$\begin{array}{c}\hfill F_{\mathrm{N}}=-2{{\displaystyle \int }}_N{\mathit{A}}_{\mathrm{p}}\times \mathit{E}·\widehat{\mathit{r}}{\mathrm{d}}^{2}x=F_{v_{\varphi },\mathrm{N}}+F_{v_{\theta },\mathrm{N}}+F_{\eta ,\mathrm{N}}.\end{array}$$(15)

The fluxes are defined similarly for the Southern hemisphere S. The net helicity fluxes in north and south need not be exactly equal and opposite, although they remain so to a good approximation. In fact, we found that it can be shown analytically that the helicity flux from supergranular diffusion is is exactly balanced between the hemispheres, and that from differential rotation is balanced to leading order (see Appendix A).

Figure 2 shows clearly that the contributions from solar rotation dominate the flux overall, followed by meridional velocity and supergranular diffusion in that order for both simulations. To quantify the difference in magnitude between the helicity fluxes, we compute the ratio

$$\begin{array}{c}\hfill f_{v_{\varphi }}(t)=\frac{1}{2}(\frac{|F_{v_{\varphi },\mathrm{N}}|}{|F_{v_{\varphi },\mathrm{N}}|+|F_{v_{\theta },\mathrm{N}}|+|F_{\eta ,\mathrm{N}}|}+\frac{|F_{v_{\varphi },\mathrm{S}}|}{|F_{v_{\varphi },\mathrm{S}}|+|F_{v_{\theta },\mathrm{S}}|+|F_{\eta ,\mathrm{S}}|}),\end{array}$$(16)

and similarly for the other two fluxes. Table 1 shows the mean, maximum and minimum values of these ratios, over the course of the W18 and J11 simulations. The mean magnitude of each term is roughly separable by an order of magnitude, consistent with the magnitudes of the shear velocities associated with each term (in the same order), as described in Sect. 2.1.

Notice that W18 and J11 agree in the ordering of contributions in Table 1, although there can be substantial differences between their actual fluxes at any given time, as seen in Fig. 2. Around Solar Maximum, differences between the helicity fluxes in W18 and J11 arise because J11 has more absolute magnetic flux at low latitudes, which tends to lead to greater helicity flux (except from diffusion, which is weaker in J11). Around Solar Minimum, the helicity fluxes from the two models are qualitatively more alike. Which is larger at Minimum varies from cycle to cycle, according to which has the larger axial dipole moment. This is because both the helicity flux and axial dipole at Minimum are dominated by the high latitude (polar) field.

3.2. Net injection per solar cycle

Table 2 show the time-integrated helicity flux during each solar cycle, for the two models. Following Whitbread et al. (2017), we define Cycle 21 as 1976 May 1–1986 March 10, Cycle 22 as 1986 March 11–1996 June 1 and Cycle 23 as 1996 June 2–2008 August 3. Concentrating on the net helicity flux (top two rows), we observe that the sign of this net flux in W18 is negative in the Northern hemisphere and positive in the Southern for all three cycles, while in J11 it has the opposite sign in Cycles 21 and 22. We show below that this difference arises because the solar-cycle helicity flux is a signed quantity that involves much cancellation between contributions from different latitudes. This makes it rather sensitive to model changes.

To illustrate this latitudinal cancellation, Fig. 3 breaks the helicity fluxes down into their contributions in time and latitude, i.e., the surface integrals (12)–(15) are performed over longitude but not latitude. Figures 3c and d show that both W18 and J11 have the same qualitative behaviour, with strong contributions of both signs from the “butterfly wings” below about 55° latitude and a weaker but systematic contribution of systematic sign from the high-latitude polar field. In the Northern hemisphere, this high-latitude helicity flux is always negative, consistent with the sense of the interplanetary Parker spiral. It arises simply because the product A_pϕB_r for a polar dipole is always positive/negative in the northern/southern hemisphere, irrespective of the polarity of B_r. (If the sign of B_r changes, so does the sign of A_p.) Evidently, the model dependence of the cycle-integrated flux in Table 2 must arise from the balance between these differently signed contributions.

Fig. 3. Latitude–time distributions of longitude-averaged radial magnetic field (a,b) alongside longitude integrated total helicity flux (c,d), and helicity flux associated individually with differential rotation (e,f), meridional velocity (g, h) and supergranular diffusion (i,j) for W18 and J11 respectively. We stress that these latitudinal distributions are not in themselves physically meaningful, since there is typically significant cancellation between different latitudes. The same colour scales are used for W18 and J11. Units are G for (a,b) and Mx2 day−1 for (c–j).

The positive and negative contributions from the butterfly wings arise not from Joy’s law (like in the usual magnetic butterfly diagram) but simply from the fact that active regions are typically bipolar in the East-West direction. To illustrate this, Fig. 4 shows B_r, A_pϕ, and the (dominant) helicity flux integrand for two active regions A and B, during W18. Region A is stronger, so its quadrupolar A_pϕ distribution dominates that from the global dipole, leading to a product A_pϕB_r whose polarity is arranged North–South. This North–South pattern is the same for any strong active region, irrespective of the sign of B_r or the hemisphere. Region B does not show a quadrupolar helicity injection because its A_pϕ is weaker than that from the global dipole. Since the helicity flux is effectively quadratic in B_r, it is dominated by the strong regions, thus leading to the latitudinal sign pattern seen in the butterfly wings.

Fig. 4. Illustration of the source of helicity flux Fvϕ at a single time during November 1977 (during the early rise phase of Cycle 21), in W18. Panels show (a) Br, (b) Aϕ and (c) the integrand in (12).

It turns out that the positive and negative helicity fluxes in each butterfly wing are not usually equal. In fact, the net contribution over a solar cycle from latitudes |λ|< 55° is always opposite in sign to that from the high latitudes, being positive/negative in the Northern/Southern Hemisphere. This imbalance arises because, later in the cycle, the equatorward polarity helicity fluxes tend to cancel across the equator, leaving an excess of the poleward polarities. Early in the cycle, there is little equatorial cancellation, and the net contribution tends to be opposite in sign because v_ϕ is slightly stronger nearer to the equator.

The net helicity flux in each hemisphere over a solar cycle is therefore a balance between oppositely signed contributions from high and low latitudes. We quantify this balance in Table 3 which shows means of the ratios

$$\begin{array}{cc}\hfill f_{\mathrm{N},|\lambda |<55}(t)& =|F_{\mathrm{N},|\lambda |<55}(t)|/(|F_{\mathrm{N},|\lambda |<55}(t)|+|F_{\mathrm{N},|\lambda |>55}(t)|),\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill f_{\mathrm{N},|\lambda |>55}(t)& =|F_{\mathrm{N},|\lambda |>55}(t)|/(|F_{\mathrm{N},|\lambda |>55}(t)||+|F_{\mathrm{N},|\lambda |>55}(t)|),\hfill \end{array}$$(18)

and similarly for the Southern hemisphere. Notice that the high latitudes dominate for W18 in all cycles, explaining the net negative/positive helicity flux in the North/South. On the other hand, the low latitudes dominate for J11 during Cycles 21 and 22, explaining the opposite signs. This tendency for the low latitude butterfly wings to contribute more in J11 simply results from the stronger magnetic flux in that simulation, because of the lower supergranular diffusivity. This higher unsigned flux is evident in Fig. 3b compared to Fig. 3a, as well as in Fig. 2a.

3.3. Comparison with earlier work

The third row of Table 2 shows the helicity fluxes estimated by Hawkes & Berger (2018) using WSO data, using the differential rotation profile of Berger & Ruzmaikin (2000). This profile is shown by the black dashed line in the right panel of Fig. 1, and is seen to be close to those used inthe present paper. The Hawkes & Berger (2018) helicity fluxes follow F_vϕ from the W18 model reasonably closely for Cycle 22, but are a factor two stronger for Cycles 21 and 23. These discrepancies can be attributed to the different spatial resolutions used for the two studies: in Hawkes & Berger (2018) the authors use only the first few degrees from the WSO spherical harmonic decomposition, and as such calculate the large scale winding without accounting for the contribution of active regions. Since active regions are important for the net helicity flux in our study, the Hawkes & Berger (2018) helicity flux is not necessarily directly comparable, but we include it for the sake of comparison.

4. Results for earlier solar cycles

Having demonstrated the basic pattern of helicity flux in surface flux-transport models, we now extend the calculation for J11 to the full time period of that simulation, years 1700–2009. Covering twenty-seven complete solar cycles, the results are shown in Fig. 5. We observe the same qualitative behaviour as described in Sect. 3 across the entire simulation period, with the terms being ordered in magnitude as before. Having multiple cycles, however, enables us to consider the variation in cycle-integrated helicity flux from one cycle to the next.

Fig. 5. Results for the full J11 model, showing unsigned magnetic flux (a), axial dipole moment (b), total helicity flux (c) and helicity flux from solar rotation (d), meridional velocity (e) and supergranular diffusion (f). In (c–f), light red shows the northern hemisphere and dark red the southern hemisphere.

Hawkes & Berger (2018) showed that the large scale helicity flux estimated from WSO data correlates more strongly with polar magnetic field strength than with sunspot number, for cycles 20–23. We find similar behaviour in our helicity fluxes for the J11 simulation when considering all 27 cycles. The polar field trend is shown (by proxy) in Fig. 6, where cycle-integrated helicity flux is plotted against the strength (absolute value) of the end-of-cycle axial dipole moment, for each cycle. Helicity fluxes above and below 55° latitude are shown separately, because they differ both in sign (as seen in Sect. 3) but also in the slope of the trend. The latter is expected because both high-latitude helicity flux and end-of-cycle axial dipole are determined by the strength of the high-latitude polar field. As is well known, the end-of-cycle polar field is not a direct function of the cycle amplitude (strength of emerging active region magnetic flux), since it depends also on the orientation and locations of active regions (Jiang et al. 2014). Thus the high-latitude contribution causes the helicity flux to correlate better with axial dipole moment than with either cycle amplitude or cycle-integrated magnetic flux.

Fig. 6. Scatter plots of integrated helicity flux per cycle in each sub-hemisphere for both W18 and J11, plotted against the axial dipole strength (absolute value) at the end of each cycle. This time is defined as the time of minimum unsigned magnetic flux in Figs.1a and 2a. Lines of best fit and Pearson’s correlation coefficient are calculated for the J11 data.

For comparison, Fig. 6 also shows the W18 data points. Most of these do not deviate greatly from the J11 trend lines, except for one low-latitude contribution in the Northern hemisphere (for Cycle 22), which has the opposite sign. We found this to be associated with the more rapid rise to maximum of this cycle compared with Cycles 21 and 23 in W18, visible in Fig. 2a. The reason for the reversed sign of helicity flux is that early active regions are further from the equator, so that the corresponding part of the butterfly wing does not make a net positive contribution (as mentioned in Sect. 3.2). Rather, the differential rotation gradient tends to win out early in the cycle and the net contribution is negative. As shown by Fig. 7, this sign-reversal phenomenon occurs in all three W18 cycles. But it only beats the later contribution in Cycle 22, and only for the Northern hemisphere, owing to the stronger early flux emergence.

Fig. 7. Time variation of the |λ|< 55° helicity flux in the Northern hemisphere for W18 (red), compared to the unsigned magnetic flux (blue). The black line denotes the zero axis for helicity flux.

5. Conclusion

To summarise, we have used surface flux transport models to estimate the amount of magnetic helicity injected into the solar corona through evolution of the large-scale magnetic field on the solar surface, on solar cycle timescales. We find a consistent pattern whereby negative/positive helicity is systematically injected at high latitudes in the northern/southern hemisphere. In the lower-latitude wings of the magnetic butterfly diagram, the net helicity injection over a solar cycle is usually opposite to this in sign, unless the flux emergence is dominated by early active regions far from the equator. The overall helicity injection rate is therefore a balance between these high and low latitude contributions, and thus quite sensitive to the details of the flux transport model.

Using the 27-cycle J11 simulation driven by statistical active region emergence (Jiang et al. 2011a), we have found that the rate of helicity injection in any given cycle correlates well with the end-of-cycle axial dipole moment. This accords with the previous results of Hawkes & Berger (2018) using lower resolution WSO magnetogram data, who found a similar relation with the polar field for solar cycles 20–23. Compared to that study – which extended the earlier work of Berger & Ruzmaikin (2000) – our models have higher spatial resolution for the magnetic field. We also include the helicity flux from meridional velocity and supergranular diffusion, although we have shown that this is essentially negligible in comparison with the flux from solar rotation.

Recently, Pipin et al. (2019) have published latitude-time maps of A_p ⋅ B for 2010–2019 using vector synoptic magnetograms from the HMI instrument on Solar Dynamics Observatory. Whilst these show only the local helicity density in the solar photosphere and not the flux of helicity out of the Sun, they are nonetheless interesting to compare with our results. In particular, we compare with their Fig. 5a showing the longitude-averaged helicity density. At low latitudes, the helicity density shows a similar tendency to our helicity flux, namely net negative earlier in the cycle and net positive later in the cycle (cf. our Fig. 7). However, at higher latitudes, the helicity density changes sign around 2014, following our sign pattern only in the first half of the cycle. Given the pioneering nature of these vector synoptic maps, and inherent uncertainties outside of active regions, it will be interesting to see whether this pattern is confirmed by other instruments.

Finally, we emphasize two limitations of our study that require further work. Firstly, we neglected the contribution of helicity due to emergence and small-scale motions in young active regions; this may be addressed in the future with modern vector magnetogram observations. Secondly, we have computed only the flux of relative magnetic helicity out of the Sun. Once in the corona, the fate of this helicity depends on whether it is ejected onto a magnetic field line whose other end is open into the heliosphere, or one whose other end closes back to the solar surface. Depending on the relative injection at their two footpoints, closed field lines can store helicity within the corona, releasing it only in the form of sporadic coronal mass ejections (Low 1994; Bieber & Rust 1995). The details of this process require time-dependent modelling of the magnetic field structure in the corona, which is beyond the scope of the present paper. See Yeates & Hornig (2016) for some preliminary calculations in this regard.

Acknowledgments

The authors would like to thank R. Cameron for providing the radial magnetic field data from the J11 simulations. We thank the UK STFC for their funding under grants ST/N504063/1 (GH) and ST/S000321/1 (ARY). This collaboration was facilitated through an ISSI International Team on Magnetic Helicity in Astrophysical Plasmas, and the NORDITA programme on Solar Helicities in Theory and Observations. GH would also like to thank Mitchell Berger for comments and supervision made during this work. We also thank the anonymous referee for useful suggestions which improved this work.

References

Appendix A: Hemispheric symmetry of magnetic helicity flux

For each of the proofs in this section, we make use of our gauge condition ∇ ⋅ A_p = 0, specifically

$$\begin{array}{c}\hfill \frac{\partial }{\partial \theta }(sin\theta A_{\mathrm{p}\theta })+\frac{\partial A_{\mathrm{p}\varphi }}{\partial \varphi }=0.\end{array}$$(A.1)

Supergranular diffusion

The net helicity flux over the whole solar surface from supergranular diffusion (following Eq. (14)) is

$$\begin{array}{c}\hfill F_{\eta }(t)=2\eta {{\displaystyle \oint }}_{r=R_{\odot }}{\mathit{A}}_{\mathrm{p}}·\mathrm{\nabla }{B}_{\mathrm{r}}\mathrm{d}S=2\eta {{\displaystyle \oint }}_{r=R_{\odot }}(\mathrm{\nabla }·({\mathit{A}}_{\mathrm{p}}{B}_{\mathrm{r}})-B_{\mathrm{r}}\mathrm{\nabla }·{\mathit{A}}_{\mathrm{p}})\mathrm{d}S.\end{array}$$(A.2)

The first term vanishes by the divergence theorem (since the surface is closed), while the second term vanishes by our gauge condition ∇ ⋅ A_p = 0. Thus F_η(t) = 0 and so the flux in northern and southern hemispheres must be equal and opposite.

Differential rotation

The helicity flux from differential rotation is

$$\begin{array}{cc}\hfill F_{v_{\varphi }}(t)& =-2{{\displaystyle \oint }}_{r=R_{\odot }}{A}_{\mathrm{p}\varphi }{v}_{\varphi }{B}_{\mathrm{r}}{\mathrm{d}}^{2}x\hfill \\ \hfill & =-2R_{\odot }^3{{\displaystyle \int }}_0^{\pi }{{\displaystyle \int }}_0^{2\pi }\mathrm{\Omega }(\theta ){sin}^2\theta A_{\mathrm{p}\varphi }{B}_{\mathrm{r}}\mathrm{d}\varphi \mathrm{d}\theta .\hfill \end{array}$$(A.3)

Now,

$$\begin{array}{cc}& R_{\odot }{sin}^2\theta A_{\mathrm{p}\varphi }{B}_r={sin}^2\theta \frac{1}{sin\theta }[\frac{\partial }{\partial \theta }(sin\theta A_{\mathrm{p}\varphi })-\frac{\partial A_{\mathrm{p}\theta }}{\partial \varphi }]A_{\mathrm{p}\varphi }\hfill \\ \hfill & =\frac{1}{2}\frac{\partial }{\partial \theta }\left({sin}^2\theta A_{\mathrm{p}\varphi }^2\right)-\frac{\partial }{\partial \varphi }(sin\theta A_{\mathrm{p}\varphi }{A}_{\mathrm{p}\theta })+sin\theta A_{\mathrm{p}\theta }\frac{\partial A_{\mathrm{p}\varphi }}{\partial \varphi }\hfill \\ \hfill & =\frac{1}{2}\frac{\partial }{\partial \theta }\left({sin}^2\theta A_{\mathrm{p}\varphi }^2\right)-\frac{\partial }{\partial \varphi }(sin\theta A_{\mathrm{p}\varphi }{A}_{\mathrm{p}\theta })-sin\theta A_{\mathrm{p}\theta }\frac{\partial }{\partial \theta }(sin\theta A_{\mathrm{p}\theta })\hfill \\ \hfill & \text{[using}(\text{A}.1)]\hfill \\ \hfill & =\frac{1}{2}\frac{\partial }{\partial \theta }\left({sin}^2\theta A_{\mathrm{p}\varphi }^2\right)-\frac{\partial }{\partial \varphi }(sin\theta A_{\mathrm{p}\varphi }{A}_{\mathrm{p}\theta })-\frac{1}{2}\frac{\partial }{\partial \theta }({sin}^2\theta A_{\mathrm{p}\theta }^2).\hfill \end{array}$$

Integration in ϕ reduces ϕ-derivative term to zero (by the periodicity of A_p), such that F_vϕ(t) is given by

$$\begin{array}{c}\hfill F_{v_{\varphi }}(t)=-R_{\odot }^2{{\displaystyle \int }}_0^{\pi }{{\displaystyle \int }}_0^{2\pi }\mathrm{\Omega }(\theta )\frac{\partial }{\partial \theta }({sin}^2\theta (A_{\mathrm{p}\varphi }^2-A_{\mathrm{p}\theta }^2))\mathrm{d}\varphi \mathrm{d}\theta .\end{array}$$(A.4)

For the case of uniform rotation Ω(θ) = Ω₀ (constant), we see that the integral vanishes. However, when Ω is a function of θ this is not the case, but we note that the flux is small because of the relative smallness of Ω₂ and Ω₄ as compared to Ω₀ (see Eq. (9)). As such, to leading order, the helicity flux associated with differential rotation is balanced between the hemispheres.

All Tables

All Figures

	Fig. 1. Velocity profiles as a function of latitude for meridional velocity vθ (left), and differential rotation vϕ (right), in W18 and J11. The black dashed curve for vϕ shows the profile used by Berger & Ruzmaikin (2000) and also by Hawkes & Berger (2018).
In the text

	Fig. 2. Time variation of unsigned magnetic flux (a), axial dipole moment (b), and hemispheric helicity fluxes (c)–(f) for the W18 and J11 models over solar cycles 21–23. The net helicity flux FN/S (c) is the sum of those due to solar rotation Fvϕ, N/S (d), meridional velocity Fvθ, N/S (e) and supergranular diffusion Fη, N/S (f) for both W18 (blue shades) and J11 (red shades).
In the text

Fig. 3. Latitude–time distributions of longitude-averaged radial magnetic field (a,b) alongside longitude integrated total helicity flux (c,d), and helicity flux associated individually with differential rotation (e,f), meridional velocity (g, h) and supergranular diffusion (i,j) for W18 and J11 respectively. We stress that these latitudinal distributions are not in themselves physically meaningful, since there is typically significant cancellation between different latitudes. The same colour scales are used for W18 and J11. Units are G for (a,b) and Mx2 day−1 for (c–j).

In the text

	Fig. 4. Illustration of the source of helicity flux Fvϕ at a single time during November 1977 (during the early rise phase of Cycle 21), in W18. Panels show (a) Br, (b) Aϕ and (c) the integrand in (12).
In the text

	Fig. 5. Results for the full J11 model, showing unsigned magnetic flux (a), axial dipole moment (b), total helicity flux (c) and helicity flux from solar rotation (d), meridional velocity (e) and supergranular diffusion (f). In (c–f), light red shows the northern hemisphere and dark red the southern hemisphere.
In the text

	Fig. 6. Scatter plots of integrated helicity flux per cycle in each sub-hemisphere for both W18 and J11, plotted against the axial dipole strength (absolute value) at the end of each cycle. This time is defined as the time of minimum unsigned magnetic flux in Figs.1a and 2a. Lines of best fit and Pearson’s correlation coefficient are calculated for the J11 data.
In the text

	Fig. 7. Time variation of the |λ|< 55° helicity flux in the Northern hemisphere for W18 (red), compared to the unsigned magnetic flux (blue). The black line denotes the zero axis for helicity flux.
In the text