{'fraction_non_alphanumeric': 0.03327725341469674, 'fraction_numerical': 0.015540791012004904, 'mean_word_length': 4.529801324503311, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computation of Relative Magnetic Helicity in Spherical Coordinates

Magnetic helicity is a quantity of great importance in solar studies because it is conserved in ideal magneto-hydrodynamics. While many methods to compute magnetic helicity in Cartesian finite volumes exist, in spherical coordinates, the natural coordinate system for solar applications, helicity is only treated approximately. We present here a method to properly compute relative magnetic helicity in spherical geometry. The volumes considered are finite, of shell or wedge shape, and the three-dimensional magnetic field is considered fully known throughout the studied domain. Testing of the method with well-known, semi-analytic, force-free magnetic-field models reveals that it has excellent accuracy. Further application to a set of nonlinear force-free reconstructions of the magnetic field of solar active regions, and comparison with an approximate method used in the past, indicates that the proposed methodology can be significantly more accurate, thus making our method a promising tool in helicity studies that employ the spherical geometry. Additionally, the range of applicability of the approximate method is determined and discussed.


Introduction
Magnetic helicity is a geometrical quantity that describes the twist, writhe, and linkage of magnetic-field lines. It is invariant in ideal magneto-hydrodynamics K. Moraitis kostas.moraitis@obspm.fr with B p = ∇ × A p the reference field, usually taken to be a potential field, and A p its vector potential. The important condition allowing gauge independence is that the studied field and the potential field must have the same normal component on the whole boundary of the system. The range of applications of relative magnetic helicity (hereafter simply helicity) in the Sun is quite broad, extending from the solar dynamo (e.g. Brandenburg and Subramanian, 2005) to the triggering of coronal mass ejections: according to Rust (1994), the existence of solar eruptions might be attributed to the need to shed the coronal magnetic helicity, constantly accumulating from a continuous injection through the photosphere.
In the solar context, different methods and approaches can be used to evaluate magnetic helicity. A review and comparison of several of these methods is presented by Valori et al. (2016). The focus of this article is on the finite volume methods (e.g. Thalmann, Inhester, and Wiegelmann, 2011;Valori, Démoulin, and Pariat, 2012;Yang et al., 2013), that is, methods that employ the definition of helicity as a volume integral and thus require the three-dimensional (3D) magnetic field in the entire volume as input.
In all of these methods, the computation of relative magnetic helicity is performed in Cartesian geometry. The natural coordinate system for the Sun however is the spherical one, and so it is necessary to be able to calculate helicity in spherical coordinates. Possible applications of such a calculation include MHD simulations in spherical geometry (e.g. Masson, Antiochos, and DeVore, 2013;Fan, 2016;Karpen et al., 2017), as well as nonlinear force-free (NLFF) field extrapolations (e.g. Gilchrist and Wheatland, 2014;Savcheva et al., 2016).
Although there are studies that compute magnetic helicity in spherical coordinates (e.g. Bobra, van Ballegooijen, and DeLuca, 2008;Fan, 2010;Savcheva et al., 2015;Karpen et al., 2017), these use various assumptions and/or simplifications. The respective helicity computations are then problem-dependent, and the methods used cannot be generalised to other datasets. In the MHD simulations of Fan (2010), for example, there is no magnetic flux penetrating the lateral boundaries of the volume, which results in the simplification of the helicity calculations. In another example, Karpen et al. (2017) drive a coronal hole jet with a rotational photospheric plasma motion that leaves approximately unchanged the initially prescribed potential magnetic field (and corresponding vector potential), thus simplifying the helicity calculations.
There are other cases, as in Bobra, van Ballegooijen, and DeLuca (2008), where a different helicity formula is used. In spherical coordinates this reads This relation stems from the original definition of relative helicity (Berger and Field, 1984) where space is divided in two regions; the volume of interest [V ], and a complementary volume [V ′ ] where B ′ = B ′ p , and thus A ′ − A ′ p = ∇χ for some scalar function χ (with primed quantities referring to V ′ ). The interface of the two volumes is denoted as S, while B r is the radial component of the field. Additionally, the boundary of V ∪V ′ is considered as a flux surface, or extending to infinity.
In the particular case examined by Bobra, van Ballegooijen, and DeLuca (2008), V is the coronal wedge-shaped volume of interest, and V ′ its sub-photospheric extension to the center of the Sun, and so neither of them is in general bounded by a flux surface. Following however Bobra, van Ballegooijen, and DeLuca (2008), V ′ is a bounded volume, and since V ∪ V ′ is bounded by a flux surface, it follows that V ′ should be bounded by a flux surface below the photosphere. This necessarily implies that the photospheric boundary S has to be flux-balanced, i.e. without net magnetic flux. These restrictions could be lifted by considering V ′ as the whole volume outside V , but this was not the choice made in the original derivation of Equation 2, as the limitation of the surface integral to the photosphere demonstrates. Indeed, the first term in Equation 2 represents the contribution to helicity from V , while the same quantity in V ′ reduces to the surface integral of the second term because of the assumption B ′ = B ′ p there. Had a different geometry of V ′ been adopted, an additional contribution to Equation 2 should be considered.
To see whether, and under which conditions, Equation 1 is equivalent to Equation 2, we apply the latter to the same volume as the former, V ∪ V ′ . The contribution to relative magnetic helicity from V ′ is however identicaly zero, since B ′ = B ′ p there. In the volume of interest, V , we can expand the terms in Equation 1 as The first term coincides with the respective term in Equation 2, while the second is a mixed term which can also be written as the surface integral after using standard vector idendities (Valori, Démoulin, and Pariat, 2012). Herê n stands for the outward-pointing unit vector on the boundary of the volume of interest ∂V .
To better compare the surface terms of the two helicity equations, we reverse the steps in the derivation of the second term of Equation 2, and get In the first equality the surface integral is extended to the whole boundary of the sub-photospheric volume, ∂V ′ , by the assumption that ∂V ′ − S is a flux surface,n ′ · B ′ = 0, withn ′ denoting the outward-pointing unit vector on ∂V ′ . The second equality follows directly from Gauss' theorem and an integration by Again, only the photospheric part survives in the surface integral since the remaining boundaries are flux surfaces. Moreover, as noted by Finn and Antonsen (1985), Equation 2 is gauge-invariant if additionally the tangential components of the two vector potentials are continuous across S, a condition that leads to Replacing in Equation 6 also thatn ′ = −n on S, we finally get The surface term in the Bobra, van Ballegooijen, and DeLuca (2008) formulation thus corresponds to the photospheric part of the H mix term. The comparison of the surface terms of the two helicity formulations, given by Equations 7 and 4, shows that they coincide if the conditionn·(A×A p ) = 0 holds on the coronal part of the boundary of the volume of interest, ∂V −S. The conditions for Equations 1 and 2 to be equivalent are thus thatn · (A × A p )| ∂V −S = 0 on the coronal part of ∂V , andn ′ · B ′ | ∂V ′ −S = 0 on the sub-photospheric part of ∂V ′ , also implying magnetic flux balance on S. The first condition is obviously gauge-dependent, unless V extents to infinity (where the vector potentials presumably vanish), or ∂V is a flux surface (so that the tangential components of A vanish there). If additionally the sub-photospheric part of ∂V ′ is treated as in Bobra, van Ballegooijen, and DeLuca (2008), i.e. it is assumed as generic and not bounded by a flux surface, with S not flux balanced, then the condition n · (A × A p ) = 0 must be valid also there, and then H mix = 0 holds throughout the boundary of V ∪ V ′ .
In deriving Equation 2 therefore, the gauge-dependent assumption H mix = 0 is made implicitly. This should be taken into account when applying this equation to finite volumes, since it is a gauge-dependent condition that is not valid in general.
In this article we extend the computationally most efficient and robust of the finite-volume methods in Cartesian coordinates (Valori, Démoulin, and Pariat, 2012;Moraitis et al., 2014) to the spherical geometry. The equations that we derive are applicable to both spherical-shell and spherical-wedge geometries.
However, we focus on the latter because it is of more common use. We also compare the results produced by this method with helicity values computed with an approximate method based on Equation 2, and we determine the applicability of the latter. In Section 2 we describe the implementation details of the method and in Section 3 we perform its validation against semi-analytic force-free magnetic field models. In Section 4 we apply it to a set of NLFF fields that are also described there, and finally in Section 5 we summarize and discuss the results of the article.

Relative Magnetic Helicity in Spherical Coordinates
Relative magnetic helicity is computed directly from its definition by Finn and Antonsen (1985), Equation 1. From a given 3D magnetic field [B] one has to estimate the other three vector fields appearing in Equation 1, namely the potential magnetic field [B p ], and the corresponding vector potentials [A and A p ] of the two magnetic fields. Note that helicity as given in Equation 1 is independent of the gauges used in the vector potentials as long as the normal components of the original and potential fields match along the whole boundary of the volume, i.e.n · B| ∂V =n · B p | ∂V .
The calculation of relative magnetic helicity can then be done in two steps, similarly to the Cartesian case (Valori, Démoulin, and Pariat, 2012;Moraitis et al., 2014), as detailed in the present section.

Calculation of the Potential Field
The potential field can always be expressed through a scalar potential Φ as B p = ∇Φ. The potential then satisfies Laplace's equation ∇ 2 Φ = 0 in V . In spherical coordinates this reads The requirement of Equation 8, ensuring gauge-invariance, leads to the following Neumann boundary conditions for the scalar potential We note that the solution of Laplace's equation under Neumann boundary conditions exists (up to an additive constant) only for flux-balanced fields, ∂V B · dS = 0, a condition that is never fully satisfied with numerical data. Laplace's equation is solved with the use of a FORTRAN routine contained in the MUDPACK library (Adams, 1989). The routine uses a multigrid iteration method to obtain the solution and is therefore efficient, quick, and computationally not demanding. The multigrid method requires a uniform computational grid that has dimensions of the special form m 2 n + 1, for n some positive integer, and m a small prime number, like 2, 3, or 5. In the case that either of these conditions is not fulfilled by the input magnetic field data, we linearly interpolate them to an appropriate grid, and we inevitably introduce numerical errors to the following helicity calculations.

Analytical Derivation
For the calculation of the vector potential of a given 3D magnetic field we follow the conceptual method initially developed by DeVore (2000), and subsequently adapted to finite volumes by Valori, Démoulin, and Pariat (2012). We start from the defining relation of the vector potential [B = ∇ × A] written in spherical coordinates We then employ the DeVore gauge (DeVore, 2000) in which the radial component of the vector potential is identically zero: A r = 0. This same gauge has also been used by Valori, Démoulin, and Pariat (2012) and respectively, and r 0 is an arbitrary radius inside the volume. These can be written more compactly as where vector that represents the vector potential on the surface r = r 0 . Substituting Equation 14 in the radial component of Equation 11, and using the divergence-freeness of the field, we get where the symbol "⊥" stands for the normal direction to the radial one. This relation can be thought of as a restriction on the possible choices for the remaining freedom in the gauge of A.

DeVore Simple Gauge
A first simple solution to Equation 15, which corresponds to a first family of gauge, is given by the relations and where θ 0 , φ 0 are arbitrary angles, and c ∈ [0, 1] a constant, typically c = 1/2. The vector potential in this gauge is given by Equations 14, 16, and 17, and its calculation requires the computation of four one-dimensional integrals. In our present implementation, the integrations are treated with the simple trapezoidal method, although more sophisticated methods can be envisioned. In the following we refer to this gauge family as the DeVore simple gauge (DVS).

DeVore Coulomb Gauge
An alternative family of gauges can be obtained when a, the vector potential on the surface r = r 0 , satisfies the Coulomb gauge, ∇ ⊥ ·a = 0. This is automatically fulfilled if a is expressed through a 2D scalar function [u], as given by By substituting Equation 18 in Equation 15 we find that the function u satisfies the 2D Poisson equation ∇ 2 ⊥ u = B r (r 0 , θ, φ). This equation can be solved easily with the additional assumption of the Dirichlet boundary condition, u = 0, along the boundary of the r = r 0 plane. The use of this boundary condition allows for a more generic B r -distribution, such as a non flux-balanced one, compared to the corresponding Neumann boundary condition, as follows from the uniqueness condition of the Poisson equation. In our particular implementation we solve the Poisson equation with the aid of a FORTRAN routine from the FISHPACK library (Swarztrauber and Sweet, 1979).
The vector potential in this gauge family is given by Equations 14, and 18, and its calculation requires the computation of two one-dimensional integrals, and of a 2D Poisson equation. While Amari et al. (2013) named this gauge the restricted DeVore gauge, we will denote it as the DeVore Coulomb gauge (DVC), because, when applied to the potential field, this satisfies the Coulomb gauge in the whole volume (see Section 2.2 in Valori, Démoulin, and Pariat, 2012).

Reference Boundary Choice
Apart from the two gauges for the integration vector (DVS and DVC) there is also the choice of the location of the reference surface in the vector potential calculation. We consider here only the two interesting cases, the bottom and the top surfaces of the volume, similarly to Equations 10 and 11 of Valori, Démoulin, and Pariat (2012), although other choices exist also (see Section 2.2.3 of Valori et al., 2016). These cases are denoted by the letters "b" and "t", respectively, and they follow the symbol of the gauge. A given vector potential can thus be in any of the four different gauges: DVSb, DVSt, DVCb, or DVCt.
The vector potential for the potential field can be obtained with the same method by just replacing B with B p in the above relations. It can thus also be in four different gauges and since in general the gauges of the two vector potentials are independent, there are 16 possible combinations in the calculation of helicity.

Test Datasets
The performance of the helicity calculation method is checked against the semianalytical, force-free field solutions of Low and Lou (1990, hereafter LL fields). The usual parameters of the LL fields are also assumed, n = m = 1, while the source is placed at a depth of ≈ 30 Mm below the AR (or l = 0.3 in LL notation), and is rotated by the angle φ = π/4 with respect to the radial direction. A plot of the LL field, denoted hereafter as B LL , with some representative field lines is shown in Figure 1. In order to check the effect of resolution we also use the same volume, but discretized by 257 × 257 × 257 grid points.
The first step in estimating the helicity of B LL is to calculate its scalar potential Φ LL as described in Section 2. The corresponding potential field is obtained from the relation B p,LL = ∇Φ LL , where the numerical derivatives used are of second order inside V , and on the boundary are given by Equation 8. The solenoidality of these magnetic fields is verified below.

Solenoidality Verification
The divergence-freeness of a magnetic field is quantified in three ways. The first is by the average fractional flux increase [ |f i | ] defined by Wheatland, Sturrock, and Roumeliotis (2000).
The second verification metric is the flux imbalance ratio [ǫ flux ]. This is defined by where Φ + and Φ − are the positively-defined fluxes entering and leaving the volume through all of its boundaries, respectively. A highly flux-balanced field, indicating also a solenoidal field, corresponds to the limit ǫ flux = 0, while the opposite happens for ǫ flux = 1. The third metric used relates to the magnetic energies. For each B LL and B p,LL at each resolution, we calculate its magnetic energy (in arbitrary units) as E = B 2 dV (cf. Table 1). By substracting the magnetic energy of the potential field from the magnetic energy of the LL field, we get the free energy: E c . Table 1 indicates that the free energy corresponds to ≈ 26 % of the total energy of the LL field, a value typical for LL fields, independently of the resolution. The solenoidality can finally be estimated thanks to the energy ratio E div /E, with E div a pseudo-energy given by This energy vanishes for perfectly solenoidal fields and quantifies the violation of Thompson's theorem by non-purely solenoidal fields (Valori et al., 2013). This quantity is more widely used to estimate non-solenoidality Valori et al., 2016;Pariat et al., 2017;Polito et al., 2017). Table 1 shows that the small values of E div /E is again an indication of the good level of solenoidality and another proof of the validity of the LL field as a test field, since nonsolenoidality was shown to be the strongest source of errors in helicity estimations (see Valori et al., 2016).

Vector-Potentials Verification
The second main step in the volume-helicity calculation is the computation of the vector potentials. It is important to check the quality of the reconstruction of the vector potential, i.e. that the curl of the computed vector potential corresponds indeed to its source magnetic field. For this, we use two methods. First, we directly compare the r, θ, and φ components of the original fields, B LL and B p,LL , with their respective vector-potential-reconstructed magnetic fields, ∇ × A LL and ∇ × A p,LL . To compare the components we compute their linear (Pearson) correlation coefficients, which are presented in Table 2.

Computation of Relative Magnetic Helicity in Spherical Coordinates
All metrics indicate the good agreement of the two fields by values close to unity. In Table 2 we present the values of these parameters for a few example reconstructions of the LL field, and separately, of the potential field. We performed different tests using: the two different gauges discussed in Section 2.2, DVS and DVC; the two different spatial resolutions, 129 3 and 257 3 ; and/or the two different reference surfaces: bottom and top. As an example, consider the first row in Table 2, where the values correspond to the reconstruction of the LL field of 129 3 grid points with a vector potential in the simple DeVore gauge, and with the top surface as the reference surface for the integration. Table 2 demonstrates that the correlation coefficients are all very close to 1. The value for the r-component presents the weakest correlation since this component involves the most numerical operations. The values for Schrijver's metrics also indicate that the reconstruction is excellent in all the cases. We further note that the reconstruction in spherical coordinates has comparable accuracy to the corresponding one for the Cartesian case (Table 8 in Valori et al., 2016).
From the values of Table 2 there can be drawn a few conclusions. First, the differences for the two different spatial resolutions are very small for both B LL , and B p,LL , with the latter exhibiting slightly better metrics in the lower-resolution field. Second, the differences between the reconstructions performed with the DVS or the DVC gauges are practically nonexistent, for both B LL and B p,LL . We note here that the vector potentials in the DVC gauge satisfy the Coulomb gauge on the surface r = r 0 to a large degree. This can be seen by taking the 129 3 -grid potential field as an example, where the absolute fractional flux increase of its vector potential in the DVCt gauge is |f i | = 1.4 × 10 −7 , a value comparable to the best-performing methods of Valori et al. (2016, Figure 7f).
The most important differences in the reconstruction metrics are found, however, between the top and bottom reference surfaces, mostly in the parameters E ′ n and E ′ m , which are the most sensitive ones. One can note that using the top boundary yields more precise results. This is a general characteristic of the helicity-computation methods based on the DeVore gauge and the properties of the data tested, as already noted by Moraitis et al. (2014); Pariat et al. (2015); Valori et al. (2016). In solar-like datasets the magnetic field at the top surface is weaker and smoother than in the bottom one, and, as a result, the integrations involved in the vector-potential computation that start from the top surface are also smoother.
Additionally, the values of the metrics for the potential fields [B p,LL ] are, in general, slightly inferior to the respective ones for B LL . This is to be expected since the potential field is derived from the LL one and thus carries any numerical errors in it.
As a final check we compute relative helicity as given by Equation 1 using all of the possible different gauge combinations for the vector potentials of B LL and its potential fields. The produced helicity values for the 129 3 -grid fields are all very similar, ranging in −145.2 ± 0.4 (in arbitrary units), indicating that helicity is indeed independent of the chosen gauge up to a factor 2 × 10 −3 .
As a conclusion of this Section we can say that the developed helicity computation method in spherical coordinates is performing extremely well. This is deduced from the high solenoidality of the constructed potential magnetic fields, and also, from the level of reproduction of the magnetic fields from the respective vector potentials. In addition, the method is computationally very efficient, as it requires only a few minutes to compute helicity for the 257 3 -grid field on a commercial laptop.

Application to Solar Active Region 3D Extrapolated Magnetic Fields
Having established that our method performs as expected on a semi-analytical magnetic field, we now apply it to typical magnetic-field data from a reconstructed solar active region. This will show the performance of our method to a more solar-like case, and also, it will give the opportunity to compare our results with more approximate helicity calculations methods.

NLFF Field Model
We use a set of data-constrained non-linear force-free field models of the NOAA AR 11060 (SOL2010-04-08T02:35:00L110C176). These have been constructed in order to topologically study the specific active region (Savcheva et al., 2016). All computations were performed in spherical coordinates. The fields were produced by the flux-rope-insertion method (van Ballegooijen, 2004), which involves various intermediate steps that we briefly summarise here. First, a global potential magnetic field is computed from a low-resolution synoptic magnetogram. Then, a modified potential-field extrapolation is performed starting from a high-resolution magnetogram centred on the active region, and with lateral boundary conditions given by the global field. The relevant calculations are done directly with the vector potential [A p,MF ], to ensure solenoidality of the respective magnetic field. This potential field [B p,MF = ∇ × A p,MF ], is then modified with two sources in the photosphere where the inserted flux rope is to be anchored.
The insertion of the flux rope is done by modifying the initial vector potential with a combination of axial and poloidal flux. The model is then relaxed towards a force-free state using the magnetofrictional (MF) method. In the MF code, it is the vector potential again that is relaxed, A MF , and not the magnetic field B MF = ∇ × A MF . During the MF relaxation,n · B MF is held constant at the photospheric boundary; the other five boundaries are very far away from the AR core so that the field there is very well approximated by the potential one. The potential magnetic field and its respective vector potential remain constant and equal to their value at the initial instant during the MF process, as this is not a temporal evolution. The full details of the method used to produce the NLFF fields are described by Savcheva et al. (2016) and references therein.
A magnetic-field model is considered as stable if, after the relaxation, it has reached a force-free equilibrium where no residual Lorentz force is present, otherwise it is considered to be marginally stable or unstable. An unstable model is obtained from a stable model by adding slightly more axial flux to the stable inserted-flux-rope parameters. The parameters used in the present study are the same as those used by Savcheva et al. (2016).
This marginally unstable NLFF model leads also to an "eruption" resembling a coronal mass ejection during the MF evolution of the field, as was demonstrated by Savcheva et al. (2016). The gradual expansion and elevation of the flux rope can also be inferred from the evolution of the field-line connectivity and of the flux-rope height that are shown in Figure 2. A data-constrained MHD simulation starting from a similar initial condition and proving that an unstable model obtained with the flux-rope insertion method is unstable as well in the MHD sense was performed by Kliem et al. (2013) and demonstrated a good correspondence with the observations. We use 38 snapshots of the magnetic-field evolution, three at very early stages at 0, 100, and 1000 iterations, and the remaining from 5000 to 175, 000 iterations with a step of 5000 iterations. The magnetic field covers an area of ≈40 • ×40 • in the θ-φ plane, and it extends 805 Mm in the radial direction. The magnetic-field datacubes have the dimensions of 385×385×385 grid points. The grid in the φdirection is uniform with a stepsize of 0.11 • , while along the r-, and θ-directions it is nonuniform. In the θ-direction the stepsize increases from 0.08 • to 0.11 • as we go to the Equator, while along r the stepsize increases from 1.4 Mm in the photosphere, to 3 Mm in the uppermost level. The grid that we use here is thus structured.

Results
In this Section we compute the instantaneous value of helicity during the relaxation of the magnetic field with the magnetofrictional method, using the method of this article, and also, the approximate method based on Equation 2. The method of this article will be denoted as the exact method since it makes no assumptions, it treats appropriately all boundaries, and it performs very well with the LL field, as was shown in the previous section.

Exact Method
To compute magnetic helicity with the method of this article we start from the magnetic field [B MF ] that is defined on a nonuniform grid, and we interpolate it to a uniform grid on the same volume. This is not a necessary step in general, but the current implementation of the method is more accurate with a uniform grid. A simple trilinear interpolation is sufficient to obtain similar levels of accuracy as in the LL case, as it will shortly be seen.
From the interpolated magnetic field [B DV ] we compute the vector fields, B p,DV , A DV , and A p,DV , using the method described in Section 2. Note that the potential field is computed at each snapshot in order to take into account possible changes ofn · B DV caused by the interpolation from B MF . The helicity computation that uses the method of this article, is thus based on the relation The evolution of magnetic helicity as computed with the exact method is shown in Figure 3. The helicity obtained is an increasing function of iteration number, after the first couple of snapshots of the flux-rope insertion, and it seems to saturate at a value a little above ≈ 1.76 × 10 42 Mx 2 . The increasing helicity evolution pattern is in contrast to the respective free energy one, which is decreasing as shown in Figure 8 of Savcheva et al. (2016), and also, to the general trend between free energy and helicity . Nevertheless, it is in line with the removal of negative helicity as the flux rope expands and subsequently erupts during the MF evolution. The level of gauge invariance is quite high, as can be deduced from the values of helicity computed with three different combinations of the DVSt and DVCt gauges for the vector potentials of the original and of the potential magnetic fields. The three different helicity curves of Figure 3 practically coincide, exhibiting differences 0.2 %.
The computed potential magnetic field is sufficiently divergence-free, as the average (over all snapshots) values for the metrics |f i | = 3.0 × 10 −4 and ǫ flux = 3.5 × 10 −4 indicate. The average values of the reconstruction metrics for the vector potential A DV (as in  a large degree, as Figure 3 shows, this leads to the conclusion that this can be considered as the actual AR helicity, and it further justifies the characterization of the method as exact.

Approximate Method
The second method that we use to compute helicity follows the reasoning of Bobra, van Ballegooijen, and DeLuca (2008), that is, it employs the definition of Equation 2. We furthermore use two different sets of vector fields in the computations. The first is the direct outcome of the magnetofrictional code. The relation that is used in this case is where all MF-related quantites are explicitly stated. The scalar quantity χ MF at a point P of the photospheric plane S is given by the line integral χ MF (P) = where O is a reference point on the same plane. Notice also that the vector fields involved in Equation 27 are the original, uninterpolated ones, and that the potential field and its respective vector potential do not change during the relaxation. The MF code evolves A MF but on the (coronal) boundary the condition A MF × A p,MF = 0 is maintained valid at all times so that Equation 27 is as accurate as possible. The second set of vector fields that we use in this approximate method is the ones from the exact method. We use again the same relation, namely with χ DV (P) = P O dl · (A DV − A p,DV ). The evolution of helicity obtained with these two approximate methods is plotted in Figure 4, along with the result of the exact method in the DVSt-DVSt gauge combination.
We note that the H R method exhibits a very similar pattern with the exact method. The helicity values are however on average 15 % higher than those of the exact method. This difference is due to the limitations of the formulation of Equation 2 with respect to the finite-volume geometry used, as was pointed out in the Introduction. We see therefore that our helicity computation method improves, by the designated percent, with respect to the approximate method of Equation 2.
For the H ′ R method, we note that the general pattern of helicity evolution is similar to the other two methods. The obtained helicity values however, are (on average) higher by 55 % than the exact ones, and by 35 % than those of the H R method. This large difference is expected since the vector potentials of Equation 28 are in the DeVore gauge and they do not assume the condition H mix = 0 that Equation 2 implies, since A DV × A p,DV = 0 on the boundary. According also to the discussion in the Introduction, the relation of Equation 2 depends on the gauges chosen for the vector potentials, unless the bounding surface is a flux surface. The use of Equation 2 with vector potentials that do not respect its validity conditions can thus lead to differences in helicity of the order of 30 − 40 %, at least in the case examined here.
We conclude that helicity calculations based on Equation 2 are problematic when using finite volumes (in whatever coordinate system). The conditions of validity of Equation 2 are rarely true in practice and so its use with finitevolume data should be avoided. Situations where it is safer to use this relation include, for example, the infinite-plane geometry originally used in Berger and Field (1984).

Discussion
A method for computing relative magnetic helicity in spherical geometry was presented in this article. The necessity for such a method stems from the fact that magnetic helicity is an important quantity with many applications in the solar context, and that the natural coordinate system for the Sun is the spherical one. The developed method treats the generic problem of computing helicity in a spherical box, and it is thus superior to previous methods that either employed simplified conditions on the boundary for the magnetic fields and/or the vector potentials, or used approximations in the computations.
Testing of the developed method with semi-analytic NLFF field models showed that it is working very well, with comparable accuracy to corresponding calculations in the Cartesian case. More specifically, the potential magnetic field produced by the method is highly divergence-free, and also, the computed vector potentials reproduce the respective magnetic fields to a high degree.
Additionally, from the application on a data-driven NLFF field model of a solar AR we already see an improvement on the values of helicity compared to approximate methods used in the past. The specific approximate method that we examined here suffers from an important limitation. It assumes a certain choice for the gauge of the vector potentials, which is hard to enforce in finite volumes. As a result, the obtained helicity values with the different gauge choices can be quite different from each other, and also, from the method developed in this article.
In all cases where the magnetic field is in the spherical geometry and it is known in the whole volume, such as in MHD simulations or NLFF field reconstructions, the method presented here is expected to give the most accurate estimation of the structure's magnetic helicity, and additionally with minimal computational effort and resources.