# Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods

Gherardo Valori, Etienne Pariat, Sergey Anfinogentov, Feng Chen, Manolis Georgoulis, Yang Guo, Yang Liu, Kostas Moraitis, Julia Thalmann, Shangbin Yang

► To cite this version:

Gherardo Valori, Etienne Pariat, Sergey Anfinogentov, Feng Chen, Manolis Georgoulis, et al. . Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods. Space Science Reviews, 2016, 201 (1-4), pp.147-200. 10.1007/s11214-016-0299-3 . hal-02458981

HAL Id: hal-02458981

https://hal.science/hal-02458981v1

Submitted on 13 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods

Gherardo Valori<fnref n="1" /> · Etienne Pariat<fnref n="2" /> · Sergey Anfinogentov<fnref n="3" /> · Feng Chen<fnref n="4" /> · Manolis K. Georgoulis<fnref n="5" /> · Yang Guo<fnref n="6" /> · Yang Liu<fnref n="7" /> · Kostas Moraitis<fnref n="5" /> · Julia K. Thalmann<fnref n="8" /> · Shangbin Yang<fnref n="9" />

Received: 10 May 2016 / Accepted: 6 October 2016 / Published online: 19 October 2016
© The Author(s) 2016. This article is published with open access at Springerlink.com

**Abstract** Magnetic helicity is a conserved quantity of ideal magneto-hydrodynamics 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 three-dimensional, force-free equilibrium, to magneto-hydrodynamical numerical simulations. Almost all methods are found to produce the same value of magnetic helicity within few percent in all tests. In the more solar-relevant 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

<figure><img src="image_3.png" /></figure>

<footnote>G. Valori
g.valori@ucl.ac.uk</footnote>

<footnote n="1">¹ Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK</footnote>

<footnote n="2">² LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 92190 Meudon, France</footnote>

<footnote n="3">³ Institute of Solar-Terrestrial Physics SB RAS 664033, P.O. Box 291, Lermontov st., 126a, Irkutsk, Russia</footnote>

<footnote n="4">⁴ Max-Plank-Institut für Sonnensystemforschung, 37077 Göttingen, Germany</footnote>

<footnote n="5">⁵ Research Center for Astronomy and Applied Mathematics of the Academy of Athens, 4 Soranou Efesiou Street, 11527 Athens, Greece</footnote>

<footnote n="6">⁶ School of Astronomy and Space Science, Nanjing University, Nanjing 210023, China</footnote>

<footnote n="7">⁷ Stanford University HEPL Annex, B210, Stanford, CA 94305-4085, USA</footnote>

<footnote n="8">⁸ Institute of Physics/IGAM, University of Graz, Universitätsplatz 5/II, 8010 Graz, Austria</footnote>

<footnote n="9">⁹ Key Laboratory of Solar Activity, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China</footnote>
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 pre-defined three-dimensional magnetic-field solution.

**Keywords** Magnetic fields · Methods: numerical · Sun: surface magnetism · Sun: corona

# 1 Introduction

The volume integral

$$ \mathcal{H}(t) = \int_{\mathcal{V}} (\mathbf{A} \cdot \mathbf{B}) \, d\mathcal{V} \qquad (1) $$

is the helicity of the vector field $\mathbf{B} = \nabla \times \mathbf{A}$ in a given volume $\mathcal{V}$, with $\mathbf{A}(\mathbf{x}, t)$ representing the corresponding space- and time-dependent vector potential. If the field consists of a discrete collection of flux tubes, $\mathcal{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 $\mathbf{A}$, the addition of the gradient of a (sufficiently regular but otherwise arbitrary) scalar function, i.e., the transformation $\mathbf{A} \rightarrow \mathbf{A} + \nabla \psi$, does not change the resulting $\mathbf{B}$. This property of the definition of $\mathbf{B}$ is called gauge-invariance. Due to this freedom in the gauge, $\mathcal{H}$ is not uniquely defined, since

$$ \mathcal{H}(\mathbf{A} + \nabla \psi) = \mathcal{H}(\mathbf{A}) + \int_{\partial \mathcal{V}} (\psi \mathbf{B}) \cdot d\mathbf{S} - \int_{\mathcal{V}} \psi (\nabla \cdot \mathbf{B}) \, d\mathcal{V}, \qquad (2) $$

where $d\mathbf{S} = d\mathbf{S}\hat{\mathbf{n}}$, with $d\mathbf{S}$ being the infinitesimal element of the bounding surface $\partial\mathcal{V}$ of the volume $\mathcal{V}$, and $\hat{\mathbf{n}}$ the outward-oriented normal to $\partial\mathcal{V}$. Hence, $\mathcal{H}$ is not gauge-invariant unless two conditions are met: First, the vector field $\mathbf{B}$ must be solenoidal, as implied by its definition as curl of $\mathbf{A}$, and, second, the volume’s bounding surface $\partial\mathcal{V}$ must be a flux surface of $\mathbf{B}$, i.e., $(\mathbf{B} \cdot \hat{\mathbf{n}})|_{\partial\mathcal{V}} = 0$. When applied to a magnetic field $\mathbf{B}$, the solenoidal requirement is satisfied by virtue of Maxwell equations, although possibly only to a finite extent in numerical experiments, and $\partial\mathcal{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 magneto-hydrodynamics (MHD, Boltzmann 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 non-ideal 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 time-evolution to helicity-conserving paths in phase-space, which, for instance, yielded the so-called 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

$$ \mathcal{H}_V \equiv \int_V (\mathbf{A} + \mathbf{A}_p) \cdot (\mathbf{B} - \mathbf{B}_p) \, dV, \qquad (3) $$

of a magnetic field $\mathbf{B}$ with respect to a reference magnetic field $\mathbf{B}_p = \nabla \times \mathbf{A}_p$. Even though Eq. (3) allows for an arbitrary reference field, here we adopt the usual choice of the electric current-free (potential) field for $\mathbf{B}_p$. This choice has the following motivation: in order for $\mathcal{H}_V$ in Eq. (3) to be gauge invariant, the input and potential fields, $\mathbf{B}$ and $\mathbf{B}_p$, respectively, must be solenoidal, and such that

$$ (\hat{\mathbf{n}} \cdot \mathbf{B})|_{\partial V} = (\hat{\mathbf{n}} \cdot \mathbf{B}_p)|_{\partial V}. \qquad (4) $$

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, $\partial V$ needs not to be a flux surface for $\mathcal{H}_V$ to be gauge-invariant, and the definition Eq. (3) can be applied to arbitrary field distributions.

Berger (2003) derived a useful decomposition of Eq. (3) as

$$ \mathcal{H}_V = \mathcal{H}_{V,J} + \mathcal{H}_{V,JP} \qquad (5) $$

where

$$ \mathcal{H}_{V,J} \equiv \int_V (\mathbf{A} - \mathbf{A}_p) \cdot (\mathbf{B} - \mathbf{B}_p) \, dV, \qquad (6) $$

$$ \mathcal{H}_{V,JP} \equiv 2 \int_V \mathbf{A}_p \cdot (\mathbf{B} - \mathbf{B}_p) \, dV. \qquad (7) $$

The two definitions in Eqs. (6), (7) have the property of being *separately* gauge-invariant under the same assumptions guaranteeing the gauge-invariance of $\mathcal{H}_V$. The first term $\mathcal{H}_{V,J}$ corresponds to the general definition of helicity equation (1), but this time of the current-carrying part of the field only, $(\mathbf{B} - \mathbf{B}_p) = \nabla \times (\mathbf{A} - \mathbf{A}_p)$. By construction, such a field has no normal component on the boundary, i.e., has $\partial V$ as a flux surface. The second term $\mathcal{H}_{V,JP}$ has no intuitive interpretation, but is a sort of mutual helicity that basically takes care of the flux threading $\partial V$ (via the transverse component of $\mathbf{A}_p$), and is gauge-invariant because only the current-carrying part of $\mathbf{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, higher-density layers
of the atmosphere, yielding basically two-dimensional 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 sub-volumes, 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)

- twist-number (TN)

- helicity-flux integration (FI)

- connectivity-based (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 force-free 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 semi-infinite 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><caption>**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</caption><thead><tr><th>Finite volume (FV)</th><th>Helicity-flux integration (FI)</th></tr></thead><tbody><tr><td><i>H</i><sub>V</sub> = ∫<sub>V</sub>(<b>A</b> + <b>A</b><sub>p</sub>)·(<b>B</b> − <b>B</b><sub>p</sub>) d<i>V</i><br>see Eq. (3)</td><td><i>dH</i><sub>V</sub> / dt = 2 ∫<sub>∂<i>V</i></sub>[(<b>A</b><sub>p</sub> · <b>B</b>)<i>v</i><sub>n</sub> − (<b>A</b><sub>p</sub> · <b>v</b><sub>t</sub>)<b>B</b><sub>n</sub>] d<i>S</i></td></tr><tr><td><ul><li>Requires <b>B</b> in <i>V</i> e.g., from MHD simulations or NLFFF</li><li>Compute <i>H</i><sub>V</sub> at one time</li><li>May employ different gauges (see Table 2)</li></ul></td><td><ul><li>Requires time evolution of vector field on ∂<i>V</i></li><li>Requires knowledge or model of flows on ∂<i>V</i></li><li>Valid for a specific set of gauge and assumptions, see Pariat et al. (2017)</li></ul></td></tr><tr><td><b>Discrete flux-tubes (DT)</b></td><td><b>Connectivity-based (CB)</b></td></tr><tr><td><i>H</i> ≈ ∑<sub>i=1</sub><sup>M</sup> <i>T</i><sub>i</sub> <i>Φ</i><sub>i</sub><sup>2</sup> + ∑<sub>i=1</sub><sup>M</sup> ∑<sub>j=1, j≠i</sub><sup>M</sup> <i>L</i><sub>i,j</sub> <i>Φ</i><sub>i</sub> <i>Φ</i><sub>j</sub>,<br>see Eq. (31)</td><td><i>H</i> = A ∑<sub>i=1</sub><sup>M</sup> α<sub>i</sub> <i>Φ</i><sub>i</sub><sup>2δ</sup> + ∑<sub>i,m=1</sub><sup>M</sup> <i>L</i><sub>lm</sub> <i>Φ</i><sub>i</sub> <i>Φ</i><sub>m</sub><br>see Eq. (35)</td></tr><tr><td><b>Twist-number (TN)</b></td><td><b>Connectivity-based (CB)</b></td></tr><tr><td><i>H</i> ≈ <i>T</i> <i>Φ</i><sup>2</sup><br>see Eq. (32)</td><td><i>H</i> = A ∑<sub>i=1</sub><sup>M</sup> α<sub>i</sub> <i>Φ</i><sub>i</sub><sup>2δ</sup> + ∑<sub>i,m=1</sub><sup>M</sup> <i>L</i><sub>lm</sub> <i>Φ</i><sub>i</sub> <i>Φ</i><sub>m</sub><br>see Eq. (35)</td></tr><tr><td><ul><li>Estimation of the twist contribution to <i>H</i></li><li>Requires <b>B</b> in <i>V</i></li><li>Requires a flux-rope-like structure for computing the twist <i>T</i></li></ul></td><td><ul><li>Requires the vector field on photosphere at one time</li><li>Models the corona connectivity as a collection of <i>M</i> force-free flux tubes</li><li>Minimal connection length principle</li></ul></td></tr></tbody></table>

The field-line 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 field-line 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 above-mentioned article.

Similarly to FV methods, the *twist-number (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).

**Helicity-flux 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 connectivity-based (CB) method (Georgoulis et al. 2012). The method is based on modeling the unknown connectivity of the coronal field with a collection of slender force-free flux tubes, each with different constant force-free parameter $\alpha \equiv (\nabla \times \mathbf{B})/\mathbf{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 flare-productive active regions. In this way, the connectivity-based 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 force-free 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 lower-limit estimation of the true helicity associated to a flux-balanced coronal field in a very fast way. In this sense, the CB method does not share the same purpose of finite volume and helicity-flux 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 finite-sized flux tubes, as opposite to a continuous three-dimensional field. We therefore categorized both methods as discrete flux-tubes (DT) methods, see Table 1.

## 1.3 Systematic Comparison of Methods

Section 1.2 briefly introduced the methods currently available for the estimation of $H_v$. 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 force-free equations, proceeding then to time-dependent 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 ISSI-Bern in the 2014–2016 period.<fnref n="1" /> As a reference for future testing, we make available the data cube of each

<footnote n="1">¹See the web page http://www.issibern.ch/teams/magnetichelicity/ for more information.</footnote>
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 $H_V$ 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 strongly-controlled-
environment, equilibrium test-cases such as the Low and Lou (1990) and Titov and Dé-
moulin (1999) solutions of the force-free 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
twist-number (TN) and the connectivity-based (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 informa-
tion 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 con-
nectivity. Also, similarly to FV methods, the connectivity-based 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 mag-
netic 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><caption>**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</caption><thead><tr><th>Method</th><th>Label</th><th>Category</th><th>Section</th><th>Reference</th></tr></thead><tbody><tr><td>Coulomb-Thalmann</td><td>Coulomb_JT</td><td>Finite volume</td><td>Sect. 2.1.1</td><td>Thalmann et al. (2011)</td></tr><tr><td>Coulomb-Yang</td><td>Coulomb_SY</td><td>Finite volume</td><td>Sect. 2.1.2</td><td>Yang et al. (2013b)</td></tr><tr><td>Coulomb-Rudenko</td><td>Coulomb_GR</td><td>Finite volume</td><td>Sect. 2.1.3</td><td>Rudenko and Anfinogentov (2014)</td></tr><tr><td>DeVore-Valori</td><td>DeVore_GV</td><td>Finite volume</td><td>Sect. 2.2.1</td><td>Valori et al. (2012)</td></tr><tr><td>DeVore-Moraitis</td><td>DeVore_KM</td><td>Finite volume</td><td>Sect. 2.2.2</td><td>Moraitis et al. (2014)</td></tr><tr><td>DeVore-Anfinogentov</td><td>DeVore_SA</td><td>Finite volume</td><td>Sect. 2.2.3</td><td>Not available</td></tr><tr><td>Twist-number</td><td>TN</td><td>Discrete flux-tubes</td><td>Sect. 2.3.1</td><td>Guo et al. (2010)</td></tr><tr><td>Connectivity-based</td><td>CB</td><td>Discrete flux-tubes</td><td>Sect. 2.3.2</td><td>Georgoulis et al. (2012)</td></tr></tbody></table>

## 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 ($\nabla \cdot \mathbf{A} = 0$) or the DeVore gauge ($A_z = 0$, DeVore 2000). Due to the gauge-invariant 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 current-carrying 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 $\mathbf{B}_p = \nabla\phi$, with $\phi$ being the scalar potential, solution of

$$ \nabla^2 \phi = 0, \qquad (8) $$

$$ (\hat{\mathbf{n}} \cdot \nabla \phi)|_{\partial \mathcal{V}} = (\hat{\mathbf{n}} \cdot \mathbf{B})|_{\partial \mathcal{V}}, \qquad (9) $$

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

$$ \nabla^2 \mathbf{A}_p = 0, \qquad (10) $$

$$ \nabla \cdot \mathbf{A}_p = 0, \qquad (11) $$

$$ \hat{\mathbf{n}} \cdot (\nabla \times \mathbf{A}_p)|_{\partial \mathcal{V}} = (\hat{\mathbf{n}} \cdot \mathbf{B})|_{\partial \mathcal{V}}, \qquad (12) $$
for the vector potential of the potential field, and

$$ \nabla^2 \mathbf{A} = -\mathbf{J}, \qquad (13) $$

$$ \nabla \cdot \mathbf{A} = 0, \qquad (14) $$

$$ \hat{\mathbf{n}} \cdot (\nabla \times \mathbf{A})|_{\partial V} = (\hat{\mathbf{n}} \cdot \mathbf{B})|_{\partial V}, \qquad (15) $$

for the vector potential of the input field, where $\mathbf{J} = \nabla \times \mathbf{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 $\mathbf{A}_p$ and $\mathbf{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 $\mathbf{A}_p$ and $\mathbf{A}$.

### 2.1.1 Coulomb_JT

In order to solve Eqs. (10)–(12) and (13)–(15), appropriate additional boundary conditions for $\mathbf{A}$ and $\mathbf{A}_p$ on the boundaries of the 3D-rectangular computational domain need to be specified. For this purpose, the method of Thalmann et al. (2011), decomposes $\mathbf{A}$ into a current-carrying and a potential (current-free) part, in the form $\mathbf{A} = \mathbf{A}_c + \mathbf{A}_p$. The reproduction of the input fields’ flux at the volume’s boundaries, $\partial V$, is entirely dedicated to $\mathbf{A}_p$ (obeying Eqs. (10)–(12)). The electric current distribution in $V$ and on $\partial V$, on the other hand, are delivered by $\mathbf{A}_c$ (Eqs. (13), (14) with $\mathbf{A}_c$ instead of $\mathbf{A}$ and Eq. (15) replaced by $\hat{\mathbf{n}} \cdot (\nabla \times \mathbf{A}_c)|_{\partial V} = 0$).

In particular, the tangential components of $\mathbf{A}_p$ on a particular face, $f$, of the 3D computational domain are specified to be the 2D stream function of a corresponding Laplacian field, $\phi_f$, in the form $\mathbf{A}_{p,t} = -\hat{\mathbf{n}} \times \nabla_t \phi_f$, where $\nabla_t$ is the 2D-gradient 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 $\nabla_t^2 \phi_f = -\hat{\mathbf{n}} \cdot \mathbf{B}$, for which boundary conditions on the four edges of each face $f$ need to be specified. This approach of defining $\mathbf{A}_p$ on $\partial 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 $\partial_n \phi_f = c_f$, where $c_f$ is a constant along a particular face and $\partial_n$ is the derivative in the direction normal to the edge of the face $f$ (see their Sect. 2.1 for details). The different $c_f$ are constructed in such a way that the total outflow through the volume’s bounding surface $\partial V$ is minimized. In this way, a vanishing tangential divergence ($\nabla_t \cdot \mathbf{A}_{p,t} = 0$) is enforced on $\partial V$, and following Gauss’ theorem, Eqs. (10)–(12) are approximately fulfilled.

Another difference of the applied Coulomb methods is how the current-carrying vector potential $\mathbf{A}_c$ is calculated and its solenoidality enforced. Thalmann et al. (2011) solves Eq. (13) for $\mathbf{A}_c$ numerically, similar to Yang et al. (2013a), just with differing boundary conditions. In the Coulomb_JT case, $\nabla_t \cdot \mathbf{A}_{c,t} = 0$ on $\partial V$ is explicitly enforced in order to fulfill the Coulomb gauge for $\mathbf{A}_c$.

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 coordinates 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 $(\hat{n} \cdot A_p)|_{\partial V} = (\hat{n} \cdot A)|_{\partial 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, $A'_p$ and $A'$. These preliminary solutions satisfies Eqs. (12) and (15), but not the Coulomb gauge condition. The Coulomb_SY method then uses a divergence-cleaning 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 $B_p = \nabla \times A_p$, rather than using Eqs. (8), (9). The corresponding boundary value problem, Eq. (10), is solved with the constraint equation (11) and the boundary condition $(A_p \cdot \hat{n})|_{\partial V} = 0$. The Laplace problem is divided into six sub-problems, 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 $B_p^m = \nabla \times A_p^m$ 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 - B_p^m$ has zero total flux on each face independently and can be correctly used as a boundary condition for the sub-problems

$$ \hat{n} \cdot (\nabla \times A_p^{fi})|_{f_j} = \delta_{ij} (\hat{n} \cdot B')|_{f_j}, \quad (16) $$

where $A_p^{fi}$ is the vector potential of sub-problem solution corresponding to the side $f_i$.

After solving all sub-problems, the full solution is then obtained by summation of the solutions of the six sub-problems $A_p^{fi}$ and the vector potential of a compensation field, $A^m$, as

$$ A_p = A_p^m + \sum_{i=1}^{6} A_p^{fi}. \quad (17) $$

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 =$
$A_c + A_p$, but the vector potential of the current-carrying part of the field is computed as

$$A_c(\mathbf{r}) = -\frac{1}{4\pi} \int_V \frac{\mathbf{r} - \mathbf{r}'}{|\mathbf{r} - \mathbf{r}'|^3} \times (\mathbf{B} - \mathbf{B}_p) \, dV. \quad (18)$$

In contrast to other methods, the solutions of the Laplace and Poisson problems for the $A^fi$ 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 $A_z = 0$ (DeVore 2000), Valori et al. (2012) derived the expression for the vector potential of the magnetic field $\mathbf{B}$ in the finite volume $V = [x_1, x_2] \times [y_1, y_2] \times [z_1, z_2]$ as

$$\mathbf{A} = \mathbf{b} + \hat{\mathbf{z}} \times \int_z^{z_2} \mathbf{B} \, d\mathbf{z}', \quad (19)$$

where the integration function $b(x, y) = A(z = z_2)$ obeys to

$$\partial_x b_y - \partial_y b_x = B_z(z = z_2), \quad (20)$$

and $b_z = 0$. The particular solution of Eq. (20) employed here is

$$b_x = -\frac{1}{2} \int_{y_1}^{y} B_z(x, y', z = z_2) \, dy', \quad (21)$$

$$b_y = \frac{1}{2} \int_{x_1}^{x} B_z(x', y, z = z_2) \, dx', \quad (22)$$

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 $\mathbf{B}$ with $\mathbf{B}_p$ everywhere in Eqs. (19)–(22). In particular, using Eqs. (21), (22) for both $A_p$ and $A$ implies $A_p = A$ at $z = z_2$, 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 $A_z = 0$, then

$$B_x = -\partial_z A_y = \partial_z \int_z^{z_2} B_x \, dz', \quad (23)$$

where Eqs. (20) and (21), (22) where used; a similar expression holds for $B_y$. Hence, the accuracy of DeVore method in reproducing $B_x$ and $B_y$ from $A$ of Eq. (19) depends only on how accurately the relation

$$\partial_z \int_z^{z_2} = \text{identity} \quad (24)$$

is verified numerically. On the other hand, even when Eq. (24) is obeyed to acceptable accuracy, one can easily show that, for a non-perfectly solenoidal field $\mathbf{B}_{ns}$, it is

$$\mathbf{B}_{ns} - \nabla \times \mathbf{A} = \hat{\mathbf{z}} \int_z^{z_2} (\nabla \cdot \mathbf{B}_{ns}) \, dz', \quad (25)$$
as derived in Eq. (B.4) of Valori et al. (2012). Hence, the accuracy in the reproduction of the z-component of the field depends on the solenoidal level of the input field (and on how accurately Eq. (20) is solved).

All DeVore-gauge 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 z-integral 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 $\phi$ 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 DeVore-Coulomb 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 $\le 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) second-order 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 $\phi$ are solved in Fourier space separately for all faces of the box. In particular, the problem is divided into six sub-problems using

$$ \phi = \phi^c + \sum_{i=1}^{6} \phi_i, \qquad (26) $$

where $\phi^c$ is the 3D scalar potential of the compensation field $\mathbf{B}^c = \nabla \phi_c$ and $\phi_i$ are 3D solutions for the potential field with the normal component given on $i$st side of $V$ and vanishing
boundary conditions on the other sides of $V$. The individual Laplace problems for each $\phi_i$ are then solved in Fourier space following the general scheme of the potential and linear force-free 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 $B_z$ from any level $z_r$ inside the data cube. In particular, by addition and subtraction to Eq. (19), one has

$$A = b + \hat{z} \times \left( \int_{z}^{z_2} B \, dz' + \int_{z_r}^{z_2} B \, dz' - \int_{z_r}^{z_2} B \, dz' \right), \quad (27)$$

which can be re-casted as

$$A = c + \hat{z} \times \left( \int_{z}^{z_2} B \, dz' - \int_{z_r}^{z_2} B \, dz' \right), \quad (28)$$

where we have defined

$$c = b + \hat{z} \times \int_{z_r}^{z_2} B \, dz'. \quad (29)$$

Taking the $x$- and $y$-derivatives of Eq. (29), and using Eq. (20) and $\nabla \cdot B = 0$, one derives

$$\partial_x c_y - \partial_y c_x = B_z(z = z_r). \quad (30)$$

The solution of Eq. (30) is then analogous to Eqs. (21), (22), where $B_z(z = z_2)$ is replaced by $B_z(z = z_r)$. Tests using the LL case of Sect. 4.1 shown that the minimal error in $A_p$ and $A$ is obtained for $z_r = (z_2 - z_1)/2$. The vector potential is finally computed following Eq. (28).

This scheme coincides with the original one of DeVore_GV if $z_r$ is taken at the top boundary of the box, i.e., for $c(z_r = z2) = b$.

## 2.3 Discrete Flux-Tubes Methods

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:

$$\mathcal{H} \approx \sum_{i=1}^{M} \mathcal{T}_i \Phi_i^2 + \sum_{i=1}^{M} \sum_{j=1, j \neq i}^{M} \mathcal{L}_{i,j} \Phi_i \Phi_j, \quad (31)$$

where $\mathcal{T}_i$ denotes the twist and writhe of magnetic flux tube $i$ with flux $\Phi_i$, and $\mathcal{L}_{i,j}$ is the linking number between two magnetic flux tubes $i$ and $j$ with fluxes $\Phi_i$ and $\Phi_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 flux-tubes methods do not provide the vector potentials and potential fields in the considered volumes like FV methods. Therefore, the comparison
between DT methods and FV methods is necessarily restricted to the helicity values only. The twist-number method and connectivity-based 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 Twist-Number

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

$$ \mathcal{H}_{\text{twist}} \approx \mathcal{T}\Phi^2, \qquad (32) $$

where $\mathcal{T}$ is the twist number of the considered magnetic flux rope with flux $\Phi$. In order to estimate $\mathcal{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 quasi-separatrix layer (QSL) that is found to wrap the flux rope (Guo et al. 2013). Then the twist density, $d\mathcal{T}/ds$, at an arc length $s$ is:

$$ \frac{d\mathcal{T}}{ds} = \frac{1}{2\pi} \mathbf{T} \cdot \mathbf{V} \times \frac{d\mathbf{V}}{ds}. \qquad (33) $$

Two unit vectors are used in Eq. (33): $\mathbf{T}(s)$, that is tangent to the axis curve, and $\mathbf{V}(s)$, that is normal to $\mathbf{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 force-free and non-force-free magnetic field models.

### 2.3.2 Connectivity-Based

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 self-consistently 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 opposite-polarity 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 field-line tracing associates connected flux with photospheric partitions, providing the magnetic connectivity matrix upon summation of individual field-line contributions. Obviously, only magnetic field lines entirely embedded in the finite volume are taken into account. In the second case, a simulated-annealing method is used to absolutely and simultaneously minimize the flux imbalance (hence achieving connections between opposite-polarity 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 simulated-annealing 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 $\Phi_i$; $i = \{1, ..., M\}$. The respective force-free parameters $\alpha_i$ are assumed constant for a given flux tube but vary between different tubes, thus implementing the nonlinear force-free (NLFF) field approximation. Force-free parameters for each flux tube are the mean values of the force-free parameters of the tubes' respective footprints, each calculated by the relation $\alpha_i = \frac{4\pi I_i}{c F_i}$; $i = \{1, ..., M'\}$ for $M'$ magnetic partitions, where $I_i$ is the total electric current of the $i$-partition and $F_i$ 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 $F_i$, $\alpha_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

$$E_{cCB} = E_{c(CB;self)} + E_{c(CB;mutual)} = A\lambda^2 \sum_{i=1}^{M} \alpha_i^2 \Phi_i^{2\delta} + \frac{1}{8\pi} \sum_{l=1}^{M} \sum_{m=1; l \neq m}^{M} \alpha_l \mathcal{L}_{lm} \Phi_l \Phi_m, \quad (34)$$

where A, $\delta$ are known fitting constants, $\lambda$ is the length element (the pixel size in observed photospheric magnetograms), and $\mathcal{L}_{lm}$ is the mutual-helicity 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 flux-tube footpoints. The locations of point-like footpoints of the slender flux tubes coincide with the flux-weighted centroids of the respective flux partitions. As included in Eq. (34), the parameter $\mathcal{L}_{lm}$ 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 $E_{cCB}$ in Eq. (34). Therefore, the above $E_{cCB}$ is already a lower limit of the actual $E_c$, assuming only “arch-like” (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 free-energy term that is due to the generation, caused by induction, of potential flux tubes around the collection of non-potential ones (Démoulin et al. 2006). Such a term would again contribute to the mutual term of the free energy.

The corresponding self-consistent relative helicity is, then,

$$H_{mCB} = H_{m(CB;self)} + H_{m(CB;mutual)} = 8\pi A\lambda^2 \sum_{i=1}^{M} \alpha_i \Phi_i^{2\delta} + \sum_{l=1}^{M} \sum_{m=1; l \neq m}^{M} \mathcal{L}_{lm} \Phi_l \Phi_m. \quad (35)$$
From Eqs. (34), (35) we identically have $E_{CCB} \equiv 0$ for potential flux tubes ($\alpha_i = 0; i \equiv \{1, \dots, M\}$). For $H_{mCB} = 0$ in this case, we further require $\sum_{l=1}^{M_p} \sum_{m=1; l \neq m}^{M_p} \mathcal{L}_{lm} \Phi_l \Phi_m = 0$, where $\mathcal{L}_{lm} \neq \mathcal{L}_{lm}$ is the mutual-helicity factor for a collection of $M_p \neq M$ collection of potential flux tubes. As Démoulin et al. (2006) discuss, this can be the case for a flux-balanced potential-field boundary condition. In practical situations of not-precisely flux-balanced magnetic configurations, however, one may approximate $H_{mCB; mutual} = 0$, in case all $\alpha_i; i \equiv \{1, \dots, M\}$ are zero within uncertainties $\delta\alpha_i$, which are fully defined in this analysis. More generally, one may use the “energy-helicity diagram” correlation of Tziotziou et al. (2012, 2014) to infer $|H_{mCB}| \propto E_{CCB}^{0.84 \pm 0.05}$ for $E_{CCB} \to 0$.

### 3 Analysis Metrics

Apart from extremely simplistic magnetic fields, the analytical computation of the relative magnetic helicity in a non-magnetically bounded system is highly non-trivial. Even with simple natural-world-relevant models relative magnetic helicity cannot be analytically estimated. Similarly, the exact value of $H_V$ 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 $\mathbf{B}$ defined by Eq. (3) involves the corresponding potential field $\mathbf{B}_p$, as well as the vector potentials for $\mathbf{B}_p$ and $\mathbf{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 $H_V$ as the helicity defined in Eq. (3) normalized to $\Phi^2$,

$$H_V \equiv \mathcal{H}_V / \Phi^2, \qquad (36)$$

where

$$\Phi(\mathbf{B}) = \frac{1}{2} \int_{z=0} |B_z(x', y', z=0)| \, dx \, dy \qquad (37)$$

is half of the unsigned flux through the bottom boundary, corresponding to the injected flux for an exactly flux-balanced 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 $H_V$ as defined by Eq. (36) is obtained by computing separately the four volume contributions of Eq. (3) and normalizing them to $\Phi^2$. We reserve the calligraphic symbol $\mathcal{H}$ for non-normalized helicities.

The numerically obtained $H_V$ 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 $H_V$ 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 research-relevant dataset, all employed tests are defined on discretized grids of moderate- to high-resolution, 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 $H_V$ 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 $H_V$ values, and a relative spread around it, defined as the standard deviation of the $H_V$ values distribution over the mean. In addition, in order to discern among different factors influencing $H_V$ 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

$$ \epsilon_N = 1 - \frac{\sum_i |\mathbf{X}_i - \mathbf{Y}_i|}{\sum_i |\mathbf{X}_i|}, \qquad (38) $$

$$ \epsilon_E = \frac{\sum_i |\mathbf{Y}_i|^2}{\sum_i |\mathbf{X}_i|^2}, \qquad (39) $$

which are, respectively, the complement of the normalized vector error and the energy ratio, introduced by Schrijver et al. (2006).<fnref n="2" /> Both are unity if $\mathbf{X}_i = \mathbf{Y}_i$ in all grid points $i$ in $V$. The metrics are applied to the pair ($\mathbf{X} = \mathbf{B}_p, \mathbf{Y} = \nabla \times \mathbf{A}_p$) or ($\mathbf{X} = \mathbf{B}, \mathbf{Y} = \nabla \times \mathbf{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 second-order, central-difference discretization scheme for points in the interior of $V$. Values on the volume-bounding surface, $\partial 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 $\nabla \cdot \mathbf{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 $E_{\text{div}}$ of the

<footnote n="2">²We use the notation $\epsilon_N$ rather than $1 - E_N$ to avoid confusion with energy symbols.</footnote>
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, $\langle |f_i(\mathbf{B})| \rangle = (\sum_i |f_i|)/n$, of the fractional flux,

$$f_i(\mathbf{B}) = \frac{\int_{\nu} \nabla \cdot \mathbf{B}_i \, dv}{\int_{\partial \nu} |\mathbf{B}_i| \, dS} \approx \frac{\nabla \cdot \mathbf{B}_i}{6|\mathbf{B}_i|/\Delta}, \quad (40)$$

through the surface $\partial 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 $\Delta$. 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 $\langle |f_i| \rangle$, the more solenoidal the field. However, the actual value of this metric depends on the number of grid points $n$ and the resolution $\Delta$, 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, $\langle |f_i(\mathbf{A})| \rangle$, 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 magneto-static and magneto-hydrodynamic 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-$\beta$ Grad-Shafranov 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, \phi = \pi/4$ in the notation of Low and Lou 1990), in the same $V = [-1, 1] \times [-1, 1] \times [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 large-scale, force-free field with large-scale smooth currents distributed
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Representative field lines of the four employed test cases: (a) low and Lou (LL) force-free equilibrium, see Sect. 4.1; (b) Titov and Démoulin (TD) force-free equilibrium for $N = 1$ and $\Delta = 0.06$, see Sect. 4.2; (c) snapshot at $t = 155$ of the stable MHD simulation (MHD-st), and (d) snapshot at $t = 140$ of the unstable MHD simulation (MHD-un), 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 semi-transparent iso-contour of the current density is shown at 15%, 35%, 11%, and 10% of its maximum in the four cases (a) to (d), respectively</figcaption></figure>

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 $E_{free} \simeq 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 $H_V$ 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 force-free 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$ end-to-end turns, bald patches, and a 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
<figure><img src="image_5.png" /><figcaption>(c) TD resolution End-to-end twist</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(e) MHD-un</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) TD twist</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) MHD-st</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) MHD-st-div(B)</figcaption></figure>

<figure><img src="image_10.png" /></figure>

**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) MHD-st; (e) MHD-un (f) divergence test MHD-st-div(B), as a function of $δ$, see Sect. 7
<table><caption>**Table 3** Parameters of the TD twist test cases, with resolution $\Delta = 0.06$; $N$ is the approximate end-to end number of turns, $a$ and $d$ are the minor radius and the depth of the center of the current ring, $q$ is the strength of the magnetic charges generating the confining field. The $N = 1$ case here is the same as the $\Delta = 0.06$ case of Table 4. We refer to Valori et al. (2010) for definitions and normalizations of the TD parameters</caption><thead><tr><th>$N$</th><th>$a/d$</th><th>$q \times 10^{12}$</th><th>Twist/$\pi$</th></tr></thead><tbody><tr><td>3</td><td>0.31</td><td>100</td><td>-6.054</td></tr><tr><td>1</td><td>0.80</td><td>100</td><td>-2.114</td></tr><tr><td>0.5</td><td>0.80</td><td>29.5</td><td>-1.004</td></tr><tr><td>0.1</td><td>0.80</td><td>5.5</td><td>-0.201</td></tr></tbody></table>

<table><caption>**Table 4** Parameters of the TD resolution test cases, with $N \simeq 1$; $\Delta$ is the spatial resolution ($n_x, n_y, n_z$) are the number of nodes in the x-, y- and z-directions, respectively. The case $\Delta = 0.06$ here is the same as the $N = 1$ case in Table 3. We refer to Valori et al. (2010) for definitions and normalizations of the TD parameters</caption><thead><tr><th>$\Delta$</th><th>$n_x$</th><th>$n_y$</th><th>$n_z$</th><th>Twist/$\pi$</th></tr></thead><tbody><tr><td>0.03</td><td>209</td><td>337</td><td>153</td><td>-2.15</td></tr><tr><td>0.06</td><td>107</td><td>171</td><td>79</td><td>-2.14</td></tr><tr><td>0.12</td><td>56</td><td>88</td><td>42</td><td>-2.09</td></tr></tbody></table>

stated, in the following we refer to twist of the TD cases as the average twist over the current ring's cross section at the flux rope apex, as computed, e.g., in Sect. 4 of Valori et al. (2010).

Since the equations in Titov and Démoulin (1999) are given in implicit form, the twist of each equilibria is obtained by trial and error. In all six TD cases considered here, the discretized volume is $[-3.18, 3.18] \times [-5.10, 5.10] \times [0.00, 4.56]$, the distance between the charges and the current ring center is $L = 0.83$, the depth of the current ring center is $d = 0.83$, and the ring's radius is $R = 1.83$. We refer to Valori et al. (2010) for definitions and normalizations of the TD parameters, where the Low\_HFT case has the same parameters as the $N = 1$ case here. Based on previous tests (see e.g., Török et al. 2004; Kliem and Török 2006), all cases presented here are stable equilibria, except for the $N = 3$ one, which is almost certainly kink-unstable.

**TD-Twist Test Cases** We consider equilibria defined by parameters given in Table 3 resulting into flux ropes with different average twist. The resolution (uniform pixel size) is $\Delta = 0.06$ for all four cases, i.e., given by a grid of $107 \times 171 \times 79$ nodes. The energy decomposition of the TD twist cases are shown in Fig. 2b. The solenoidal errors vary between 0.1 % and 2 %, monotonically increasing with twist except for the $N = 3$ case, which corresponds to a much thinner flux tube that satisfies the local cylindrical approximation better. Their contribution in relative energy is always one order of magnitude smaller than the free energy. The only exception is the $N \simeq 0.05$ case where the twist is lower and the field is almost potential ($E_{\text{free}} \simeq 0.1$ %). In this case solenoidal errors are slightly larger than the free energy, but anyway extremely small ($E_{\text{div}} = 0.2$ %).

**TD-Resolution Test Cases** The three equilibria for this test are the N1 case in Table 4 at resolution $\Delta = 0.06$, and two additional equilibria with exactly the same parameters but with $\Delta = 0.03$ and 0.12, respectively. The energy decomposition of the TD resolution tests summarized in Fig. 2c shows that the free energy is independent of resolution (at about 2 % of the total energy), and that the solenoidal errors depend only very weakly on it (cf. Table 7).
While the former is expected, the latter confirms that solenoidal errors in the (non-relaxed) numerical implementations of TD equilibria have a stronger dependence on twist than on resolution that is mostly generated by the match at the interface between flux rope and ambient potential field.

In the construction of TD equilibria, the local matching between current ring and poten-
tial environment field at the flux rope’s boundary is done in locally cylindrical coordinates,
hence some spurious Lorentz forces and solenoidal errors are present at the interface be-
tween the two flux systems. For such reasons, when employed as an initial state of numeri-
cal simulations, an MHD relaxation is normally applied beforehand to remove such residual
forces and errors, see e.g., Török and Kliem (2003). For the purpose of computing the rel-
ative magnetic helicity, however, the relaxation step is unnecessary as long as the errors in
the solenoidal property of the field are small enough, which is the case here.

## 4.3 MHD-Emergence Stable and MHD-Emergence Unstable

The MHD cases are simulations of stable (hereafter, MHD-st) and unstable (hereafter, MHD-un) evolutions of flux emergence published in by Leake et al. (2013) and (2014), respectively. These simulations offer the possibility of studying the time evolution of the helicity as the flux breaks through the photospheric layers, slowly accumulates in the coronal volume, and either reaches stability or erupts.

More in detail, the MHD-st case is a sub-domain of the Strong Dipole case by Leake
et al. (2013), which was obtained using a stretched grid. Here we consider five snap-
shots of that simulation, and for each one the magnetic field is interpolated using a grid
of 233 × 233 × 174 nodes with a uniform mesh size Δ = 0.86 discretizing the volume
[−100, 100] × [−100, 100] × [0, 149] (in units of L₀, see Leake et al. 2013 for more de-
tails). The interpolation is expected to introduce spurious solenoidal errors, but is required
to accommodate for the requirements of all helicity computation methods. Compared to
the original data, the bottom boundary of the interpolated mesh corresponds to the photo-
spheric level, the top boundary to z = 150, and the domain is reduced in the lateral extension.
The time series includes snapshots at t = [50, 85, 120, 155, 190]. The MHD-un case, corre-
sponding to the Medium Dipole case of Leake et al. (2014), is prepared in a similar way,
but the time series is t = [50, 80, 110, 140, 150, 160, 190]. Leake et al. (2014) identify the
starting time of the eruption around t = 120 and the erupting structure is leaving the domain
between t = 150 and t = 160, see, e.g., their Fig. 11.

In the MHD simulation, the stable and unstable cases are obtained only by changing
strength and orientation of the coronal field, but nothing of the emerging flux rope. In this
sense, the two simulations are very similar, albeit resulting in a very different end state. Note
that in the MHD cases, *E*<sub>free</sub>, *E*<sub>p</sub>, and hence *E* are all dynamically evolving in time. In the
unstable case, the small variation of the ratio *E*<sub>free</sub>(*t*)/*E*(*t*) (of about 0.05 points in Table 7)
during the eruption, actually corresponds to a variation of 93 units of free energy with the
maximum total energy before the eruption being equal to 464 units.

Figure 2d, e show that, in both MHD-st and MHD-un cases, the solenoidal errors slightly
decrease during the MHD evolution from at most 1.6% at the beginning of the simulations,
to 0.5% at their end. Due to the interpolation, the actual values of the solenoidal errors
are larger than for the original simulations, but the trend in time is presumably the same. In
addition to the helicity estimation discussed in Sect. 4.3, the snapshot at *t* = 50 of MHD-st is
also used as the reference case for the test on the effect of nonzero divergence independently
of resolution, as described in Sects. 4.4 and 7.
<table><caption>**Table 5** Parameters $\delta$ and corresponding values of relative nonsolenoidal energy for the MHD-st-div($\mathbf{B}$) test, see Eq. (44)</caption><thead><tr><th>$\delta$</th><th>0.48</th><th>-0.45</th><th>-1.28</th><th>-2.87</th><th>-6.04</th><th>-16.00</th></tr></thead><tbody><tr><th>$E_{\text{div}}(\%)$</th><td>0.2</td><td>1.1</td><td>2.0</td><td>4.0</td><td>8.2</td><td>14.4</td></tr><tr><th>$E_{\text{ns}}(\%)$</th><td>-0.1</td><td>-1.0</td><td>-1.9</td><td>-3.6</td><td>-6.3</td><td>-8.3</td></tr></tbody></table>

With respect to TD and LL cases, the MHD-st and MHD-un tests confront the finite volume methods with a higher complexity in the field, with the formation of small scales, and with the coronal re-organization of the connectivity in time, see Fig. 1c, d. Moreover, these cases contain large currents and large free energies, of the order of 50 % to 60 % of the total magnetic energy, depending on the type and stage of the evolution, see Table 7 and Fig. 2d, e. Naturally, such cases are of interest for our discussion since they are supposed to mimic more realistically the difficulties that helicity estimation methods face in solar applications.

## 4.4 Parametric Nonsolenoidal Test Case

The application of Eq. (3) to a field with nonzero divergence makes simply no sense in terms of helicity because Eq. (3) becomes gauge-dependent. However, in the routine situation of numerical studies, a finite value of divergence is always present. The question arises, what is the value of divergence that is tolerable in computing helicity, i.e., that is producing a helicity value close to the one obtained for the solenoidal field? A complication of the problem is that the value of helicity is in general not known, not even for the test cases presented here. Hence, we are forced to reformulate our question in terms of variation of $H_\nu$ obtained by each method as a function of the increasing divergence. In practice, we check how each method behaves for increasing violation of the solenoidal property, but we are not in the position of stating which method gives the more “correct” value as divergence increases.

To this purpose, we designed a dedicated MHD-emergence stable div($\mathbf{B}$) test (hereafter, MHD-st-div($\mathbf{B}$)) by considering a numerically solenoidal field to which divergence is added in a controlled way and without changing the resolution, as done in Eqs. (14), (15) of Valori et al. (2013), to whom we refer for the details. In brief, starting from the snapshot of the magnetic field $\mathbf{B}$ at $t = 50$ of the MHD-st case, a solenoidal field $\mathbf{B}_s$ is produced by removing the divergence part $\mathbf{B}_{ns}$ of $\mathbf{B}$. A parametric, generally nonsolenoidal field,

$$ \mathbf{B}_{\delta} = \mathbf{B}_{s} + \delta \mathbf{B}_{ns}, \qquad (41) $$

is then constructed by adding the divergence part back, multiplied by a scalar amplitude $\delta$. In this way, the original spatial structure of $\mathbf{B}_{ns}$ is kept in $\mathbf{B}_{\delta}$ but its amplitude is modified according to the chosen value of $\delta$. The case $\delta = 0$ corresponds the numerically solenoidal field $\mathbf{B}_s$, whereas $\delta = 1$ reproduces the original test field $\mathbf{B}$, i.e., the snapshot of MHD-st at $t = 50$. The method for building the vector field $\mathbf{B}_{ns}$ is detailed in Sect. 7.1 of Valori et al. (2013).

By increasing the value of $\delta$, progressively more divergence can be added, yielding a larger amplitude of the nonsolenoidal component. A trial-and-error tuning of $\delta$ resulted into six MHD-st-div($\mathbf{B}$) test fields with nonsolenoidal contributions as a function of the parameter $\delta$ as summarized in Table 5 and Fig. 2f. In the following, different MHD-st-div($\mathbf{B}$) test cases are identified by the correspondent fraction of $E_{\text{div}}$.

The dependence of FV methods on solenoidal errors as modeled by Eq. (41) is discussed in Sect. 7.
<figure><img src="image_11.png" /><figcaption>Fig. 3 TD-twist test: (a) normalized helicity $H_V$; (b) complement of the normalized vector for the test field $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$</figcaption></figure>

# 5 Results: Finite Volume Methods

In this section we discuss the accuracy of FV methods presented in Sect. 2 using the test cases introduced in Sect. 4. As mentioned above, such methods require the full knowledge of the 3D magnetic field. Obviously, a non-perfectly solenoidal input field cannot be reproduced by a vector potential, and will affect each method in a different way. Therefore, the comparison between methods presented below should be read against the properties of the input field as described in Sect. 4. The purpose of these comparisons is not only to assess absolute accuracy in specific tests, but also to address the sensitivity of the methods towards certain studied parameters, such as resolution, twist, and topological complexity.

As anticipated in Sect. 3.1, the vector potentials obtained by FV methods are judged on their ability of reproducing the test fields and their corresponding potential fields. Indirectly, this is also a measure of the accuracy of the methods in the computation of $H_V$. The full metrics' values are provided in Tables 8–9, together with the computed values of helicity, and a partial but representative selection is reproduced in the plots of the next sections.

## 5.1 Dependence on Twist in the Titov and Démoulin Case

Figure 3a shows the dependence on the twist of $H_V$ computed by the different FV methods. With the exception of the Coulomb_GR method, all other methods are basically producing the same value of $H_V$ at all twists. More quantitatively, the spread in $H_V$ computed including all twist values and all methods except for the Coulomb_GR method, is only 2 %. In fact, DeVore methods are practically indistinguishable from each other. For a given twist value, for instance at $N = 1$, the spread of $H_V$ values around the average −0.079 is only 2.3 %, whereas average and spread become −0.084 and 16 %, respectively, if Coulomb_GR is included. The $H_V$ from the Coulomb_GR method, on the other hand, follows the same trend as the other methods, but has values about a factor two larger. In addition, there is no apparent correlation between the spread of $H_V$ and the value of the twist. All methods seem to be unaffected by this particular aspect of complexity of the field, at least at the resolution considered here. On the other hand, recalling the dependence on twist of the solenoidal errors discussed in Sect. 4.2 and Fig. 2b, it is found that larger spreads in $H_V$ correlates to larger values of divergence.
In order to address specifically the methods’ accuracy, we show in Fig. 3b the comple-
ment of the normalized vector error, ε<sub>N</sub>(**B**, **∇** × **A**) (cf. Eq. (38)), between the test field **B**
and the rotation of the corresponding vector potential **A** computed by each method. The De-
Vore methods sport extremely high accuracy metric of ε<sub>N</sub> = 0.994 or higher, and are indeed
indistinguishable from each other. The Coulomb_SY accuracy correlates inversely with the
solenoidal errors of that test (cf. Fig. 2b), in the sense that the accuracy is lower for propor-
tionally larger values of E<sub>div</sub>, as could be expected. Given the high accuracy in all twist cases,
such a dependence in not very clearly visible in Fig. 3b. Among the Coulomb methods, the
Coulomb_SY method has accuracy values (e.g., ε<sub>N</sub> = 0.976 at N = 1) that are only slightly
worse than the DeVore methods. Similarly to DeVore-gauge methods, the Coulomb_SY ac-
curacy correlates inversely with the solenoidal errors, but in this case with more pronounced
variations. The vector potential computed by the Coulomb_GR method has the largest error
in reproducing the input field, with values of the complement of the error vector as small as
ε<sub>N</sub> = 0.88. A trend similar to Coulomb_GR’s one is found for the Coulomb_JT method but
with a slightly smaller error (ε<sub>N</sub> = 0.903). Apparently, both Coulomb_GR and Coulomb_JT
method show an accuracy in this test that directly correlates with the solenoidal errors, i.e.,
the accuracy is lower (smaller values of ε<sub>N</sub>) for smaller values of E<sub>div</sub>. This counterintuitive
trend, however, is not confirmed by the tests in Sect. 7, and has to be regarded as insignifi-
cant. On the other hand, one could equally conjecture that Coulomb_GR and Coulomb_JT
methods show a direct correlation between ε<sub>N</sub> and the free energy E<sub>free</sub>, see again Fig. 2b. In
this sense, Coulomb_GR and Coulomb_JT could be said to obtain more accurate results for
fields with higher currents, which can be understood in terms of a stronger source term in
Eq. (13), and is not contradicted by the tests presented in the following sections. We notice
that, even for similar values of ε<sub>N</sub>, the helicity values obtained by the Coulomb_GR and
Coulomb_JT methods are quite different, the latter aligning with the DeVore ones within
2 %-variation.

As for the average $H_V$ values, the TD-twist test cases in Fig. 3a demonstrate that more twist does not necessarily translate into larger helicity values. Indeed, the $N=3$ case has a smaller helicity than $N=1$. The higher twist of the $N=3$ case is obtained by reducing the radius of the current channel ($a$ in Table 3) without changing the ambient field (i.e., the magnetic “charges”, $q$). The dependence of the current on $a$ is weak (see Eq. (6) of Titov and Démoulin 1999), and the difference in $H_V$ is mostly due to the difference in the—dominant—mutual helicity between ambient field and current channel.

In conclusion, in the TD-twist test cases all methods, with the exception of the
Coulomb_GR one, provide the same value of helicity within 2 % of variation. Furthermore,
in a range of different twist ranging from practically zero to 3 full turns, no evidence of direct
influence of the amount of twist on the accuracy of the methods is found. The variation of
the accuracy with twist clearly correlates for most of the methods with the solenoidal errors,
but an inverse correlation is also found for the Coulomb_JT and Coulomb_GR methods.

## 5.2 Time Evolution in the MHD-Emergence Stable and MHD-Emergence Unstable Cases

Figure 4 shows the helicity evolution for the MHD-st and MHD-un cases. No solution from the Coulomb_GR method is included for the MHD-st and MHD-un cases due to the long computation time that is required to obtain the two time series. The spread in $H_V$ computed including all available methods and time snapshots is basically negligible in the MHD-st case, amounting to just 0.2 %, see Fig. 4a. Hence, in the MHD-st case, even more than in the TD-twist case of Sect. 5.1, all considered methods yield essentially the same value
<figure><img src="image_12.png" /><figcaption>**Fig. 4** Normalized helicity $H_V$ for the (a) MHD-st and (b) MHD-un, as a function of time. No result from the Coulomb_GR method is included here due to computing-time limitations</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 5** Complement of the normalized vector for the test field, $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$, for the MHD-st (a) and the MHD-un (b) cases, as a function of time</figcaption></figure>

of $H_V$, which is a very encouraging result in view of future applications. The spread in $H_V$ is very small also in the MHD-un case, being 3% at the end of the simulation. Even though this represents a modestly larger spread in $H_V$ values, one has to recall that the MHD-un case corresponds to a time evolution where an eruption rearranges drastically the magnetic field, with ejection of field and currents from the top boundary. The challenge posed to numerical accuracy in such cases is indeed very high.

Concerning the accuracy metrics for the magnetic field, Fig. 5 shows the $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$ (cf. Eq. (38)), for the MHD-st and MHD-un cases, respectively, as a function of time. The accuracy metric shows that Coulomb_SY and all DeVore methods have comparable accuracy, of about 0.96 and 0.97 on average, respectively. On the other hand, the Coulomb_JT method has an accuracy for the complement of the vector error of $\epsilon_N = 0.54$ at the beginning of the simulation, increasing up to $\epsilon_N = 0.77$ at the end. Comparing with the evolution in the MHD-un case, we notice that the metric $\epsilon_N$ drops to worse values in correspondence of
the eruption, which correlates in time with the passage of current carrying field through the top boundary (without necessarily contradicting the direct correlation between $\epsilon_N$ and $E_{free}$ noticed for the Coulomb_JT method in Sect. 5.1, cf. Fig. 5b with Fig. 2e).

Indeed, the simulation MHD-st and MHD-un are similar in many respects, but one essential difference from the point of view of the $H_\nu$ computation is that the eruption in the MHD-un case generates the transit of a plasmoid carrying significant field through the top boundary. This seems to have hardly any consequence for DeVore methods, but may have more severe consequences for Coulomb methods that use the normal component of the field to specify the boundary conditions for $A$ (cf. Eq. (15)). The reason of the increase in the error might then be due to the fact that the Coulomb_JT method practically enforces flux balance by concentrating to the top boundary any possible error deriving from a nonsolenoidal input field (see Sect. 2.1). As significant field transits through the boundary as a consequence of the eruption, larger solenoidal errors affect the value of the field there and, as a consequence, the solution of $A_p$ that is computed form the boundary values. The correlation between the drop in the metrics of Coulomb_JT in correspondence with the plasmoid passage through the top boundary supports this speculation.

As for the average $H_\nu$ values, it is worth noticing that, for the simulation in exam, “larger $H_\nu$ value” does not immediately translate into “more likely to erupt”. Indeed, Fig. 5 shows that, according to all methods, the MHD-un case has less helicity than the stable one, MHD-st.

In conclusion, the helicity values $H_\nu$ in this most relevant test of MHD evolutions, both stable and unstable, show a very good agreement between different methods, namely within 3 % in the most challenging MHD-un case. It is only using very sensitive metrics such as $\epsilon_N$ that differences between methods of vector potentials’ computations can be disclosed. The results presented in this section are very encouraging for applications of helicity estimation methods to numerical simulations. The obtained values are largely independent of the specific methods employed, provided that the solenoidal property is sufficiently fulfilled, which allows for accurate and reliable studies of the properties of helicity as tracer of magnetic field evolution.

# 6 Dependence on Resolution

Numerical resolution, taken here to be the voxel size as customary in simulations, affects the solenoidal property of the input field as well as the accuracy of finite volume (FV) methods in solving for the vector potentials. Here we address the issue directly, trying to separate each contribution.

## 6.1 TD Resolution Test: Major Field Complexity with Little Flux Through Boundaries

Figure 2c shows that the solenoidal error in the TD-resolution cases is practically independent of resolution, and amounts to $E_{div} \simeq 3 \%$ at most in the examined resolution interval, $\Delta = [0.03, 0.12]$. Similarly, the $H_\nu$ values for each method separately, see Fig. 6a, are all essentially independent of resolution, except for a small variation of about 1 % for the Coulomb_JT case. On the other hand, the spread in values between methods over all resolutions is definitely more significant, being equal to 20 % if all methods are included, and 4 % if Coulomb_GR is excluded. Hence, this is a first indication that differences between
<figure><img src="image_14.png" /><figcaption>**Fig. 6** TD resolution test: (a) normalized helicity $H_V$; (b) complement of the normalized vector $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$</figcaption></figure>

methods (4 % at least) are more important than resolution (1 % at most) in determining the value of $H_V$, for similar levels of solenoidal errors.

Figure 6b shows the complement of the normalized vector for the test field $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$ (cf. Eq. (38)). The analysis of the complement of the vector error again separates the methods more markedly. There are no appreciable differences between the three DeVore methods, all around $\epsilon_N = 0.99$. Errors are larger on average going from Coulomb_SY ($\epsilon_N = 0.97$), to Coulomb_JT ($\epsilon_N = 0.90$), to the Coulomb_GR ($\epsilon_N = 0.89$).

Since there is no strong dependence on resolution of the solenoidal error in the TD-resolution case, then any such dependence that is found in the value of $H_V$ should be due to specific sensitivity of the different method to resolution.

## 6.2 LL Resolution Test: Minor Field Complexity with Significant Flux Through the Boundaries

In the tests analyzed so far, differences in $H_V$ are limited, except for the Coulomb_GR method and, to a slightly lesser extent, the Coulomb_JT method. In order to progress further in the analysis we consider the LL case described in Sect. 4.1.

Figure 7 summarizes the detailed analysis of the resolution effects on the computation of $H_V$ in that case. The main result is that a very clear dependence of $H_V$ on resolution is found when FV methods are applied to the LL case, see in particular Fig. 7a. In more detail, excluding again Coulomb_GR, the spread in $H_V$ at $n = 64$ is only 0.8 % (11 % if Coulomb_GR is included). This is even smaller than the spread in $H_V$ for the TD case of Sect. 6.1, at a comparable level of $E_{div}$ (cf. the end-to-end twist = 2.11 data point in Fig. 2a with the $n = 64$ data point of Fig. 2b). In the LL case we can decrease the resolution even further, thereby following the trends over a longer interval.

All methods show an increase of $H_V$ to more negative values for larger pixel size, with the exception of the $n = 32$ data point of the Coulomb_GR method. This is markedly different from the trend in Sect. 6.1 where, for a test field with basically the same $E_{div}$ across different resolutions, the $H_V$ values found by different methods were also basically independent of resolution, albeit not the same. Therefore, the dependence of $H_V$ on resolution in the LL case can be directly related to the presence of larger $\nabla \cdot \mathbf{B}$ at lower resolutions, as
<figure><img src="image_15.png" /><figcaption>**Fig. 7** LL resolution test: (a) normalized helicity, $H_\gamma$; complement of the normalized vector error for the test field for the (b) potential, $\epsilon_N(\mathbf{B}_p, \nabla \times \mathbf{A}_p)$, and (c) test field, $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$; energy ratios for the (d) potential, $\epsilon_E(\mathbf{B}_p, \nabla \times \mathbf{A}_p)$, and (e) test field, $\epsilon_E(\mathbf{B}, \nabla \times \mathbf{A})$; (f) $\langle |f_i \cdot \mathbf{A}_p| \rangle$ and $\langle |f_i \cdot \mathbf{A}| \rangle$ for the Coulomb methods only; as a function of resolution</figcaption></figure>
Fig. 2a shows. On the other hand, for a given resolution, the spread due to different methods ranges from 0.3 % at higher resolution, to 2.3 % at the lower one (5.4 % if *Coulomb_GR* is included). Among the different methods, *DeVore_GV* give less variable values with $\Delta$, whereas *Coulomb_JT* is more sensitive to it (and *Coulomb_GR* even more).

The complement of the normalized vector error and the energy metrics show as usual more differences between the methods. The DeVore methods are all very accurate in the computation of $\mathbf{A}_p$, with DeVore_SA slightly less sensitive than DeVore_KM and DeVore_GV (see Fig. 7b). Differences become entirely negligible in the computation of $\mathbf{A}$ (cf. Fig. 7c). The energy metrics $\epsilon_E$ (cf. Eq. (39)) provide similar information, with more accentuated spreads (see Fig. 7d and e). Therefore, we conclude that the most important difference between the three implementations of the DeVore method (DeVore_SA, DeVore_KM, DeVore_GV) is how the potential field is computed.

Among the Coulomb methods, *Coulomb_SY* is the more accurate one on average, and is only marginally less accurate than the DeVore methods, having slightly larger inaccuracies in computing $\mathbf{A}_p$ than $\mathbf{A}$. The *Coulomb_GR* method is even more accurate in this case than most of the DeVore methods in the computation of $\mathbf{A}_p$ (in both $\epsilon_N$ and $\epsilon_E$ metrics), but yields contradictory metrics for $\mathbf{A}$ (poor $\epsilon_N$ but good $\epsilon_E$). This results into $H_V$ values that stand more clearly apart from the trend. The *Coulomb_GR* improvement of $\epsilon_N$ with coarser pixel size in Fig. 7c is difficult to interpret without additional testing. The *Coulomb_JT* method lies somewhere in between the other two Coulomb methods, but overall yields $H_V$ values in line with the DeVore methods. The computation of $\mathbf{A}_p$ (respectively, $\mathbf{A}$) with the *Coulomb_JT* method yields $\epsilon_N \simeq 0.88$ (respectively, $\epsilon_N \simeq 0.90$) but is not particularly sensitive to resolution. Interestingly, Fig. 7d shows that the spread in energy metric for the potential field has an overall value including all methods of just 1.2 %, with a maximum value of 2 % at the lowest resolution. The same metrics for the input fields are 1.6 % and 2.4 %, respectively. In summary, the energy metric $\epsilon_E$, even though not exactly reproducing the accuracy of the methods as quantified by $\epsilon_N$, indicates a small spread of energy in the field recomputed from the vector potentials, even at low resolution and for all methods.

DeVore methods can impose the gauge $A_z = 0$ and $A_{p,z} = 0$ exactly also in numerical implementations. On the contrary, the gauge conditions $\nabla \cdot \mathbf{A} = 0$ and $\nabla \cdot \mathbf{A}_p = 0$ in Coulomb methods are numerically fulfilled only up to a finite precision. Hence, errors in fulfilling the gauge requirements might become a source of inaccuracy by generating spurious, nonsolenoidal components of the vector potentials. In order to quantify how well the solenoidal property of the vector potentials is satisfied for each Coulomb-gauge-based FV method, we show in Fig. 7f the fractional fluxes, $\langle |f_i(\mathbf{A})| \rangle$ and $\langle |f_i(\mathbf{A}_p)| \rangle$, as defined in Eq. (40), for the $\mathbf{A}_p$ and $\mathbf{A}$ solution of the LL resolution test cases, respectively. All Coulomb methods satisfy the gauge better (i.e., have lower $\langle |f_i| \rangle$) at higher resolutions, with similar rates of change, except for the *Coulomb_SY* method that has a minimum in correspondence of the $n = 64$ data point. The *Coulomb_SY* method satisfies the Coulomb gauge with almost identical accuracy for both $\mathbf{A}_p$ and $\mathbf{A}$, to a relatively low degree ($\langle |f_i| \rangle \simeq 10^{-4}$). The other two Coulomb methods show larger differences in the fulfillment of the solenoidal property of $\mathbf{A}_p$ and $\mathbf{A}$. In particular, both the *Coulomb_GR* method and the *Coulomb_JT* method respect the gauge condition better for $\mathbf{A}_p$ than for $\mathbf{A}$. The difference is minimal for *Coulomb_GR*—about $3 \times 10^{-5}$—and about one order of magnitude for *Coulomb_JT*. At the same time, *Coulomb_JT* provides the most solenoidal $\mathbf{A}_p$ and $\mathbf{A}$ of all Coulomb methods, at all resolutions, with values of $\langle |f_i| \rangle$ below $10^{-8}$ for both potentials at the highest resolution.

In summary, the proper fulfillment of the Coulomb gauge is attained at different level of accuracy by the different methods, especially as far as $\mathbf{A}$ is concerned, but no obvious correlation between violation of the solenoidal property of the vector potentials and accuracy of the vector potential is found.
# 7 Dependence on Divergence

In Sects. 6.1 and 6.2 we show that the pixel size affects not only the accuracy of methods in solving for the vector potentials, but also the degree of divergence that is caused by the resolution of the test fields. More generally, in practical cases, when trying to estimate the helicity of a 3D dataset, a nonzero divergence of **B** is always present. It is thus of practical importance to determine the level of confidence of an helicity measurement given its level of nonsolenoidality, here expressed in function of $E_{\text{div}}$. The present section aims at giving a first test of the impact of the nonsolenoidal effect on the degree of precision of helicity measurements.

Similarly to classical helicity (cf. Eq. (2)), assuming that **B** and **B**<sub>p</sub> effectively have the same distribution of the normal component on $\partial V$, a gauge transformation ($A, A_p$) —→ ($A + \nabla \psi, A_p + \nabla \psi_p$) induce the following difference on the relative magnetic helicity:

$$ \mathcal{H}_V(\mathbf{A} + \nabla \psi, \mathbf{A}_p + \nabla \psi_p) = \mathcal{H}_V(\mathbf{A}, \mathbf{A}_p) - \int_V (\psi + \psi_p)(\nabla \cdot \mathbf{B} - \nabla \cdot \mathbf{B}_p) dV. \quad (42) $$

Comparing one FV method against another is theoretically equivalent to performing a gauge transformation: for a purely solenoidal field they should provide identical values of helicity within the numerical precision of each method. While for a finite nonsolenoidal field it is always theoretically possible to find a gauge transformation that will lead to a significantly large difference, in practice most methods are solving for **A** and **A**<sub>p</sub> in such a way that the amplitude of the vector potentials turns out to be dominant over the contributions from the gauges $\psi$ and $\psi_p$. Hence, for a given dataset, the difference of amplitude of the vector potentials computed by each method is relatively small. Thus, for a given $E_{\text{div}}$, the different methods provide helicity values that remain within a certain range of each other. As $E_{\text{div}}$ is infinitely small, the methods are expected to provide helicity values which are close to each other, while for large $E_{\text{div}}$, the methods should provide helicity values presenting a larger spread.

In the present section, we consider the MHD-st-div(**B**) test described in Sect. 4.4, which allows us to estimate the dispersion of the helicity values obtained with the different methods for different level of nonsolenoidality. In the present test the Coulomb_GR method is not included because of computational limitations.

The amplitude of the nonsolenoidal component in the MHD-st-div(**B**) test field as quantified by $E_{\text{div}}$ for different values of $\delta$ in Eq. (41) is reported in Table 5 and shown in Fig. 2f, where $E_{\text{div}}$ is shown growing from 0.2 % up to 14 % for the considered range of $\delta$ values (the corresponding values of $|f_i|$ can be found in Table 7). The contributions $E_{1,ns}$ and $E_{\text{mix}}$ to $E_{\text{div}}$ in Eq. (44) and their variation with $\delta$ are shown in Fig. 8a. The energy $E_{\text{div}}$ turns out to be a nonlinear function of $\delta$, with a minimum close, but not at, $\delta = 0$. This apparent contradiction is due to the definition equation (44) that forbids cancellation between different contributions, and practically shifts the zero of $E_{\text{div}}$ to slightly higher values. For completeness, the nonsolenoidal error $E_{\text{ns}}$ given by Eq. (44) without the absolute value, i.e., allowing for cancellation between terms, is also plotted in Fig. 8a.

The values of $H_V$ obtained in the MHD-st-div(**B**) test cases are plotted in Fig. 8b. We confirmed that for low values of $E_{\text{div}}$, the methods are providing $H_V$ which are close to each other, while the dispersion of the values obtained are spreading as $E_{\text{div}}$ is increasing. All method follow a qualitatively similar exponential trend as a function of $E_{\text{div}}$, both in $H_V$ (Fig. 8b) as well as in the accuracy metric $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$ (Fig. 8c). The only exception to
<figure><img src="image_16.png" /><figcaption>**Fig. 8** MHD-st-div(B) test (a) nonsolenoidal contributions to the energy, as a function of the divergence parameter, $\delta$, see Eqs. (41), (44). As a function of $E_{\text{div}}$: (b) helicity $H_V$; (c) complement of the normalized vector error for the test field, $\epsilon_N(\mathbf{B}, \nabla \times \mathbf{A})$; (d) $\langle |f_i \mathbf{A}| \rangle$ for the Coulomb methods only</figcaption></figure>

a smooth trend is the $E_{\text{div}} = 2\%$-case of the Coulomb_JT method, for which no particular reason was identified.

In the most solenoidal case, for $E_{\text{div}} = 0.2\%$, the relative dispersion of $H_V$ obtained is of $0.8\%$. Excluding the $E_{\text{div}} = 2\%-$case of the Coulomb_JT method, for each case with $E_{\text{div}} \le 4\%$, the relative spread in $H_V$ remains lower than $1\%$. At $E_{\text{div}} = 2\%$, taking all the methods into account, the relative spread in the $H_V$ is of the order of $6\%$.

Excluding the Coulomb_GR method, we note that in all the previous tests, the maximum spread of $H_V$ observed was of $4\%$ (in the TD-resolution test of Sect. 6.1). This inclines us to state that within that range of $E_{\text{div}}$, the differences between the methods are not related to the nonsolenoidality but rather to the intrinsic numerical error within each method.

For $E_{\text{div}} = 8.2\%$, the relative dispersion of $H_V$ is also relatively low, equal to $1.9\%$. However, in the least solenoidal case, for $E_{\text{div}} = 14.4\%$, the dispersion in helicity estimation obtained for the different methods reaches $18\%$ of the average value. This is much higher than most of the errors that we have encountered so far. We believe that in this range of $E_{\text{div}}$, the gauge dependence is directly impacting the estimation of $H_V$.
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 $\epsilon_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 $H_\nu$ and $\epsilon_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 Coulomb-gauge condition extremely well for $A_p$ ($\langle|f_i|\rangle \simeq 10^{-9}$), and still excellently for $A$ ($\langle|f_i|\rangle \simeq 10^{-6}$ or lower). Interestingly, $\langle|f_i(A_p)|\rangle$ is practically independent of $E_{\text{div}}$, 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 $\mathcal{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 $\langle|f_i(A_p)|\rangle \simeq \langle|f_i(A)|\rangle \simeq 10^{-4}$ than those by the Coulomb_JT method. In this sense, the divergence-cleaning strategy adopted by Coulomb_SY is less efficient in imposing the Coulomb gauge. However, such values of $\langle|f_i|\rangle$ seem still low enough to guarantee a relative high accuracy, as testified by $\epsilon_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 well-controlled 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 $E_{\text{div}}$ is below 1%, but become abruptly more important above that threshold. For a dataset with $E_{\text{div}}$ comprised between 1 and 8%, using on FV method or the other would lead up to 6% difference in the estimation of $H_\nu$ (excluding Coulomb_GR method). For higher $E_{\text{div}}$, 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 $E_{\text{div}}$ (and hence its influence on $H_\nu$) 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 $H_\nu$ 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 solar-like separation between a current-carrying 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 $H_V$, 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 $H_V$. This is in a way confirmed by the remarkable absence of differences in the MHD-st case, in which $H_V$ 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 self-helicity $H_J$ is almost a factor 100 smaller than the total helicity $H_V$, 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, $E_{div}$.

We find a rapid increase of $H_V$ fluctuations as $E_{div}$ 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 under-resolved 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 $H_V$ 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 $c_f$ 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 $A_p$, in a way contrary to the DeVore methods. In particular, the *Coulomb_JT* method seems to be sensitive to current on the boundaries (as MHD-un 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 $H_V$ 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 $\langle |f_i(\mathbf{A})| \rangle$ below $10^{-3}$ yields helicity values comparable with those from DeVore methods (e.g., within 0.8 % in the moderate-resolution 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 $H_V$. 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 Flux-Tubes 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 Twist-Number Method

The twist-number method needs identifiable flux-rope 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, MHD-st and MHD-un cases, as well as some nonlinear force-free 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 $H_{\text{twist}}$ 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 self-helicity only, Fig. 10 compares the not-normalized value of $H_{V,J}$ of Eq. (6) from the DeVore_GV method with the $H_{\text{twist}}$ 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, $H_{\text{twist}}$ has the same value of $H_{V,J}$ within errors. In the other cases, larger differences between the TN method and DeVore_GV method are present.
<figure><img src="image_17.png" /><figcaption>**Fig. 9** Application of the TN method to the TD *N* = 3 case. Left: a vertical slice of the Q map in the *xz*-plane 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</figcaption></figure>

<table><caption>**Table 6** Helicity in the twist (upper) and resolution (lower) TD tests, computed with the FV-DeVore_GV method (2nd to 4th column), and with the TN method (5th column). Note that $\mathcal{H}_V$, $\mathcal{H}_V,J$, and $\mathcal{H}_{\text{twist}}$ are not normalized</caption><thead><tr><th>TD case</th><th>H<sub>V</sub><br>(FV-DeVore_GV)</th><th>H<sub>V</sub><br>(FV-DeVore_GV)</th><th>H<sub>V,J</sub><br>(TN)</th><th>H<sub>twist</sub><br>(TN)</th></tr></thead><tbody><tr><td>N = 0.1</td><td>-0.0057</td><td>-2.23147</td><td>0.00084</td><td>-0.16 ± 0.06</td></tr><tr><td>N = 0.5</td><td>-0.0290</td><td>-8.33340</td><td>-0.08374</td><td>-0.66 ± 0.13</td></tr><tr><td>N = 1</td><td>-0.0782</td><td>-7.21022</td><td>-0.55240</td><td>-0.50 ± 0.05</td></tr><tr><td>N = 3</td><td>-0.0527</td><td>-1.78752</td><td>-0.08976</td><td>-0.087 ± 0.020</td></tr><tr><td>Δ = 0.03</td><td>-0.0782</td><td>-7.20592</td><td>-0.54921</td><td>-0.56 ± 0.06</td></tr><tr><td>Δ = 0.06</td><td>-0.0782</td><td>-7.21022</td><td>-0.55240</td><td>-0.50 ± 0.05</td></tr><tr><td>Δ = 0.12</td><td>-0.0782</td><td>-7.21039</td><td>-0.55444</td><td>-0.49 ± 0.05</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Fig. 10** Application of the TN method to the TD-twist (a) and resolution (b) tests. The two curves represent the values of $\mathcal{H}_{\text{twist}}$ of Eq. (32) for the TN method, and of the (nonnormalized) value of $\mathcal{H}_{V,J}$ for the DeVore_GV method</figcaption></figure>
By contrast, the not-normalized values of the total helicity $\mathcal{H}_V = -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 $\mathcal{H}_{\text{twist}}$ provided by the TN method is an estimation of $\mathcal{H}_{V,J}$ rather than of $\mathcal{H}$ (or equivalently $H_V$), 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 $\mathcal{H}_{V,J}$ is reproduced.

## 9.2 Connectivity-Based Method

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 (force-freeness and minimal connectivity-length 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 connectivity-based 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 force-free field approach of Georgoulis and LaBonte (2007), in which a single value of the force-free parameter $\alpha$ 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 LFF-mode tends to pick up mostly the potential field component. On the other hand, when applied to the LL cases, where a large-scale $\alpha$ is present, the CB method yields values of $H_V$ 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 force-free theory. Hence, we do not report further on such applications here.

In the MHD-st and MHD-un 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 MHD-st and MHD-un cases. Figure 11a shows an example of the flux partition and of the resulting connectivity matrix for the MHD-st case at $t = 150$.

Figure 12 compare the CB method with the DeVore_GV finite volume method.

In the MHD-st case (Fig. 12a), between $t = 50$ and $t = 95$, the of $H_V$ 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 $H_V$ values obtained by the CB method settle on a lower, roughly constant value $H_V = 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 MHD-un 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
<figure><img src="image_19.png" /><figcaption>**Fig. 11** Application of the CB method to the MHD-st case. (a) Magnetic connectivity of the MHD-emergence 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 gray-scale, with *black/white contours* outlining positive/negative-polarity partitions. The flux-weighted centroids of the partitions are denoted by *crosses*. Magnetic connections are projected on the lower boundary by *colored lines*, with *different colors* denoting different connected-flux contents. The connected flux includes approx. 93.6% of the total flux present in the field of view. (b) The force-freeness metric *σ*<sub>J</sub> as a function of the height of the bottom boundary for MHD-st 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</figcaption></figure>

present at the time of the eruption, very much the same as for the DeVore_GV method. How-
ever, 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 H<sub>V</sub> in the CB time se-
ries. As a matter of fact, the examination of the time evolution of the magnetic field at z = 0
in the MHD-st and MHD-un 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 MHD-st from MHD-un 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 H<sub>V</sub> for both the MHD-st (equal to 0.015, from t = 95 onwards) and the
MHD-un case (equal to 0.070 on the same time interval).

Besides the similarity of the MHD-st and MHD-un field distributions at *z* = 0, differ-
ences with the DeVore_GV method may have several origins. First of all, the CB-method
approach yields a *minimal* value of helicity for a given photospheric configuration. In this
sense, it is expected that the CB *H*<sub>V</sub>-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 force-free. Figure 11b shows at a representative time, however, that
in large part of the volume of the MHD-st case the field is not force-free, as quantified
by the relative ratio of the current that is perpendicular to the field in the volume, σ<sub>J</sub> ≡
<figure><img src="image_20.png" /><figcaption>**Fig. 12** Comparison of $H_V$ between the CB and the DeVore_GV method applied to the MHD-st (a, c) and MHD-un cases (b, d) in the full domain (a, b) and in the reduced, more force-free 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</figcaption></figure>

($\int_V |J_\perp| dV$)/($\int_V |J| dV$). In particular, the value of $\sigma_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 force-free assumption in this case, we repeat the $H_V$ calculation with the CB and DeVore_GV methods in a reduced volume starting at $z = 8.9$, where $\sigma_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 MHD-st and MHD-un, respectively. Since the volume is now changed, also the corresponding DeVore_GV estimations are recalculated for this reduced, more force-free volume. In the MHD-st case, the average CB $H_V$ 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 full-volume case). In the MHD-un case in Fig. 12d, the curves obtained by the two methods are slightly closer than in the full-volume case of Fig. 12b, although only marginally (the ratio of the end values being 1.5). Hence, the application to the more force-free, upper part of the volume improves the match between the
CB and the DeVore_GV results, markedly so in the MHD-st case, which is a clear indication that the fulfillment of the force-free requirement may help to partially compensate for the $H_\nu$ 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 force-free 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 MHD-st case improves the matching between CB and DeVore_GV methods. The CB-average (after $t = 95$) $H_\nu$ value is slightly above 0.10, which results in just a factor two between the corresponding end values. In the MHD-un case, differences with the standard application of the CB method that does not used the three-dimensional 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 $H_\nu$ 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 MHD-un 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 MHD-st case the helicity is significantly under-estimated, by a factor eight at the end of the simulation. The deficit in $H_\nu$ that the CB method shows in the MHD-st case, can be partially ascribed to the exiguous differences between the MHD-st and MHD-un cases in terms of flux distributions at $z = 0$. Additionally, the field at that plane is not quite force-free, and the CB results are shown to improve for a more force-free test volume. Furthermore, the minimization of the connectivity matrix seems to represent connectivity fairly well, in the sense that the implementation of the full three-dimensional information, although improving the CB estimation in the MHD-st case, does not entirely remove the under-estimation of $H_\nu$.

# 10 Conclusions

In this work we review, benchmark, and compare the currently available methods for the computation of the relative magnetic helicity, $H_\nu$, in finite volumes. Given that the three-dimensional magnetic-field 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 ($\nabla \cdot \mathbf{A} = 0$) and DeVore ($A_z = 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 128³ 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 solar-physics applications, from 3D force-free equilibria (LL and TD of Sects. 4.1 and 4.2, respectively), to snapshots of time-dependent non-force-free MHD simulations of flux emergence (Sect. 4.3). Depending on details of the test field being studied, the accuracy in the computation of $H_V$ by different methods is found to vary to some extent, especially for Coulomb methods. We can, however, definitely conclude that in solar-like cases practically all FV methods converge to the same helicity value within few percent.<fnref n="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 $H_V$ in the most relevant test of time-dependent MHD evolution in a coronal model volume (i.e., the MHD-st and MHD-un 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 TD-twist case are within 2 %. In other words, when helicity computations are applied to numerical volumes as in this article, differences in the way $H_V$ 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 helicity-flux integration methods applied to MHD-st and MHD-un 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, MHD-st, and MHD-un cases. For the TD cases, we find that the TN method yields relatively accurate estimations of twist and possibly of $H_{V,j}$ in high-twist cases, but not of the total helicity $H_V$. A report on the application of the TN method to other cases, including the MHD-st and MHD-un ones, is in preparation by Guo et al. (2017).

<footnote n="3">³From this statistic, the Coulomb_GR method is excluded.</footnote>
The CB method is designed to be applied to complex photospheric flux distributions with an unknown coronal magnetic field, where a multi-polar partition of the flux yields a coronal field approximated by a collection of flux tubes of constant $J_z/B_z$. 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 flux-tube, purely linear approximation of the (force-free) field. In the more complex MHD-un 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 MHD-st 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 $H_v$ values. The impact on helicity values of the employed coronal model, being it from a nonlinear force-free extrapolation or from a data-driven 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 2014-051. E. Pariat acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-15-CE31-0001 (project HeliSol). S. Anfinogentov acknowledges support of the Russian Foundation of Basic Research under grants 15-02-01089-a, 15-02-03835-a, 16-32-00315-mol-a, and 15-32-20504-mol-a-ved, and the Federal Agency for Scientific Organisations base project II.16.3.2 “Non-stationary 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), NKBRFS (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): P25383-N27. F. Chen acknowledge the support by the International Max-Planck Research School (IMPRS) for Solar System Science at the University of Göttingen. S. Yang is supported by grants KJCX2-EWT07 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 quadri-core bi-Xeon 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 = E_{p,s} + E_{J,s} + E_{p,ns} + E_{J,ns} + E_{mix}, \qquad (43)$$

where $E_{p,s}$ and $E_{J,s}$ are the energies associated to the potential and current-carrying solenoidal contributions, $E_{p,ns}$ and $E_{J,ns}$ are those of the nonsolenoidal contributions, and $E_{mix}$ 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 $E_{mix}$. For a perfectly solenoidal field, it is $E_{p,s} = E_p$, $E_{J,s} = E - E_p$, and $E_{p,ns} = E_{J,ns} = E_{mix} = 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

$$E_{\text{div}} = E_{p,ns} + E_{J,ns} + |E_{\text{mix}}|, \qquad (44)$$

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 $E_{ns}$ is considered, corresponding to Eq. (44) with $E_{mix}$ instead of $|E_{mix}|$. 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 current-carrying part, i.e., $E_{free} = E_{J,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

$$C_{\text{vec}} = \frac{\sum_i \mathbf{X}_i \cdot \mathbf{Y}_i}{(\sum_i |\mathbf{X}_i|^2 \sum_i |\mathbf{Y}_i|^2)^{1/2}} \qquad (45)$$

$$C_{\text{CV}} = \frac{1}{N} \sum_{i} \frac{\mathbf{X}_{i} \cdot \mathbf{Y}_{i}}{|\mathbf{X}_{i}| |\mathbf{Y}_{i}|} \qquad (46)$$

$$\epsilon_M = 1 - \frac{1}{N} \sum_i \frac{|\mathbf{Y}_i - \mathbf{X}_i|}{|\mathbf{X}_i|} \qquad (47)$$
<table><caption>**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 [$E_{p,s}$], current-carrying/solenoidal [$E_{J,s}$], potential/nonsolenoidal [$E_{p,ns}$], current-carrying/nonsolenoidal [$E_{J,ns}$], nonsolenoidal mixed [$E_{mix}$] contributions normalized to the corresponding total magnetic energy, E. The last column contains the average fractional flux [$\langle |f_i(\mathbf{B})| \rangle$]. In the text, any contribution $E_{xx}$ is intended to be normalized to E, i.e., it is here indicated as $E_{xx}/E$, for each test case separately</caption><thead><tr><th>Label</th><th>E</th><th>E<sub>p,s</sub>/E</th><th>E<sub>J,s</sub>/E</th><th>E<sub>p,ns</sub>/E</th><th>E<sub>J,ns</sub>/E</th><th>E<sub>mix</sub>/E</th><th>⟨|f<sub>i</sub>(B)|⟩ × 10<sup>5</sup></th></tr></thead><tbody><tr><td>LL n = 32</td><td>5.48e+02</td><td>0.78</td><td>0.25</td><td>2.34e-03</td><td>6.52e-04</td><td>-3.82e-02</td><td>26.43</td></tr><tr><td>LL n = 64</td><td>5.30e+02</td><td>0.76</td><td>0.25</td><td>2.60e-04</td><td>8.22e-05</td><td>-1.33e-02</td><td>3.34</td></tr><tr><td>LL n = 128</td><td>5.25e+02</td><td>0.75</td><td>0.26</td><td>2.37e-05</td><td>9.15e-06</td><td>-3.90e-03</td><td>0.62</td></tr><tr><td>LL n = 256</td><td>5.24e+02</td><td>0.75</td><td>0.26</td><td>2.38e-06</td><td>1.32e-06</td><td>-1.07e-03</td><td>0.27</td></tr><tr><td>TD N = 0.1 Δ = 0.06</td><td>4.69e+01</td><td>1.00</td><td>0.00</td><td>1.18e-05</td><td>1.02e-05</td><td>-1.20e-03</td><td>0.39</td></tr><tr><td>TD N = 0.5 Δ = 0.06</td><td>3.35e+01</td><td>0.98</td><td>0.02</td><td>1.18e-05</td><td>2.68e-05</td><td>3.02e-03</td><td>0.71</td></tr><tr><td>TD N = 1 Δ = 0.03</td><td>1.08e+01</td><td>0.81</td><td>0.17</td><td>3.08e-06</td><td>5.94e-04</td><td>2.53e-02</td><td>0.75</td></tr><tr><td>TD N = 1 Δ = 0.06</td><td>1.09e+01</td><td>0.81</td><td>0.17</td><td>1.21e-05</td><td>6.12e-04</td><td>2.44e-02</td><td>1.89</td></tr><tr><td>TD N = 1 Δ = 0.12</td><td>1.09e+01</td><td>0.81</td><td>0.17</td><td>6.33e-05</td><td>6.78e-04</td><td>2.14e-02</td><td>6.86</td></tr><tr><td>TD N = 3 Δ = 0.06</td><td>3.80e+00</td><td>0.85</td><td>0.14</td><td>1.35e-05</td><td>1.16e-04</td><td>8.07e-03</td><td>1.35</td></tr><tr><td>MHD-st t = 50</td><td>4.04e+02</td><td>0.52</td><td>0.50</td><td>3.86e-05</td><td>3.34e-04</td><td>-1.67e-02</td><td>9.99</td></tr><tr><td>MHD-st t = 85</td><td>5.31e+02</td><td>0.53</td><td>0.48</td><td>3.49e-05</td><td>3.00e-04</td><td>-1.15e-02</td><td>13.08</td></tr><tr><td>MHD-st t = 120</td><td>6.03e+02</td><td>0.46</td><td>0.55</td><td>2.08e-05</td><td>2.62e-04</td><td>-7.06e-03</td><td>13.66</td></tr><tr><td>MHD-st t = 155</td><td>6.42e+02</td><td>0.41</td><td>0.59</td><td>1.49e-05</td><td>2.31e-04</td><td>-5.46e-03</td><td>12.27</td></tr><tr><td>MHD-st t = 190</td><td>6.67e+02</td><td>0.38</td><td>0.62</td><td>1.34e-05</td><td>2.22e-04</td><td>-4.49e-03</td><td>12.10</td></tr><tr><td>MHD-un t = 50</td><td>3.67e+02</td><td>0.39</td><td>0.62</td><td>3.92e-05</td><td>2.40e-04</td><td>-1.10e-02</td><td>11.25</td></tr><tr><td>MHD-un t = 80</td><td>4.43e+02</td><td>0.42</td><td>0.58</td><td>4.36e-05</td><td>2.55e-04</td><td>-8.34e-03</td><td>22.83</td></tr><tr><td>MHD-un t = 110</td><td>4.64e+02</td><td>0.40</td><td>0.60</td><td>3.09e-05</td><td>2.93e-04</td><td>-5.11e-03</td><td>27.10</td></tr><tr><td>MHD-un t = 140</td><td>4.46e+02</td><td>0.41</td><td>0.59</td><td>2.16e-05</td><td>2.50e-04</td><td>-3.63e-03</td><td>32.29</td></tr><tr><td>MHD-un t = 150</td><td>4.07e+02</td><td>0.44</td><td>0.56</td><td>2.23e-05</td><td>2.57e-04</td><td>-3.26e-03</td><td>46.16</td></tr><tr><td>MHD-un t = 160</td><td>3.81e+02</td><td>0.47</td><td>0.53</td><td>2.32e-05</td><td>3.02e-04</td><td>-4.17e-03</td><td>43.88</td></tr><tr><td>MHD-un t = 190</td><td>3.62e+02</td><td>0.49</td><td>0.51</td><td>2.33e-05</td><td>3.35e-04</td><td>-5.14e-03</td><td>32.42</td></tr><tr><td>MHD-st-div(B) 0.2 %</td><td>4.14e+02</td><td>0.50</td><td>0.50</td><td>3.92e-05</td><td>3.14e-04</td><td>-1.45e-03</td><td>14.95</td></tr><tr><td>MHD-st-div(B) 1 %</td><td>4.06e+02</td><td>0.51</td><td>0.50</td><td>3.85e-05</td><td>1.70e-04</td><td>-1.08e-02</td><td>9.77</td></tr><tr><td>MHD-st-div(B) 2 %</td><td>4.03e+02</td><td>0.52</td><td>0.50</td><td>3.88e-05</td><td>4.90e-04</td><td>-1.98e-02</td><td>16.36</td></tr><tr><td>MHD-st-div(B) 4 %</td><td>4.09e+02</td><td>0.53</td><td>0.51</td><td>4.16e-05</td><td>2.25e-03</td><td>-3.82e-02</td><td>47.11</td></tr><tr><td>MHD-st-div(B) 8 %</td><td>4.64e+02</td><td>0.49</td><td>0.57</td><td>5.43e-05</td><td>9.19e-03</td><td>-7.25e-02</td><td>85.67</td></tr><tr><td>MHD-st-div(B) 14 %</td><td>1.01e+03</td><td>0.31</td><td>0.77</td><td>9.64e-05</td><td>3.06e-02</td><td>-1.12e-01</td><td>165.15</td></tr></tbody></table>

respectively the vector correlation, the Cauchy-Schwartz 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 MHD-st and MHD-un, and Table 10 those for the divergence and gauge tests of Sect. 7.
<table><caption>**Table 8** Metrics of accuracy of vector potentials as defined in Eqs. (38), (39) and Eqs. (45)–(47), and helicity $H_\gamma$ as defined in Eq. (36), for the LL (first page) and TD (second and third pages) test cases</caption><thead><tr><th rowspan="2">Label</th><th colspan="4">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="4">B vs. &nabla; &times; A</th><th rowspan="2">H<sub>&gamma;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th></tr></thead><tbody><tr><td>Coulomb_SY 32</td><td>0.9997</td><td>1.0000</td><td>0.9885</td><td>0.9938</td><td>1.0033</td><td>0.9996</td><td>0.9924</td><td>0.9742</td><td>0.9130</td><td>0.9965</td><td>-0.1462</td></tr><tr><td>Coulomb_SY 64</td><td>1.0000</td><td>1.0000</td><td>0.9963</td><td>0.9981</td><td>1.0014</td><td>0.9999</td><td>0.9970</td><td>0.9880</td><td>0.9496</td><td>0.9994</td><td>-0.1427</td></tr><tr><td>Coulomb_SY 128</td><td>1.0000</td><td>1.0000</td><td>0.9989</td><td>0.9992</td><td>1.0005</td><td>1.0000</td><td>0.9990</td><td>0.9943</td><td>0.9726</td><td>1.0001</td><td>-0.1418</td></tr><tr><td>Coulomb_SY 256</td><td>1.0000</td><td>1.0000</td><td>0.9990</td><td>0.9978</td><td>1.0001</td><td>1.0000</td><td>0.9996</td><td>0.9965</td><td>0.9836</td><td>1.0003</td><td>-0.1415</td></tr><tr><td>DeVore_KM 32</td><td>0.9992</td><td>0.9998</td><td>0.9799</td><td>0.9854</td><td>0.9915</td><td>0.9997</td><td>1.0000</td><td>0.9891</td><td>0.9958</td><td>0.9978</td><td>-0.1483</td></tr><tr><td>DeVore_KM 64</td><td>0.9999</td><td>0.9999</td><td>0.9923</td><td>0.9929</td><td>0.9963</td><td>1.0000</td><td>1.0000</td><td>0.9968</td><td>0.9989</td><td>0.9988</td><td>-0.1441</td></tr><tr><td>DeVore_KM 128</td><td>1.0000</td><td>0.9999</td><td>0.9970</td><td>0.9965</td><td>0.9988</td><td>1.0000</td><td>1.0000</td><td>0.9991</td><td>0.9997</td><td>0.9996</td><td>-0.1426</td></tr><tr><td>DeVore_KM 256</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9982</td><td>0.9996</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>0.9999</td><td>-0.1419</td></tr><tr><td>Coulomb_JT 32</td><td>0.9974</td><td>0.8772</td><td>0.8841</td><td>0.6307</td><td>1.0320</td><td>0.9978</td><td>0.8583</td><td>0.8856</td><td>0.5295</td><td>1.0500</td><td>-0.1499</td></tr><tr><td>Coulomb_JT 64</td><td>0.9974</td><td>0.8658</td><td>0.8843</td><td>0.5771</td><td>1.0242</td><td>0.9980</td><td>0.8409</td><td>0.8963</td><td>0.4688</td><td>1.0410</td><td>-0.1437</td></tr><tr><td>Coulomb_JT 128</td><td>0.9967</td><td>0.8595</td><td>0.8810</td><td>0.5420</td><td>1.0166</td><td>0.9975</td><td>0.8313</td><td>0.9001</td><td>0.4272</td><td>1.0275</td><td>-0.1414</td></tr><tr><td>Coulomb_JT 256</td><td>0.9958</td><td>0.8562</td><td>0.8778</td><td>0.5207</td><td>1.0121</td><td>0.9968</td><td>0.8262</td><td>0.9009</td><td>0.4011</td><td>1.0177</td><td>-0.1404</td></tr><tr><td>DeVore_GV 32</td><td>0.9988</td><td>0.9997</td><td>0.9781</td><td>0.9797</td><td>0.9677</td><td>0.9995</td><td>1.0000</td><td>0.9884</td><td>0.9958</td><td>0.9798</td><td>-0.1427</td></tr><tr><td>DeVore_GV 64</td><td>0.9998</td><td>0.9999</td><td>0.9915</td><td>0.9892</td><td>0.9851</td><td>0.9999</td><td>1.0000</td><td>0.9965</td><td>0.9989</td><td>0.9918</td><td>-0.1420</td></tr><tr><td>DeVore_GV 128</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9942</td><td>0.9949</td><td>1.0000</td><td>1.0000</td><td>0.9990</td><td>0.9997</td><td>0.9974</td><td>-0.1416</td></tr><tr><td>DeVore_GV 256</td><td>1.0000</td><td>1.0000</td><td>0.9985</td><td>0.9968</td><td>0.9985</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>0.9993</td><td>-0.1415</td></tr><tr><td>DeVore_SA 32</td><td>0.9996</td><td>1.0000</td><td>0.9887</td><td>0.9959</td><td>1.0044</td><td>0.9997</td><td>1.0000</td><td>0.9894</td><td>0.9959</td><td>0.9960</td><td>-0.1471</td></tr><tr><td>DeVore_SA 64</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9989</td><td>1.0015</td><td>1.0000</td><td>1.0000</td><td>0.9969</td><td>0.9989</td><td>0.9977</td><td>-0.1428</td></tr><tr><td>DeVore_SA 128</td><td>1.0000</td><td>1.0000</td><td>0.9991</td><td>0.9997</td><td>1.0005</td><td>1.0000</td><td>1.0000</td><td>0.9991</td><td>0.9997</td><td>0.9989</td><td>-0.1418</td></tr><tr><td>DeVore_SA 256</td><td>1.0000</td><td>1.0000</td><td>0.9998</td><td>0.9999</td><td>1.0001</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>0.9995</td><td>-0.1415</td></tr><tr><td>Coulomb_GR 32</td><td>0.9997</td><td>1.0000</td><td>0.9882</td><td>0.9950</td><td>1.0042</td><td>0.9992</td><td>0.9525</td><td>0.9247</td><td>0.7829</td><td>1.0036</td><td>-0.1655</td></tr><tr><td>Coulomb_GR 64</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9989</td><td>1.0015</td><td>0.9977</td><td>0.8661</td><td>0.8471</td><td>0.5270</td><td>1.0048</td><td>-0.1861</td></tr><tr><td>Coulomb_GR 128</td><td>1.0000</td><td>1.0000</td><td>0.9991</td><td>0.9997</td><td>1.0005</td><td>0.9974</td><td>0.8573</td><td>0.8388</td><td>0.4970</td><td>1.0013</td><td>-0.1846</td></tr></tbody></table>
<table><caption>**Table 8** (Continued)</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="3">B vs. &nabla; &times; A</th><th colspan="3">H<sub>&lambda;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th></tr></thead><tbody><tr><td>Coulomb_SY N01 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9994</td><td>0.9991</td><td>1.0002</td><td>1.0000</td><td>0.9997</td><td>0.9939</td><td>0.9816</td><td>0.9999</td><td>-0.0057</td></tr><tr><td>Coulomb_SY N05 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9993</td><td>0.9991</td><td>1.0002</td><td>1.0000</td><td>0.9997</td><td>0.9921</td><td>0.9792</td><td>0.9959</td><td>-0.0289</td></tr><tr><td>Coulomb_SY N1 0.03</td><td>1.0000</td><td>1.0000</td><td>0.9987</td><td>0.9960</td><td>1.0001</td><td>0.9998</td><td>0.9995</td><td>0.9771</td><td>0.9683</td><td>0.9745</td><td>-0.0768</td></tr><tr><td>Coulomb_SY N1 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9992</td><td>0.9989</td><td>1.0003</td><td>0.9998</td><td>0.9991</td><td>0.9763</td><td>0.9628</td><td>0.9751</td><td>-0.0769</td></tr><tr><td>Coulomb_SY N1 0.12</td><td>1.0000</td><td>1.0000</td><td>0.9977</td><td>0.9977</td><td>1.0009</td><td>0.9998</td><td>0.9982</td><td>0.9739</td><td>0.9510</td><td>0.9755</td><td>-0.0769</td></tr><tr><td>Coulomb_SY N3 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9992</td><td>0.9989</td><td>1.0003</td><td>0.9999</td><td>0.9995</td><td>0.9884</td><td>0.9759</td><td>0.9889</td><td>-0.0524</td></tr><tr><td>DeVore_KM N01 0.06</td><td>1.0000</td><td>0.9999</td><td>0.9920</td><td>0.9910</td><td>0.9889</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>1.0003</td><td>-0.0057</td></tr><tr><td>DeVore_KM N05 0.06</td><td>1.0000</td><td>0.9999</td><td>0.9919</td><td>0.9910</td><td>0.9890</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9996</td><td>1.0025</td><td>-0.0290</td></tr><tr><td>DeVore_KM N1 0.03</td><td>1.0000</td><td>1.0000</td><td>0.9959</td><td>0.9954</td><td>0.9943</td><td>0.9996</td><td>1.0000</td><td>0.9949</td><td>0.9992</td><td>1.0166</td><td>-0.0783</td></tr><tr><td>DeVore_KM N1 0.06</td><td>1.0000</td><td>0.9999</td><td>0.9920</td><td>0.9911</td><td>0.9892</td><td>0.9996</td><td>1.0000</td><td>0.9948</td><td>0.9990</td><td>1.0151</td><td>-0.0784</td></tr><tr><td>DeVore_KM N1 0.12</td><td>0.9999</td><td>0.9998</td><td>0.9844</td><td>0.9832</td><td>0.9806</td><td>0.9996</td><td>1.0000</td><td>0.9936</td><td>0.9981</td><td>1.0118</td><td>-0.0786</td></tr><tr><td>DeVore_KM N3 0.06</td><td>1.0000</td><td>0.9999</td><td>0.9919</td><td>0.9911</td><td>0.9892</td><td>0.9998</td><td>1.0000</td><td>0.9976</td><td>0.9994</td><td>1.0032</td><td>-0.0529</td></tr><tr><td>Coulomb_JT N01 0.06</td><td>0.9973</td><td>0.9446</td><td>0.8899</td><td>0.7118</td><td>1.0238</td><td>0.9975</td><td>0.9461</td><td>0.8913</td><td>0.7198</td><td>1.0361</td><td>-0.0060</td></tr><tr><td>Coulomb_JT N05 0.06</td><td>0.9973</td><td>0.9444</td><td>0.8901</td><td>0.7114</td><td>1.0237</td><td>0.9976</td><td>0.9470</td><td>0.8939</td><td>0.7240</td><td>1.0319</td><td>-0.0302</td></tr><tr><td>Coulomb_JT N1 0.03</td><td>0.9975</td><td>0.9444</td><td>0.8943</td><td>0.7087</td><td>1.0172</td><td>0.9981</td><td>0.9548</td><td>0.9074</td><td>0.7491</td><td>0.9971</td><td>-0.0810</td></tr><tr><td>Coulomb_JT N1 0.06</td><td>0.9974</td><td>0.9429</td><td>0.8918</td><td>0.7069</td><td>1.0231</td><td>0.9980</td><td>0.9543</td><td>0.9068</td><td>0.7506</td><td>1.0097</td><td>-0.0818</td></tr><tr><td>Coulomb_JT N1 0.12</td><td>0.9972</td><td>0.9403</td><td>0.8864</td><td>0.7043</td><td>1.0330</td><td>0.9979</td><td>0.9534</td><td>0.9028</td><td>0.7529</td><td>1.0290</td><td>-0.0833</td></tr><tr><td>Coulomb_JT N3 0.06</td><td>0.9975</td><td>0.9443</td><td>0.8924</td><td>0.7108</td><td>1.0232</td><td>0.9979</td><td>0.9518</td><td>0.9018</td><td>0.7391</td><td>1.0231</td><td>-0.0551</td></tr></tbody></table>
<table><caption>Table 8 (Continued)</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A</th><th colspan="3">H<sub>&nu;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th></tr></thead><tbody><tr><td>DeVore_GV N01 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9990</td><td>0.9988</td><td>0.9985</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>0.9991</td><td>-0.0057</td></tr><tr><td>DeVore_GV N05 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9981</td><td>0.9985</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9996</td><td>1.0014</td><td>-0.0290</td></tr><tr><td>DeVore_GV N1 0.03</td><td>1.0000</td><td>1.0000</td><td>0.9991</td><td>0.9979</td><td>0.9995</td><td>0.9996</td><td>1.0000</td><td>0.9950</td><td>0.9992</td><td>1.0164</td><td>-0.0782</td></tr><tr><td>DeVore_GV N1 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9980</td><td>0.9960</td><td>0.9983</td><td>0.9996</td><td>1.0000</td><td>0.9948</td><td>0.9990</td><td>1.0146</td><td>-0.0782</td></tr><tr><td>DeVore_GV N1 0.12</td><td>1.0000</td><td>0.9999</td><td>0.9950</td><td>0.9924</td><td>0.9943</td><td>0.9996</td><td>1.0000</td><td>0.9939</td><td>0.9983</td><td>1.0099</td><td>-0.0782</td></tr><tr><td>DeVore_GV N3 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9981</td><td>0.9965</td><td>0.9983</td><td>0.9998</td><td>1.0000</td><td>0.9977</td><td>0.9995</td><td>1.0035</td><td>-0.0527</td></tr><tr><td>DeVore_SA N01 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>1.0001</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9997</td><td>0.9998</td><td>-0.0057</td></tr><tr><td>DeVore_SA N05 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>1.0001</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9996</td><td>1.0021</td><td>-0.0290</td></tr><tr><td>DeVore_SA N1 0.03</td><td>1.0000</td><td>1.0000</td><td>0.9999</td><td>0.9999</td><td>1.0000</td><td>0.9996</td><td>1.0000</td><td>0.9950</td><td>0.9992</td><td>1.0164</td><td>-0.0785</td></tr><tr><td>DeVore_SA N1 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9994</td><td>0.9997</td><td>1.0001</td><td>0.9996</td><td>1.0000</td><td>0.9948</td><td>0.9990</td><td>1.0147</td><td>-0.0785</td></tr><tr><td>DeVore_SA N1 0.12</td><td>1.0000</td><td>1.0000</td><td>0.9980</td><td>0.9987</td><td>1.0005</td><td>0.9996</td><td>1.0000</td><td>0.9938</td><td>0.9982</td><td>1.0111</td><td>-0.0787</td></tr><tr><td>DeVore_SA N3 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9994</td><td>0.9997</td><td>1.0001</td><td>0.9998</td><td>1.0000</td><td>0.9977</td><td>0.9995</td><td>1.0027</td><td>-0.0528</td></tr><tr><td>Coulomb_GR N01 0.0</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9997</td><td>1.0001</td><td>0.9967</td><td>0.9523</td><td>0.8760</td><td>0.6852</td><td>1.0088</td><td>-0.0543</td></tr><tr><td>Coulomb_GR N05 0.0</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9997</td><td>1.0001</td><td>0.9969</td><td>0.9578</td><td>0.8797</td><td>0.6924</td><td>1.0071</td><td>-0.0731</td></tr><tr><td>Coulomb_GR N1 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>1.0001</td><td>0.9976</td><td>0.9706</td><td>0.8939</td><td>0.7218</td><td>0.9970</td><td>-0.1113</td></tr><tr><td>Coulomb_GR N1 0.12</td><td>1.0000</td><td>1.0000</td><td>0.9982</td><td>0.9988</td><td>1.0004</td><td>0.9978</td><td>0.9733</td><td>0.8999</td><td>0.7438</td><td>0.9974</td><td>-0.1119</td></tr><tr><td>Coulomb_GR N3 0.06</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9997</td><td>1.0001</td><td>0.9974</td><td>0.9682</td><td>0.8912</td><td>0.7141</td><td>1.0015</td><td>-0.0843</td></tr></tbody></table>
<table><caption>**Table 9** Same as Table 8 for the MHD-st (first page) and MHD-un (second and third pages) test cases</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="6">B vs. &nabla; &times; A</th><th rowspan="2">H<sub>&nu;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;N</th><th>&epsilon;M</th><th>&epsilon;E</th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;N</th><th>&epsilon;M</th><th>&epsilon;E</th></tr></thead><tbody><tr><td>Coulomb_SY 50</td><td>1.0000</td><td>1.0000</td><td>0.9986</td><td>0.9982</td><td>1.0005</td><td>0.9996</td><td>0.9994</td><td>0.9573</td><td>0.9432</td><td>1.0026</td><td>0.0360</td></tr><tr><td>Coulomb_SY 85</td><td>1.0000</td><td>1.0000</td><td>0.9987</td><td>0.9983</td><td>1.0006</td><td>0.9995</td><td>0.9994</td><td>0.9652</td><td>0.9502</td><td>1.0061</td><td>0.1138</td></tr><tr><td>Coulomb_SY 120</td><td>1.0000</td><td>1.0000</td><td>0.9987</td><td>0.9983</td><td>1.0005</td><td>0.9994</td><td>0.9994</td><td>0.9690</td><td>0.9545</td><td>1.0047</td><td>0.1664</td></tr><tr><td>Coulomb_SY 155</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9983</td><td>1.0004</td><td>0.9995</td><td>0.9995</td><td>0.9726</td><td>0.9585</td><td>1.0035</td><td>0.2112</td></tr><tr><td>Coulomb_SY 190</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9983</td><td>1.0004</td><td>0.9996</td><td>0.9995</td><td>0.9746</td><td>0.9612</td><td>1.0022</td><td>0.2463</td></tr><tr><td>DeVore_KM 50</td><td>1.0000</td><td>1.0000</td><td>0.9972</td><td>0.9974</td><td>0.9969</td><td>0.9995</td><td>0.9996</td><td>0.9707</td><td>0.9694</td><td>0.9838</td><td>0.0360</td></tr><tr><td>DeVore_KM 85</td><td>1.0000</td><td>1.0000</td><td>0.9973</td><td>0.9974</td><td>0.9969</td><td>0.9993</td><td>0.9996</td><td>0.9712</td><td>0.9678</td><td>0.9894</td><td>0.1141</td></tr><tr><td>DeVore_KM 120</td><td>1.0000</td><td>1.0000</td><td>0.9972</td><td>0.9973</td><td>0.9970</td><td>0.9992</td><td>0.9996</td><td>0.9720</td><td>0.9673</td><td>0.9896</td><td>0.1670</td></tr><tr><td>DeVore_KM 155</td><td>1.0000</td><td>1.0000</td><td>0.9972</td><td>0.9972</td><td>0.9970</td><td>0.9993</td><td>0.9996</td><td>0.9740</td><td>0.9674</td><td>0.9892</td><td>0.2121</td></tr><tr><td>DeVore_KM 190</td><td>1.0000</td><td>1.0000</td><td>0.9972</td><td>0.9971</td><td>0.9970</td><td>0.9993</td><td>0.9996</td><td>0.9748</td><td>0.9674</td><td>0.9880</td><td>0.2475</td></tr><tr><td>Coulomb_JT 50</td><td>0.9498</td><td>0.7208</td><td>0.5794</td><td>0.2706</td><td>1.0077</td><td>0.9769</td><td>0.6899</td><td>0.5553</td><td>0.1634</td><td>1.0097</td><td>0.0343</td></tr><tr><td>Coulomb_JT 85</td><td>0.9635</td><td>0.7396</td><td>0.6200</td><td>0.3131</td><td>1.0058</td><td>0.9824</td><td>0.7419</td><td>0.6611</td><td>0.2918</td><td>1.0127</td><td>0.1102</td></tr><tr><td>Coulomb_JT 120</td><td>0.9657</td><td>0.7541</td><td>0.6389</td><td>0.3434</td><td>1.0068</td><td>0.9851</td><td>0.7921</td><td>0.7239</td><td>0.4082</td><td>1.0132</td><td>0.1630</td></tr><tr><td>Coulomb_JT 155</td><td>0.9660</td><td>0.7616</td><td>0.6490</td><td>0.3566</td><td>1.0095</td><td>0.9866</td><td>0.8288</td><td>0.7601</td><td>0.4850</td><td>1.0146</td><td>0.2085</td></tr><tr><td>Coulomb_JT 190</td><td>0.9652</td><td>0.7641</td><td>0.6513</td><td>0.3565</td><td>1.0139</td><td>0.9874</td><td>0.8575</td><td>0.7834</td><td>0.5426</td><td>1.0159</td><td>0.2447</td></tr><tr><td>DeVore_GV 50</td><td>1.0000</td><td>1.0000</td><td>0.9963</td><td>0.9956</td><td>0.9951</td><td>0.9997</td><td>0.9996</td><td>0.9721</td><td>0.9696</td><td>0.9927</td><td>0.0361</td></tr><tr><td>DeVore_GV 85</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9959</td><td>0.9955</td><td>0.9997</td><td>0.9996</td><td>0.9747</td><td>0.9687</td><td>0.9930</td><td>0.1138</td></tr><tr><td>DeVore_GV 120</td><td>1.0000</td><td>1.0000</td><td>0.9967</td><td>0.9959</td><td>0.9963</td><td>0.9997</td><td>0.9996</td><td>0.9760</td><td>0.9685</td><td>0.9920</td><td>0.1665</td></tr><tr><td>DeVore_GV 155</td><td>1.0000</td><td>1.0000</td><td>0.9967</td><td>0.9957</td><td>0.9968</td><td>0.9997</td><td>0.9997</td><td>0.9774</td><td>0.9685</td><td>0.9915</td><td>0.2116</td></tr><tr><td>DeVore_GV 190</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9954</td><td>0.9971</td><td>0.9997</td><td>0.9997</td><td>0.9777</td><td>0.9684</td><td>0.9908</td><td>0.2470</td></tr><tr><td>DeVore_SA 50</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>1.0004</td><td>0.9997</td><td>0.9996</td><td>0.9839</td><td>0.9788</td><td>0.9878</td><td>0.0360</td></tr><tr><td>DeVore_SA 85</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9999</td><td>1.0005</td><td>0.9995</td><td>0.9995</td><td>0.9830</td><td>0.9774</td><td>0.9949</td><td>0.1140</td></tr><tr><td>DeVore_SA 120</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9999</td><td>1.0004</td><td>0.9994</td><td>0.9995</td><td>0.9835</td><td>0.9770</td><td>0.9964</td><td>0.1668</td></tr><tr><td>DeVore_SA 155</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>1.0003</td><td>0.9995</td><td>0.9996</td><td>0.9857</td><td>0.9778</td><td>0.9969</td><td>0.2119</td></tr><tr><td>DeVore_SA 190</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>1.0003</td><td>0.9996</td><td>0.9996</td><td>0.9868</td><td>0.9787</td><td>0.9966</td><td>0.2473</td></tr></tbody></table>

<figure><img src="image_21.png" /></figure>
<table><caption>**Table 9** (Continued)</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="6">B vs. &nabla; &times; A</th><th rowspan="2">H<sub>&lambda;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th></tr></thead><tbody><tr><td>Coulomb_SY 50</td><td>1.0000</td><td>1.0000</td><td>0.9983</td><td>0.9979</td><td>1.0006</td><td>0.9996</td><td>0.9994</td><td>0.9571</td><td>0.9415</td><td>0.9954</td><td>0.0114</td></tr><tr><td>Coulomb_SY 80</td><td>1.0000</td><td>1.0000</td><td>0.9983</td><td>0.9978</td><td>1.0008</td><td>0.9994</td><td>0.9992</td><td>0.9634</td><td>0.9471</td><td>0.9998</td><td>0.0356</td></tr><tr><td>Coulomb_SY 110</td><td>1.0000</td><td>1.0000</td><td>0.9983</td><td>0.9979</td><td>1.0006</td><td>0.9992</td><td>0.9992</td><td>0.9623</td><td>0.9410</td><td>0.9998</td><td>0.0606</td></tr><tr><td>Coulomb_SY 140</td><td>1.0000</td><td>1.0000</td><td>0.9985</td><td>0.9980</td><td>1.0005</td><td>0.9992</td><td>0.9992</td><td>0.9670</td><td>0.9414</td><td>0.9975</td><td>0.0837</td></tr><tr><td>Coulomb_SY 150</td><td>1.0000</td><td>1.0000</td><td>0.9985</td><td>0.9981</td><td>1.0005</td><td>0.9991</td><td>0.9990</td><td>0.9690</td><td>0.9542</td><td>0.9963</td><td>0.0913</td></tr><tr><td>Coulomb_SY 160</td><td>1.0000</td><td>1.0000</td><td>0.9987</td><td>0.9983</td><td>1.0005</td><td>0.9991</td><td>0.9991</td><td>0.9660</td><td>0.9532</td><td>0.9970</td><td>0.0809</td></tr><tr><td>Coulomb_SY 190</td><td>1.0000</td><td>1.0000</td><td>0.9988</td><td>0.9985</td><td>1.0005</td><td>0.9991</td><td>0.9992</td><td>0.9650</td><td>0.9517</td><td>0.9967</td><td>0.0862</td></tr><tr><td>DeVore_KM 50</td><td>1.0000</td><td>1.0000</td><td>0.9969</td><td>0.9973</td><td>0.9964</td><td>0.9995</td><td>0.9995</td><td>0.9713</td><td>0.9697</td><td>0.9840</td><td>0.0114</td></tr><tr><td>DeVore_KM 80</td><td>1.0000</td><td>1.0000</td><td>0.9968</td><td>0.9971</td><td>0.9964</td><td>0.9992</td><td>0.9994</td><td>0.9690</td><td>0.9661</td><td>0.9887</td><td>0.0355</td></tr><tr><td>DeVore_KM 110</td><td>1.0000</td><td>1.0000</td><td>0.9968</td><td>0.9970</td><td>0.9966</td><td>0.9989</td><td>0.9992</td><td>0.9662</td><td>0.9654</td><td>0.9885</td><td>0.0601</td></tr><tr><td>DeVore_KM 140</td><td>1.0000</td><td>1.0000</td><td>0.9969</td><td>0.9969</td><td>0.9968</td><td>0.9988</td><td>0.9992</td><td>0.9699</td><td>0.9671</td><td>0.9877</td><td>0.0820</td></tr><tr><td>DeVore_KM 150</td><td>1.0000</td><td>1.0000</td><td>0.9965</td><td>0.9964</td><td>0.9968</td><td>0.9987</td><td>0.9990</td><td>0.9654</td><td>0.9612</td><td>0.9832</td><td>0.0880</td></tr><tr><td>DeVore_KM 160</td><td>1.0000</td><td>1.0000</td><td>0.9963</td><td>0.9959</td><td>0.9968</td><td>0.9985</td><td>0.9986</td><td>0.9585</td><td>0.9533</td><td>0.9817</td><td>0.0782</td></tr><tr><td>DeVore_KM 190</td><td>1.0000</td><td>1.0000</td><td>0.9966</td><td>0.9963</td><td>0.9969</td><td>0.9985</td><td>0.9989</td><td>0.9626</td><td>0.9582</td><td>0.9811</td><td>0.0826</td></tr><tr><td>Coulomb_JT 50</td><td>0.9577</td><td>0.6606</td><td>0.5260</td><td>0.1438</td><td>1.0201</td><td>0.9844</td><td>0.6659</td><td>0.5489</td><td>0.0916</td><td>1.0049</td><td>0.0135</td></tr><tr><td>Coulomb_JT 80</td><td>0.9694</td><td>0.6525</td><td>0.5468</td><td>0.1176</td><td>1.0225</td><td>0.9873</td><td>0.7210</td><td>0.6565</td><td>0.2344</td><td>1.0111</td><td>0.0406</td></tr><tr><td>Coulomb_JT 110</td><td>0.9666</td><td>0.6500</td><td>0.5440</td><td>0.1262</td><td>1.0229</td><td>0.9863</td><td>0.6805</td><td>0.6754</td><td>0.1751</td><td>1.0113</td><td>0.0661</td></tr><tr><td>Coulomb_JT 140</td><td>0.9794</td><td>0.7475</td><td>0.6665</td><td>0.3743</td><td>1.0027</td><td>0.9908</td><td>0.7064</td><td>0.7771</td><td>0.3502</td><td>1.0017</td><td>0.0848</td></tr></tbody></table>

<figure><img src="image_22.png" /></figure>
<table><caption>Table 9 (Continued)</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. &nabla; &times; A<sub>p</sub></th><th colspan="3">B<sub>vs. &nabla; &times; A</sub></th><th rowspan="2">H<sub>&nu;</sub></th></tr><tr><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th><th>C<sub>Vec</sub></th><th>C<sub>CS</sub></th><th>&epsilon;<sub>N</sub></th><th>&epsilon;<sub>M</sub></th><th>&epsilon;<sub>E</sub></th></tr></thead><tbody><tr><td>Coulomb_JT 150</td><td>0.9486</td><td>0.6605</td><td>0.4911</td><td>0.0177</td><td>1.0546</td><td>0.9768</td><td>0.7915</td><td>0.6735</td><td>0.2263</td><td>1.0220</td><td>0.0799</td></tr><tr><td>Coulomb_JT 160</td><td>0.9513</td><td>0.7466</td><td>0.5273</td><td>0.1055</td><td>1.0873</td><td>0.9760</td><td>0.8024</td><td>0.6717</td><td>0.3748</td><td>1.0391</td><td>0.0821</td></tr><tr><td>Coulomb_JT 190</td><td>0.9955</td><td>0.9477</td><td>0.8599</td><td>0.7440</td><td>1.0001</td><td>0.9964</td><td>0.9748</td><td>0.9002</td><td>0.8177</td><td>0.9985</td><td>0.0829</td></tr><tr><td>DeVore_GV 50</td><td>0.9999</td><td>1.0000</td><td>0.9957</td><td>0.9949</td><td>0.9934</td><td>0.9997</td><td>0.9995</td><td>0.9728</td><td>0.9701</td><td>0.9958</td><td>0.0116</td></tr><tr><td>DeVore_GV 80</td><td>0.9999</td><td>1.0000</td><td>0.9958</td><td>0.9948</td><td>0.9933</td><td>0.9997</td><td>0.9995</td><td>0.9710</td><td>0.9667</td><td>0.9958</td><td>0.0354</td></tr><tr><td>DeVore_GV 110</td><td>1.0000</td><td>1.0000</td><td>0.9958</td><td>0.9948</td><td>0.9943</td><td>0.9996</td><td>0.9994</td><td>0.9714</td><td>0.9672</td><td>0.9945</td><td>0.0599</td></tr><tr><td>DeVore_GV 140</td><td>1.0000</td><td>1.0000</td><td>0.9958</td><td>0.9946</td><td>0.9953</td><td>0.9996</td><td>0.9993</td><td>0.9753</td><td>0.9693</td><td>0.9943</td><td>0.0817</td></tr><tr><td>DeVore_GV 150</td><td>1.0000</td><td>1.0000</td><td>0.9954</td><td>0.9941</td><td>0.9955</td><td>0.9995</td><td>0.9992</td><td>0.9703</td><td>0.9632</td><td>0.9910</td><td>0.0877</td></tr><tr><td>DeVore_GV 160</td><td>1.0000</td><td>1.0000</td><td>0.9955</td><td>0.9940</td><td>0.9957</td><td>0.9994</td><td>0.9987</td><td>0.9635</td><td>0.9552</td><td>0.9899</td><td>0.0780</td></tr><tr><td>DeVore_GV 190</td><td>1.0000</td><td>1.0000</td><td>0.9957</td><td>0.9944</td><td>0.9960</td><td>0.9994</td><td>0.9990</td><td>0.9676</td><td>0.9600</td><td>0.9900</td><td>0.0825</td></tr><tr><td>DeVore_SA 50</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9999</td><td>1.0006</td><td>0.9996</td><td>0.9994</td><td>0.9830</td><td>0.9782</td><td>0.9859</td><td>0.0115</td></tr><tr><td>DeVore_SA 80</td><td>1.0000</td><td>1.0000</td><td>0.9994</td><td>0.9999</td><td>1.0008</td><td>0.9993</td><td>0.9993</td><td>0.9795</td><td>0.9742</td><td>0.9922</td><td>0.0355</td></tr><tr><td>DeVore_SA 110</td><td>1.0000</td><td>1.0000</td><td>0.9994</td><td>0.9999</td><td>1.0006</td><td>0.9991</td><td>0.9991</td><td>0.9748</td><td>0.9701</td><td>0.9927</td><td>0.0601</td></tr><tr><td>DeVore_SA 140</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9999</td><td>1.0005</td><td>0.9990</td><td>0.9991</td><td>0.9773</td><td>0.9701</td><td>0.9915</td><td>0.0822</td></tr><tr><td>DeVore_SA 150</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9998</td><td>1.0005</td><td>0.9989</td><td>0.9992</td><td>0.9756</td><td>0.9698</td><td>0.9903</td><td>0.0892</td></tr><tr><td>DeVore_SA 160</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9998</td><td>1.0005</td><td>0.9989</td><td>0.9980</td><td>0.9731</td><td>0.9654</td><td>0.9902</td><td>0.0786</td></tr><tr><td>DeVore_SA 190</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9998</td><td>1.0004</td><td>0.9989</td><td>0.9990</td><td>0.9759</td><td>0.9714</td><td>0.9987</td><td>0.0831</td></tr></tbody></table>
<table><caption>**Table 10** Same as Table 8 for the MHD-st-div(B) test cases</caption><thead><tr><th rowspan="2">Label</th><th colspan="3">B<sub>p</sub> vs. ∇ × A<sub>p</sub></th><th colspan="3">B<sub>vs. ∇ × A</sub></th><th rowspan="2">H<sub>γ</sub></th></tr><tr><th>C<sub>vec</sub></th><th>C<sub>CS</sub></th><th>1 − E<sub>n</sub></th><th>ε<sub>E</sub></th><th>C<sub>vec</sub></th><th>C<sub>CS</sub></th><th>1 − E<sub>n</sub></th><th>1 − E<sub>m</sub></th><th>ε<sub>E</sub></th></tr></thead><tbody><tr><td>Coulomb_SY divB at 0.2 %</td><td>1.0000</td><td>1.0000</td><td>0.9985</td><td>0.9981</td><td>1.0005</td><td>0.9990</td><td>0.9979</td><td>0.9697</td><td>0.9694</td><td>0.9854</td><td>0.0356</td></tr><tr><td>Coulomb_SY divB at 1.1 %</td><td>1.0000</td><td>1.0000</td><td>0.9986</td><td>0.9981</td><td>1.0005</td><td>0.9996</td><td>0.9991</td><td>0.9727</td><td>0.9656</td><td>0.9970</td><td>0.0359</td></tr><tr><td>Coulomb_SY divB at 2.0 %</td><td>1.0000</td><td>1.0000</td><td>0.9986</td><td>0.9982</td><td>1.0005</td><td>0.9995</td><td>0.9990</td><td>0.9455</td><td>0.9272</td><td>1.0051</td><td>0.0361</td></tr><tr><td>Coulomb_SY divB at 4.0 %</td><td>1.0000</td><td>1.0000</td><td>0.9986</td><td>0.9984</td><td>1.0005</td><td>0.9976</td><td>0.9951</td><td>0.8747</td><td>0.8375</td><td>1.0137</td><td>0.0367</td></tr><tr><td>Coulomb_SY divB at 8.2 %</td><td>1.0000</td><td>1.0000</td><td>0.9987</td><td>0.9986</td><td>1.0005</td><td>0.9983</td><td>0.9877</td><td>0.7532</td><td>0.6784</td><td>1.0048</td><td>0.0380</td></tr><tr><td>Coulomb_SY divB at 14.4 %</td><td>1.0000</td><td>1.0000</td><td>0.9989</td><td>0.9991</td><td>1.0006</td><td>0.9529</td><td>0.9663</td><td>0.5348</td><td>0.2345</td><td>0.8928</td><td>0.0440</td></tr><tr><td>DeVore_KM divB at 0.2 %</td><td>0.9994</td><td>0.9985</td><td>0.9892</td><td>0.9879</td><td>0.9984</td><td>0.9986</td><td>0.9964</td><td>0.9708</td><td>0.9725</td><td>0.9936</td><td>0.0358</td></tr><tr><td>DeVore_KM divB at 1.1 %</td><td>0.9999</td><td>0.9997</td><td>0.9945</td><td>0.9942</td><td>0.9971</td><td>0.9996</td><td>0.9988</td><td>0.9805</td><td>0.9805</td><td>0.9873</td><td>0.0359</td></tr><tr><td>DeVore_KM divB at 2.0 %</td><td>0.9999</td><td>0.9999</td><td>0.9959</td><td>0.9960</td><td>0.9968</td><td>0.9993</td><td>0.9987</td><td>0.9606</td><td>0.9567</td><td>0.9821</td><td>0.0360</td></tr><tr><td>DeVore_KM divB at 4.0 %</td><td>0.9991</td><td>0.9983</td><td>0.9862</td><td>0.9852</td><td>0.9983</td><td>0.9958</td><td>0.9785</td><td>0.9031</td><td>0.8819</td><td>0.9737</td><td>0.0365</td></tr><tr><td>DeVore_KM divB at 8.2 %</td><td>0.9941</td><td>0.9946</td><td>0.9676</td><td>0.9662</td><td>1.0083</td><td>0.9794</td><td>0.9107</td><td>0.8044</td><td>0.7545</td><td>0.9642</td><td>0.0380</td></tr><tr><td>DeVore_KM divB at 14.4 %</td><td>0.9635</td><td>0.9865</td><td>0.9277</td><td>0.9249</td><td>1.0765</td><td>0.9196</td><td>0.8895</td><td>0.6264</td><td>0.4642</td><td>0.9778</td><td>0.0475</td></tr><tr><td>Coulomb_JT divB at 0.2 %</td><td>0.9996</td><td>0.9973</td><td>0.9643</td><td>0.9381</td><td>1.0102</td><td>0.9984</td><td>0.9956</td><td>0.9524</td><td>0.9375</td><td>0.9925</td><td>0.0357</td></tr><tr><td>Coulomb_JT divB at 1.1 %</td><td>0.9990</td><td>0.9927</td><td>0.9413</td><td>0.8966</td><td>0.9945</td><td>0.9989</td><td>0.9915</td><td>0.9384</td><td>0.8945</td><td>0.9975</td><td>0.0356</td></tr><tr><td>Coulomb_JT divB at 2.0 %</td><td>0.9512</td><td>0.8492</td><td>0.5392</td><td>0.2158</td><td>1.2367</td><td>0.9707</td><td>0.8307</td><td>0.4435</td><td>−0.0145</td><td>1.1233</td><td>0.0382</td></tr><tr><td>Coulomb_JT divB at 4.0 %</td><td>0.9985</td><td>0.9847</td><td>0.9554</td><td>0.9168</td><td>0.9944</td><td>0.9967</td><td>0.9617</td><td>0.8717</td><td>0.8182</td><td>1.0125</td><td>0.0367</td></tr><tr><td>Coulomb_JT divB at 8.2 %</td><td>0.9866</td><td>0.9378</td><td>0.8266</td><td>0.7542</td><td>0.9362</td><td>0.9838</td><td>0.9029</td><td>0.7600</td><td>0.6830</td><td>0.9685</td><td>0.0373</td></tr><tr><td>Coulomb_JT divB at 14.4 %</td><td>0.9039</td><td>0.8344</td><td>0.5738</td><td>0.5227</td><td>0.7202</td><td>0.9331</td><td>0.7435</td><td>0.5131</td><td>0.2749</td><td>0.7689</td><td>0.0395</td></tr><tr><td>DeVore_GV divB at 0.2 %</td><td>0.9976</td><td>0.9963</td><td>0.9834</td><td>0.9812</td><td>1.0007</td><td>0.9989</td><td>0.9964</td><td>0.9756</td><td>0.9752</td><td>1.0047</td><td>0.0359</td></tr><tr><td>DeVore_GV divB at 1.1 %</td><td>0.9996</td><td>0.9991</td><td>0.9937</td><td>0.9933</td><td>0.9961</td><td>0.9998</td><td>0.9988</td><td>0.9835</td><td>0.9822</td><td>0.9965</td><td>0.0360</td></tr><tr><td>DeVore_GV divB at 2.0 %</td><td>0.9999</td><td>0.9997</td><td>0.9963</td><td>0.9961</td><td>0.9952</td><td>0.9995</td><td>0.9987</td><td>0.9624</td><td>0.9574</td><td>0.9911</td><td>0.0361</td></tr><tr><td>DeVore_GV divB at 4.0 %</td><td>0.9964</td><td>0.9955</td><td>0.9797</td><td>0.9777</td><td>1.0017</td><td>0.9959</td><td>0.9780</td><td>0.9077</td><td>0.8835</td><td>0.9861</td><td>0.0366</td></tr><tr><td>DeVore_GV divB at 8.2 %</td><td>0.9760</td><td>0.9888</td><td>0.9501</td><td>0.9484</td><td>1.0438</td><td>0.9792</td><td>0.9110</td><td>0.8143</td><td>0.7566</td><td>0.9963</td><td>0.0380</td></tr><tr><td>DeVore_GV divB at 14.4 %</td><td>0.8687</td><td>0.9772</td><td>0.8864</td><td>0.8885</td><td>1.3247</td><td>0.9232</td><td>0.8897</td><td>0.6463</td><td>0.4670</td><td>1.0873</td><td>0.0475</td></tr><tr><td>DeVore_SA divB at 0.2 %</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9999</td><td>1.0004</td><td>0.9990</td><td>0.9974</td><td>0.9759</td><td>0.9780</td><td>0.9872</td><td>0.0356</td></tr><tr><td>DeVore_SA divB at 1.1 %</td><td>1.0000</td><td>1.0000</td><td>0.9997</td><td>0.9999</td><td>1.0004</td><td>0.9996</td><td>0.9992</td><td>0.9861</td><td>0.9853</td><td>0.9885</td><td>0.0358</td></tr><tr><td>DeVore_SA divB at 2.0 %</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9999</td><td>1.0004</td><td>0.9998</td><td>0.9986</td><td>0.9773</td><td>0.9676</td><td>0.9869</td><td>0.0361</td></tr><tr><td>DeVore_SA divB at 4.0 %</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9999</td><td>1.0004</td><td>0.9978</td><td>0.9842</td><td>0.9389</td><td>0.9174</td><td>0.9777</td><td>0.0366</td></tr><tr><td>DeVore_SA divB at 8.2 %</td><td>1.0000</td><td>1.0000</td><td>0.9996</td><td>0.9998</td><td>1.0004</td><td>0.9886</td><td>0.9694</td><td>0.8724</td><td>0.8678</td><td>0.9393</td><td>0.0378</td></tr><tr><td>DeVore_SA divB at 14.4 %</td><td>1.0000</td><td>1.0000</td><td>0.9995</td><td>0.9998</td><td>1.0005</td><td>0.9542</td><td>0.9256</td><td>0.7519</td><td>0.7180</td><td>0.8068</td><td>0.0429</td></tr></tbody></table>

<figure><img src="image_23.png" /></figure>
References

C.E. Alissandrakis, On the computation of constant alpha force-free 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/20122078

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 force-free 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 active-region solar corona. II. Nonlinear force-free 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/138

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 X-class 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 self-helicity in coronal magnetic fields. Astrophys. J. **674**, 1130–1143 (2008). doi:10.1086/524011

B.C. Low, Y.Q. Lou, Modeling solar force-free 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 ejec-
tions. 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 solar-like 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, Gauge-invariant helicity for force-free 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. Wilmot-Smith, 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 force-free 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 energy-helicity 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 relative-helicity 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 Titov-Dé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 force-free 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 force-free 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

<figure><img src="image_24.png" /></figure>