Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods
S
Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods
Gherardo Valori 0 1 2 3 4 5 6 7
Etienne Pariat 0 1 2 3 4 5 6 7
Sergey Anfinogentov 0 1 2 3 4 5 6 7
Feng Chen 0 1 2 3 4 5 6 7
Manolis K. Georgoulis 0 1 2 3 4 5 6 7
Yang Guo 0 1 2 3 4 5 6 7
Yang Liu 0 1 2 3 4 5 6 7
Kostas Moraitis 0 1 2 3 4 5 6 7
Julia K. Thalmann 0 1 2 3 4 5 6 7
Shangbin Yang 0 1 2 3 4 5 6 7
B G. Valori 0 1 2 3 4 5 6 7
0 Institute of SolarTerrestrial Physics SB RAS 664033 , P.O. Box 291, Lermontov st., 126a, Irkutsk , Russia
1 LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, UPMC Univ. Paris 06, Univ. Paris Diderot , Sorbonne Paris Cité, 92190 Meudon , France
2 Mullard Space Science Laboratory, University College London , Holmbury St. Mary, Dorking, Surrey, RH5 6NT , UK
3 Key Laboratory of Solar Activity, National Astronomical Observatories, Chinese Academy of Sciences , Beijing 100012 , China
4 Institute of Physics/IGAM, University of Graz , Universitätsplatz 5/II, 8010 Graz , Austria
5 Stanford University HEPL Annex , B210, Stanford, CA 943054085 , USA
6 School of Astronomy and Space Science, Nanjing University , Nanjing 210023 , China
7 Research Center for Astronomy and Applied Mathematics of the Academy of Athens , 4 Soranou Efesiou Street, 11527 Athens , Greece
Magnetic helicity is a conserved quantity of ideal magnetohydrodynamics characterized by an inverse turbulent cascade. Accordingly, it is often invoked as one of the basic physical quantities driving the generation and structuring of magnetic fields in a variety of astrophysical and laboratory plasmas. We provide here the first systematic comparison of six existing methods for the estimation of the helicity of magnetic fields known in a finite volume. All such methods are reviewed, benchmarked, and compared with each other, and specifically tested for accuracy and sensitivity to errors. To that purpose, we consider four groups of numerical tests, ranging from solutions of the threedimensional, forcefree equilibrium, to magnetohydrodynamical numerical simulations. Almost all methods are found to produce the same value of magnetic helicity within few percent in all tests. In the more solarrelevant and realistic of the tests employed here, the simulation of an eruptive flux rope, the spread in the computed values obtained by all but one method is only 3 %, indicating the reliability and mutual consistency of such methods in appropriate parameter
MaxPlankInstitut für Sonnensystemforschung; 37077 Göttingen; Germany

ranges. However, methods show differences in the sensitivity to numerical resolution and to
errors in the solenoidal property of the input fields. In addition to finite volume methods, we
also briefly discuss a method that estimates helicity from the field lines’ twist, and one that
exploits the field’s value at one boundary and a coronal minimal connectivity instead of a
predefined threedimensional magneticfield solution.
1 Introduction
H (t ) ≡
V
is the helicity of the vector field B = ∇ × A in a given volume V , with A(x, t )
representing the corresponding space and timedependent vector potential. If the field consists of a
discrete collection of flux tubes, H is the winding number (Moffatt 1969) expressing their
degree of mutual linkage. By extension, Eq. (1) is a measure of the entangling
(“knottedness”) of the field’s streamlines.
For a given vector potential A, the addition of the gradient of a (sufficiently regular but
otherwise arbitrary) scalar function, i.e., the transformation A −→ A + ∇ψ , does not change
the resulting B. This property of the definition of B is called gaugeinvariance. Due to this
freedom in the gauge, H is not uniquely defined, since
∂V
where dS = dS nˆ, with dS being the infinitesimal element of the bounding surface ∂ V of
the volume V , and nˆ the outwardoriented normal to ∂ V . Hence, H is not gaugeinvariant
unless two conditions are met: First, the vector field B must be solenoidal, as implied by
its definition as curl of A, and, second, the volume’s bounding surface ∂ V must be a flux
surface of B, i.e., (B · nˆ)∂V = 0. When applied to a magnetic field B, the solenoidal
requirement is satisfied by virtue of Maxwell equations, although possibly only to a finite extent
in numerical experiments, and ∂ V is a flux surface if no magnetic field line is threading the
boundary. This latter requirement is rarely satisfied in natural systems, and makes Eq. (1) of
little interest for practical use.
On the other hand, helicity has the fundamental property of being strictly conserved in
ideal magnetohydrodynamics (MHD, Woltjer 1958). Since MHD evolution in the absence
of dissipation preserves the topology of the magnetic field, the field lines’ winding
numbers, or, more generally, magnetic helicity, cannot be changed during the evolution. Even
more appealing is the fact that magnetic helicity, contrary to magnetic energy, is very well
conserved in nonideal dynamics as well (Berger 1984), as expected theoretically because
it cascades to large scales rather than to the small, dissipative scales (see, e.g., Frisch et al.
1975). Thanks to these properties, helicity has the possibility of being used as a constraint
for the magnetic field evolution. For isolated systems, conservation of helicity effectively
restricts the allowed timeevolution to helicityconserving paths in phasespace, which, for
instance, yielded the socalled Taylor hypothesis on magnetic field relaxation (see, e.g.,
Taylor 1986). In the solar context, helicity conservation is involved in magnetic field dynamos
(see, e.g., Brandenburg and Subramanian 2005), as well as a potential trigger of CMEs (see,
e.g., Rust 1994)
1.1 Relative Magnetic Helicity
In order to overcome the limitations of Eq. (1), Berger and Field (1984) and Finn and
Antonsen (1985) introduced the relative magnetic helicity
HV ≡
(A + Ap) · (B − Bp) dV,
of a magnetic field B with respect to a reference magnetic field Bp = ∇ × Ap. Even though
Eq. (3) allows for an arbitrary reference field, here we adopt the usual choice of the electric
currentfree (potential) field for Bp. This choice has the following motivation: in order for
HV in Eq. (3) to be gauge invariant, the input and potential fields, B and Bp, respectively,
must be solenoidal, and such that
(nˆ · B)∂V = (nˆ · Bp)∂V .
With such a prescription, the potential field that is chosen as a reference is uniquely defined
and represents the minimal energy state for a given distribution of the normal component of
the field on the boundaries (see, e.g., Valori et al. 2013). Moreover, ∂ V needs not to be a flux
surface for HV to be gaugeinvariant, and the definition Eq. (3) can be applied to arbitrary
field distributions.
Berger (2003) derived a useful decomposition of Eq. (3) as
HV = HV,J + HV,J P
HV,J ≡
(A − Ap) · (B − Bp) dV,
HV,J P ≡ 2
Ap · (B − Bp) dV.
The two definitions in Eqs. (6), (7) have the property of being separately gaugeinvariant
under the same assumptions guaranteeing the gaugeinvariance of HV . The first term HV,J
corresponds to the general definition of helicity equation (1), but this time of the
currentcarrying part of the field only, (B − Bp) = ∇ × (A − Ap). By construction, such a field has no
normal component on the boundary, i.e., has ∂ V as a flux surface. The second term HV,J P
has no intuitive interpretation, but is a sort of mutual helicity that basically takes care of the
flux threading ∂ V (via the transverse component of Ap), and is gaugeinvariant because only
the currentcarrying part of B appears.
The conservation properties of the relative magnetic helicity were numerically tested
by Pariat et al. (2015), confirming that helicity is a very well conserved quantity even in
presence of very strong dissipation. The equation regulating the relative magnetic helicity
variation rate due to dissipation and flux through the boundaries is also derived by Pariat
et al. (2015). In the particular case of the simulation of a jet eruption examined in that article,
relative helicity is conserved more than one order of magnitude better than free energy.
In the particular case of applications to the solar atmosphere, one additional
complication is that the magnetic field cannot be measured in the solar corona with the resolution
necessary for the computation of magnetic helicity. The magnetic field is instead inferred by
inverting spectropolarimetric measurements of emission from lower, higherdensity layers
of the atmosphere, yielding basically twodimensional maps of the field vector mostly at
photospheric heights, see e.g., Lites (2000). Therefore, in order to compute the helicity, it is
first necessary to introduce a model of the solar corona based on the observed photospheric
field values. In this work we exclude addressing this problem directly by using numerical
models of magnetic fields, in this way testing the different helicity methods in a strictly
controlled environment.
In summary, magnetic helicity is a fundamental quantity of plasma physics that is almost
exactly conserved in most conditions. This can be relevant in fusion plasmas, as well as in
astrophysical ones. In this article we focus on solar applications, but the conclusions derived
are general enough to be extended to other fields.
In the following we refer to the relative magnetic helicity equation (3) computed with
respect to the potential field defined by the boundary condition Eq. (4) simply as helicity,
and we consider V to be a single, rectangular, 3D, finite volume. More general formulations
are possible, e.g., Longcope and Malanushenko (2008) introduces a procedure for defining
reference fields on multiple subvolumes, possibly covering the entire open space. We also
refer to the discussion and references in Prior and Yeates (2014) for a physical interpretation
of H in Eq. (1) under several different gauges in the presence of open field lines threading
opposite faces of V .
1.2 Overview of the Methods for the Estimation of Helicity
Several methods of helicity estimation are currently available. A practical categorization,
according to decreasing levels of required input information, results into
– finite volume (FV)
– twistnumber (TN)
– helicityflux integration (FI)
– connectivitybased (CB)
methods. In practical applications, some assumption about the unknown coronal magnetic
field need to be made, and the above groups of methods essentially differ in the nature of
this assumption and in the correspondent helicity definition. A synoptic view of the available
methods for the estimation of magnetic helicity in finite volumes is presented in Table 1.
Finite volume (FV) methods entirely rely on external techniques, such as nonlinear
forcefree field extrapolations or MHD simulations, to produce numerical models of the coronal
magnetic field, (see, e.g., Wiegelmann et al. 2015). The “finite volume” characterization
indicates that the methods are designed to provide the helicity value in a bounded volume,
typically one employed in a 3D numerical simulation, as opposed to methods that estimate
the helicity in a semiinfinite domain. The helicity in a given volume at a given time can be
directly computed if the magnetic field is known at each location in V at that time.
Therefore, FV methods are direct implementation of Eq. (3) which requires only the computation
of the vector potentials for a given discretized magnetic field B in V , see e.g., Thalmann
et al. (2011), Valori et al. (2012), Yang et al. (2013b), Amari et al. (2013), Rudenko and
Anfinogentov (2014), Moraitis et al. (2014). Despite the apparent straightforward task that
such methods have, differences in the gauge, in the implementations, and in the sensitivity
to the input discretized magnetic field may impact on the accuracy of the helicity
estimation. To test the accuracy of finite volume methods is the main focus of this article, and such
methods are extensively described in Sect. 2.
Table 1 Synoptic view of helicity computation methods, their properties and formulation, as described in
Sect. 1.2. The subset of methods actually tested in this paper is listed in Table 2
d HdtV = 2 ∂V [(Ap · B)vn − (Ap · vt)Bn] dS
– Requires B in V e.g., from MHD simulations or – Requires time evolution of vector field on ∂V
NLFFF – Requires knowledge or model of flows on ∂V
–– CMoamy peumteplHoyVdaiftfeorneenttigmaeuges (see Table 2) – sVeaeliPdafroiartaestpael.ci(fi2c01s7e)t of gauge and assumptions,
H T Φ2
H = A
– Estimation of the twist contribution to H – Requires the vector field on photosphere at one
– Requires B in V time
– Requires a fluxropelike structure for computing – Models the corona connectivity as a collection of
the twist T M forcefree flux tubes
– Minimal connection length principle
The fieldline helicity method by Russell et al. (2015) is also using the full 3D
vector magnetic field in a volume as an input to the method. Rather than producing a single
number for the value of Eq. (3) in V , the fieldline helicity method provides the value of
helicity associated to a single flux tube, and follows its evolution in time. In this sense, the
FL method is a powerful investigating tool for studying the distribution of helicity in
numerical simulations, especially those involving reconnection processes. Given its peculiar
and focussed applications, we do not discuss this method further, and refer the reader to the
abovementioned article.
Similarly to FV methods, the twistnumber (TN) method (see, e.g., Guo et al. 2010)
requires in input the 3D discretized magnetic field vector. The method also assumes the
presence of a flux rope in the coronal volume, and proceeds by relating the twist of that structure
with helicity. The level of approximation of the true helicity value that is implied by such a
technique is assessed here for the first time. Application of this method to observations can
be found in Guo et al. (2013).
Helicityflux integration (FI) methods, do not make any assumption about the coronal
field, but rather assume that the helicity accumulated in a given volume is the result of the
helicity flux through the volume boundaries, from a given point in time onward, see, e.g.,
Eq. (5) in Berger and Field (1984) for negligible dissipation. Such an estimation requires the
knowledge of the time evolution of magnetic and velocity fields on the bounding surface of
the considered volume (see e.g., Eq. (16) in Berger 1999. The vector potential of the potential
field appearing there can be derived from the field distribution on the boundaries). Under
these assumptions, in the case of negligible dissipation, no information on the magnetic
field inside the volume is necessary.
In practical applications, such methods follow the time evolution of the photospheric field
and assume that the flux of helicity through that boundary accumulates in the coronal field,
see e.g., Chae et al. (2001). Since only the flux is computed, FI methods can only estimate
the variation of accumulated helicity with respect to an unknown initial state. Methods that
exploit this approach appear in Nindos et al. (2003), Pariat et al. (2005), LaBonte et al.
(2007), Liu and Schuck (2012), among others, but direct comparisons between them do not
yet exist.
A method of computing the helicity that also uses only the distribution of magnetic field
on the bottom boundary is the connectivitybased (CB) method (Georgoulis et al. 2012).
The method is based on modeling the unknown connectivity of the coronal field with a
collection of slender forcefree flux tubes, each with different constant forcefree parameter
α ≡ (∇ × B)/B. More specifically, the CB method takes in input the photospheric
observations and models the coronal field as a single (linear) or a collection of (nonlinear) flux
tube(s) as an integral step of the helicity computation itself. The set of flux tubes is
obtained assuming that the line connectivity is, globally, the shortest possible, mimicking the
compact character of the more flareproductive active regions. In this way, the
connectivitybased method requires as input only the knowledge of the magnetic field distribution at
the photosphere, at each time. Different flux tubes have a different value of the forcefree
parameter, hence the characterization of the method as “nonlinear”, despite the
simplification of neglecting the braiding between different flux tubes that is used in summing up
the helicity and energy contributions of individual flux tubes. Therefore, the CB method is
an approximate, nonlinear method that is meant to produce a lowerlimit estimation of the
true helicity associated to a fluxbalanced coronal field in a very fast way. In this sense, the
CB method does not share the same purpose of finite volume and helicityflux integration
methods, which, in the ideal situation, are in principle capable of obtaining the true value of
helicity in a volume, at the price of requiring more information in input.
Both the TN and the CB methods are based on the representation of the magnetic field
as a collection of a discrete number of finitesized flux tubes, as opposite to a continuous
threedimensional field. We therefore categorized both methods as discrete fluxtubes (DT)
methods, see Table 1.
1.3 Systematic Comparison of Methods
Section 1.2 briefly introduced the methods currently available for the estimation of HV .
Many of those have been independently tested and already applied to observations.
However, the accuracy, mutual consistency, and sensitivity of these methods are not sufficiently
tested, while a systematic comparison of methods in both their theoretical as well as practical
aspects is necessary. This work is the first one of a series of papers undertaking such a task.
Beside benchmarking, the ultimate goal of this series is no less than to assess if and how can
helicity be meaningfully used as a tracer of the evolution of the magnetic field in magnetized
plasmas. To that purpose, we designed a collection of progressively more complex and
realistic discretized test fields, starting with 3D solutions of the forcefree equations, proceeding
then to timedependent MHD simulations of flux emergence, and finally to applications to
real solar observations.
The present article is focussed on FV methods, whereas in Pariat et al. (2017) we mostly
focus on FI methods and the comparison with the CB method. Subsequent papers are
dedicated to the TN method (Guo et al. 2017), and to applications to observed solar active
regions. The results of this and of subsequent articles presented in this section are the direct
outcome of the International Team on magnetic helicity hosted at ISSIBern in the 2014–
2016 period.1 As a reference for future testing, we make available the data cube of each
1See the web page http://www.issibern.ch/teams/magnetichelicity/ for more information.
employed test field and the complete results of the analysis in tabular form, of which the
plots presented here are a subset. That material can be found in the section Publications of
the ISSI website given in the footnote.
1.4 This Article
In view of the large scope of the project outlined in Sect. 1.3, it is necessary in the first
place to determine the respective limits, field of applications, and precision of each of the
existing FV methods, and check if FV methods can be reliably used to benchmark other,
more approximate methods.
More in detail, we wish to properly quantify the reliability of HV estimations when the
field is known in the volume V . Such a reliability can be tested by quantifying the sensitivity
and robustness of the estimations with respect to resolution, energy and helicity properties
of the input field, and sensitivity to violation of the solenoidal constraint by the discretized
field.
In order to compare and benchmark existing methods against the above properties in a
representative variety of relevant setting, we consider test cases that confront the methods
with increasing complexity, uncertainties, and realism. We consider
stronglycontrolledenvironment, equilibrium testcases such as the Low and Lou (1990) and Titov and
Démoulin (1999) solutions of the forcefree equations. Such tests provide basic benchmarking
as they differ for helicity content, resolution, and accuracy of the solenoidal property. Then,
two series of snapshots from MHD numerical simulations of flux emergence resulting in
stable (Leake et al. 2013) and unstable (Leake et al. 2014) configurations are considered.
The flux emergence test cases are also meant to build a bridge toward FI methods studied in
Pariat et al. (2017), which is more focussed in understanding how the helicity information is
modified when the knowledge of the magnetic field is limited (typically to the photosphere
only).
In addition to FV methods, in this article we also consider two DT methods, namely the
twistnumber (TN) and the connectivitybased (CB) methods. These methods were already
used to obtain estimates of the magnetic helicity in several observational studies (Guo et al.
2013; Georgoulis and LaBonte 2007; Georgoulis et al. 2012; Tziotziou et al. 2012, 2013,
2014; Moraitis et al. 2014). Benchmarking DT methods together with the FV ones enables
the reader to better interpret results of past and future studies applying such methods. From
a different point of view, the TN method is included because it requires the same
information as the other finite volume methods, i.e., the full knowledge of the magnetic field in the
entire considered volume. The CB method is included because, despite requiring only the
photospheric vector magnetogram, it can use any available information of the coronal
connectivity. Also, similarly to FV methods, the connectivitybased method can be applied to a
single time snapshot, rather than requiring a series of temporal snapshots, which is the case
for the FI methods. A list of the methods tested in this article is given in Table 2.
The methods applied in this article are described in Sect. 2, whereas the numerical
magnetic fields used as tests are discussed in Sect. 4. Section 3 summarizes the diagnostic tools
used for the comparing the results. The main results of the comparison are then discussed in
Sect. 5, with specific discussion of the dependence on resolution presented in Sect. 6 and of
the dependence on the solenoidal property given in Sect. 7. An overview of the FV methods
results is given in Sect. 8, whereas the comparative tests with DT methods are presented in
Sect. 9. Finally our conclusions are summarized in Sect. 10.
Table 2 Summary of the methods employed in this article, their short label, their categorization according
to Table 1, the section of this article where they are described, and their main bibliographic reference
CoulombThalmann Coulomb_JT Finite volume Sect. 2.1.1 Thalmann et al. (2011)
CoulombYang Coulomb_SY Finite volume Sect. 2.1.2 Yang et al. (2013b)
CoulombRudenko Coulomb_GR Finite volume Sect. 2.1.3 Rudenko and Anfinogentov (2014)
DeVoreValori DeVore_GV Finite volume Sect. 2.2.1 Valori et al. (2012)
DeVoreMoraitis DeVore_KM Finite volume Sect. 2.2.2 Moraitis et al. (2014)
DeVoreAnfinogentov DeVore_SA Finite volume Sect. 2.2.3 Not available
Twistnumber TN Discrete fluxtubes Sect. 2.3.1 Guo et al. (2010)
Connectivitybased CB Discrete fluxtubes Sect. 2.3.2 Georgoulis et al. (2012)
2 Helicity Computation Methods
Finite volume methods require the knowledge of the magnetic field B in the entire volume V ,
and differ essentially in the way in which the vector potentials are computed. The methods
presented here compute vector potentials employing either the Coulomb gauge (∇ · A = 0)
or the DeVore gauge (Az = 0, DeVore 2000). Due to the gaugeinvariant property of Eq. (3),
the employed gauge should be irrelevant for the helicity value. It may have, however,
consequences on the number and type of equations to be solved for that purpose. Methods using
the Coulomb gauge differ in the way in which the magnetic fields and the corresponding
vector potentials are split into potential and currentcarrying parts. Hence, they differ to
some extent in the equations that they solve. Methods applying the DeVore gauge are
applications of the method in Valori et al. (2012) that differ only in the details of the numerical
implementation. In the following, different FV methods are identified by the gauge they
employ (DeVore or Coulomb), followed by the initial of the author of the reference article
describing its implementation (e.g., Coulomb_GR labels the Coulomb method described in
Rudenko and Myshyakov 2011), see Table 2.
All the FV methods considered in this article, except for the Coulomb_GR method, define
the reference potential field as Bp = ∇φ, with φ being the scalar potential, solution of
such that the constraint equation (4) is satisfied. Errors in solving equations (8), (9) are a
first source of inaccuracy for the methods.
2.1 Methods Employing the Coulomb Gauge
Vector potentials in the Coulomb gauge satisfy
∇2Ap = 0,
∇ · Ap = 0,
nˆ · (∇ × Ap)∂V = (nˆ · B)∂V ,
for the vector potential of the potential field, and
∇2A = −J,
∇ · A = 0,
nˆ · (∇ × A)∂V = (nˆ · B)∂V ,
for the vector potential of the input field, where J = ∇ × B. The conditions Eqs. (12) and (15)
are the translation into vector potential representation of Eq. (4). The accuracy of Coulomb
methods depend on the accuracy in solving numerically the above Laplace and Poisson
problems. This includes the accuracy in fulfilling numerically the gauge condition, i.e., the
solenoidal property of the vector potentials Ap and A.
From the computational point of view, the numerical effort required to solve for the vector
potentials consists, in general, in the solutions of Eqs. (10)–(12) and (13)–(15), i.e., of six
3D Poisson/Laplace problems, one for each Cartesian component of the vector potentials
Ap and A.
2.1.1 Coulomb_JT
In order to solve Eqs. (10)–(12) and (13)–(15), appropriate additional boundary conditions
for A and Ap on the boundaries of the 3Drectangular computational domain need to be
specified. For this purpose, the method of Thalmann et al. (2011), decomposes A into a
currentcarrying and a potential (currentfree) part, in the form A = Ac + Ap. The
reproduction of the input fields’ flux at the volume’s boundaries, ∂ V , is entirely dedicated to Ap
(obeying Eqs. (10)–(12)). The electric current distribution in V and on ∂ V , on the other
hand, are delivered by Ac (Eqs. (13), (14) with Ac instead of A and Eq. (15) replaced by
nˆ · (∇ × Ac)∂V = 0).
In particular, the tangential components of Ap on a particular face, f , of the 3D
computational domain are specified to be the 2D stream function of a corresponding Laplacian
field, φf , in the form Ap,t = −nˆ × ∇t φf , where ∇t is the 2Dgradient tangential operator
on the face f . The Laplacian field itself is gained by substituting in Eq. (12) and seeking the
solution of the derived 2D Laplace problem ∇t2φf = −nˆ · B, for which boundary conditions
on the four edges of each face f need to be specified. This approach of defining Ap on ∂ V
is in principle used by all Coulomb methods considered in the present study, but the specific
way in which the 2D Laplace problems are formulated is different.
Thalmann et al. (2011) use Neumann conditions in the form ∂nφf = cf , where cf is a
constant along a particular face and ∂n is the derivative in the direction normal to the edge
of the face f (see their Sect. 2.1 for details). The different cf are constructed in such a way
that the total outflow through the volume’s bounding surface ∂ V is minimized. In this way,
a vanishing tangential divergence (∇t · Ap,t = 0) is enforced on ∂ V , and following Gauss’
theorem, Eqs. (10)–(12) are approximately fulfilled.
Another difference of the applied Coulomb methods is how the currentcarrying
vector potential Ac is calculated and its solenoidality enforced. Thalmann et al. (2011) solves
Eq. (13) for Ac numerically, similar to Yang et al. (2013a), just with differing boundary
conditions. In the Coulomb_JT case, ∇t · Ac,t = 0 on ∂ V is explicitly enforced in order to fulfill
the Coulomb gauge for Ac.
The method discussed in Thalmann et al. (2011) is implemented in C. The Poisson and
Laplace problems are solved numerically using the Helmholtz solver in Cartesian
coordi
nates of the Intel® Mathematical Kernel Library.
2.1.2 Coulomb_SY
The Coulomb_SY method is described in Yang et al. (2013a,b). In the original
formulation of Yang et al. (2013a), the method requires a balanced magnetic flux through each of
the side boundaries of the volume. This restriction has been further removed in Yang et al.
(2013b). In order to solve Eqs. (10), (12) and (13), (15), the Coulomb_SY method
additionally enforces the boundary condition (nˆ · Ap)∂V = (nˆ · A)∂V = 0 at all boundaries. Then,
the transverse vector potential at the boundaries and the vector potential at the edges is
obtained by using Gauss’ theorem. After obtaining the boundary values, Yang et al. (2013a,b)
firstly resolve the Laplace equation (10) and the Poisson equation (13) to obtain an initial
guess of the solution, Ap and A . These preliminary solutions satisfies Eqs. (12) and (15),
but not the Coulomb gauge condition. The Coulomb_SY method then uses a
divergencecleaning technique based on the Helmholtz vector decomposition to iteratively impose the
Coulomb constraint to the vector potentials, without modifying their values at the
boundaries. Comparing with the Coulomb_JT method, in Coulomb_SY are the vector potentials
that are decomposed, rather than the boundary contributions. The method is implemented in
Fortran; Poisson and Laplace problems are solved numerically using the Helmholtz solver
in Cartesian coordinates of the IMSL® (International Mathematics and Statistics Library).
2.1.3 Coulomb_GR
The Coulomb_GR method is described in Rudenko and Myshyakov (2011). A distinctive
feature of the algorithm is that the Coulomb_GR method defines the reference potential
field in terms of vector potential Bp = ∇ × Ap, rather than using Eqs. (8), (9). The
corresponding boundary value problem, Eq. (10), is solved with the constraint equation (11)
and the boundary condition (Ap · nˆ)∂V = 0. The Laplace problem is divided into six
subproblems, one for each side of V . Such a splitting of the Laplace problem is correct only if
the total magnetic flux is zero (balanced) on each side of the box independently. To satisfy
this requirement, a compensation field Bm
p = ∇ × Apm is introduced. It is built as a field of
5 magnetic monopoles located outside of the box. Positions and charge of the monopoles are
selected such as to compensate unbalanced flux on each side of the volume independently.
The modified magnetic field B = B − Bpm has zero total flux on each face independently and
can be correctly used as a boundary condition for the subproblems
where Apfi is the vector potential of subproblem solution corresponding to the side fi .
After solving all subproblems, the full solution is then obtained by summation of the
solutions of the six subproblems Afi and the vector potential of a compensation field, Am,
as
i=1
The field described by the first term of Eq. (17), instead, is flux balanced on each side of the
box independently.
Instead of solving numerically the Poisson problem equations (13)–(15), the Coulomb_
GR methods adopts a decomposition similar to the one in Coulomb_JT method, i.e., A =
Ac + Ap, but the vector potential of the currentcarrying part of the field is computed as
In contrast to other methods, the solutions of the Laplace and Poison problems for the Afi
components are derived analytically as decompositions into a set of orthonormal basis
functions. The detailed description of the strategy for solving these equations can be found in
the original paper by Rudenko and Myshyakov (2011).
In the current implementation the method is relatively demanding in terms of running
time. Therefore, it is applied here only to a subset of test cases.
2.2 Methods Employing the DeVore Gauge
Using DeVore gauge Az = 0 (DeVore 2000), Valori et al. (2012) derived the expression for
the vector potential of the magnetic field B in the finite volume V = [x1, x2] × [y1, y2] ×
[z1, z2] as
where the integration function b(x, y) = A(z = z2) obeys to
A = b + zˆ ×
∂x by − ∂y bx = Bz(z = z2),
and bz = 0. The particular solution of Eq. (20) employed here is
but see Valori et al. (2012) for alternative options. The above equations are applied in the
computation of the vector potential of the potential field too by substituting B with Bp
everywhere in Eqs. (19)–(22). In particular, using Eqs. (21), (22) for both Ap and A implies
Ap = A at z = z2, although this is not necessarily required by the method.
The DeVore gauge can be exactly imposed also in numerical applications, which is
generally not the case for the Coulomb gauge. On the other hand, since Az = 0, then
Bx = −∂zAy = ∂z
where Eqs. (20) and (21), (22) where used; a similar expression holds for By. Hence, the
accuracy of DeVore method in reproducing Bx and By from A of Eq. (19) depends only on
how accurately the relation
is verified numerically. On the other hand, even when Eq. (24) is obeyed to acceptable
accuracy, one can easily show that, for a nonperfectly solenoidal field Bns, it is
∂z
= identity
Bns − ∇ × A = zˆ
(∇ · Bns) dz ,
as derived in Eq. (B.4) of Valori et al. (2012). Hence, the accuracy in the reproduction of
the zcomponent of the field depends on the solenoidal level of the input field (and on how
accurately Eq. (20) is solved).
All DeVoregauge methods discussed in this study employs Eqs. (8), (9) and (19)–(22),
but they differ in the way integrals are defined, and in the way the solution to Eq. (8) is
implemented. The computationally most demanding part of the method is the solution of the
3D scalar Laplace equation for the computation of the potential field, Eq. (8). This makes
DeVore methods computationally appealing since they require very little computation time.
2.2.1 DeVore_GV
DeVore_GV is the original implementation described in Valori et al. (2012), where the
requirement equation (24) is enforced by defining the zintegral operator as the numerical
inverse operation of the second order central differences operator, see Sect. 4.2 in Valori
et al. (2012). The Poisson problem for the determination of the scalar potential φ in Eqs. (8),
(9) is solved numerically using the Helmholtz solver in the proprietary Intel® Mathematical
Kernel Library (MKL).
Following Eq. (39) in Valori et al. (2012), the DeVore gauge for the potential field can be
reduced to the Coulomb gauge. We checked the effect of this gauge choice in the tests below,
and found no significant difference with the standard DeVore gauge. The DeVoreCoulomb
gauge for the potential field is thus no further discussed here.
2.2.2 DeVore_KM
DeVore_KM is described in Moraitis et al. (2014). This implementation has two main
differences with the one of DeVore_GV. The first is in the solver of Laplace’s equation.
DeVore_KM uses the routine HW3CRT that is included in the freely available FISHPACK
library (Swartztrauber and Sweet 1979). A test, however, with the corresponding Intel MKL
solver revealed minor differences in the solutions obtained with the two routines, and a
factor of ≤ 2 more computational time required by the FISHPACK solver. The second and
most important difference with the DeVore_GV method is in the numerical calculation of
integrals and derivatives in Eqs. (19)–(22). In DeVore_KM integrations are made with the
modified Simpson’s rule of error estimate 1/N 4 (Press et al. 1992), with N being the
number of integration points, and, in the special case N = 2, with the trapezoidal rule instead.
In addition, differentiations are made using the appropriate (centered, forward or backward)
secondorder numerical derivative, without trying to numerically realize Eq. (24). Finally,
Eqs. (19)–(22) in DeVore_KM are used in the same way for both the potential and the
reference fields.
2.2.3 DeVore_SA
DeVore_SA follows the general scheme of the DeVore_GV method with two differences.
The first one is that Eqs. (8), (9) for the potential φ are solved in Fourier space separately
for all faces of the box. In particular, the problem is divided into six subproblems using
i=1
where φc is the 3D scalar potential of the compensation field Bc = ∇φc and φi are 3D
solutions for the potential field with the normal component given on ist side of V and vanishing
boundary conditions on the other sides of V . The individual Laplace problems for each φi
are then solved in Fourier space following the general scheme of the potential and linear
forcefree field extrapolation employing the fast Fourier transform by Alissandrakis (1981).
For the application here, the original extrapolation algorithm is modified to take into
account the imposed boundary conditions. This method of solving equations (8), (9) will be
described in a dedicated forthcoming paper.
The second difference with the DeVore_GV method is that Eq. (19) is modified by
introducing a new integration function c that is computed using Bz from any level zr inside the
data cube. In particular, by addition and subtraction to Eq. (19), one has
A = b + zˆ ×
B dz −
which can be recasted as
where we have defined
A = c + zˆ ×
B dz −
Berger and Field (1984) and Démoulin et al. (2006) have shown that the relative magnetic
helicity can be approximated as the summation of the helicity of M flux tubes:
i=1
i=1 j=1,j=i
where Ti denotes the twist and writhe of magnetic flux tube i with flux Φi , and Li,j is the
linking number between two magnetic flux tubes i and j with fluxes Φi and Φj , respectively.
The first and second term on the right hand side of Eq. (31) represents the self and mutual
helicity, respectively. With the approximation of discrete magnetic flux tubes, the physical
quantity of the magnetic helicity is related to the topological concept of the writhe, twist and
linking number of curves, and the magnetic flux associated with those curves. The formulae
of computing these topological quantities for both close and open curves have been derived
in Berger and Prior (2006) and Démoulin et al. (2006). For the purpose of our
comparison it must be noticed that discrete fluxtubes methods do not provide the vector potentials
and potential fields in the considered volumes like FV methods. Therefore, the comparison
c = b + zˆ ×
∂x cy − ∂y cx = Bz(z = zr ).
The solution of Eq. (30) is then analogous to Eqs. (21), (22), where Bz(z = z2) is replaced
by Bz(z = zr ). Tests using the LL case of Sect. 4.1 shown that the minimal error in Ap and A
is obtained for zr = (z2 − z1)/2. The vector potential is finally computed following Eq. (28).
This scheme coincides with the original one of DeVore_GV if zr is taken at the top
boundary of the box, i.e., for c(zr = z2) = b.
2.3 Discrete FluxTubes Methods
between DT methods and FV methods is necessarily restricted to the helicity values only.
The twistnumber method and connectivitybased method presented in this section adopt
different assumptions in the helicity formulae and the magnetic field models to compute the
magnetic helicity.
2.3.1 TwistNumber
The TN method is described in Guo et al. (2010, 2013). This method is aimed at computing
the helicity of a highly twisted magnetic structure, such as a magnetic flux rope. A magnetic
flux rope is considered as an isolated, single flux tube such that only the self magnetic
helicity is computed. The helicity contributed by the writhe is also omitted assuming that
the flux rope is not highly kinked. With these two assumptions, the magnetic helicity of a
single highly twisted structure is simplified as
where T is the twist number of the considered magnetic flux rope with flux Φ . In order to
estimate T , the formula derived in Berger and Prior (2006) to compute the twist number of
a sample curve referred to an axis is employed. Practically, the axis can be determined by
the symmetry of a magnetic configuration or by other assumptions, such as requiring it to
be horizontal and to follow the polarity inversion line (Guo et al. 2010). The boundary of
the flux rope is determined by the quasiseparatrix layer (QSL) that is found to wrap the flux
rope (Guo et al. 2013). Then the twist density, dT /ds, at an arc length s is:
Two unit vectors are used in Eq. (33): T(s), that is tangent to the axis curve, and V(s), that
is normal to T and pointing from the axis curve to the sample curve. By integrating the
equation along the axis curve the total twist number is derived. Equation (33) is suitable
for smooth curves in arbitrary geometries without self intersection. Since it makes no
assumption about the magnetic field, it can be applied to both forcefree and nonforcefree
magnetic field models.
2.3.2 ConnectivityBased
The CB method was introduced by Georgoulis et al. (2012) and was used by a number of
studies thereafter. In principle, the method requires only the full (vector magnetic field)
photospheric boundary condition to selfconsistently estimate a lower limit of the free energy
and the corresponding relative helicity.
A key element of the method is the discretization of a given, continuous photospheric
flux distribution into a set of partitions with known spatial extent and flux content. Each
partition is then treated as the collective footprint of one or more flux tubes. To map the
relative locations of these footprints, one either infers or calculates the coronal magnetic
connectivity that distributes the partitioned magnetic flux into oppositepolarity connections,
treated thereafter as discrete magnetic flux tubes. The flux content of these connections, with
both ends within the photospheric field of view (FOV), constitutes the magnetic connectivity
matrix corresponding to the given photospheric boundary condition.
The unknown coronal connectivity is either inferred by any explicit solution of the
volume magnetic field or calculated with respect to the existing photospheric boundary
condition. In the first case, individual fieldline tracing associates connected flux with
photospheric partitions, providing the magnetic connectivity matrix upon summation of individual
fieldline contributions. Obviously, only magnetic field lines entirely embedded in the finite
volume are taken into account. In the second case, a simulatedannealing method is used to
absolutely and simultaneously minimize the flux imbalance (hence achieving connections
between oppositepolarity partitions) and the (photospheric) connection length. This
criterion is designed to emphasize photospheric magnetic polarity inversion lines by assigning
higher priority to connections alongside them. The converged simulatedannealing solution,
that provides the connectivity matrix, is unique for a given photospheric partition map. More
information and examples are provided in Georgoulis et al. (2012) and Tziotziou et al. (2012,
2013).
The connectivity matrix in a collection of partitions of both polarities will reveal a
number of M discrete, assumed slender, flux tubes with flux contents Φi ; i ≡ {1, . . . , M }. The
respective forcefree parameters αi are assumed constant for a given flux tube but vary
between different tubes, thus implementing the nonlinear forcefree (NLFF) field
approximation. Forcefree parameters for each flux tube are the mean values of the forcefree
4π Ii ;
parameters of the tubes’ respective footprints, each calculated by the relation αi = c Fi
i ≡ {1, . . . , M } for M magnetic partitions, where Ii is the total electric current of the
ipartition and Fi its flux content. The total current is calculated by applying the integral form
of Ampére’s law along the outlining contour of the partition.
Knowing Fi , αi , and the relative positions of each flux tube’s footpoints, Georgoulis et al.
(2012) showed that a lower limit of the free magnetic energy for a collection of M flux tubes
is
EcCB ≡ Ec(CB;self ) + Ec(CB;mutual) = Aλ2
i=1
l=1 m=1;l=m
where A, δ are known fitting constants, λ is the length element (the pixel size in observed
photospheric magnetograms), and Llm is the mutualhelicity parameter for a pair (l, m) of
flux tubes. This parameter is inferred geometrically, by means of trigonometric interior
angles for the relative positions of the two pairs of fluxtube footpoints. The locations of
pointlike footpoints of the slender flux tubes coincide with the fluxweighted centroids of the
respective flux partitions. As included in Eq. (34), the parameter Llm does not include
braiding between the two flux tubes, that can be found only by the explicit knowledge of the
coronal connectivity. Additional complexity via braiding will only add to the free energy
EcCB in Eq. (34). Therefore, the above EcCB is already a lower limit of the actual Ec,
assuming only “archlike” (i.e., one above or below the other) flux tubes that do not intertwine
around each other. In addition, Eq. (34) does not include an unknown freeenergy term that
is due to the generation, caused by induction, of potential flux tubes around the collection of
nonpotential ones (Démoulin et al. 2006). Such a term would again contribute to the mutual
term of the free energy.
The corresponding selfconsistent relative helicity is, then,
HmCB ≡ Hm(CB;self ) + Hm(CB;mutual) = 8π Aλ2
From Eqs. (34), (35) we identically have EcCB ≡ 0 for potential flux tubes (αi = 0; i ≡
{1, . . . , M }). For HmCB = 0 in this case, we further require lM=p1 mM=p1;l=m LlmP Φl Φm = 0,
where LlmP = Llm is the mutualhelicity factor for a collection of Mp = M collection of
potential flux tubes. As Démoulin et al. (2006) discuss, this can be the case for a
fluxbalanced potentialfield boundary condition. In practical situations of notprecisely
fluxbalanced magnetic configurations, however, one may approximate HmCB;mutual = 0, in case all
αi ; i ≡ {1, . . . , M } are zero within uncertainties δαi , which are fully defined in this analysis.
More generally, one may use the “energyhelicity diagram” correlation of Tziotziou et al.
(2012, 2014) to infer HmCB  ∝ E0.84±0.05 for EcCB −→ 0.
cCB
3 Analysis Metrics
Apart from extremely simplistic magnetic fields, the analytical computation of the relative
magnetic helicity in a nonmagnetically bounded system is highly nontrivial. Even with
simple naturalworldrelevant models relative magnetic helicity cannot be analytically
estimated. Similarly, the exact value of HV in the finite volume of the 3D discretized magnetic
fields used here as tests is, in general, not known. Hence, we need to provide indirect
accuracy metrics to judge the examined methods.
The main goal of the analysis presented below is to compare the helicity values that are
obtained employing the potential field and vector potentials computed with the methods
described in Sects. 2.1 and 2.2. Since the helicity of B defined by Eq. (3) involves the
corresponding potential field Bp, as well as the vector potentials for Bp and B, this basically
implies providing a quantitative estimation of the accuracy of such fields. To that purpose,
we introduce normalized quantities and metrics as follows:
For each discretized magnetic field, we define HV as the helicity defined in Eq. (3)
normalized to Φ2,
is half of the unsigned flux through the bottom boundary, corresponding to the injected flux
for an exactly fluxbalanced configuration. In that normalization, a uniformly twisted flux
rope with field lines having N turns has an helicity equal to N (see e.g., Démoulin and
Pariat 2009). In computing the helicity values with different methods we refrain from using
simplifications of Eq. (3) coming from the specific gauge in use, in this way keeping the
comparison as general as possible. Hence, for each FV method and for each test case, the
value of HV as defined by Eq. (36) is obtained by computing separately the four volume
contributions of Eq. (3) and normalizing them to Φ2. We reserve the calligraphic symbol H
for nonnormalized helicities.
The numerically obtained HV values depend in principle on many factors. In the first
place, different methods may have a different level of accuracy in computing the vector
potentials of the test and potential fields, depending on the strategy applied to solve the
relevant equations, see Sect. 2. Second, the reference potential field is uniquely defined by
the requirement that HV is gauge invariant, yielding to Eq. (4). However, without
violating that requirement, the potential field can equivalently be computed both as the curl of
the vector potential, as in some methods of Sect. 2.1, or as the gradient of the scalar
potential. Numerically, the derived potential fields may not be identical for different methods.
Finally, since helicity estimation methods are developed for application to researchrelevant
dataset, all employed tests are defined on discretized grids of moderate to highresolution,
and therefore violate the solenoidal property to some extent (cf. Valori et al. 2013). Different
methods might be affected differently by small violations of the solenoidal property of the
test field. Substantial violations of the solenoidal property are not considered here since the
very definition of HV is devoid of meaning in that case.
We report in Tables 8–10 the complete listing of all employed metrics for all FV methods
and test cases considered in this study. On the other hand, in the next sections we provide
concise summaries of the tables’ values for subsets of test cases and/or methods. To this
purpose, we compute the mean of the relevant HV values, and a relative spread around it,
defined as the standard deviation of the HV values distribution over the mean. In addition, in
order to discern among different factors influencing HV values, several diagnostic metrics
are here introduced.
3.1 Accuracy of Vector Potentials
The vector potentials required in the helicity computation of Eq. (36) must reproduce the
correspondent magnetic fields as accurately as possible. In order to compare two vector
fields X and Y in V we employ the metrics
N = 1 −
i Xi − Yi  ,
E =
which are, respectively, the complement of the normalized vector error and the energy ratio,
introduced by Schrijver et al. (2006).2 Both are unity if Xi = Yi in all grid points i in V .
The metrics are applied to the pair (X = Bp, Y = ∇ × Ap) or (X = B, Y = ∇ × A), to
quantitatively assess the accuracy of a vector potential in reproducing the corresponding
magnetic field. Additional metrics defined in Schrijver et al. (2006) are either not particularly
sensitive, or not providing essential additional information in the cases examined below. For
the interested reader, they are listed in Appendix B.
For the metrics defined here, as well as for the integral in Eq. (36), we use standard
numerical prescriptions, as those in Appendix A of Valori et al. (2013). In particular, we
compute the curl and divergence operators using a secondorder, centraldifference discretization
scheme for points in the interior of V . Values on the volumebounding surface, ∂ V , are taken
from the input magnetic fields.
3.2 Quantification of the Solenoidal Property
The test cases used in this article must have a value of ∇ · B small enough to be considered
numerically solenoidal. In order to quantify the level of solenoidality, we apply here the
decomposition of the energy of the magnetic field in V into solenoidal and nonsolenoidal
contributions, as in Valori et al. (2013). Using that decomposition, a fraction Ediv of the
2We use the notation N rather than 1 − EN to avoid confusion with energy symbols.
total magnetic energy can be associated to the nonsolenoidal component of the field. In this
article, all energy contributions are normalized to the total energy, E, of the test case in
exam. More details on the decomposition can be found in Appendix A.
For reference, we occasionally include the divergence metric proposed in Wheatland
et al. (2000) and often used in the literature to test the solenoidal property of discretized
fields, defined as the average over all n grid nodes, fi (B) = ( i fi )/n, of the fractional
flux,
fi (B) =
v ∇ · Bi dv
∂v Bi  dS
through the surface ∂ v of an elementary volume v including the node i. The rightmost
expression in Eq. (40) holds for a grid of uniform and homogeneous resolution . Therefore,
it may be appropriately used as a metric for the methods analyzed in this study, since they all
are based on uniform Cartesian grids. The smaller the value of fi  , the more solenoidal the
field. However, the actual value of this metric depends on the number of grid points n and
the resolution , therefore it makes most sense to apply it to identical discretized volumes.
In addition, this metric is used when the energy one is not applicable, for instance, when
checking the solenoidal property of vector potentials in Coulomb gauge methods, fi (A) ,
as in Sect. 7.
4 Test Fields
In order to be able to critically test the different helicity computation methods of Sect. 2,
the test fields used in this study are chosen such that they represent challenging tests, at
the same time including aspects of relevance for solar physics. From the point of view of
the fields’ structure, we include both compact structures with concentrated currents as in
Fig. 1b and c, as well as more extended structures with currents threading also the lateral
and top boundaries, as in Fig. 1a and d. Similarly, we consider both static (Fig. 1a and b)
and time evolving (Fig. 1c and d) fields, representing typical magnetostatic and
magnetohydrodynamic applications. In the following, the most relevant properties of the test fields
used in this study are discussed in some detail. In particular, the solenoidal properties of
the discretized magnetic fields are quantified using the method described in Sect. 3.1. The
results of the analysis of the input fields is given in full in Table 7, and a selection thereof
is graphically presented in the following sections. In particular, Table 7 also reports the
fractional flux defined in Eq. (40) for all test fields considered here.
4.1 Low and Lou
The Low and Lou model (Low and Lou 1990; Wiegelmann and Neukirch 2006, hereafter
LL) is a 2.5D solution of the zeroβ GradShafranov equation, analytical except for the
numerical solution of an ordinary differential equation. In particular, the four LL cases
considered here are different discretizations of the same solution (corresponding to n = 1, l = 0.3,
φ = π/4 in the notation of Low and Lou 1990), in the same V = [−1, 1] × [−1, 1] × [0, 2]
volume. From this solution, four test cases are constructed where the above volume is
discretized using 32, 64, 128, and 256 nodes per side. Since the solution of the ordinary
differential equation that defines the LL field is always the same in the four cases, the only
factor changing among the different LL cases is the resolution. As a test for helicity
methods, LL represents a largescale, forcefree field with largescale smooth currents distributed
Fig. 1 Representative field lines of the four employed test cases: (a) low and Lou (LL) forcefree equilibrium,
see Sect. 4.1; (b) Titov and Démoulin (TD) forcefree equilibrium for N = 1 and = 0.06, see Sect. 4.2;
(c) snapshot at t = 155 of the stable MHD simulation (MHDst), and (d) snapshot at t = 140 of the unstable
MHD simulation (MHDun), see Sect. 4.3. Yellow field lines depict the flux rope, red field lines belong to the
ambient magnetic field, except for the LL case where no flux rope is present and the field lines’ color scheme
does not apply. A cyan semitransparent isocontour of the current density is shown at 15 %, 35 %, 11 %,
and 10 % of its maximum in the four cases (a) to (d), respectively
in the entire volume. This latter aspect is not supported by observations of solar active
regions. However, being almost analytic, the LL equilibrium offers a tightly controlled test
case. Here, this test field is used mostly in Sect. 6.2 for exploring the dependence of the
finite volume methods on spatial resolution.
Figure 2a shows that for the LL test cases, while Efree 26 % and hardly changes with
resolution, the solenoidal error of the input field varies from 0.1 % to 4 % as pixels become
coarser (see also Table 7). Hence, the value of HV for the LL cases computed by the
different methods will be affected simultaneously by the combined effects of resolution in the
computation of the vector potentials on one hand, and by the different degree of violation of
the solenoidal property by the test field on the other.
4.2 Titov and Démoulin
The Titov and Démoulin model (Titov and Démoulin 1999, hereafter TD) is a parametric
solution of the 3D forcefree equations constituted by a current ring embedded in a confining
potential field. By considering the portion of the ring above a given (photospheric) plane, the
TD equilibrium is possibly the simplest 3D model of a bipolar active region with localized
direct currents, see Fig. 1b. Differently from the LL case, the current is tethering only the
bottom boundary, while the field is potential on the lateral and top boundaries of the volume
considered here.
The TD model has significant topological complexity, in the sense that, for different
values of its defining parameters, can exhibit a finite twist of N endtoend turns, bald patches,
and an hyperbolic flux tube (Titov et al. 2002). In this article we consider six realizations
of that solution, in different combinations of twist and spatial resolution. Unless differently
Fig. 2 Normalized free (blue connected triangles) and nonsolenoidal (green connected crosses) energies
for the test cases (a) LL; (b) TD as a function of twist; (c) TD as a function of resolution; (d) MHDst;
(e) MHDun (f) divergence test MHDstdiv(B), as a function of δ, see Sect. 7
The implementations of DeVore’s method (DeVore_KM, DeVore_GV and DeVore_SA)
yields very similar results, both in values and in trends. The trend already noticed in Fig. 7b
is confirmed by the accuracy of the vector potential quantified by N in Fig. 8c. From the
analysis in Sect. 6.2, we tend to attribute this difference to the accuracy in the solution of
Eq. (8) being higher for DeVore_SA with respect to DeVore_KM and DeVore_GV.
Similarly, the Coulomb_SY and Coulomb_JT methods deliver analogous HV and N curves (see
Fig. 8b and c). However, they strongly differ in the accuracy with which the Coulomb gauge
condition is respected, as Fig. 8d shows. Confirming the result of the LL case of Sect. 4.1,
the Coulomb_JT method respect the Coulombgauge condition extremely well for Ap
( fi  10−9), and still excellently for A ( fi  10−6 or lower). Interestingly, fi (Ap)
is practically independent of Ediv, confirming that the strategy adopted by the Coulomb_JT
method for solving equation (10) is very much able to handle flux unbalance resulting
from solenoidal errors in V . However, even though the solenoidal property of the vector
potentials is definitely better verified by the Coulomb_JT method than by the Coulomb_SY,
the accuracy of the vector potential does not reflect this trend. The Coulomb_SY vector
potentials, indeed, have much higher values of fi (Ap) fi (A) 10−4 than those
by the Coulomb_JT method. In this sense, the divergencecleaning strategy adopted by
Coulomb_SY is less efficient in imposing the Coulomb gauge. However, such values of
fi  seem still low enough to guarantee a relative high accuracy, as testified by N of
Fig. 8c.
The test discussed here, of course, does not pretend to be general in assessing the
influence of the nonsolenoidality on helicity estimations, and further tests are likely to be
required. However, it represents one wellcontrolled example that enables an estimate of the
degree of confidence of a helicity estimations for given a finite nonsolenoidality level.
In summary, according to our tests, errors in respecting the solenoidal constraint might
be still ignorable as long as Ediv is below 1 %, but become abruptly more important above
that threshold. For a dataset with Ediv comprised between 1 and 8 %, using on FV method or
the other would lead up to 6 % difference in the estimation of HV (excluding Coulomb_GR
method). For higher Ediv, the gauge invariance starts to have significant effects.
According to the above discussion, most of the tests employed in the remaining sections
of this work have nonsolenoidal contributions that are mostly ignorable, with few data points
where Ediv (and hence its influence on HV ) is of the order of few percent (cf. Fig. 2 and
Table 7).
8 Discussion of FV Methods
The previous sections discuss the results of FV methods when applied to a variety of test
fields that differ for topological complexity, importance and distribution of currents in the
volume, stability properties. Factors that seem to influence the accuracy in the computation
of HV in FV methods range from the distribution of strong currents at the top boundary (see
Coulomb_JT in Sect. 5.2 in particular), to the solver employed for the construction of the
reference potential field, as for the DeVore methods in Sect. 6.1. In this sense, the solarlike
separation between a currentcarrying and a more potential part of the coronal field seems
to favor accuracy in helicity computations, possibly because in this way currents do not
affect the boundary conditions for Coulomb methods. However, as TD and LL tests show,
the accuracy of the helicity computed by different methods is not always directly related
to the accuracy of the vector potentials in reproducing the corresponding fields. Two of the
Coulomb methods, Coulomb_JT and Coulomb_SY, despite their lesser accuracy in solving
for the vector potentials (when compared to the accuracy of the DeVore methods), they still
deliver a helicity in line with that obtained by the DeVore methods. On the other hand,
while, e.g., Coulomb_GR is of similar (in)accuracy as the other two Coulomb methods,
it delivers slightly different helicity values. Likely, given the nonlocal nature of HV , such
details depend on the particular spatial distribution of the solenoidal errors on a case by case
basis, which may affect how the different distributions of values combine into HV . This is
in a way confirmed by the remarkable absence of differences in the MHDst case, in which
HV is dominated by the large contribution of the potential field. Similarly, we find no strong
influence on the methods’ accuracy by the amount of twist in the field of the TD case, where
the selfhelicity HJ is almost a factor 100 smaller than the total helicity HV , Sect. 5.1.
A strong effect on helicity values is found from errors in the solenoidal property of
the input field. In Sect. 7 a test field is considered that has increasing value of solenoidal
error, as measured by the normalized fraction of the energy associated with magnetic
monopoles, Ediv.
We find a rapid increase of HV fluctuations as Ediv grows above 1 %, see e.g., Fig. 8d. We
also extensively test the effect of resolution, which is found to affect helicity values basically
in two ways: Directly, by affecting the solution of the Poisson solvers that compute the scalar
or vector potentials, and indirectly, by increasing the solenoidal error in the input field and
consequently weakening the consistency between potentials and boundary conditions. By
increasing the solenoidal error, resolution may influence the helicity value significantly in
heavily underresolved cases, as in Fig. 7, for which systematic quantifications of solenoidal
errors must be put in place. On the other hand, when resolution is not pushed to limits,
differences between methods account for larger spread in HV values than resolution, as
Fig. 6a and Sect. 6.1 show.
From the point of view of the accuracy, vector potentials computed with the DeVore
methods reproduce the input field more accurately in most of the cases, see e.g., Fig. 6.
Among the DeVore methods, accuracy is improved mostly for better solutions of the Laplace
equation defining the potential field, Eq. (8). Other details of implementation, such as the
definition of derivative and integrals, play a minor role.
Coulomb methods, need to impose the solenoidality of the vector potential in the entire
volume, and may suffer more from the inconsistent boundary values for Laplace/Poisson
equations that the methods solve. In this respect, the Coulomb_SY strategy of divergence
cleaning is not as efficient as the parametric tuning of the integration constants (the cf
constants of Sect. 2) employed by Coulomb_JT (see Figs. 7f and 8d). On the other hand,
the accuracy of the Coulomb_SY method in the construction of the vector potentials is
found in all tests performed here to be always better than that of the Coulomb_JT method.
Therefore, the accuracy of the vector potential is not directly influenced by the accuracy
with which the gauge condition is satisfied. The computation of A poses in general more
difficulties for Coulomb methods than that of Ap, in a way contrary to the DeVore methods.
In particular, the Coulomb_JT method seems to be sensitive to current on the boundaries (as
MHDun show, see Fig. 5), possibly, because they yield inconsistent boundary values for
Eqs. (10), (13). The Coulomb_GR method, finally, has mostly issues in computing an A that
reproduces the field accurately enough, which result in the largest departures of HV from
average values of all methods.
In the tests presented here it is clearly found that the fulfillment of the Coulomb gauge is
quite variable among the Coulomb methods. However, using the LL test case as a reference,
we find that an average value of the fractional flux of fi (A) below 10−3 yields helicity
values comparable with those from DeVore methods (e.g., within 0.8 % in the
moderateresolution case of n = 64, see Sect. 4.1).
Hence, in a balance between accuracy and applicability, a clear advantage of DeVore
methods over Coulomb methods is that the gauge can be imposed exactly in the former,
whereas it need to be insured numerically on the latter. This translates into simpler and
more efficient implementations of the method, and to more accurate estimations of HV .
On the other hand, the Coulomb gauge yields generally simpler analytical expression of
more straightforward interpretation by eliminating those terms that depend on the divergence
of the vector potentials. This offers, for instance, a possibility of a better comparison and
integration with FI methods. Also, the Coulomb gauge allows a natural interpretation of
helicity in terms of Gauss linking number, see e.g., Berger (1999) and references therein,
which is an interpretation tool often used in helicity studies.
9 Comparison with Discrete FluxTubes Methods
On the grounds of the discussion of the previous sections, and with the limitations there
specified, we consider here one of the FV methods, namely the DeVore_GV one, to give
the correct value of helicity, and we compare the DT methods against it. As discussed in
Sect. 2.3, the comparison between DT and FV methods is by necessity limited to the
estimated helicity value, given the different level of approximation and required information of
the two groups of methods.
9.1 TwistNumber Method
The twistnumber method needs identifiable fluxrope structures in order to be applied. The
dependence of the method on some of its parameters (e.g., on the choice of the QSL surface
defining the flux rope and on the location of the axes of the flux rope, see Sect. 2.3.1) is
specific to the method and requires a testing strategy that is different from the one applicable
to the other FV methods. Therefore, we defer that discussion to a separate dedicated paper,
in preparation at the time of writing, where the TN method is tested using TD, MHDst and
MHDun cases, as well as some nonlinear forcefree field models (Guo et al. 2017). Here,
we only report the results of the application of the TN method to the TD cases.
The Q map (as defined by Eq. (24) in Titov et al. 2002) defining the QSL used to identify
the flux rope volume is shown in the left panel of Fig. 9 in the N = 3 case. We select 100
sample field lines, which are randomly distributed within the QSL, for the magnetic flux
rope. We compute the twist of each field line referred to the axis as described in Sect. 2.3.
The (nonnormalized) twist of the magnetic flux rope Htwist defined in Eq. (32) is computed
as the average of the twist of the sample field lines, and the uncertainties are computed by the
standard deviation of their distribution (see Table 6). An example of the twist distribution of
the sample magnetic field lines as a function of the distance from the flux rope axis is given
in the right panel of Fig. 9.
Since the TN method approximates the total helicity by the selfhelicity only, Fig. 10
compares the notnormalized value of HV,J of Eq. (6) from the DeVore_GV method with
the Htwist from the TN method. A more extensive view of the comparison is reported in
Table 6.
We find that the accuracy of the estimation increases for higher resolutions and higher
twist. In particular, in the N = 3 case, Htwist has the same value of HV,J within errors. In the
other cases, larger differences between the TN method and DeVore_GV method are present.
Fig. 9 Application of the TN method to the TD N = 3 case. Left: a vertical slice of the Q map in the xzplane
at y = 0 corresponding to the middle of the flux rope along its axis. Red (blue) plus sign indicates the position
of the axis (sample field lines). Right: twist of the sample magnetic field lines along the distance, r , from the
origin
Table 6 Helicity in the twist
(upper) and resolution (lower)
TD tests, computed with the
FVDeVore_GV method (2nd to
4th column), and with the TN
method (5th column). Note that
HV , HV,J , and Htwist are not
normalized
N = 0.1
N = 0.5
N = 1
N = 3
= 0.03
= 0.06
= 0.12
−0.0057
−0.0290
−0.0782
−0.0527
−0.0782
−0.0782
−0.0782
−2.23147
−8.33340
−7.21022
−1.78752
−7.20592
−7.21022
−7.21039
−0.08374
−0.55240
−0.08976
−0.54921
−0.55240
−0.55444
−0.16 ± 0.06
−0.66 ± 0.13
−0.50 ± 0.05
−0.087 ± 0.020
−0.56 ± 0.06
−0.50 ± 0.05
−0.49 ± 0.05
Fig. 10 Application of the TN method to the TDtwist (a) and resolution (b) tests. The two curves represent
the values of Htwist of Eq. (32) for the TN method, and of the (nonnormalized) value of HV,J for the
DeVore_GV method
By contrast, the notnormalized values of the total helicity HV = −7.2 for DeVore_GV for
the case N = 1, i.e., the total helicity is two orders of magnitude larger than the current
helicity in the N = 1 case.
We conclude that the quantity Htwist provided by the TN method is an estimation of HV,J
rather than of H (or equivalently HV ), likely because no modeling of the mutual helicity
part is included. Moreover, the accuracy of the method is higher for highly twisted structures,
above N = 1 in out tests. For N = 3, the correct value of HV,J is reproduced.
As for the TN method, the CB method can be compared to FV methods only regarding the
estimated helicity values. However, we test here for the first time some of the assumptions
made in the derivation of the CB method (forcefreeness and minimal connectivitylength
principle, see Sect. 2.3.2), and their impact on the obtained helicity values. These novel
comparisons are made possible by the reliability assessment of FV methods discussed in
the previous sections. Moreover, the comparison of the CB methods with the FI methods is
discussed in detail by Pariat et al. (2017).
The connectivitybased method is designed for applications to solar magnetograms with
a complex flux distribution and an unknown coronal magnetic field. If the lower boundary
includes only two connected partitions the code automatically uses the linear forcefree field
approach of Georgoulis and LaBonte (2007), in which a single value of the forcefree
parameter α is used for the entire volume. When the NLFF mode switches on, the CB method
has the possibility to replaces the total photospheric magnetic flux by the connected
magnetic flux, namely the one included in the magnetic connectivity matrix, see Sect. 2.3. The
CB code was tuned to use almost the entire flux of the magnetograms in all case presented
here, except when differently explicitly stated.
When applied to the TD cases, where currents are localized, the CB in LFFmode tends
to pick up mostly the potential field component. On the other hand, when applied to the LL
cases, where a largescale α is present, the CB method yields values of HV that are three
to four times larger than FV methods. This overestimation effect in case of the LFF field
approximation was also reported by Georgoulis et al. (2012). For these two sets of tests,
where a single dipole appears, the CB method is hampered by the limitations of the linear
forcefree theory. Hence, we do not report further on such applications here.
In the MHDst and MHDun cases, instead, the magnetogram is complex enough to have
more than one connectivity domain. The CB method worked in the NLFF mode for all
snapshots at all times, in both MHDst and MHDun cases. Figure 11a shows an example of
the flux partition and of the resulting connectivity matrix for the MHDst case at t = 150.
Figure 12 compare the CB method with the DeVore_GV finite volume method.
In the MHDst case (Fig. 12a), between t = 50 and t = 95, the of HV values obtained
by the CB method match reasonably well those from the DeVore_GV one, within a factor 2
at most. Starting from t = 95 onwards, however, the HV values obtained by the CB method
settle on a lower, roughly constant value HV = 0.016, on average. At the end of the
simulations this average is about eight times lower than the value reached by the DeVore_GV
method.
On the other hand, in the MHDun case (Fig. 12b), the agreement between the
DeVore_GV and the CB methods is very good (within 9 % on average, from t = 95 onwards),
and the two curves overlap for most of that phase. A local maxima in the CB curve is even
Fig. 11 Application of the CB method to the MHDst case. (a) Magnetic connectivity of the
MHDemergence stable model at t = 150 used for the calculation of magnetic helicity in the CB method. The partitioned
normal magnetic field component on the lower (assumed photospheric) boundary is shown in grayscale, with
black/white contours outlining positive/negativepolarity partitions. The fluxweighted centroids of the
partitions are denoted by crosses. Magnetic connections are projected on the lower boundary by colored lines,
with different colors denoting different connectedflux contents. The connected flux includes approx. 93.6 %
of the total flux present in the field of view. (b) The forcefreeness metric σJ as a function of the height of
the bottom boundary for MHDst at t = 155. The vertical orange dotted line represents the location of the
bottom boundary in the tests of Sect. 9.2, corresponding to z = 8.9
present at the time of the eruption, very much the same as for the DeVore_GV method.
However, this is not distinguishable from previous, even more pronounced ones, and it would be
challenging to identify the time of the eruption only as a decrease of HV in the CB time
series. As a matter of fact, the examination of the time evolution of the magnetic field at z = 0
in the MHDst and MHDun simulations shows that there is very little differences between
the two cases. Since the CB method aims at an approximate estimation of helicity that is
based only on the flux distribution at the (photospheric) bottom boundary at a given time,
the task of distinguishing MHDst from MHDun based only on the field on that one plane is
a very arduous one, indeed. This is probably the reason why the CB method provides similar
average values of HV for both the MHDst (equal to 0.015, from t = 95 onwards) and the
MHDun case (equal to 0.070 on the same time interval).
Besides the similarity of the MHDst and MHDun field distributions at z = 0,
differences with the DeVore_GV method may have several origins. First of all, the CBmethod
approach yields a minimal value of helicity for a given photospheric configuration. In this
sense, it is expected that the CB HV curves lay, on average, below the DeVore_GV ones, as
they do at varying levels in all panels of Fig. 12.
Secondly, the CB method models the coronal field as a discrete collection of a finite
number of constantα flux tubes. Hence, it can be expected to have better chances of success
if the magnetic field is forcefree. Figure 11b shows at a representative time, however, that
in large part of the volume of the MHDst case the field is not forcefree, as quantified
by the relative ratio of the current that is perpendicular to the field in the volume, σJ ≡
Fig. 12 Comparison of HV between the CB and the DeVore_GV method applied to the MHDst (a, c) and
MHDun cases (b, d) in the full domain (a, b) and in the reduced, more forcefree domain (c, d). In panels
(c), (d) CB method (blue crosses), the CB method with 3D connectivity information included (green crosses),
and the DeVore_GV method (red squares). Note that the DeVore_GV values (c, d) are not the same as (a, b)
and in Fig. 4 because the considered volume is different
( V J⊥ dV)/( V J dV). In particular, the value of σJ computed for increasing heights of
the bottom boundary decreases from 0.7 to 0.25 in the first 20 pixels above the bottom
boundary.
In order to test how important is the forcefree assumption in this case, we repeat the HV
calculation with the CB and DeVore_GV methods in a reduced volume starting at z = 8.9,
where σJ at this height has dropped to the value 0.29, see Fig. 11b. The corresponding
curves are shown by the blue crosses in Fig. 12c and d for MHDst and MHDun,
respectively. Since the volume is now changed, also the corresponding DeVore_GV estimations
are recalculated for this reduced, more forcefree volume. In the MHDst case, the average
CB HV after t = 95 is 0.07, which is 3.5 times smaller than the final value obtained by the
DeVore_GV method (against a factor 8 of the fullvolume case). In the MHDun case in
Fig. 12d, the curves obtained by the two methods are slightly closer than in the fullvolume
case of Fig. 12b, although only marginally (the ratio of the end values being 1.5). Hence, the
application to the more forcefree, upper part of the volume improves the match between the
CB and the DeVore_GV results, markedly so in the MHDst case, which is a clear indication
that the fulfillment of the forcefree requirement may help to partially compensate for the
HV underestimation in Fig. 12a.
Thirdly, the connectivity map between flux partitions is obtained by the CB method as
part of the minimization method discussed in Sect. 2.3. However, when the 3D coronal field
is available, the true connectivity map can be constructed from the numerical simulation,
and the influence of the minimization tested. The result of such a test are represented by
the green symbols in Fig. 12c and d, for the stable and unstable cases, respectively (again
in the reduced, more forcefree volume). It must be noticed that using the 3D information
implies a decrease of the amount of flux included in the connectivity matrix to 60–80 %,
compared to the 95 % or more employed in the two cases above, since tracking of field
lines intersecting the lateral and top boundaries cannot be completed in the CB method,
even in case these lines return to the simulation volume by intersecting a different boundary
location. The contribution of these lines is, then, ignored.
In general, Fig. 12c demonstrates that the knowledge of the true connectivity in
the MHDst case improves the matching between CB and DeVore_GV methods. The
CBaverage (after t = 95) HV value is slightly above 0.10, which results in just a factor
two between the corresponding end values. In the MHDun case, differences with the
standard application of the CB method that does not used the threedimensional connectivity
information are less significant, except for a slight increase of the mismatch with the
DeVore_GV method (e.g., the ratio between end values of HV obtained by the two methods
is down to 1.9). Such relatively small variations can be explained in terms of reduced
connected flux. We also recall that, for the CB method, an error analysis by Moraitis et al. (2014)
is available that was not included in the discussion presented in this article. It is likely that
some of the fluctuations discussed above fall within the error estimation provided by that
analysis.
In conclusion, in the MHDun case the agreement between CB and DeVore_GV (and, by
extension, FV) methods is within 10 %, which is very good considering the much more
limited amount of information that the CB method requires. On the other hand, in the MHDst
case the helicity is significantly underestimated, by a factor eight at the end of the
simulation. The deficit in HV that the CB method shows in the MHDst case, can be partially
ascribed to the exiguous differences between the MHDst and MHDun cases in terms of
flux distributions at z = 0. Additionally, the field at that plane is not quite forcefree, and the
CB results are shown to improve for a more forcefree test volume. Furthermore, the
minimization of the connectivity matrix seems to represent connectivity fairly well, in the sense
that the implementation of the full threedimensional information, although improving the
CB estimation in the MHDst case, does not entirely remove the underestimation of HV .
10 Conclusions
In this work we review, benchmark, and compare the currently available methods for the
computation of the relative magnetic helicity, HV , in finite volumes. Given that the
threedimensional magneticfield solutions we use are common for all tested methods, the problem
essentially reduces to computing the vector potential of a given discretized input magnetic
field. The considered methods group into Coulomb (∇ · A = 0) and DeVore (Az = 0)
methods, according to the gauge in which the vector potentials are written. A total of six
different implementations including three Coulomb methods (Coulomb_JT, Coulomb_SY, and
Coulomb_GR) and three DeVore methods (DeVore_SA, DeVore_KM, and DeVore_GV) are
included which differ in the equations they solve and/or in their numerical implementations.
Details of these implementations can be found in Thalmann et al. (2011) (Coulomb_JT),
Yang et al. (2013b) (Coulomb_SY), Rudenko and Anfinogentov (2014) (Coulomb_GR),
Valori et al. (2012) (DeVore_GV), Moraitis et al. (2014) (DeVore_KM), and in Sect. 2.2
(DeVore_SA). Accordingly, a different level of numerical complexity and computational
effort is required to solve for the vector potentials, with the Coulomb methods being
essentially far more demanding than DeVore methods (see Sect. 2). As a case in point, the
Coulomb_GR method could not be tested on cases above 1283 pixels due to the large
running time required by its current implementation.
The tested methods are put under severe strain by choosing a variety of numerical test
input fields that are considered to be relevant for helicity studies in solarphysics applications,
from 3D forcefree equilibria (LL and TD of Sects. 4.1 and 4.2, respectively), to snapshots of
timedependent nonforcefree MHD simulations of flux emergence (Sect. 4.3). Depending
on details of the test field being studied, the accuracy in the computation of HV by different
methods is found to vary to some extent, especially for Coulomb methods. We can, however,
definitely conclude that in solarlike cases practically all FV methods converge to the same
helicity value within few percent.3 Such a spread is likely to be overrun by other sources of
errors in applications to observed—and reconstructed—coronal fields, but it is definitely to
be considered in helicity estimations of numerical simulations.
More in detail, the helicity values HV in the most relevant test of timedependent MHD
evolution in a coronal model volume (i.e., the MHDst and MHDun tests of Sect. 4.3) show
a very good agreement of few percent between different methods, somewhat independently
of the details in the vector potentials computation, see e.g., Fig. 4. Such errors are as small as
0.2 % in the case where the field is slowly evolving, and always below 3 % even in the highly
dynamical eruptive phase. Similarly, excluding the Coulomb_GR case, despite differences
in the accuracy of the vector potential computation, helicity values computed by different
methods for the TDtwist case are within 2 %. In other words, when helicity computations
are applied to numerical volumes as in this article, differences in the way HV is computed
can amount to 3 % at most.
Such an agreement is tested and verified to hold independently of the dynamical evolution
of the field and consistently throughout the MHD evolution of an eruption, thereby justifying
the application of FV methods to the study of helicity in numerical simulations, and for
benchmarking other helicity computation methods. For instance, Pariat et al. (2017) employs
finite volume methods to benchmark helicityflux integration methods applied to MHDst
and MHDun evolution, where the flux of helicity is estimated by the photospheric evolution
of the field.
In addition to finite volume methods, we also include other two methods that use (Guo
et al. 2010, TN) or optionally make use of (Georgoulis et al. 2012, CB) the 3D information of
the magnetic field in the volume, see Sect. 9. The TN method estimates the helicity content
of a field by parametrically fitting a flux rope to it. Therefore, it is applicable to tests that
include an identifiable flux rope, namely to TD, MHDst, and MHDun cases. For the TD
cases, we find that the TN method yields relatively accurate estimations of twist and possibly
of HV,J in hightwist cases, but not of the total helicity HV . A report on the application of
the TN method to other cases, including the MHDst and MHDun ones, is in preparation
by Guo et al. (2017).
3From this statistic, the Coulomb_GR method is excluded.
The CB method is designed to be applied to complex photospheric flux distributions
with an unknown coronal magnetic field, where a multipolar partition of the flux yields
a coronal field approximated by a collection of flux tubes of constant Jz/Bz. In addition,
a minimal free energy and the corresponding relative helicity are sought. Cases like the
LL and TD have a too simple connectivity for the CB method, which then falls back to a
single fluxtube, purely linear approximation of the (forcefree) field. In the more complex
MHDun case, the CB method provides an estimation of the helicity that is, on average,
within 10 % of the one obtained by FV methods. This is a positive result given that the CB
method employs only the photospheric information, whereas FV methods use the full 3D
information about the coronal field. The MHDst case poses more difficulties to the coronal
field approximation within the CB method, which, in this case, underestimate significantly
the helicity content of the field.
Finite volume methods can be used when the magnetic field is known in the entire
volume of interest. However, in applications to solar observations, the magnetic field is typically
known only on a surface at photospheric heights, and only with limited accuracy. Hence, in
order to know the magnetic helicity in a given coronal volume, a model need to be computed
that approximates the coronal field on the base of its photospheric values, which introduce
an additional dependence to the estimated HV values. The impact on helicity values of the
employed coronal model, being it from a nonlinear forcefree extrapolation or from a
datadriven simulation, is yet to be tested. Alternatively, the CB method can be used, as it is
designed for such cases. A further alternative would be to compute the flux of helicity
passing through the “photospheric plane” in time. Reviewing and benchmarking FI methods for
the estimation of the helicity flux is the subject of Pariat et al. (2017).
Acknowledgements This article results from the work of the ISSI International Team on Magnetic
Helicity estimations in models and observations of the solar magnetic field. The authors are pleased to thank Marc
deRosa, James Leake, Bernhard Kliem, and Tibor Török for making their numerical data available. G. Valori
acknowledges the support of the Leverhulme Trust Research Project Grant 2014051. E. Pariat
acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR15CE310001
(project HeliSol). S. Anfinogentov acknowledges support of the Russian Foundation of Basic Research
under grants 150201089a, 150203835a, 163200315mola, and 153220504molaved, and the Federal
Agency for Scientific Organisations base project II.16.3.2 “Nonstationary and wave processes in the solar
atmosphere” and thanks George Rudenko for providing the source code of Coulomb_GR method and
contributing to development of DeVore_SA method. Y. Guo is supported by NSFC (11203014 and 11533005),
NKBRSF (2011CB811402 and 2014CB744203), and the mobility grant from the Belgian Federal Science
Policy Office (BELSPO). K. Moraitis and M.K. Georgoulis were partially supported by the European Union
(European Social Fund—ESF) and Greek national funds through the Operational Program “Education and
Lifelong Learning” of the National Strategic Reference Framework (NSRF)—Research Funding Program:
Thales. Investing in knowledge society through the European Social Fund. J.K. Thalmann acknowledges
support from the Austrian Science Fund (FWF): P25383N27. F. Chen acknowledge the support by the
International MaxPlanck Research School (IMPRS) for Solar System Science at the University of Göttingen.
S. Yang is supported by grants KJCX2EWT07 and XDB09040200 from the Chinese Academy of Sciences;
grants 11221063, 11078012, 11178016, 11173033, 11125314, 10733020, 10921303, 11303053, 11573037,
and 10673016 of National Natural Science Foundation of China; grant 2011CB811400 of National Basic
Research Program of China; and KLSA2015 of the Collaborating Research Program of National Astronomical
Observatories. Part of the calculations presented here were done on the quadricore biXeon computers of the
Cluster of the Division Informatique de l’Observatoire de Paris (DIO).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0
International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source,
provide a link to the Creative Commons license, and indicate if changes were made.
Appendix A: Decomposition of the Magnetic Energy
The method that we employ to quantify the error in the solenoidal property (Valori et al.
2013) is basically a numerical verification of Thomson’s theorem, and allows to quantify
the effect of a (numerical) finite divergence of the magnetic field in terms of associated
energies. In order to obtain such a decomposition, the magnetic field is firstly separated
into a potential and a current carrying part. Secondly, each part is additionally split into a
solenoidal and a nonsolenoidal contribution, using a Helmholtz decomposition. As a result,
the magnetic energy E is split into
E = Ep,s + EJ,s + Ep,ns + EJ,ns + Emix,
where Ep,s and EJ,s are the energies associated to the potential and currentcarrying
solenoidal contributions, Ep,ns and EJ,ns are those of the nonsolenoidal contributions, and
Emix is a nonsolenoidal mixed term, see Eqs. (7), (8) in Valori et al. (2013) for the
corresponding expressions. All terms in Eq. (43) are positively defined, except for Emix. For a
perfectly solenoidal field, it is Ep,s = Ep, EJ,s = E − Ep, and Ep,ns = EJ,ns = Emix = 0, i.e.,
Thomson’s theorem is recovered.
In most of the analysis in this article we consider a single number for characterizing the
energy associated to nonsolenoidal components of the field, given by
Ediv = Ep,ns + EJ,ns + Emix,
which generally overestimates such errors by hindering possible cancellations. However,
since we consider numerically accurate models, we neglect such overestimations, unless
explicitly stated. In that case, also the sum with sign Ens is considered, corresponding to
Eq. (44) with Emix instead of Emix. The full decomposition for each test case considered
in the article can be found in Table 7. Also, in the article we associate the free energy of the
field with the solenoidal component of its currentcarrying part, i.e., Efree = EJ,s.
For a given discretized field, the accuracy of the decomposition can be easily quantified
by checking to what precision the sum of the right hand side of Eq. (43) reproduces the total
energy E, normalized to E. In relative terms, such an error is smaller than 10−7 in all cases
discussed in this article. In this article, all energy contributions are normalized to the total
energy, E, of the test case in exam.
Appendix B: Accuracy of Vector Potentials
In addition to the metrics Eqs. (38), (39), we include here also the remaining
CVec = (
1
CCV = N
1
M = 1 − N
Yi − Xi 
Table 7 Energy metrics of test cases, see Appendix A. The first two columns indicate the test field
identification label and the total magnetic energy in model units [E]. The following columns indicate
the potential/solenoidal [Ep,s], currentcarrying/solenoidal [EJ,s], potential/nonsolenoidal [Ep,ns],
currentcarrying/nonsolenoidal [EJ,ns], nonsolenoidal mixed [Emix] contributions normalized to the corresponding
total magnetic energy, E. The last column contains the average fractional flux [ fi (B) ]. In the text, any
contribution Exx is intended to be normalized to E, i.e., it is here indicated as Exx/E, for each test case
separately
Ep,s/E EJ,s/E Ep,ns/E
respectively the vector correlation, the CauchySchwartz metric, and the complement of
the mean vector error, from Schrijver et al. (2006). We provide the full set of values for
each case considered in the paper in the following tables. In particular, Table 8 reports the
metrics for the LL and TD, Table 9 those for MHDst and MHDun, and Table 10 those for
the divergence and gauge tests of Sect. 7.
7 0 2 2 2 7 7 0 5 5 7 8 3 1 3 9 3
V 0050. 0290. 0780. 0780. 0780. 0520. 0050. 0290. 0708. 0078. 0078. 0052. 0054. 0073. 0111. 0111. 0084.
H
− − − − − − − − − − − − − − − − −
0 8 4 2 3 0 1 0 1 5 3 2 0 5 7 1 8 5 6 0 0 0 8 9 3
63 31 66 11 64 63 41 76 21 74 34 10 63 08 44 36 13 66 11 4 3 1 6 1 4
7 6 4 6 1 7
l_bouoSYm l_bouoSYm l_bouoSYm l_bouoSYm l_bouoSYm l_bouoSYm l_buooSYm ree_oKVM ree_oKVM ree_oKVM ree_oKVM ree_oKVM ree_oKVM ree_oKVM lJb_uooTm lJb_ouoTm lJbuo_oTm lJbou_oTm
C C C C C C C D D D D D D D C C C C
C.E. Alissandrakis , On the computation of constant alpha forcefree magnetic field . Astron. Astrophys . 100 , 197  200 ( 1981 )
T. Amari , J.J. Aly , A. Canou , Z. Mikic , Reconstruction of the solar coronal magnetic field in spherical geometry . Astron. Astrophys . 553 , 43 ( 2013 ). doi:10.1051/ 0004  6361 /201220787
M.A. Berger , Rigorous new limits on magnetic helicity dissipation in the solar corona . Geophys. Astrophys. Fluid Dyn . 30 , 79  104 ( 1984 ). doi:10.1080/03091928408210078
M.A. Berger , Introduction to magnetic helicity . Plasma Phys. Control. Fusion 41 , 167  175 ( 1999 )
M.A. Berger , in Topological Quantities in Magnetohydrodynamics, ed. by K. Zhang, A. Soward , C. Jones ( 2003 ), pp. 345  374 . doi:10.1201/9780203493137.ch10
M.A. Berger , G.B. Field , The topological properties of magnetic helicity . J. Fluid Mech . 147 , 133  148 ( 1984 ). doi:10.1017/S0022112084002019
M.A. Berger , C. Prior , The writhe of open and closed curves . J. Phys. A, Math. Gen . 39 , 8321  8348 ( 2006 ). doi:10.1088/ 0305  4470 /39/26/005
A. Brandenburg , K. Subramanian , Astrophysical magnetic fields and nonlinear dynamo theory . Phys. Rep . 417 , 1  209 ( 2005 ). doi:10.1016/j.physrep. 2005 .06.005
J. Chae , H. Wang , J. Qiu , P.R. Goode , L. Strous , H.S. Yun , The formation of a prominence in active region NOAA 8668 . I. SOHO/MDI observations of magnetic field evolution . Astrophys. J . 560 , 476  489 ( 2001 ). doi:10.1086/322491
P. Démoulin , E. Pariat , Modelling and observations of photospheric magnetic helicity . Adv. Space Res . 43 , 1013  1031 ( 2009 ). doi:10.1016/j.asr. 2008 .12.004
P. Démoulin , E. Pariat , M.A. Berger , Basic properties of mutual magnetic helicity . Sol. Phys . 233 , 3  27 ( 2006 ). doi:10.1007/s11207 006  0010 z
C.R. DeVore , Magnetic helicity generation by solar differential rotation . Astrophys. J . 539 , 944  953 ( 2000 ). doi:10.1086/309274
J.H. Finn , T.M.J. Antonsen , Magnetic helicity: what is it, and what is it good for? Comments Plasma Phys. Control. Fusion 9 , 111  126 ( 1985 )
U. Frisch , A. Pouquet , J. Leorat , A. Mazure , Possibility of an inverse cascade of magnetic helicity in magnetohydrodynamic turbulence . J. Fluid Mech . 68 , 769  778 ( 1975 ). doi:10.1017/S002211207500122X
M.K. Georgoulis , B.J. LaBonte , Magnetic energy and helicity budgets in the active region solar corona . I. Linear forcefree approximation . Astrophys. J . 671 , 1034  1050 ( 2007 ). doi:10.1086/521417
M.K. Georgoulis , K. Tziotziou , N.E. Raouafi , Magnetic energy and helicity budgets in the activeregion solar corona . II. Nonlinear forcefree approximation . Astrophys. J . 759 , 1 ( 2012 ). doi:10.1088/ 0004  637X/759/1/1
Y. Guo , M.D. Ding , B. Schmieder , H. Li , T. Török , T. Wiegelmann , Driving mechanism and onset condition of a confined eruption . Astrophys. J . 725 ( 1 ), 38  42 ( 2010 ). doi:10.1088/ 2041  8205 /725/1/l38
Y. Guo , M.D. Ding , X. Cheng , J. Zhao , E. Pariat , Twist accumulation and topology structure of a solar magnetic flux rope . Astrophys. J . 779 , 157 ( 2013 ). doi:10.1088/ 0004 637X/779/2/157
Y. Guo , E. Pariat , G. Valori , S. Anfinogentov , F. Chen , M. Georgoulis , Y. Liu , K. Moraitis , J.K. Thalmann , S. Yang , Magnetic helicity estimations in models and observations of the solar magnetic field . Part III: the twist number method. Astron. Astrophys . ( 2017 ), submitted
B. Kliem , T. Török , Torus instability . Phys. Rev. Lett . 96 (25), 255002 ( 2006 ). doi:10.1103/PhysRevLett.96. 255002
B.J. LaBonte , M.K. Georgoulis , D.M. Rust , Survey of magnetic helicity injection in regions producing Xclass flares . Astrophys. J . 671 , 955  963 ( 2007 ). doi:10.1086/522682
J.E. Leake , M.G. Linton , T. Török , Simulations of emerging magnetic flux . I. The formation of stable coronal flux ropes . Astrophys. J . 778 , 99 ( 2013 ). doi:10.1088/ 0004 637X/778/2/99
J.E. Leake , M.G. Linton , S.K. Antiochos , Simulations of emerging magnetic flux . II. The formation of unstable coronal flux ropes and the initiation of coronal mass ejections . Astrophys. J . 787 , 46 ( 2014 ). doi:10.1088/ 0004 637X/787/1/46
B.W. Lites , Remote sensing of solar magnetic fields . Rev. Geophys . 38 , 1  36 ( 2000 )
Y. Liu , P.W. Schuck , Magnetic energy and helicity in two emerging active regions in the Sun . Astrophys. J . 761 , 105 ( 2012 ). doi:10.1088/ 0004 637X/761/2/105
D.W. Longcope , A. Malanushenko , Defining and calculating selfhelicity in coronal magnetic fields . Astrophys. J . 674 , 1130  1143 ( 2008 ). doi:10.1086/524011
B.C. Low , Y.Q. Lou , Modeling solar forcefree magnetic fields . Astrophys. J . 352 , 343  352 ( 1990 )
H.K. Moffatt , The degree of knottedness of tangled vortex lines . J. Fluid Mech . 35 (01), 117  129 ( 1969 )
K. Moraitis , K. Tziotziou , M.K. Georgoulis , V. Archontis , Validation and benchmarking of a practical free magnetic energy and relative magnetic helicity budget calculation in solar magnetic structures . Sol. Phys . 289 , 4453  4480 ( 2014 ). doi:10.1007/s11207 014  0590 y
A. Nindos , J. Zhang , H. Zhang , The magnetic helicity budget of solar active regions and coronal mass ejections . Astrophys. J . 594 , 1033  1048 ( 2003 ). doi:10.1086/377126
E. Pariat , P. Démoulin , M.A. Berger , Photospheric flux density of magnetic helicity . Astron. Astrophys . 439 , 1191  1203 ( 2005 ). doi:10.1051/ 0004  6361 : 20052663
E. Pariat , G. Valori , P. Démoulin , K. Dalmasse , Testing magnetic helicity conservation in a solarlike active event . Astron. Astrophys . 580 , 128 ( 2015 ). doi:10.1051/ 0004  6361 /201525811
E. Pariat , G. Valori , S. Anfinogentov , F. Chen , M. Georgoulis , Y. Guo , Y. Liu , K. Moraitis , J.K. Thalmann , S. Yang , Magnetic helicity estimations in models and observations of the solar magnetic field. Part II: Flux injection methods . Space Sci. Rev . ( 2017 ), submitted
W.H. Press , S.A. Teukolsky , W.T. Vetterling , B.P. Flannery , Numerical Recipes in FORTRAN. The Art of Scientific Computing ( 1992 )
C. Prior , A.R. Yeates , On the helicity of open magnetic fields . Astrophys. J . 787 , 100 ( 2014 ). doi:10.1088/ 0004 637X/787/2/100
G.V. Rudenko , S.A. Anfinogentov , Very fast and accurate azimuth disambiguation of vector magnetograms . Sol. Phys . 289 , 1499  1516 ( 2014 ). doi:10.1007/s11207 013  0437 y
G.V. Rudenko , I.I. Myshyakov , Gaugeinvariant helicity for forcefree magnetic fields in a rectangular box . Sol. Phys . 270 , 165  173 ( 2011 ). doi:10.1007/s11207 011  9743 4
A.J.B. Russell , A.R. Yeates , G. Hornig , A.L. WilmotSmith , Evolution of field line helicity during magnetic reconnection . Phys. Plasmas 22 ( 3 ), 032106 ( 2015 ). doi:10.1063/1.4913489
D.M. Rust , Spawning and shedding helical magnetic fields in the solar atmosphere . Geophys. Res. Lett . 21 , 241  244 ( 1994 ). doi:10.1029/94GL00003
C.J. Schrijver , M.L. DeRosa , T.R. Metcalf , Y. Liu , J. McTiernan , S. Régnier , G. Valori , M.S. Wheatland , T. Wiegelmann , Nonlinear forcefree modeling of coronal magnetic fields part I: a quantitative comparison of methods . Sol. Phys . 235 , 161  190 ( 2006 ). doi:10.1007/s11207 006  0068 7
P.N. Swartztrauber , R.A. Sweet , Algorithm 541: efficient Fortran subprograms for the solution of separable elliptic partial differential equations [D3] . ACM Trans. Math. Softw . 5 ( 3 ), 352  364 ( 1979 ). doi:10. 1145/355841.355850
J.B. Taylor , Relaxation and magnetic reconnection in plasmas . Rev. Mod. Phys . 58 , 741  763 ( 1986 ). doi:10. 1103/RevModPhys.58.741
J.K. Thalmann , B. Inhester , T. Wiegelmann , Estimating the relative helicity of coronal magnetic fields . Sol. Phys . 272 , 243  255 ( 2011 ). doi:10.1007/s11207 011  9826 2
V.S. Titov , P. Démoulin , Basic topology of twisted magnetic configurations in solar flares . Astron. Astrophys . 351 , 707  720 ( 1999 )
V.S. Titov , G. Hornig , P. Démoulin , Theory of magnetic connectivity in the solar corona . J. Geophys. Res. Space Phys . 107 , 1164 ( 2002 ). doi:10.1029/2001JA000278
T. Török , B. Kliem , The evolution of twisting coronal magnetic flux tubes . Astron. Astrophys . 406 , 1043  1059 ( 2003 ). doi:10.1051/ 0004  6361 : 20030692
T. Török , B. Kliem , V.S. Titov , Ideal kink instability of a magnetic loop equilibrium . Astron. Astrophys . 413 , 27  30 ( 2004 ). doi:10.1051/ 0004  6361 : 20031691
K. Tziotziou , M.K. Georgoulis , N.E. Raouafi , The magnetic energyhelicity diagram of solar active regions . Astrophys. J. Lett . 759 , 4 ( 2012 ). doi:10.1088/ 2041  8205 /759/1/L4
K. Tziotziou , M.K. Georgoulis , Y. Liu , Interpreting eruptive behavior in NOAA AR 11158 via the region's magnetic energy and relativehelicity budgets . Astrophys. J . 772 , 115 ( 2013 ). doi:10.1088/ 0004 637X/ 772/2/115
K. Tziotziou , K. Moraitis , M.K. Georgoulis , V. Archontis , Validation of the magnetic energy vs. helicity scaling in solar magnetic structures . Astron. Astrophys . 570 , 1 ( 2014 ). doi:10.1051/ 0004  6361 /201424864
G. Valori , B. Kliem , T. Török , V.S. Titov , Testing magnetofrictional extrapolation with the TitovDémoulin model of solar active regions . Astron. Astrophys . 519 , 44 ( 2010 ). doi:10.1051/ 0004  6361 /201014416
G. Valori , P. Démoulin , E. Pariat , Comparing values of the relative magnetic helicity in finite volumes . Sol. Phys . 278 , 347  366 ( 2012 ). doi:10.1007/s11207 012  9951 6
G. Valori , P. Démoulin , E. Pariat , S. Masson , Accuracy of magnetic energy computations . Astron. Astrophys . 553 , 38 ( 2013 ). doi:10.1051/ 0004  6361 /201220982
M.S. Wheatland , P.A. Sturrock , G. Roumeliotis , An optimization approach to reconstructing forcefree fields . Astrophys. J . 540 , 1150  1155 ( 2000 ). doi:10.1086/309355
T. Wiegelmann , T. Neukirch , An optimization principle for the computation of MHD equilibria in the solar corona . Astron. Astrophys . 457 , 1053  1058 ( 2006 ). doi:10.1051/ 0004  6361 : 20065281
T. Wiegelmann , G.J.D. Petrie , P. Riley , Coronal magnetic field models . Space Sci. Rev . ( 2015 ). doi:10.1007/ s11214 015  0178 3
L. Woltjer , The stability of forcefree magnetic fields . Astrophys. J . 128 , 384 ( 1958 ). doi:10.1086/146551
S. Yang , J. Büchner , J.C. Santos , H. Zhang , Evolution of relative magnetic helicity: method of computation and its application to a simulated solar corona above an active region . Sol. Phys . 283 , 369  382 (2013a). doi:10.1007/s11207 013  0236 5
S. Yang , J. Büchner , J.C. Santos , H. Zhang , Method of relative magnetic helicity computation ii: boundary conditions for the vector potentials (2013b) . arXiv:1304.3526