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

Space Science Reviews, Oct 2016

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11214-016-0299-3.pdf

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

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


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11214-016-0299-3.pdf

Gherardo Valori, Etienne Pariat, Sergey Anfinogentov, Feng Chen, Manolis K. Georgoulis, Yang Guo, Yang Liu, Kostas Moraitis, Julia K. Thalmann, Shangbin Yang. Magnetic Helicity Estimations in Models and Observations of the Solar Magnetic Field. Part I: Finite Volume Methods, Space Science Reviews, 2016, 147-200, DOI: 10.1007/s11214-016-0299-3