© ESO, 2016

1. Introduction

Magnetic helicity is well-known to be invariant in ideal magnetohydrodynamics (MHD), and in highly conducting plasmas it is almost conserved even for finite resistivity (Berger 1984; Pariat et al. 2015). Helicity may be interpreted as a net linking or winding of magnetic field lines around one another (Moffatt 1969). This linking can put a lower bound on the magnetic energy (Moffatt 1990; Freedman & He 1991; Berger 1993), reflecting the physical barrier that, in ideal MHD, magnetic field lines are unable to pass through one another or to reconnect.

In the corona, which is highly conducting and close to ideal, there are two primary sources of helicity, both coming from the solar interior: shearing of magnetic field lines by footpoint motions, and emergence of twisted magnetic fields. In this paper, we will consider only the production of helicity by footpoint motions. These footpoint motions may arise either from large-scale flows (primarily differential rotation), or from small-scale convection. Here we model the net effect of convection on the large-scale magnetic field with an isotropic diffusion on the solar surface. We neglect the additional helicity injection that could arise if the convective motions had a net sign of vorticity (Antiochos 2013; Mackay et al. 2014; Knizhnik et al. 2015).

Perhaps the most important practical consequence of helicity in the corona is the formation of twisted magnetic flux ropes, and their eruption as coronal mass ejections (e.g., Chen 2011). These remove helicity from the corona and send it out into the heliosphere. However, a difficulty in quantifying this notion arises because helicity is a volume integral, and is not conserved on an arbitrary sub-volume of the corona. Previous authors have quantified the helicity generated by solar rotation in two extreme cases: entire hemispheres (Berger & Ruzmaikin 2000), and a single isolated active region (DeVore 2000). The goal of this paper is to show how one can meaningfully study the spatial distribution of helicity within the corona, with the ultimate aim of understanding the origin of solar eruptions. The basic idea is to decompose the corona into infinitesimal tubular volumes around each magnetic field line. The helicity of each of these sub-volumes is an ideal invariant (provided that the field line endpoints are fixed), called the field line helicity.

This idea of field line helicity is not a new one, and goes back to Taylor (1974). Subsequently, Berger (1988) derived lower energy bounds based on field line helicity (see also Aly 2014; Yeates et al. 2014), but the concept was not significantly developed for a number of years. More recently, field line helicity has been found to be an invaluable tool for understanding the turbulent relaxation of braided magnetic fields. For cylindrical domains, Yeates & Hornig (2013, 2014) proved that knowing the field line helicity on each field line uniquely determines the field line mapping from one end of the cylinder to the other. And in sufficiently complex magnetic fields, Russell et al. (2015) showed that field line helicity is efficiently redistributed by reconnection, but is not destroyed on dynamical timescales. It therefore acts as a constraint on magnetic relaxation, demonstrated in the numerical experiments of Pontin et al. (2011).

The aim of this paper is to apply the tool of field line helicity to the Sun’s corona, in which the magnetic field has a rather more complex topology than the cylinder. Accordingly, the “completeness” proof of Yeates & Hornig (2014) no longer applies, though field line helicity retains its importance as a topological invariant.

Since the global magnetic field in the coronal volume can not be measured directly, a numerical model is required. However, for the study of helicity, a model more sophisticated than potential field extrapolation is needed. Primarily, this is because potential extrapolations lack the free magnetic energy that is associated with helicity. But it is also because they do not evolve continuously over time, so do not preserve the connectivity of magnetic field lines associated with an ideal evolution. In other words, a sequence of potential field extrapolations could “undo” the field line entanglement imposed in reality by the footpoint motions. Instead, in order to model the gradual injection of helicity over time, a time dependent model is required. Here we apply the magneto-frictional model (van Ballegooijen et al. 2000), as a compromise that retains sufficient physics but is less computationally expensive than full-MHD simulations. The importance of retaining a continuous time dependence has been demonstrated before, but will be shown rather clearly by the field line helicity.

The paper is organised as follows. Section 2 explains the physical interpretation of field line helicity, and gives a practical definition, then Sect. 3 describes the magneto-frictional model. We then study three situations of increasing complexity: a dipolar field (Sect. 4), a quadrupolar field (Sect. 5), and finally a more realistic, non-axisymmetric configuration (Sect. 6). We conclude in Sect. 7.

2. Field line helicity

We model the solar corona by a spherical shell D = { (r,θ,φ) | r₀<r<r₁ }. The field line helicity of a magnetic field line L ⊂ D is defined as (1)where e_B = B/ | B | is the unit vector aligned with the local direction of the magnetic field B = ∇ × A. It follows that is undefined on ergodic magnetic field lines, which have infinite length. In generic coronal magnetic fields, this situation does not usually arise since field lines are typically finite in length, ending on one or more of the domain boundaries r = r₀ and r = r₁. The choice of A will be discussed below.

Since there is a unique field line through each point (except for magnetic null points where B = 0), we can also assign values of to points x ∈ D, and think of as a function on D. This can be useful for visualization. This function is evidently constant along magnetic field lines, and will, furthermore, be continuous in regions of continuous field line mapping. In the presence of magnetic null points, will generally be discontinuous across their separatrix surfaces, like any field-line integrated quantity.

Although we have defined as a line integral, it may also be written as the limit (2)where D_ϵ is the magnetic flux tube of radius ϵ around the field line L, with Φ_ϵ being the flux of this tube. This motivates the name “field line helicity” (Berger 1988). Integrating Eq. (2)over all field lines, weighted by their flux, will recover the total helicity H = ∫_DA·B dV. In this sense, is a meaningful density for H, describing how topological sub-structure is distributed within D.

2.1. Physical interpretation as magnetic flux

The physical meaning of is clear when L is a closed curve such as L₁ in Fig. 1. In that case, Stokes’ theorem implies that is simply the magnetic flux that links through L. It follows that must be an ideal invariant when L is a closed curve. This conclusion does not depend on the chosen gauge of A, and indeed the value of is independent of this gauge (provided that A is single-valued).

Fig. 1 Physical interpretation of field line helicity as the magnetic flux linking a closed (L1) or open (L2) magnetic field line.

In the coronal situation, closed magnetic field lines are rare, and we must consider field lines ending on one or more of the domain boundaries (r = r₀ and r = r₁). For example, consider the field line L₂ in Fig. 1, linking two points x₁ and x₂ on r = r₀. Under a gauge transformation from A to A^′ = A + ∇χ, the field line helicity changes from to . However, whichever gauge is used, we can always find a corresponding surface S whose magnetic flux is exactly . This is equivalent to the existence of some curve γ between x₁ and x₂ (as in Fig. 1) such that ∫_γA· dl = 0, which is demonstrated in Appendix A.

In different gauges, the possible curves γ will differ, so that the physical fluxes represented by will depend on the chosen gauge. But, whichever gauge is chosen, we may interpret in terms of the physical linking of fluxes. Note that it is only the gauge on the boundary ∂D that matters; this is analogous to the situation with relative helicity (Berger & Field 1984; Prior & Yeates 2014). For field lines such as L₂, Antiochos (1987) defined a “flux-per-field-line” which is equivalent to if we take γ to be a geodesic on ∂D between x₁ and x₂. However, rather than specifying the curve γ explicitly like this, we will specify the gauge of A explicitly, as described in the next section. This makes easier to compute, although it will lead (in general) to more complex γ.

The interpretation in terms of linked fluxes shows that remains an ideal invariant in the case of non-closed field lines, provided that the gauge is fixed and, in addition, that there are no motions of field line footpoints on ∂D. Of course, an important feature of the corona is the injection of helicity through footpoint motions, which will necessarily lead to a change in . This will be demonstrated in Sects. 4–6.

2.2. Gauge choice

For practical application, we must choose a gauge for A. In this paper, we use the so-called DeVore gauge, chosen since it is straightforward to compute and leads to a clear physical interpretation for . The gauge was introduced by DeVore (2000) in Cartesian coordinates, and has been used by a number of authors in both Cartesian geometry (Valori et al. 2012; Moraitis et al. 2014) and spherical geometry (Amari et al. 2013).

The DeVore gauge condition is that A_r ≡ 0. From B = ∇ × A, we get (3)which may be integrated in r to give (4)Here A₀ is the vector potential on the initial surface r = r₀. Under our assumption that A_r ≡ 0, it follows that A₀ has only θ and φ components, which must satisfy (5)but are otherwise arbitrary. We will follow Amari et al. (2013) and fix A₀ with the condition ∇·A₀ = 0, so that it may be written as (6)Note that ψ is (up to a constant) the poloidal flux function from a poloidal-toroidal decomposition of B. Equation (5)then requires that (7)Solving this Poisson equation for ψ on the sphere determines A₀. The function ψ is determined only up to an additive constant, which we may fix by requiring that ∫_{r = r0}ψ dΩ = 0. In spherical geometry, the Green’s function for Eq. (7)is known, and the solution may be expressed analytically (Kimura & Okamoto 1987) as (8)where (9)In the special case of a potential field B = B_p (i.e., ∇ × B_p = 0), it follows from Eq. (4)that ∇·A_p = 0 everywhere in V. For a potential field in this (Coulomb) gauge, we have ∫_DA_p·B_p dV = 0 (Berger 1984), although need not vanish for any individual field line – we will see an example of this in Sect. 6.

2.3. Physical interpretation of the gauge choice

One advantage of the DeVore gauge Eq. (4)is its explicit physical interpretation. We consider the contributions to from both the integral term in Eq. (4)and the boundary term rA₀.

The integral term contributes to A_θ(r,θ,φ) when B_φ(r^′,θ,φ) is non-zero at some radius r^′ between r₀ and r. Similarly, it contributes to A_φ(r,θ,φ) when B_θ(r^′,θ,φ) is non-zero. So the contribution to from this term represents the net perpendicular magnetic flux beneath the field line concerned (Fig. 2). For a field line with both footpoints on r = r₀, this is rather like choosing the curve γ from Sect. 2.1 to be the radial projection γ^′ of the field line on r = r₀. Accordingly, this term will measure the twisting of magnetic field lines with height in the corona, and the net linking of flux beneath magnetic arcades. But it should be borne in mind that the projected curve γ^′ will generally have , and the additional contribution from the boundary term rA₀ is needed to ensure that gives an ideal-invariant flux.

Fig. 2 Physical interpretation of field line helicity in the DeVore gauge. The red circles show the direction of A0 (contours of ψ) arising from a strong magnetic source Br> 0. The shaded surfaces are radial projections of the field lines L1 and L2. Both field lines have a contribution to from any flux linking through these surfaces (owing to the second term of Eq. (4)), but only L1 has a contribution from A0, since the projection of L2 is perpendicular to A0.

With our choice of A₀ in Eq. (6), the boundary term contribution to represents the winding of coronal field lines around strong sources of magnetic flux on r = r₀. From Eq. (6), we see that the integral curves of A₀ are the curves of constant ψ (cf. Hornig 2006). Since ψ solves the Poisson equation with source term B_r(r₀,θ,φ), these curves are analogous to the surfaces of equal temperature in a solution to the heat equation where B_r(r₀,θ,φ) corresponds to a distribution of heat sources and sinks. Larger contributions to arise when field lines (in projection) are aligned with these curves, which encircle the sources of (locally) strongest | B_r |.

So, overall, in our gauge represents the net effect of two contributions: twisting of magnetic field lines with height, and winding around centres of strong flux on the boundary r = r₀. The examples of Sect. 6 suggest that both terms are generally significant.

Fig. 3 Illustration of the dipolar simulation with ν0 = 0.36 × 10-5 s-1 on days 0, 1, 2, and 20. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 0.5), and projected coronal magnetic field lines traced from height r = R⊙ are coloured (red/blue) according to , saturated at ± 0.01 Mx with white indicating .

3. Magneto-frictional model

To approximate the evolution of non-potential magnetic fields in the corona, we use the magneto-frictional model introduced by van Ballegooijen et al. (2000) and subsequently applied to the global corona by Yeates et al. (2008). In this model, the coronal magnetic field evolves through a continuous quasi-static sequence of approximately force-free equilibria, in response to continual shearing by photospheric footpoint motions. In this paper we use a uniform (but stretched) grid to cover the domain { r₀<r<r₁,θ₀<θ<θ₁,0 <φ< 2π } with a resolution 28 × 160 × 192. Here we take r₀ = R_⊙ (the photosphere), r₁ = 2.5 R_⊙, θ₀ = 0.05π, and θ₁ = 0.95π. Omitting the poles from our domain does not significantly affect the results presented here; the solution Eq. (8)for ψ remains valid provided that , which we impose in our initial condition.

3.1. Coronal evolution

In the magneto-frictional model, the coronal vector potential evolves according to the induction equation (10)where j = ∇ × B and η represents a turbulent resistivity arising from the cumulative effect of small-scale coronal flows. For simplicity in this paper, we follow Mackay & van Ballegooijen (2006) and set (11)where η₀ is a constant background value and the second term acts only in regions of strong current density | j | to limit the formation of unresolved gradients in B. An alternative would be to consider higher-order hyperdiffusion (as in Yeates 2014), but the simpler form suffices here.

The gauge Φ in Eq. (10)is, of course, arbitrary. For computation itself we use the Weyl gauge Φ ≡ 0, but for calculating the field line helicity we subsequently recompute the DeVore gauge A from B, as defined in Sect. 2.2.

The main simplification in the magneto-frictional method is to forego solving the full MHD equations and instead approximate the plasma velocity by (12)Here the first term is a friction-like term that enforces relaxation towards a force-free equilibrium. The factor | B | ² prevents relaxation from being inhibited in weak-field regions, although it must be limited away from zero near null points where | B | = 0. The coefficient ν has the same dimensions as η, and is set to ν = ν₀r²sin²θ. The second term in Eq. (12)is a radial outflow imposed only near the outer boundary. This term simulates (crudely) the effect of the solar wind in radially opening out the magnetic field lines, while allowing horizontal field to pass through the upper boundary if necessary.

Equations (10)and (12)are solved on a staggered grid (Yee 1966) using finite differences. Zero-gradient boundary conditions are imposed at r = r₁, and B_θ = 0 is imposed at θ = θ₀,θ₁. At r = r₀, we do not prescribe v according to Eq. (12), but rather determine ∂A_θ/∂t and ∂A_φ/∂t from our imposed photospheric driving. (No boundary condition on A_r is needed, owing to the staggered grid.) For a given photospheric driver, the coronal model is then determined by three parameters: ν₀, η₀ and v_out. The friction coefficient ν₀ controls the speed of coronal relaxation relative to the surface evolution, while η₀ controls the rate of coronal diffusion. Rather than η₀, we vary the dimensionless number , which measures the relative importance of diffusion compared to relaxation in the corona (cf. Cheung & DeRosa 2012). For this paper, we fix the radial outflow speed v_out = 100 km s^-1.

3.2. Photospheric driving

The magneto-frictional method simulates the evolution of the coronal magnetic field in response to shearing by surface motions. In this paper, we consider the effect of these motions on three different initial magnetic fields. The motions are modelled by a simple surface flux transport model (Sheeley 2005; Mackay & Yeates 2012; Jiang et al. 2014) in which, at r = r₀, we impose The first term in Eq. (13)represents differential rotation. The simulations are carried out in the carrington frame, and we choose the Snodgrass (1983) angular velocity (in degrees per day) (15)This implies that the coronal magnetic field is relaxing relative to the carrington frame, rather than to the background stars. The coefficient D = 600 km s^-1 represents “supergranular diffusion” of B_r, namely the net large-scale effect of the random walk of magnetic elements under supergranular convection on the solar surface. For illustrative purposes, we neglect other flux transport effects such as meridional flow, as well as the emergence of new magnetic flux.

4. Dipolar field

Our first simulation is initialized with a potential field extrapolation from the photospheric boundary condition (16)along with B_θ = 0 on the boundaries θ = θ₀,θ₁ and B_θ = B_φ = 0 on the outer boundary r = r₁. The potential field is computed using the eigenfunction method of van Ballegooijen et al. (2000). The coronal magnetic field is then evolved with magneto-friction, as described in Sect. 3. For this example the coronal field remains close to potential, with low electric currents, so the results are insensitive to η₀. Accordingly we will illustrate only the effect of varying ν₀, while holding the dimensionless ratio fixed at 2.89 × 10^-5 (a typical value from previous simulations). For the first day of evolution, no photospheric motions are applied, so as to illustrate the effect of switching on differential rotation from day 1 onwards.

Fig. 4 Various integrated quantities as a function of time, for the dipolar simulations with different ν0 (indicated by line styles). Panel a) shows the open magnetic flux ∫r = r1 | Br | dΩ, panel b) shows ∫D | j | dV, and panel c) shows HN (asterisks) and HS (circles). Panel d) shows the terms in Eq. (18)for the northern hemisphere, with asterisks denoting S0, circles S1, squares Seq, and diamonds SV.

The evolution of the magnetic field structure for one of the runs is shown in Fig. 3. The evolution is straightforward: firstly there is an opening out of the magnetic field, due to the radial outflow at the upper boundary. This expansion takes approximately 1 day, and creates electric currents near the outer boundary associated with the extended “streamer” structure at the equator. Once the surface motions are switched on, the field then relaxes to a dynamical equilibrium between the footpoint shearing and the magneto-frictional relaxation. As we will discuss below, there is non-zero field line helicity associated with this dynamical equilibrium (shown by the colours in Fig. 3). The time taken to reach equilibrium depends on ν₀.

Indeed, it is instructive to consider the effect of ν₀ (the rate of frictional relaxation) on the evolution. Figure 4 shows a number of integrated quantities, as a function of time for four runs with different ν₀. As seen in Fig. 4a, weaker friction allows the radial outflow to open out the field further, leading to more open magnetic field lines, although this is not a particularly strong effect. More striking, in this example, is the increased coronal electric current that weaker friction allows (Fig. 4b). A small part of this difference in current arises from the greater initial expansion of the field, but most arises from the character of the dynamical equilibrium. When friction is weaker, the field lines relax back less in response to shearing of their footpoints by differential rotation, so more current is stored in the corona. The shearing of field lines in the equilibrium is actually rather small and hard to discern in Fig. 3, although it is just visible at the south pole on day 20.

Fig. 5 Latitudinal distribution of field line helicity for the dipolar simulations, showing how the peak value tends to zero as the friction parameter ν0 is successively doubled. A log-log fit shows that .

Next we consider the evolution of magnetic helicity. In this unusually symmetric situation, it is helpful to consider the net helicity in each hemisphere, (17)By symmetry these are equal and opposite (so that the total helicity vanishes). They are shown in Fig. 4c. Before the surface motions are switched on there is no helicity, since B_φ ≡ 0 and A_θ ≡ 0. After the motions are switched on, the helicity increases to a steady value in each hemisphere. It is clear from Fig. 4c that this steady value is larger when the friction is weaker, in accordance with the greater shear of the equilibrium field lines.

The equilibrium distribution of field line helicity is shown both by the colours in Fig. 3 and, as a function of latitude, in Fig. 5. It is clear that the helicity in each hemisphere is not distributed uniformly among all field lines, but is stored only on open field lines. This is due to the symmetry of the configuration: closed field lines cross the equator and pick up equal and opposite contributions to from each hemisphere.

It is also interesting to consider the helicity flux through the boundaries. To calculate this, it is most convenient to use the form (18)which is easily derived using Faraday’s law (19)This is valid for the helicity in any subdomain V, whether magnetically closed or not. Here we use A computed in our gauge to estimate ∂A/∂t, and we also record the electric field E during the simulation. When we apply this formula to H_N (or H_S), we obtain four contributions: the volume dissipation term S_V, and three contributions to the surface integral from different boundaries, namely S₀ (r = r₀), S₁ (r = r₁), and S_eq (θ = π/ 2), as in Fig. 6. There is no contribution from the latitudinal boundaries θ = θ₀ and θ = θ₁ owing to our boundary conditions in the simulation. Figure 4d shows these four contributions for each of the dipole simulations, for the northern hemisphere. (The southern hemisphere contributions are equal and opposite.) Firstly, the volume dissipation term is small compared to the surface terms. The main contributions are an injection S₀ of helicity through r = r₀, by differential rotation, and an output S₁ through the upper boundary. The latter would correspond to winding up of the solar wind (the Parker spiral). However, the helicity output is rather less than the input (0.0026 Mx² day^-1 compared to 0.0032 Mx² day^-1 for ν₀ = 0.36 × 10^-5 s^-1). The difference is accounted for by S_eq, which represents a net transfer of helicity across the equator on closed field lines. During the relaxation phase, there is a slight imbalance between these terms, allowing the equilibrium helicity to build up in each hemisphere. The overall flow of helicity is summarised in Fig. 6.

Note that, as the friction parameter ν₀ is increased, the stored field line helicity in each hemisphere, along with H_N and H_S, tends to zero approximately as . However, helicity is injected by differential rotation through the photosphere at a constant rate S₀ that is independent of ν₀. Figure 4d shows that the lack of stored helicity is compensated by the other surface terms S₁ and S_eq during the relaxation phase. Even for the finite values of ν₀ considered here, the stored helicity in the corona is little more than the helicity injected in a single day by differential rotation. However, we will see in the subsequent examples that much more helicity can be stored if we break the symmetry of the magnetic configuration.

Fig. 6 Schematic of helicity flow in the dipolar example, for the northern hemisphere. Arrows show the direction of positive helicity transfer as measured by the surface terms S0, S1, and Seq in Eq. (18).

Fig. 7 Illustration of the quadrupolar simulation with ν0 = 0.36 × 10-5 s-1 on days 0, 20, 66, and 68. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 2 G), and projected coronal magnetic field lines traced from height r = 1.2 R⊙ are coloured (red/blue) according to , saturated at ± 0.2 Mx with white indicating .

It is interesting to note that the sign of helicity injected into the solar wind is opposite to that of Berger & Ruzmaikin (2000), who estimated the injection of helicity into the volume r>r₀ by solar rotation. Figure 6 shows that the outward helicity flux in the northern hemisphere is positive in our example, since S₁ is negative (Fig. 4d). This sign is opposite to the direction of winding of the Parker spiral. However, this is an apparent difference caused by our use of the carrington frame. If the constant 27-day rotation rate were added back in, the sign would reverse.

5. Quadrupolar field

In more realistic configurations, differential rotation is able to build up field line helicity on closed field lines. Our second axisymmetric example gives a simple demonstration of this process, starting from a potential field extrapolated from the photospheric distribution (20)With B₁ = B₂ = 100, the bipolar rings each contain the same unsigned magnetic flux 2πd², provided that they overlap neither each other nor the poles. We take d = 0.1, and locate them at θ₁ = 0.3π and θ₂ = 0.55π (illustrated in Fig. 7). Since the rings are asymmetrically placed with respect to the equator, we expect differential rotation above each PIL to build up helicity at different rates, even though both rings contain the same magnetic flux. The same photospheric motions are imposed as in Sect. 4, except that they are switched on immediately. A slightly different value is used, although we will also consider the effect of varying this parameter below.

Figure 7 shows how the magnetic field evolves over 68 days, while Fig. 8 shows various integrated quantities, analogous to Fig. 4. The most striking difference from the dipolar case is that the quadrupolar system does not reach a dynamical equilibrium, in spite of the fact that the rate S₀ of helicity injection by differential rotation remains fairly constant, albeit higher than before owing to the greater magnetic flux on r = r₀. (The slight decay in S₀ over time arises from diffusive decay of the more concentrated photospheric field, visible in Fig. 7.) Instead, current and helicity continue to be injected into the corona. The open flux does initially level off, but then increases as the magnetic arcades are sheared and energised.

Fig. 8 Various integrated quantities as a function of time, for the quadrupolar simulations with different ν0 (indicated by line styles). The format is the same as Fig. 4. For clarity, panel d) shows only the run with ν0 = 0.36 × 10-5 s-1, and only for the northern hemisphere, although the hemispheres are no longer symmetric. The peak value of S1 during the flux rope eruption is not shown, and is much larger, about − 1.16 Mx2 day-1.

It is clear from Fig. 7 that the additional helicity is primarily stored along closed field lines, particularly those that do not cross the equator. This arises because the footpoints are no longer symmetric about the equator, so that differential rotation shears the magnetic loops. The sign injected is opposite for the arcades in each hemisphere. More helicity is injected in the northern hemisphere, simply because the two bipolar rings are asymmetrically placed and the northern ring lies at a latitude with greater shear in the differential rotation. This asymmetry also leads, eventually, to negative helicity in the arcade straddling the equator (Fig. 7). For finite ν₀, the open field lines near the poles actually store a similar amount of field line helicity as in the dipolar example, but this is insignificant compared to that stored at lower latitudes. Moreover, this lower latitude helicity is almost independent of ν₀, since it is enforced topologically by the footpoint motions and cannot be removed by ideal relaxation, however rapid. As the open flux increases gradually with energisation of the field, the helicity output S₁ increases, consistent with the dipolar example where open field lines act as continuous “conduits” of helicity from the photosphere out to the solar wind. The increasing proportion of open field lines leads to the levelling off of hemispheric helicity from about day 30 onwards. In this example, the cross-equatorial helicity flux S_eq is negligible.

We remark that the sign of helicity in each hemisphere, in this example, is opposite to the typical hemispheric pattern of helicity on the Sun, which is negative in the northern hemisphere and positive in the south (Pevtsov & Balasubramaniam 2003). This arises from the East-West orientation of the polarity inversion line in our axisymmetric model. On the real Sun, polarity inversion lines at active latitudes are often aligned North-South, so that differential rotation injects helicity of the observed majority sign. This was illustrated by the simulations of DeVore (2000) in Cartesian geometry, and Yeates & Mackay (2009b) in spherical geometry.

As is evident in Fig. 8, the amount of electric current and helicity in the corona does not build up indefinitely, but is suddenly reduced on about day 67 of the simulation. This sudden reduction results from ejection of the magnetic flux rope that forms above the northern polarity inversion line. The flux rope is visible on day 66 in Fig. 7, but has been ejected through the outer boundary r = r₁ by day 68, leaving only a weakly sheared arcade behind it. The mechanism by which the flux rope forms is well understood (van Ballegooijen & Martens 1989); essentially, it is a combination of reconnection of sheared magnetic loops accompanied by flux cancellation due to supergranular diffusion on r = r₀, which leaves horizontal magnetic field in the corona above polarity inversion lines. Due to the symmetry in this rather artificial example, the flux rope that forms is detached from the photosphere, encircling the whole Sun. Since the corresponding magnetic field lines are either closed or ergodic (infinite length), the field line helicity in the rope is undefined. (In Fig. 7, the colour scale is saturated at ± 0.2 Mx.) Nevertheless, integrating for a finite length clearly indicates the location of the rope. As is evident in Fig. 8d, the eruption causes a very high, sudden, spike in the helicity output S₁, and a consequent sudden reduction in the total helicity H_N in the northern hemisphere (Fig. 8c). After the eruption, the helicity begins to build up again since the footpoint shearing continues.

Finally we consider the effect of the simulation parameters ν₀ and η₀. We have already seen that the hemispheric helicity and open flux depend only weakly on ν₀. We have also run simulations with different but the same value of ν₀ = 0.36 × 10^-5 s^-1. We find that increasing from 1.73 × 10^-5 to 6.94 × 10^-5 delays the flux rope eruption by three days, but has little impact on the magnitude of open flux or hemispheric helicity overall. This is consistent with the findings of Yeates & Mackay (2009a), who showed that higher diffusion limits the speed at which highly twisted flux ropes are able to form, by dissipating the concentrated electric currents in the ropes. This led to a lower eruption rate.

This example has shown how field line helicity reveals the storage of helicity on closed magnetic field lines in the corona, as well as the sudden expulsion of this helicity in the form of flux rope eruptions. In the next section, we see these processes at work in a more realistic global configuration.

6. Non-axisymmetric field

Fig. 9 ADAPT maps of Br on r = r1, used to generate the initial potential field extrapolations for periods A and B (white positive, black negative, saturated at ± 20 G).

Our final example is a more realistic global magnetic configuration. For simplicity, we still consider only driving by large-scale surface motions, and continue to neglect the emergence of new magnetic flux.

Two simulations are presented: period A and period B. Each starts from a potential field extrapolation, as before, but now these are computed from full-surface B_r maps modelling the real Sun on two dates: 2011-Jan.-01 and 2015-Mar.-10 (Fig. 9). The maps are taken from the Air Force Data-Assimilative Photospheric Flux Transport (ADAPT) model (Arge et al. 2010; Henney et al. 2012; Hickmann et al. 2015), which assimilates observed magnetograms for the visible side of the Sun into a surface flux transport model. Here we simply take the two ADAPT maps as initial conditions for our two simulations. The maps in question come from ADAPT runs based on GONG magnetograms, and have been remapped to our simulated grid. A multiplicative flux correction has been applied to ensure flux balance. The choice of ADAPT maps, as opposed to any other model, is not particularly important; we simply wanted a realistic distribution of magnetic flux on the full solar surface. Period B represents a more active time in the solar cycle than period A, with considerably larger total flux. By this time in early 2015, the polar fields visible in period A have been almost completely removed by cancellation with magnetic flux from Cycle 24 active regions.

The magneto-frictional simulations use the same parameters as Sect. 5, driven by the same differential rotation and supergranular diffusion, with the parameters ν₀ = 0.36 × 10^-5 s^-1 and . The evolution is followed for much longer, up to 180 days. Whilst it is unrealistic to evolve the global magnetic field for so long without any new flux emergence, our purpose is to explore how the field line helicity responds to photospheric motions. On the Sun, the large-scale magnetic fields at higher latitudes do indeed result from many months of evolution of old active region fields (Petrie 2015).

Figure 10 shows four snapshots of the magnetic field during each period, along with the field line helicity on a grid of magnetic field lines. The effects of both differential rotation and of supergranular diffusion are apparent on r = r₀, where the pattern of B_r is both sheared and significantly smoothed out, removing smaller features. Once again, we see how field line helicity is injected and stored in closed magnetic arcades. However, the distribution of among closed field lines is far from uniform. The amount of field line helicity stored in any particular magnetic arcade is dependent on the degree of shearing of the arcade, which depends both on the orientation and the B_r pattern on r = r₀ (see also Yeates & Mackay 2009b). In fact, the non-uniform distribution of magnetic flux across the solar surface means that some field lines have non-zero even in the initial potential field (day 0 for each period in Fig. 10). But much stronger field line helicity builds up at particular locations where the field orientation is favourable to shearing by differential rotation. This often reverses the initial sign of at a particular location (e.g., around 300^° longitude in the southern hemisphere between days 0 and 60 of period B). As in the quadrupolar example (Sect. 5), the largest values of lie in twisted magnetic flux ropes.

Fig. 10 Projected magnetic field lines in the period A (left column) and B (right column) simulations, on days 0, 60, 120, and 180. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 10 G), and projected coronal magnetic field lines traced from height r = r0 are coloured (red/blue) according to , saturated at ± 0.5 Mx. Animated versions of this figure for periods A and B are available online.

Fig. 11 Various integrated quantities as a function of time, for the non-axisymmetric simulations (periods A and B). Panel a) shows the total photospheric magnetic flux ∫r = r0 | Br | dΩ, panel b) shows the total open flux ∫r = r1 | Br | dΩ, panel c) shows ∫D | j | dV, panel d) shows HN (asterisks) and HS (circles), and panel e) shows the root-mean-square field line helicity . The vertical grey lines indicate times of strong flux rope ejections, as explained in the text.

Figure 11 shows how global quantities evolve in the simulations for periods A and B. In the initial map there is about twice as much flux in period B as in period A, leading to correspondingly higher open flux, total current, and hemispheric helicity throughout the simulation. The total current and open flux increase gradually over about the first 60 days as the coronal field is energised, before reaching saturation and then gradually decaying over the rest of the simulation (owing to the decaying photospheric flux). The root-mean-square also takes about 2 months to reach its maximum value, although this does not seem to decay over the remainder of the simulation. This indicates how helicity is stored in the coronal magnetic field through memory of the footpoint motions. The hemispheric helicities are harder to interpret, highlighting the greater utility of as a diagnostic when helicity is non-uniformly distributed through the corona. Nevertheless, it is generally true that the helicity has greater magnitude in period B than in period A, commensurate with the greater flux.

A significant feature of the evolution are the multiple flux rope ejections that occur over the 180-day simulations. These are visible as transient peaks in the open flux (Fig. 11b), similar to the quadrupolar example (Fig. 8) although less pronounced owing to their more localised nature. There are many more ejections in period B than in period A, due to the more complex corona in period B. The vertical grey lines in Fig. 11 indicate the times of significant flux rope ejections. These have been determined not from the open flux, but from monitoring the horizontal magnetic field at the outer boundary r = r₁. This is enhanced significantly during the ejection of flux ropes, as is shown by the left column of Fig. 12. This shows a running difference of at r = r₁, for a particular ejection during period A. The ejection times shown in Fig. 11 were found by computing the number of grid points on each day with B_⊥> 0.05 G day^-1, then identifying local maxima in this time series.

The second column of Fig. 12 shows the distribution of on the solar surface r = r₀. The black circles are the footpoints of field lines with B_⊥> 0.05 G day^-1, so represent the footpoints of the ejected flux rope. It is clear that the region of strongest is the erupting flux rope (and its overlying arcade). In fact, this flux rope is clearly seen in Fig. 10. Further evidence that these are the footpoints of the erupting rope comes from the significant weakening of in this region following ejection of the rope. Similar behaviour is found for all of the ejections in periods A and B.

Finally, consider the right-most column of Fig. 12. This shows the dimensionless quantity (21)which Liu et al. (2016) call the twist number. In a force-free field, which is approximately the case in our model, we have j = αB with α constant along each field line, so Tw is simply α × length(L). This measure is also an indicator of where twisted structures are located within the magnetic field, clearly identifying the erupting flux rope in Fig. 12. For this rope, and Tw agree that this is the most significant twisted structure present, and agree on the sign of twist. But in general, Tw and have different relative magnitude and sign. Partly, the difference in relative magnitude can be explained by the fact that has units of magnetic flux while Tw is dimensionless. Put simply, a flux rope with the same field line curves but lower field strength would have the same Tw, but weaker . This accounts for the lower values of at high latitudes in Fig. 12, because this is a weak field region. However, significant differences in and Tw also arise because Tw depends only on the local twist around a single field line, whereas is a more global quantity. This tends to give a smoother distribution within each magnetic subdomain, as is evident in Fig. 12. Overall, there is a significant correlation between and Tw, although the (rank) correlation coefficient is only about 0.6. (This value remains steady after an initial transient phase of about 24 days where the correlation is lower.) Perhaps the most compelling reason to use rather than Tw is that Tw is not an ideal invariant (see Moffatt & Ricca 1992; Berger & Prior 2006).

7. Conclusion

We have shown how field line helicity is an invaluable tool for quantifying the distribution of topological structure within the Sun’s corona. It is straightforward to compute from a 3D magnetic field, by first computing an appropriate vector potential. It is a physically meaningful measure representing the linkage of magnetic flux around each magnetic field line in the domain. In particular, it is invariant under ideal motions within the domain, provided that the field line footpoints on the boundary remain fixed.

Fig. 12 Example of a flux rope ejection from period A. From top to bottom, the rows show days 122, 124, 126, and 128. The left column shows the (absolute) running daily difference of horizontal field at the outer boundary r = r1. The middle column shows the distribution of on r = r0 (saturated at ± 2 Mx), and the right column shows the distribution of Tw at r = r0 (saturated at ± 40). The dashed lines in the left column show the neutral line where Br(r1,θ,φ) = 0. Black circles in the other columns identify footpoints of field lines traced down from locations at r = r1 where the running difference of B⊥ exceeds 0.05 G day-1. Animated versions of this figure for both periods A and B are available online.

Although the value of for a given field line is computed by integrating A along that single field line, the vector potential A being integrated is fundamentally a non-local quantity. This enables to measure the linking with other magnetic field lines, but it does mean that knowledge of the wider magnetic field is required even to compute on a single field line.

We mentioned, in Sect. 2, that is a meaningful density for the total magnetic helicity, being the limiting helicity on an infinitesimal tubular domain around each magnetic field line. In fact, this is the finest possible decomposition of magnetic helicity into subdomains that will remain ideal invariants. Any finer decomposition would necessarily have interfacial surfaces in the corona with B·n ≠ 0, across which there would be helicity fluxes even in an ideal evolution.

Although the decomposition into field lines has an infinite number of subdomains, we have seen (e.g., Fig. 12) that the distribution of tends to be rather smooth, on account of its non-local definition. One could therefore give a first-order characterisation of the magnetic structure by integrating over discrete topological subdomains, following decomposition of the magnetic skeleton (e.g., Haynes & Parnell 2010). We have not pursued this idea here, as identifying the skeleton is computationally challenging in non-potential fields (cf. Edwards et al. 2015). However, a similar idea was proposed by Longcope & Malanushenko (2008), who defined the “additive self-helicity” of a sub-domain. Computations for simulations of a twisted magnetic flux tube were able to relate this quantity to the stability of the flux tube (Malanushenko et al. 2009).

The gauge dependence of arises purely from the fact that coronal magnetic field lines end on the boundaries rather than being closed loops. This gauge dependence is unavoidable; however, we have shown that every gauge is physically meaningful, corresponding to a different definition of what it means for flux to be linked with a magnetic loop. We have suggested that the DeVore gauge is a practical choice where not only is A is easy to compute, but the resulting field line helicity is appropriate for measuring twisted structures forming in the lower corona. An alternative way to choose a gauge would be to fix a vector potential where A × n matches some chosen reference field on the boundary, as in the commonly-used relative helicity (Berger & Field 1984). But really the choice of reference field is just another way of viewing the choice of gauge (Prior & Yeates 2014).

Having shown that field line helicity is a useful tool for coronal simulations, there are many possible future applications. An obvious one is to try to identify the locations where flux rope eruptions will occur, but there are many others. For example, we have, in this paper, neglected the direct emergence of already-twisted structures from the solar interior, and we have also neglected the net injection of helicity by small-scale convective motions. The relative importance of these two effects compared to surface shearing is important to establish, particularly for explaining the hemispheric pattern of helical structures in the corona (Pevtsov & Balasubramaniam 2003). It will also be needed in order to make improved estimates of the Sun’s helicity output over the solar cycle (cf. DeVore 2000). Another application will be to compare different methods of simulating the coronal magnetic field evolution – for example, how accurate is the magneto-frictional approximation? What is the importance of including thermodynamics? What is the effect of different parametrizations of turbulent diffusion in the corona? Or how best do we drive coronal simulations based on limited photospheric data (e.g., Kazachenko et al. 2014)? We hope to address some of these questions in future research.

Acknowledgments

This work was supported by STFC consortium grant ST/K001043/1 to the universities of Dundee and Durham. A.R.Y. also thanks the US Air Force Office of Scientific Research for support through a grant from the Basic Research Initiative “Understanding the interaction of CMEs with the solar-terrestrial environment”. We are grateful to Carl Henney for supplying the ADAPT maps for Sect. 6, and thank Alexander Russell, Christopher Prior, and Mitchell Berger for useful discussions.

References

Appendix A: Existence of a surface with flux

Consider the magnetic field line L₂ in Fig. 1, whose endpoints x₁ and x₂ both lie on the boundary r = r₀. We will show that there exists a curve γ from x₁ to x₂, lying on the surface r = r₀, such that (A.1)This means that the surface bounded by L and γ has flux equal to .

To see that such a curve exists, suppose that we continuously deform the original curve γ into either γ₊ or γ₋, as shown in Fig. A.1. Since x₂ is a field line footpoint, we must have B_r(x₂) ≠ 0. If B_r(x₂) > 0, then the curve γ₊ will have a larger value of ∫A· dl than γ, and the curve γ₋ will have a smaller value. By further deforming these curves to encircle x₂ more than once, we may ensure that and . By continuity, there must exist some intermediate curve with vanishing integral.

It is easy to see that there are many such curves γ with the required property, for any pair of footpoints x₁, x₂, and even if the gauge of A is fixed. But all of the corresponding surfaces will have the same flux , and this will be an ideal invariant if footpoint motions are disallowed.

Clearly this argument applies equally if both footpoints lie on r = r₁ (a rarer situation in the corona). But what about an

open field line, where x₁ lies on r = r₀ and x₂ on r = r₁? Now the curve γ that completes the loop must pass through D, rather than lying on the boundary. But, provided this portion of γ is chosen to be a magnetic field line, the resulting surface will again have an ideal-invariant flux. Again, this can be made equal to (the field line helicity of the original field line) by appropriately choosing the portions of γ on the two boundaries. So still represents an ideal-invariant flux, even if the field line is open.

Fig. A.1 The original and deformed curves, all on the boundary r = r0.

Online material

Movie of Fig. 10, left column Access here

Movie of Fig. 10, right column Access here

Movie of Fig. 12, period A Access here

Movie of Fig. 12, period B Access here

All Figures

	Fig. 1 Physical interpretation of field line helicity as the magnetic flux linking a closed (L1) or open (L2) magnetic field line.
In the text

Fig. 2 Physical interpretation of field line helicity in the DeVore gauge. The red circles show the direction of A0 (contours of ψ) arising from a strong magnetic source Br> 0. The shaded surfaces are radial projections of the field lines L1 and L2. Both field lines have a contribution to from any flux linking through these surfaces (owing to the second term of Eq. (4)), but only L1 has a contribution from A0, since the projection of L2 is perpendicular to A0.

In the text

	Fig. 3 Illustration of the dipolar simulation with ν0 = 0.36 × 10-5 s-1 on days 0, 1, 2, and 20. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 0.5), and projected coronal magnetic field lines traced from height r = R⊙ are coloured (red/blue) according to , saturated at ± 0.01 Mx with white indicating .
In the text

Fig. 4 Various integrated quantities as a function of time, for the dipolar simulations with different ν0 (indicated by line styles). Panel a) shows the open magnetic flux ∫r = r1 | Br | dΩ, panel b) shows ∫D | j | dV, and panel c) shows HN (asterisks) and HS (circles). Panel d) shows the terms in Eq. (18)for the northern hemisphere, with asterisks denoting S0, circles S1, squares Seq, and diamonds SV.

In the text

	Fig. 5 Latitudinal distribution of field line helicity for the dipolar simulations, showing how the peak value tends to zero as the friction parameter ν0 is successively doubled. A log-log fit shows that .
In the text

	Fig. 6 Schematic of helicity flow in the dipolar example, for the northern hemisphere. Arrows show the direction of positive helicity transfer as measured by the surface terms S0, S1, and Seq in Eq. (18).
In the text

	Fig. 7 Illustration of the quadrupolar simulation with ν0 = 0.36 × 10-5 s-1 on days 0, 20, 66, and 68. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 2 G), and projected coronal magnetic field lines traced from height r = 1.2 R⊙ are coloured (red/blue) according to , saturated at ± 0.2 Mx with white indicating .
In the text

Fig. 8 Various integrated quantities as a function of time, for the quadrupolar simulations with different ν0 (indicated by line styles). The format is the same as Fig. 4. For clarity, panel d) shows only the run with ν0 = 0.36 × 10-5 s-1, and only for the northern hemisphere, although the hemispheres are no longer symmetric. The peak value of S1 during the flux rope eruption is not shown, and is much larger, about − 1.16 Mx2 day-1.

In the text

	Fig. 9 ADAPT maps of Br on r = r1, used to generate the initial potential field extrapolations for periods A and B (white positive, black negative, saturated at ± 20 G).
In the text

Fig. 10 Projected magnetic field lines in the period A (left column) and B (right column) simulations, on days 0, 60, 120, and 180. Greyscale shading on r = r0 shows Br (white positive, black negative, saturated at ± 10 G), and projected coronal magnetic field lines traced from height r = r0 are coloured (red/blue) according to , saturated at ± 0.5 Mx. Animated versions of this figure for periods A and B are available online.

In the text

Fig. 11 Various integrated quantities as a function of time, for the non-axisymmetric simulations (periods A and B). Panel a) shows the total photospheric magnetic flux ∫r = r0 | Br | dΩ, panel b) shows the total open flux ∫r = r1 | Br | dΩ, panel c) shows ∫D | j | dV, panel d) shows HN (asterisks) and HS (circles), and panel e) shows the root-mean-square field line helicity . The vertical grey lines indicate times of strong flux rope ejections, as explained in the text.

In the text

Fig. 12 Example of a flux rope ejection from period A. From top to bottom, the rows show days 122, 124, 126, and 128. The left column shows the (absolute) running daily difference of horizontal field at the outer boundary r = r1. The middle column shows the distribution of on r = r0 (saturated at ± 2 Mx), and the right column shows the distribution of Tw at r = r0 (saturated at ± 40). The dashed lines in the left column show the neutral line where Br(r1,θ,φ) = 0. Black circles in the other columns identify footpoints of field lines traced down from locations at r = r1 where the running difference of B⊥ exceeds 0.05 G day-1. Animated versions of this figure for both periods A and B are available online.

In the text

	Fig. A.1 The original and deformed curves, all on the boundary r = r0.
In the text