The global distribution of magnetic helicity in the solar corona^{⋆}
^{1} Department of Mathematical Sciences, Durham University, Durham, DH1 3LE, UK
email:
anthony.yeates@durham.ac.uk
^{2} Division of Mathematics, University of Dundee, Dundee, DD1 4HN, UK
Received: 15 June 2016
Accepted: 5 August 2016
By defining an appropriate field line helicity, we apply the powerful concept of magnetic helicity to the problem of global magnetic field evolution in the Sun’s corona. As an idealmagnetohydrodynamic invariant, the field line helicity is a meaningful measure of how magnetic helicity is distributed within the coronal volume. It may be interpreted, for each magnetic field line, as a magnetic flux linking with that field line. Using magnetofrictional simulations, we investigate how field line helicity evolves in the nonpotential corona as a result of shearing by largescale motions on the solar surface. On open magnetic field lines, the helicity injected by the Sun is largely output to the solar wind, provided that the coronal relaxation is sufficiently fast. But on closed magnetic field lines, helicity is able to build up. We find that the field line helicity is nonuniformly distributed, and is highly concentrated in twisted magnetic flux ropes. Eruption of these flux ropes is shown to lead to sudden bursts of helicity output, in contrast to the steady flux along the open magnetic field lines.
Key words: magnetic fields / magnetohydrodynamics (MHD) / Sun: magnetic fields / Sun: coronal mass ejections (CMEs) / Sun: corona
Movies are available at http://www.aanda.org
© ESO, 2016
1. Introduction
Magnetic helicity is wellknown 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 largescale flows (primarily differential rotation), or from smallscale convection. Here we model the net effect of convection on the largescale 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 subvolume 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 subvolumes 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 magnetofrictional model (van Ballegooijen et al. 2000), as a compromise that retains sufficient physics but is less computationally expensive than fullMHD 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 magnetofrictional model. We then study three situations of increasing complexity: a dipolar field (Sect. 4), a quadrupolar field (Sect. 5), and finally a more realistic, nonaxisymmetric configuration (Sect. 6). We conclude in Sect. 7.
2. Field line helicity
We model the solar corona by a spherical shell D = { (r,θ,φ)  r_{0}<r<r_{1} }. 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_{0} and r = r_{1}. 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 fieldline 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 = ∫_{D}A·B dV. In this sense, is a meaningful density for H, describing how topological substructure 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_{1} 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 singlevalued).
Fig. 1 Physical interpretation of field line helicity as the magnetic flux linking a closed (L_{1}) or open (L_{2}) magnetic field line. 

Open with DEXTER 
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_{0} and r = r_{1}). For example, consider the field line L_{2} in Fig. 1, linking two points x_{1} and x_{2} on r = r_{0}. 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_{1} and x_{2} (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_{2}, Antiochos (1987) defined a “fluxperfieldline” which is equivalent to if we take γ to be a geodesic on ∂D between x_{1} and x_{2}. 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 nonclosed 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 socalled 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_{0} is the vector potential on the initial surface r = r_{0}. Under our assumption that A_{r} ≡ 0, it follows that A_{0} has only θ and φ components, which must satisfy (5)but are otherwise arbitrary. We will follow Amari et al. (2013) and fix A_{0} with the condition ∇·A_{0} = 0, so that it may be written as (6)Note that ψ is (up to a constant) the poloidal flux function from a poloidaltoroidal decomposition of B. Equation (5)then requires that (7)Solving this Poisson equation for ψ on the sphere determines A_{0}. 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 ∫_{D}A_{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_{0}.
The integral term contributes to A_{θ}(r,θ,φ) when B_{φ}(r^{′},θ,φ) is nonzero at some radius r^{′} between r_{0} and r. Similarly, it contributes to A_{φ}(r,θ,φ) when B_{θ}(r^{′},θ,φ) is nonzero. 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_{0}, this is rather like choosing the curve γ from Sect. 2.1 to be the radial projection γ^{′} of the field line on r = r_{0}. 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_{0} is needed to ensure that gives an idealinvariant flux.
Fig. 2 Physical interpretation of field line helicity in the DeVore gauge. The red circles show the direction of A_{0} (contours of ψ) arising from a strong magnetic source B_{r}> 0. The shaded surfaces are radial projections of the field lines L_{1} and L_{2}. Both field lines have a contribution to from any flux linking through these surfaces (owing to the second term of Eq. (4)), but only L_{1} has a contribution from A_{0}, since the projection of L_{2} is perpendicular to A_{0}. 

Open with DEXTER 
With our choice of A_{0} in Eq. (6), the boundary term contribution to represents the winding of coronal field lines around strong sources of magnetic flux on r = r_{0}. From Eq. (6), we see that the integral curves of A_{0} are the curves of constant ψ (cf. Hornig 2006). Since ψ solves the Poisson equation with source term B_{r}(r_{0},θ,φ), these curves are analogous to the surfaces of equal temperature in a solution to the heat equation where B_{r}(r_{0},θ,φ) 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_{0}. 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 = r_{0} shows B_{r} (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 . 

Open with DEXTER 
3. Magnetofrictional model
To approximate the evolution of nonpotential magnetic fields in the corona, we use the magnetofrictional 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 quasistatic sequence of approximately forcefree 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_{0}<r<r_{1},θ_{0}<θ<θ_{1},0 <φ< 2π } with a resolution 28 × 160 × 192. Here we take r_{0} = R_{⊙} (the photosphere), r_{1} = 2.5 R_{⊙}, θ_{0} = 0.05π, and θ_{1} = 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 magnetofrictional 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 smallscale coronal flows. For simplicity in this paper, we follow Mackay & van Ballegooijen (2006) and set (11)where η_{0} 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 higherorder 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 magnetofrictional method is to forego solving the full MHD equations and instead approximate the plasma velocity by (12)Here the first term is a frictionlike term that enforces relaxation towards a forcefree equilibrium. The factor  B  ^{2} prevents relaxation from being inhibited in weakfield 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 ν = ν_{0}r^{2}sin^{2}θ. 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. Zerogradient boundary conditions are imposed at r = r_{1}, and B_{θ} = 0 is imposed at θ = θ_{0},θ_{1}. At r = r_{0}, 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: ν_{0}, η_{0} and v_{out}. The friction coefficient ν_{0} controls the speed of coronal relaxation relative to the surface evolution, while η_{0} controls the rate of coronal diffusion. Rather than η_{0}, 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 magnetofrictional 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_{0}, 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 largescale 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 θ = θ_{0},θ_{1} and B_{θ} = B_{φ} = 0 on the outer boundary r = r_{1}. The potential field is computed using the eigenfunction method of van Ballegooijen et al. (2000). The coronal magnetic field is then evolved with magnetofriction, 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 η_{0}. Accordingly we will illustrate only the effect of varying ν_{0}, 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}  B_{r}  dΩ, panel b) shows ∫_{D}  j  dV, and panel c) shows H_{N} (asterisks) and H_{S} (circles). Panel d) shows the terms in Eq. (18)for the northern hemisphere, with asterisks denoting S_{0}, circles S_{1}, squares S_{eq}, and diamonds S_{V}. 

Open with DEXTER 
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 magnetofrictional relaxation. As we will discuss below, there is nonzero field line helicity associated with this dynamical equilibrium (shown by the colours in Fig. 3). The time taken to reach equilibrium depends on ν_{0}.
Indeed, it is instructive to consider the effect of ν_{0} (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 ν_{0}. 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 loglog fit shows that . 

Open with DEXTER 
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_{0} (r = r_{0}), S_{1} (r = r_{1}), and S_{eq} (θ = π/ 2), as in Fig. 6. There is no contribution from the latitudinal boundaries θ = θ_{0} and θ = θ_{1} 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_{0} of helicity through r = r_{0}, by differential rotation, and an output S_{1} 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^{2} day^{1} compared to 0.0032 Mx^{2} day^{1} for ν_{0} = 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 ν_{0} 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_{0} that is independent of ν_{0}. Figure 4d shows that the lack of stored helicity is compensated by the other surface terms S_{1} and S_{eq} during the relaxation phase. Even for the finite values of ν_{0} 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 S_{0}, S_{1}, and S_{eq} in Eq. (18). 

Open with DEXTER 
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 = r_{0} shows B_{r} (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 . 

Open with DEXTER 
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_{0} by solar rotation. Figure 6 shows that the outward helicity flux in the northern hemisphere is positive in our example, since S_{1} 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 27day 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_{1} = B_{2} = 100, the bipolar rings each contain the same unsigned magnetic flux 2πd^{2}, provided that they overlap neither each other nor the poles. We take d = 0.1, and locate them at θ_{1} = 0.3π and θ_{2} = 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_{0} of helicity injection by differential rotation remains fairly constant, albeit higher than before owing to the greater magnetic flux on r = r_{0}. (The slight decay in S_{0} 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 S_{1} during the flux rope eruption is not shown, and is much larger, about − 1.16 Mx^{2} day^{1}. 

Open with DEXTER 
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 ν_{0}, 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 ν_{0}, 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_{1} 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 crossequatorial 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 EastWest orientation of the polarity inversion line in our axisymmetric model. On the real Sun, polarity inversion lines at active latitudes are often aligned NorthSouth, 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_{1} 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_{0}, 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_{1}, 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 ν_{0} and η_{0}. We have already seen that the hemispheric helicity and open flux depend only weakly on ν_{0}. We have also run simulations with different but the same value of ν_{0} = 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. Nonaxisymmetric field
Fig. 9 ADAPT maps of B_{r} on r = r_{1}, used to generate the initial potential field extrapolations for periods A and B (white positive, black negative, saturated at ± 20 G). 

Open with DEXTER 
Our final example is a more realistic global magnetic configuration. For simplicity, we still consider only driving by largescale 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 fullsurface B_{r} maps modelling the real Sun on two dates: 2011Jan.01 and 2015Mar.10 (Fig. 9). The maps are taken from the Air Force DataAssimilative 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 magnetofrictional simulations use the same parameters as Sect. 5, driven by the same differential rotation and supergranular diffusion, with the parameters ν_{0} = 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 largescale 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_{0}, 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_{0} (see also Yeates & Mackay 2009b). In fact, the nonuniform distribution of magnetic flux across the solar surface means that some field lines have nonzero 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 = r_{0} shows B_{r} (white positive, black negative, saturated at ± 10 G), and projected coronal magnetic field lines traced from height r = r_{0} are coloured (red/blue) according to , saturated at ± 0.5 Mx. Animated versions of this figure for periods A and B are available online. 

Open with DEXTER 
Fig. 11 Various integrated quantities as a function of time, for the nonaxisymmetric simulations (periods A and B). Panel a) shows the total photospheric magnetic flux ∫_{r = r0}  B_{r}  dΩ, panel b) shows the total open flux ∫_{r = r1}  B_{r}  dΩ, panel c) shows ∫_{D}  j  dV, panel d) shows H_{N} (asterisks) and H_{S} (circles), and panel e) shows the rootmeansquare field line helicity . The vertical grey lines indicate times of strong flux rope ejections, as explained in the text. 

Open with DEXTER 
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 rootmeansquare 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 nonuniformly 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 180day 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_{1}. 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_{1}, 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_{0}. 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 rightmost column of Fig. 12. This shows the dimensionless quantity (21)which Liu et al. (2016) call the twist number. In a forcefree 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 = r_{1}. The middle column shows the distribution of on r = r_{0} (saturated at ± 2 Mx), and the right column shows the distribution of Tw at r = r_{0} (saturated at ± 40). The dashed lines in the left column show the neutral line where B_{r}(r_{1},θ,φ) = 0. Black circles in the other columns identify footpoints of field lines traced down from locations at r = r_{1} 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. 

Open with DEXTER 
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 nonlocal 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 nonlocal definition. One could therefore give a firstorder 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 nonpotential fields (cf. Edwards et al. 2015). However, a similar idea was proposed by Longcope & Malanushenko (2008), who defined the “additive selfhelicity” of a subdomain. 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 commonlyused 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 alreadytwisted structures from the solar interior, and we have also neglected the net injection of helicity by smallscale 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 magnetofrictional 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 solarterrestrial 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
 Aly, J.J. 2014, J. Phys. Conf. Ser., 544, 012003 [NASA ADS] [CrossRef] [Google Scholar]
 Amari, T.,Aly, J.J.,Canou, A., &Mikic, Z. 2013, A&A, 553, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Antiochos, S. K. 1987, ApJ, 312, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Antiochos, S. K. 2013, ApJ, 772, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Arge, C. N.,Henney, C. J.,Koller, J., et al. 2010, Twelfth International Solar Wind Conference, 1216, 343 [NASA ADS] [Google Scholar]
 Berger, M. A. 1984, Geophys. Astrophys. Fluid Dyn., 30, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A. 1988, A&A, 201, 355 [NASA ADS] [Google Scholar]
 Berger, M. A. 1993, Phys. Rev. Lett., 70, 705 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Berger, M. A., &Field, G. B. 1984, J. Fluid Mech., 147, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A., &Prior, C. 2006, J. Phys. A Mathematical General, 39, 8321 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, M. A., &Ruzmaikin, A. 2000, J. Geophys. Res., 105, 10481 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, P. F. 2011, Liv. Rev. Sol. Phys., 8 [Google Scholar]
 Cheung, M. C. M., &DeRosa, M. L. 2012, ApJ, 757, 147 [NASA ADS] [CrossRef] [Google Scholar]
 DeVore, C. R. 2000, ApJ, 539, 944 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, S. J.,Yeates, A. R.,Bocquet, F.X., &Mackay, D. H. 2015, Sol. Phys., 290, 2791 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, M. H., &He, Z. X. 1991, Ann. Math., 134, 189 [CrossRef] [Google Scholar]
 Haynes, A. L., &Parnell, C. E. 2010, Physics of Plasmas, 17, 092903 [NASA ADS] [CrossRef] [Google Scholar]
 Henney, C. J.,Toussaint, W. A.,White, S. M., &Arge, C. N. 2012, Space Weather, 10, S02011 [NASA ADS] [CrossRef] [Google Scholar]
 Hickmann, K. S.,Godinez, H. C.,Henney, C. J., &Arge, C. N. 2015, Sol. Phys., 290, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Hornig, G. 2006, ArXiv eprints [arXiv:astroph/0606694] [Google Scholar]
 Jiang, J.,Hathaway, D. H.,Cameron, R. H., et al. 2014, Space Sci. Rev., 186, 491 [NASA ADS] [CrossRef] [Google Scholar]
 Kazachenko, M. D.,Fisher, G. H., &Welsch, B. T. 2014, ApJ, 795, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Kimura, Y., &Okamoto, H. 1987, J. Phys. Soc. Japan, 56, 4203 [NASA ADS] [CrossRef] [Google Scholar]
 Knizhnik, K. J.,Antiochos, S. K., &DeVore, C. R. 2015, ApJ, 809, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, R.,Kliem, B.,Titov, V. S., et al. 2016, ApJ, 818, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Longcope, D. W., &Malanushenko, A. 2008, ApJ, 674, 1130 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D., &Yeates, A. 2012, Liv. Rev. Sol. Phys., 9, 6 [Google Scholar]
 Mackay, D. H., & van Ballegooijen, A. A. 2006, ApJ, 641, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D. H.,DeVore, C. R., &Antiochos, S. K. 2014, ApJ, 784, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Malanushenko, A.,Longcope, D. W.,Fan, Y., &Gibson, S. E. 2009, ApJ, 702, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, H. K. 1969, J. Fluid Mech., 35, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, H. K. 1990, Nature, 347, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, H. K., &Ricca, R. L. 1992, Proc. R. Soc. London Ser. A, 439, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Moraitis, K.,Tziotziou, K.,Georgoulis, M. K., &Archontis, V. 2014, Sol. Phys., 289, 4453 [NASA ADS] [CrossRef] [Google Scholar]
 Pariat, E.,Valori, G.,Démoulin, P., &Dalmasse, K. 2015, A&A, 580, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petrie, G. J. D. 2015, Liv. Rev. Sol. Phys., 12, 5 [CrossRef] [Google Scholar]
 Pevtsov, A. A., &Balasubramaniam, K. S. 2003, Adv. Space Res., 32, 1867 [NASA ADS] [CrossRef] [Google Scholar]
 Pontin, D. I.,WilmotSmith, A. L.,Hornig, G., &Galsgaard, K. 2011, A&A, 525, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prior, C., &Yeates, A. R. 2014, ApJ, 787, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Russell, A. J. B.,Yeates, A. R.,Hornig, G., &WilmotSmith, A. L. 2015, Physics of Plasmas, 22, 032106 [NASA ADS] [CrossRef] [Google Scholar]
 Sheeley, Jr., N. R. 2005, Liv. Rev. Sol. Phys., 2, 5 [Google Scholar]
 Snodgrass, H. B. 1983, ApJ, 270, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, J. B. 1974, Phys. Rev. Lett., 33, 1139 [NASA ADS] [CrossRef] [Google Scholar]
 Valori, G.,Démoulin, P., &Pariat, E. 2012, Sol. Phys., 278, 347 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A. A., &Martens, P. C. H. 1989, ApJ, 343, 971 [NASA ADS] [CrossRef] [Google Scholar]
 van Ballegooijen, A. A.,Priest, E. R., &Mackay, D. H. 2000, ApJ, 539, 983 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R. 2014, Sol. Phys., 289, 631 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., &Hornig, G. 2013, Physics of Plasmas, 20, 012102 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., &Hornig, G. 2014, J. Phys. Conf. Ser., 544, 012002 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., &Mackay, D. H. 2009a, ApJ, 699, 1024 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R., &Mackay, D. H. 2009b, Sol. Phys., 254, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R.,Mackay, D. H., & van Ballegooijen, A. A. 2008, Sol. Phys., 247, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Yeates, A. R.,Bianchi, F.,Welsch, B. T., &Bushby, P. J. 2014, A&A, 564, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yee, K. 1966, IEEE Transactions on Antennas and Propagation, 14, 302 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Existence of a surface with flux
Consider the magnetic field line L_{2} in Fig. 1, whose endpoints x_{1} and x_{2} both lie on the boundary r = r_{0}. We will show that there exists a curve γ from x_{1} to x_{2}, lying on the surface r = r_{0}, 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_{2} is a field line footpoint, we must have B_{r}(x_{2}) ≠ 0. If B_{r}(x_{2}) > 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_{2} 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_{1}, x_{2}, 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_{1} (a rarer situation in the corona). But what about an
open field line, where x_{1} lies on r = r_{0} and x_{2} on r = r_{1}? 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 idealinvariant 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 idealinvariant flux, even if the field line is open.
Fig. A.1 The original and deformed curves, all on the boundary r = r_{0}. 

Open with DEXTER 
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 (L_{1}) or open (L_{2}) magnetic field line. 

Open with DEXTER  
In the text 
Fig. 2 Physical interpretation of field line helicity in the DeVore gauge. The red circles show the direction of A_{0} (contours of ψ) arising from a strong magnetic source B_{r}> 0. The shaded surfaces are radial projections of the field lines L_{1} and L_{2}. Both field lines have a contribution to from any flux linking through these surfaces (owing to the second term of Eq. (4)), but only L_{1} has a contribution from A_{0}, since the projection of L_{2} is perpendicular to A_{0}. 

Open with DEXTER  
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 = r_{0} shows B_{r} (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 . 

Open with DEXTER  
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}  B_{r}  dΩ, panel b) shows ∫_{D}  j  dV, and panel c) shows H_{N} (asterisks) and H_{S} (circles). Panel d) shows the terms in Eq. (18)for the northern hemisphere, with asterisks denoting S_{0}, circles S_{1}, squares S_{eq}, and diamonds S_{V}. 

Open with DEXTER  
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 loglog fit shows that . 

Open with DEXTER  
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 S_{0}, S_{1}, and S_{eq} in Eq. (18). 

Open with DEXTER  
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 = r_{0} shows B_{r} (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 . 

Open with DEXTER  
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 S_{1} during the flux rope eruption is not shown, and is much larger, about − 1.16 Mx^{2} day^{1}. 

Open with DEXTER  
In the text 
Fig. 9 ADAPT maps of B_{r} on r = r_{1}, used to generate the initial potential field extrapolations for periods A and B (white positive, black negative, saturated at ± 20 G). 

Open with DEXTER  
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 = r_{0} shows B_{r} (white positive, black negative, saturated at ± 10 G), and projected coronal magnetic field lines traced from height r = r_{0} are coloured (red/blue) according to , saturated at ± 0.5 Mx. Animated versions of this figure for periods A and B are available online. 

Open with DEXTER  
In the text 
Fig. 11 Various integrated quantities as a function of time, for the nonaxisymmetric simulations (periods A and B). Panel a) shows the total photospheric magnetic flux ∫_{r = r0}  B_{r}  dΩ, panel b) shows the total open flux ∫_{r = r1}  B_{r}  dΩ, panel c) shows ∫_{D}  j  dV, panel d) shows H_{N} (asterisks) and H_{S} (circles), and panel e) shows the rootmeansquare field line helicity . The vertical grey lines indicate times of strong flux rope ejections, as explained in the text. 

Open with DEXTER  
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 = r_{1}. The middle column shows the distribution of on r = r_{0} (saturated at ± 2 Mx), and the right column shows the distribution of Tw at r = r_{0} (saturated at ± 40). The dashed lines in the left column show the neutral line where B_{r}(r_{1},θ,φ) = 0. Black circles in the other columns identify footpoints of field lines traced down from locations at r = r_{1} 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. 

Open with DEXTER  
In the text 
Fig. A.1 The original and deformed curves, all on the boundary r = r_{0}. 

Open with DEXTER  
In the text 