Magnetic helicity estimations in models and observations of the solar magnetic field. Part I: Finite volume methods

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 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.

is the helicity of the vector field B = ∇ × A in a given volume V, with A(x, t) representing the corresponding space-and time-dependent vector potential. If the field consists of a discrete collection of flux tubes, H is the winding number (Moffatt 1969) expressing their degree of mutual linkage. By extension, Eq. (1) is a measure of the entangling ("knottedness") of the field's streamlines.
For a given vector potential A, the addition of the gradient of a (sufficiently regular but otherwise arbitrary) scalar function, i.e., the transformation A −→ A + ∇ψ, does not change the resulting B. This property of the definition of B is called gaugeinvariance. Due to this freedom in the gauge, H is not uniquely defined, since where dS = dSn, with dS being the infinitesimal element of the bounding surface ∂V of the volume V, andn the outward-oriented normal to ∂V. Hence, H is not gaugeinvariant unless two conditions are met: First, the vector field B must be solenoidal, as implied by its definition as curl of A, and, second, the volume's bounding surface ∂V must be a flux surface of B, i.e., (B ·n)| ∂V = 0. When applied to a magnetic field B, the solenoidal requirement is satisfied by virtue of Maxwell equations, although possibly only to a finite extent in numerical experiments, and ∂V is a flux surface if no magnetic field line is threading the boundary. This latter requirement is rarely satisfied in natural systems, and makes Eq. (1) of little interest for practical use.
On the other hand, helicity has the fundamental property of being strictly conserved in ideal magneto-hydrodynamics (MHD, Woltjer 1958). Since MHD evolution in the absence of dissipation preserves the topology of the magnetic field, the field lines' winding numbers, or, more generally, magnetic helicity, cannot be changed during the evolution. Even more appealing is the fact that magnetic helicity, contrary to magnetic energy, is very well conserved in 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)

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 of a magnetic field B with respect to a reference magnetic field B p = ∇ × 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 B p . This choice has the following motivation: in order for H V in Eq.
(3) to be gauge invariant, the input and potential fields, B and B p , respectively, must be solenoidal, and such that (n · B) ∂V = (n · B p ) ∂V . (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, ∂V needs not to be a flux surface for 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 where The two definitions in Eqs. (6, 7) have the property of being separately gaugeinvariant under the same assumptions guaranteeing the gauge-invariance of H V . The first term H V,J corresponds to the general definition of helicity Eq. (1), but this time of the current-carrying part of the field only, (B−B p ) = ∇×(A−A p ). By construction, such a field has no normal component on the boundary, i.e., has ∂V as a flux surface. The second term H V,JP has no intuitive interpretation, but is a sort of mutual helicity that basically takes care of the flux threading ∂V (via the transverse component of A p ), and is gauge-invariant because only the current-carrying part of B appears. The conservation properties of the relative magnetic helicity were numerically tested by Pariat et al. (2015), confirming that helicity is a very well conserved quantity even in presence of very strong dissipation. The equation regulating the relative magnetic helicity variation rate due to dissipation and flux through the boundaries is also derived by Pariat et al. (2015). In the particular case of the simulation of a jet eruption examined in that article, relative helicity is conserved more than one order of magnitude better than free energy.
In the particular case of applications to the solar atmosphere, one additional complication is that the magnetic field cannot be measured in the solar corona with the resolution necessary for the computation of magnetic helicity. The magnetic field is instead inferred by inverting spectropolarimetric measurements of emission from lower, 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. 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.

Finite volume (FV)
-Requires B in V e.g., from MHD simulations or NLFFF -Compute H V at one time -May employ different gauges (see Table 2)

Helicity-flux integration (FI)
-Requires time evolution of vector field on ∂V -Requires knowledge or model of flows on ∂V -Valid for a specific set of gauge and assumptions, see  Discrete flux-tubes (DT)

Connectivity-based (CB)
-Requires the vector field on photosphere at one time -Models the corona connectivity as a collection of M force-free flux tubes -Minimal connection length principle 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 Eq.
(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.

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); . 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.
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);; 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 . 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 α ≡ (∇ × B)/B. More specifically, the CB method takes in input the photospheric observations and models the coronal field as a single (linear) or a collection of (nonlinear) flux tube(s) as an integral step of the helicity computation itself. The set of flux tubes is obtained assuming that the line connectivity is, globally, the shortest possible, mimicking the compact character of the more 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  we mostly focus on FI methods and the comparison with the CB method. Subsequent papers are dedicated to the TN method , 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-16 period 1 . As a reference for future testing, we make available the data cube of each 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.

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 , 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 mehtods, 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;Tziotziou et al. , 2013. Benchmarking DT methods together with the FV ones enables the reader to better interpret results of past and future studies applying such methods. From a different point of view, the TN method is included because it requires the same information as the other finite volume methods, i.e., the full knowledge of the magnetic field in the entire considered volume. The CB method is included because, despite requiring only the photospheric vector magnetogram, it can use any available information of the coronal connectivity. Also, similarly to FV methods, the 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 magnetic fields used as tests are discussed in Sect. 4. Section 3 summarizes the diagnostic tools used for the comparing the results. The main results of the comparison are then discussed in Sect. 5, with specific discussion of the dependence on resolution presented in Sect. 6 and of the dependence on the solenoidal property given in Sect. 7. An overview of the FV methods results is given in Sect. 8, whereas the the comparative tests with DT methods are presented in Sect. 9. Finally our conclusions are summarized in Sect. 10.

Helicity computation methods
Finite volume methods require the knowledge of the magnetic field B in the entire volume V, and differ essentially in the way in which the vector potentials are computed. The methods presented here compute vector potentials employing either the Coulomb gauge (∇ · A = 0) or the DeVore gauge (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 B p = ∇φ, with φ being the scalar potential, solution of such that the constraint Eq. (4) is satisfied. Errors in solving Eqs. (8,9) are a first source of inaccuracy for the methods.

Methods employing the Coulomb gauge
Vector potentials in the Coulomb gauge satisfy for the vector potential of the potential field, and for the vector potential of the input field, where J = ∇ × B. The conditions Eq. (12) and Eq. (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 A p and A.
From the computational point of view, the numerical effort required to solve for the vector potentials consists, in general, in the solutions of Eqs. (10 -12) and Eqs. (13 -15), i.e., of six 3D Poisson/Laplace problems, one for each Cartesian component of the vector potentials A p and A.

Coulomb JT
In order to solve Eqs. (10 -12) and Eqs. (13 -15), appropriate additional boundary conditions for A and 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 A into a current-carrying and a potential (current-free) part, in the form A = A c + A p . The reproduction of the input fields' flux at the volume's boundaries, ∂V, is entirely dedicated to A p (obeying Eqs. (10 -12)). The electric current distribution in V and on ∂V, on the other hand, are delivered by A c (Eqs. (13,14) with A c instead of A and Eq. (15) replaced byn · (∇ × A c )| ∂V = 0).
In particular, the tangential components of 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, φ f , in the form A p,t = −n × ∇ t φ f , where ∇ 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 ∇ 2 t φ f = −n · B, for which boundary conditions on the four edges of each face f need to be specified. This approach of defining A p on ∂V is in principle used by all Coulomb methods considered in the present study, but the specific way in which the 2D Laplace problems are formulated is different. Thalmann et al. (2011) use Neumann conditions in the form ∂ n φ f = c f , where c f is a constant along a particular face and ∂ n is the derivative in the direction normal to the edge of the face f (see their Sect. 2.1 for details). The different c f are constructed in such a way that the total outflow through the volume's bounding surface ∂V is minimized. In this way, a vanishing tangential divergence (∇ t · A p,t = 0) is enforced on ∂V, and following Gauss' theorem, Eqs. (10 -12) are approximately fulfilled.
Another difference of the applied Coulomb methods is how the current-carrying vector potential A c is calculated and its solenoidality enforced. Thalmann et al. (2011) solves Eq. (13) for A c numerically, similar to Yang et al. (2013a), just with differing boundary conditions. In the Coulomb JT case, ∇ t · A c,t = 0 on ∂V is explicitly enforced in order to fulfill the Coulomb gauge for 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 R Mathematical Kernel Library.

Coulomb SY
The Coulomb SY method is described in Yang et al. (2013a); Yang et al. (2013b). 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 Eqs. (13, 15), the Coulomb SY method additionally enforces the boundary condition (n · A p )| ∂V = (n · A)| ∂V = 0 at all boundaries. Then, the transverse vector potential at the boundaries and the vector potential at the edges is obtained by using Gauss theorem. After obtaining the boundary values, Yang et al. (2013a); Yang et al. (2013b) firstly resolve the Laplace Eq. (10) and the Poisson Eq. (13) to obtain an initial guess of the solution, A p and A. These preliminary solutions satisfies Eq. (12) and Eq. (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 R (International Mathematics and Statistics Library).

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 = ∇ × A p , rather than using Eqs. (8,9). The corresponding boundary value problem, Eq. (10), is solved with the constraint Eq. (11) and the boundary condition (A p ·n)| ∂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 m p = ∇ × A m p 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 m p has zero total flux on each face independently and can be correctly used as a boundary condition for the sub-problemŝ where A f i p 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 f i and the vector potential of a compensation field, A m , as 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 Eqs. (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 In contrast to other methods, the solutions of the Laplace and Poison problems for the A f i 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.

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 B in the finite volume V = where the integration function b(x, y) = A(z = z 2 ) obeys to and b z = 0. The particular solution of Eq. (20) employed here is but see Valori et al. (2012) for alternative options. The above equations are applied in the computation of the vector potential of the potential field too by substituting B with 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 where Eq. (20) 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 Eqs. (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.

DeVore GV
DeVore GV is the original implementation described in Valori et al. (2012), where the requirement Eq. (24) is enforced by defining the z-integral operator as the numerical inverse operation of the second order central differences operator, see Section 4.2 in Valori et al. (2012). The Poisson problem for the determination of the scalar potential φ in Eqs. (8,9) is solved numerically using the Helmholtz solver in the proprietary Intel R 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.

DeVore KM
DeVore KM is described in . This implementation has two main differences with the one of DeVore GV. The first is in the solver of Laplace's equation. DeVore KM uses the routine HW3CRT that is included in the freely available FISHPACK library (Swartztrauber and Sweet 1979). A test, however, with the corresponding Intel MKL solver revealed minor differences in the solutions obtained with the two routines, and a factor of ≤ 2 more computational time required by the FISHPACK solver. The second and most important difference with the DeVore GV method is in the numerical calculation of integrals and derivatives in Eqs. (19 -22). In DeVore KM integrations are made with the modified Simpson's rule of error estimate 1/N 4 (Press et al. 1992), with N being the number of integration points, and, in the special case N = 2, with the trapezoidal rule instead. In addition, differentiations are made using the appropriate (centered, forward or backward) 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.

DeVore SA
DeVore SA follows the general scheme of the DeVore GV method with two differences. The first one is that the Eqs. (8,9) for the potential φ are solved in Fourier space separately for all faces of the box. In particular, the problem is divided into six sub-problems using where φ c is the 3D scalar potential of the compensation field B c = ∇φ c and φ 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 φ 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 Eqs. (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 which can be re-casted as where we have defined Taking the x-and y-derivatives of Eq. (29), and using Eq. (20) and 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. 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:

Discrete flux-tubes methods
where T i denotes the twist and writhe of magnetic flux tube i with flux Φ i , and L i, j is the linking number between two magnetic flux tubes i and j with fluxes Φ i and Φ j , respectively. The first and second term on the right hand side of Eq. (31) represents the self and mutual helicity, respectively. With the approximation of discrete magnetic flux tubes, the physical quantity of the magnetic helicity is related to the topological concept of the writhe, twist and linking number of curves, and the magnetic flux associated with those curves. The formulae of computing these topological quantities for both close and open curves have been derived in Berger and Prior (2006) and Démoulin et al. (2006). For the purpose of our comparison it must be noticed that discrete 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.

Twist-number
The TN method is described in Guo et al. (2010) and Guo et al. (2013). This method is aimed at computing the helicity of a highly twisted magnetic structure, such as a magnetic flux rope. A magnetic flux rope is considered as an isolated, single flux tube such that only the self magnetic helicity is computed. The helicity contributed by the writhe is also omitted assuming that the flux rope is not highly kinked. With these two assumptions, the magnetic helicity of a single highly twisted structure is simplified as where T is the twist number of the considered magnetic flux rope with flux Φ. In order to estimate T , the formula derived in Berger and Prior (2006) to compute the twist number of a sample curve referred to an axis is employed. Practically, the axis can be determined by the symmetry of a magnetic configuration or by other assumptions, such as requiring it to be horizontal and to follow the polarity inversion line (Guo et al. 2010). The boundary of the flux rope is determined by the quasi-separatrix layer (QSL) that is found to wrap the flux rope (Guo et al. 2013). Then the twist density, dT /ds, at an arc length s is: Two unit vectors are used in Eq. (33): T(s), that is tangent to the axis curve, and V(s), that is normal to T and pointing from the axis curve to the sample curve. By integrating the equation along the axis curve the total twist number is derived. Equation (33) is suitable for smooth curves in arbitrary geometries without self intersection. Since it makes no assumption about the magnetic field, it can be applied to both force-free and non-force-free magnetic field models.

Connectivity-based
The CB method was introduced by  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  and Tziotziou et al. ( , 2013. The connectivity matrix in a collection of partitions of both polarities will reveal a number of M discrete, assumed slender, flux tubes with flux contents Φ i ; i ≡ {1, ..., M}. The respective force-free parameters α 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 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 , α i , and the relative positions of each flux tube's footpoints,  showed that a lower limit of the free magnetic energy for a collection of M flux tubes is where A, δ are known fitting constants, λ is the length element (the pixel size in observed photospheric magnetograms), and 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 fluxtube 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 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 c CB in Eq. (34). Therefore, the above E c CB 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, L lm is the mutual-helicity factor for a collection of M p 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 m CB;mutual = 0, in case all α i ; i ≡ {1, ..., M} are zero within uncertainties δα i , which are fully defined in this analysis. More generally, one may use the "energy-helicity diagram" correlation of

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 Sect. 2.1 and Sect. 2.2. Since the helicity of B defined by Eq. (3) involves the corresponding potential field B p , as well as the vector potentials for B p and B, this basically implies providing a quantitative estimation of the accuracy of such fields. To that purpose, we introduce normalized quantities and metrics as follows: For each discretized magnetic field, we define H V as the helicity defined in Eq. (3) normalized to Φ 2 , where 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 Φ 2 . We reserve the calligraphic symbol 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 discretised 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.

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 which are, respectively, the complement of the normalized vector error and the energy ratio, introduced by Schrijver et al. (2006) , 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, centraldifference discretization scheme for points in the interior of V. Values on the volumebounding surface, ∂V, are taken from the input magnetic fields.

Quantification of the solenoidal property
The test cases used in this article must have a value of ∇ · B small enough to be considered numerically solenoidal. In order to quantify the level of solenoidality, we apply here the decomposition of the energy of the magnetic field in V into solenoidal and nonsolenoidal contributions, as in Valori et al. (2013). Using that decomposition, a fraction E div of the 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, through the surface ∂v of an elementary volume v including the node i. The rightmost expression in Eq. (40) holds for a grid of uniform and homogeneous resolution ∆. Therefore, it may be appropriately used as a metric for the methods analysed in this study, since they all are based on uniform Cartesian grids. The smaller the value of | f i | , the more solenoidal the field. However, the actual value of this metric depends on the number of grid points n and the resolution ∆, therefore it makes most sense to apply it to identical discretized volumes. In addition, this metric is used when the energy one is not applicable, for instance, when checking the solenoidal property of vector potentials in Coulomb gauge methods, | f i (A)| , as in Sect. 7.

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 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. 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.

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-β 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, φ = π/4 in the notation of Low and Lou (1990)), in the same V =[-1,1]×[-1,1]×[0,2] volume. From this solution, four test cases are constructed where the above volume is discretized using 32, 64, 128, and 256 nodes per side. Since the solution of the ordinary differential equation that defines the LL field is always the same in the four cases, the only factor changing among the different LL cases is the resolution. As a test for helicity methods, LL represents a large-scale, force-free field with large-scale smooth currents distributed 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 ≃ 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.

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 an hyperbolic flux tube (Titov et al. 2002). In this article we consider six realizations of that solution, in different combinations of twist and spatial resolution. Unless differently 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]×[-5.10, 5.10]×[0.00, 4.56], the distance between the charges and the current ring centre is L = 0.83, the depth of the current ring centre 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 ∆ = 0.06 for all four cases, i.e., given by a grid of 107×171×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 ≃ 0.05 Table 3 Parameters of the TD twist test cases, with resolution ∆ = 0.06; N is the approximate end-to end number of turns, a and d are the minor radius and the depth of the centre 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 ∆ = 0.06 case of Table 4. We refer to Valori et al. (2010) for definitions and normalizations of the TD parameters.
N a/d q × 10 12 Twist/π 3 0.31 100 -6.054 1 0.80 100 -2.114 0.5 0.80 29.5 -1.004 0.1 0.80 5.5 -0.201 case where the twist is lower and the field is almost potential (E free ≃ 0.1%). In this case solenoidal errors are slightly larger than the free energy, but anyway extremely small (E div = 0.2%). Table 4 Parameters of the TD resolution test cases, with N ≃ 1; ∆ 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 ∆ = 0.06 here is the same the N = 1 case in Table 3. We refer to Valori et al. (2010)  TD-resolution test cases. The three equilibria for this test are the N1 case in Table 4 at resolution ∆ = 0.06, and two additional equilibria with exactly the same parameters but with ∆ = 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 potential 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 between the two flux systems. For such reasons, when employed as an initial state of numerical 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 relative 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.

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 Leake et al. (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 snapshots 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 0 , see Leake et al. (2013) for more details). 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 photospheric 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, corresponding 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 free , E p , and hence E are all dynamically evolving in time. In the unstable case, the small variation of the ratio E free (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 Sect. 4.4 and Sect. 7.
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.

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 V 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(B) test (hereafter, MHD-st-div(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 B at t = 50 of the MHD-st case, a solenoidal field B s is produced by removing the divergence part B ns of B. A parametric, generally nonsolenoidal field, is then constructed by adding the divergence part back, multiplied by a scalar amplitude δ. In this way, the original spatial structure of B ns is kept in B δ but its amplitude is modified according to the chosen value of δ. The case δ = 0 corresponds the numerically solenoidal field B s , whereas δ = 1 reproduces the original test field B, i.e., the snapshot of MHD-st at t = 50. The method for building the vector field B ns is detailed in Sect. 7.1 of Valori et al. (2013). By increasing the value of δ, progressively more divergence can be added, yielding a larger amplitude of the nonsolenoidal component. A trial-and-error tuning of δ resulted into six MHD-st-div(B) test fields with nonsolenoidal contributions as a function of the parameter δ as summarized in Table 5 and Fig. 2f. In the following, different MHD-st-div(B) test cases are identified by the correspondent fraction of E div . The dependence of FV methods on solenoidal errors as modeled by Eq. (41) is discussed in Sect. 7.

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. 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.

Dependence on twist in the Titov and Démoulin case
In order to address specifically the methods' accuracy, we show in Fig. 3b the complement of the normalized vector error, ǫ N (B, ∇ × A) (cf. Eq. 38), between the test field B and the rotation of the corresponding vector potential A computed by each method. The DeVore methods sport extremely high accuracy metric of ǫ N = 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 proportionally larger values of E div , 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., ǫ N = 0.976 at N = 1) that are only slightly worse than the DeVore methods. Similarly to DeVore-gauge methods, the Coulomb SY accuracy 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 ǫ N = 0.88. A trend similar to Coulomb GR's one is found for the Coulomb JT method but with a slightly smaller error (ǫ N = 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 ǫ N ) for smaller values of E div . This counterintuitive trend, however, is not confirmed by the tests in Sect. 7, and has to be regarded as insignificant. On the other hand, one could equally conjecture that Coulomb GR and Coulomb JT methods show a direct correlation between ǫ N and the free energy E free , 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 ǫ N , 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. Fur- thermore, 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. 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 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.

Time evolution in the MHD-emergence stable and MHD-emergence unstable cases
Concerning the accuracy metrics for the magnetic field, Fig. 5 shows the ǫ N (B, ∇× 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 ǫ N = 0.54 at the beginning of the simulation, increasing up to ǫ N = 0.77 at the end. Comparing with the evolution in the MHD-un case, we notice that the metric ǫ 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 ǫ 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 V 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 De-Vore 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 V values, it is worth noticing that, for the simulation in exam, "larger H V 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 V 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 ǫ 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.

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 ≃ 3% at most in the examined resolution interval, ∆ = [0.03, 0.12]. Similarly, the H V 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 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 ǫ N (B, ∇× A) (cf. Eq. 38). The analysis of the complement of the vector error again sepa-rates the methods more markedly. There are no appreciable differences between the three DeVore methods, all around ǫ N =0.99. Errors are larger on average going from Coulomb SY (ǫ N =0.97), to Coulomb JT (ǫ N =0.90), to the Coulomb GR (ǫ 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 ∇ · B at lower resolutions, as 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 ∆, 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 A p , with DeVore SA slightly less sensitive than De-Vore KM and DeVore GV (see Fig. 7b). Differences become entirely negligible in the computation of A (cf Fig. 7c). The energy metrics ǫ E (cf. Eq 39) provide similar information, with more accentuated spreads (see Fig. 7d and Fig. 7e). 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 A p than A. The Coulomb GR method is even more accu- rate in this case than most of the DeVore methods in the computation of A p (in both ǫ N and ǫ E metrics), but yields contradictory metrics for A (poor ǫ N but good ǫ E ). This results into H V values that stand more clearly apart form the trend. The Coulomb GR improvement of ǫ 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 A p (respectively, A) with the Coulomb JT method yields ǫ N ≃ 0.88 (respectively, ǫ N ≃ 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 ǫ E , even though not exactly reproducing the accuracy of the methods as quantified by ǫ 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 ∇ · A = 0 and ∇ · 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, | f i (A)| and | f i (A p )| , as defined in Eq. (40), for the A p and A solution of the LL resolution test cases, respectively. All Coulomb methods satisfy the gauge better (i.e., have lower | f i | ) 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 A p and A, to a relatively low degree ( | f i | ≃ 10 −4 ). The other two Coulomb methods show larger differences in the fulfillment of the solenoidal property of A p and A. In particular, both the Coulomb GR method and the Coulomb JT method respect the gauge condition better for A p than for A. The difference is minimal for Coulomb GR -about 3 × 10 −5 -and about one order of magnitude for Coulomb JT. At the same time, Coulomb JT provides the most solenoidal A p and A of all Coulomb methods, at all resolutions, with values of | f i | 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 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.

Dependence on divergence
In Sect. 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 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 p effectively have the same distribution of the normal component on ∂V, a gauge transformation (A, A p ) −→ (A + ∇ψ, A p + ∇ψ p ) induce the following difference on the relative magnetic helicity: 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 p in such a way that the amplitude of the vector potentials turns out to be dominant over the contributions from the gauges ψ and ψ 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 div , the different methods provide helicity values that remain within a certain range of each other. As E div is infinitely small, the methods are expected to provide helicity values which are close to each other, while for large E 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 div for different values of δ in Eq. (41) is reported in Table 5 and shown in Fig. 2f, where E div is shown growing from 0.2% up to 14% for the considered range of δ values (the corresponding values of | f i | can be find in Table 7). The contributions E J,ns and E mix to E div in Eq. (44) and their variation with δ are shown in Fig. 8a. The energy E div turns out to be a nonlinear function of δ, with a minimum close, but not at, δ = 0. This apparent contradiction is due to the definition Eq. (44) that forbids cancellation between different contributions, and practically shifts the zero of E div to slightly higher values. For completeness, the nonsolenoidal error E 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 div , the methods are providing H V which are close to each other, while the dispersion of the values obtained are spreading as E div is increasing. All method follow a qualitatively similar exponential trend as a function of E div , both in H V (Fig. 8b) as well as in the accuracy metric ǫ N (B, ∇ × A) (Fig. 8c). The only exception to a smooth trend is the E div =2%-case of the Coulomb JT method, for which no particular reason was identified.
In the most solenoidal case, for E div = 0.2%, the relative dispersion of H V obtained is of 0.8%. Excluding the E div =2%-case of the Coulomb JT method, for each case with E div ≤ 4%, the relative spread in H V remains lower than 1%. At E 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 incline us to state that within that range of E div , the differences between the methods are not related to the nonsolenoidality but rather to the intrinsic numerical error within each method.
For E 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 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 div , the gauge dependence is directly impacting the estimation of H V .
The implementations of DeVore's method (DeVore KM, DeVore GV and De-Vore SA) yields very similar results, both in values and in trends. The trend already noticed in Fig. 7b is confirmed by the accuracy of the vector potential quantified by ǫ N in Fig. 8c. From the analysis in Sect. 6.2, we tend to attribute this difference to the accuracy in the solution of Eq. (8) being higher for DeVore SA with respect to DeVore KM and DeVore GV. Similarly, the Coulomb SY and Coulomb JT methods deliver analogous H V and ǫ N curves (see Fig. 8b and c). However, they strongly differ in the accuracy with which the Coulomb gauge condition is respected, as Fig. 8d shows. Confirming the result of the LL case of Sect. 4.1, the Coulomb JT method respect the Coulomb-gauge condition extremely well for A p ( | f i | ≃ 10 −9 ), and still excellently for A ( | f i | ≃ 10 −6 or lower). Interestingly, | f i (A p )| is practically independent of E div , confirming that the strategy adopted by the Coulomb JT method for solving Eq. (10) is very much able to handle flux unbalance resulting from solenoidal errors in V. However, even though the solenoidal property of the vector potentials is definitely better verified by the Coulomb JT method than by the Coulomb SY, the accuracy of the vector potential does not reflect this trend. The Coulomb SY vector potentials, indeed, have much higher values of | f i (A p )| ≃ | f i (A)| ≃ 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 | f i | seem still low enough to guarantee a relative high accuracy, as testified by ǫ N of Fig. 8c.
The test discussed here, of course, does not pretend to be general in assessing the influence of the nonsolenoidality on helicity estimations, and further tests are likely to be required. However, it represents one 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 div is below 1%, but become abruptly more important above that threshold. For a dataset with E div comprised between 1 and 8%, using on FV method or the other would lead up to 6% difference in the estimation of H V (excluding Coulomb GR method). For higher E 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 div (and hence its influence on H V ) is of the order of few percent (cf. Fig. 2 and Table 7).

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 V in FV methods range from the distribution of strong cur-rents 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 De-Vore 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 Fig. 7f and Fig. 8d). On the other hand, the accuracy of the Coulomb SY method in the con-struction 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 MHDun show, see Fig. 5), possibly, because they yield inconsistent boundary values for Eqs. (10, 13). The Coulomb GR method, finally, has mostly issues in computing an A that reproduces the field accurately enough, which result in the largest departures of 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 | f i (A)| 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 De-Vore 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.

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.

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 . 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 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 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 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. By contrast, the not-normalized values of the total helicity 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 H twist provided by the TN method is an estimation of H V,J rather than of 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 H V,J is reproduced.

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 .
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 , in which a single value of the force-free parameter α is used for the entire volume. When the NLFF mode switches on, the CB method has the possibility to replaces the total photo- spheric magnetic flux by the connected magnetic flux, namely the one included in the magnetic connectivity matrix, see Sect. 2.3. The CB code was tuned to use almost the entire flux of the magnetograms in all case presented here, except when differently explicitly stated.
When applied to the TD cases, where currents are localized, the CB in LFFmode tends to pick up mostly the potential field component. On the other hand, when applied to the LL cases, where a large-scale α 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 . 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 (Figure 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 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.
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 (Figure 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 present at the time of the eruption, very much the same as for the DeVore GV method. However, this is not distinguishable from previous, even more pronounced ones, and it would be challenging to identify the time of the eruption only as a decrease of H V in the CB time series. As a matter of fact, the examination of the time evolution of the magnetic field at z = 0 in the 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 V 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, differences 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 V -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, σ J ≡ ( V |J ⊥ | dV)/( V |J| dV). In particular, the value of σ J computed for increasing heights of the bottom boundary decreases from 0.7 to 0.25 in the first 20 pixels above the bottom boundary.
In order to test how important is the 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 σ 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 Fig. 12d 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 V 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 V 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 V 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  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 V 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 V .

Conclusions
In this work we review, benchmark, and compare the currently available methods for the computation of the relative magnetic helicity, H V , 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 (∇ · 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)  , and in Sect. 2.2 (De-Vore 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 3 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 Sect. 4. 1 and Sect. 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 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,  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 , 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 .
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 (forcefree) 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 nolinear 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 . 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 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 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, normalised 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.

B Accuracy of vector potentials
In addition to the metrics Eqs. (38, 39), we include here also the remaining 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.