Issue 
A&A
Volume 553, May 2013



Article Number  A38  
Number of page(s)  14  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201220982  
Published online  29 April 2013 
Accuracy of magnetic energy computations^{⋆}
^{1}
LESIA, Observatoire de Paris, CNRS, Université Pierre et Marie Curie,
Université Denis Diderot,
92195
Meudon,
France
email:
gherardo.valori@obspm.fr
^{2}
University of Potsdam, Institute of Physics and Astronomy,
14476
Potsdam,
Germany
^{3}
Space Weather Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA
Received:
21
December
2012
Accepted:
26
February
2013
Context. For magnetically driven events, the magnetic energy of the system is the prime energy reservoir that fuels the dynamical evolution. In the solar context, the free energy (i.e., the energy in excess of the potential field energy) is one of the main indicators used in space weather forecasts to predict the eruptivity of active regions. A trustworthy estimation of the magnetic energy is therefore needed in threedimensional (3D) models of the solar atmosphere, e.g., in coronal fields reconstructions or numerical simulations.
Aims. The expression of the energy of a system as the sum of its potential energy and its free energy (Thomson’s theorem) is strictly valid when the magnetic field is exactly solenoidal. For numerical realizations on a discrete grid, this property may be only approximately fulfilled. We show that the imperfect solenoidality induces terms in the energy that can lead to misinterpreting the amount of free energy present in a magnetic configuration.
Methods. We consider a decomposition of the energy in solenoidal and nonsolenoidal parts which allows the unambiguous estimation of the nonsolenoidal contribution to the energy. We apply this decomposition to six typical cases broadly used in solar physics. We quantify to what extent the Thomson theorem is not satisfied when approximately solenoidal fields are used.
Results. The quantified errors on energy vary from negligible to significant errors, depending on the extent of the nonsolenoidal component of the field. We identify the main source of errors and analyze the implications of adding a variable amount of divergence to various solenoidal fields. Finally, we present pathological unphysical situations where the estimated free energy would appear to be negative, as found in some previous works, and we identify the source of this error to be the presence of a finite divergence.
Conclusions. We provide a method of quantifying the effect of a finite divergence in numerical fields, together with detailed diagnostics of its sources. We also compare the efficiency of two divergencecleaning techniques. These results are applicable to a broad range of numerical realizations of magnetic fields.
Key words: magnetic fields / methods: numerical / Sun: surface magnetism / Sun: corona
Appendices are available in electronic form at http://www.aanda.org
© ESO, 2013
Numerical tests of Thomson’s theorem.
1. Introduction
Many astrophysical phenomena, such as stellar and solar jets, flares, and coronal mass ejections, are driven magnetically (e.g., Tajima & Shibata 2002; Schrijver & Zwaan 2008, and references therein). Magnetically dominated plasmas are systems where the longrange, magnetic interaction dominates other forces, e.g., plasma pressure and gravitational forces. A typical example is the lowcorona (e.g., Priest 2003; Golub & Pasachoff 2009). There, the amount of energy associated with the magnetic field is much larger than other energy sources, and the dynamics of the coronal configuration is determined by the evolution of its magnetic field (e.g., Forbes 2000). This includes solar flares, where large currents develop in relatively small volumes (e.g., Shibata & Magara 2011; Aulanier et al. 2012), and coronal mass ejections (CMEs), which are powerful expulsions of coronal material that change the local configuration of the magnetic field drastically (e.g., Forbes 2000; Amari et al. 2003; Fan 2010). In the coronal plasma, the magnetic energy is therefore the prime energy reservoir that fuels the dynamical evolution of these events.
However, not all the magnetic energy is available for conversion into other forms of energy. Without changing the field significantly at the boundaries of the considered volume, the energy that can be converted into kinetic and thermal energies is given by the free energy, i.e., by the difference between the total magnetic energy and the energy of the corresponding currentfree (potential) field. This very general result is known as Thomson’s theorem, and it is based on the decomposition of the field into the sum of a currentcarrying and a potential part. It does not depend on the presence of other forces, and is valid at any instant in time.
The separation in the potential and free energies of Thomson’s theorem is especially relevant for systems like the lowcoronal field, that have different evolution time scales, as follows. The time scale of the coronal potential field is determined by the underlying photosphere, which is an inertiadominated plasma, unlike the corona. This implies that the magnetic field at the photosphere has an evolution time scale that is much longer than the coronal one and that it is relatively insensitive to coronal changes. Since the magnetic field at the photosphere largely determines the coronal field’s currentfree component, the latter also evolves on the long photospheric time scale. As a consequence of Thomson’s theorem, relatively fast events, such as flares and CMEs, can only be powered by converting part of the magnetic free energy (e.g. Aulanier et al. 2010; Karpen et al. 2012).
In other words, the magnetic free energy is a sufficient condition for triggering active events, and it is considered in the forecast of eruptions in the space weather context (see, e.g., Forbes et al. 2006; Chen 2011). Therefore, in this and similar applications, an accurate estimation of the free energy is paramount for understanding the observed magnetic field dynamics and the maximum energy that can be released in a flare or in a CME (Emslie et al. 2012; Aulanier et al. 2013, and references therein).
On the other hand, the free energy only provides an upper limit to the energy available for coronal dynamics. For instance, in the case of a flare/eruption, the postevent magnetic field configuration does not need to be potential (see, e.g., Berger 1985; Taylor 1986; Low 2001). Indeed, flare (reconnected) loops are frequently observed to be sheared after a flare/eruption (see e.g., Asai et al. 2003; Lin et al. 2010; Savage et al. 2012, and references therein), a feature that is also reproduced in numerical simulations (Aulanier et al. 2012). This is an indication that postevent configurations have finite free energy, and the actual energy removed by the event is given by the difference between the free energy of the pre and postevent configurations. An assessment of the true energy budget related to a flare/eruptive event requires an accurate and reliable estimation of the magnetic energy.
Another motivation of this study is to address the occurrence of unphysical magnetic configurations. This is the case in some nonlinear forcefree field (NLFFF) extrapolations when nonpreprocessed, observed vector magnetograms are used as boundary conditions. The most obvious evidence of the nonphysical nature of some solutions is when the energy of the extrapolated field is lower than the potential field energy. This happens, for instance, in some of the solutions given in Table 3 of Metcalf et al. (2008) and Table 1 of Schrijver et al. (2008) for three of the considered extrapolation methods, including one used in the present manuscript (Valori et al. 2010). More generally, for all methods, the estimated coronal energy depends on the manipulations performed on the observed data prior to their use in the actual extrapolation. (This step is called preprocessing, Wiegelmann et al.2006; Fuhrmann et al.2007.) A significant part of the energy difference can eventually result from the details of the undergone preprocessing.
As a result, the understanding of basic physical processes in the solar atmosphere requires an accurate estimations of the magnetic free energy. On the other hand, coronal models like NLFFF extrapolations, have shown that such accurate estimations are not easily obtained. In such cases, Thomson’s theorem can be exploited to address the accuracy of (free) energy estimations. The fundamental assumption in Thomson’s theorem is that the magnetic field is solenoidal. Such a property is only approximately fulfilled in numerical simulations and, more generally, in magnetic fields that are discretized on a mesh. A quantitative estimation of the effects caused by nonvanishing field divergence is complicated by its nonlocal nature.
The main aim of this article is to quantify the effect of the presence of a nonsolenoidal component on the energy of a discretized magnetic field. This is studied using six different test magnetic fields that are a sample of the typical and characteristic examples used in the context of coronal solar physics. In the first part of the article, the energy of each test field is decomposed and interpreted using an extension of Thomson’s theorem that can be applied to nonsolenoidal fields. In the second part we study how the energy changes, starting from a solenoidal version of each test field and adding a parametric divergent component. The method and results of this study are of interest when working with any discretization of magnetic fields, e.g., for 3D coronal magnetic field extrapolations, as well as for magnetohydrodynamic (MHD) simulations.
In Sect. 2 the Thomson theorem for the energy of a magnetic field is summarized. The extension to nonsolenoidal, discretized fields is presented in Sect. 3. Section 4 introduces the six discretized fields together with their corresponding solenoidal versions that are used as test cases for applying our analysis, whose results are given in Sect. 5. Possible source of errors in our analysis are sort out in Sect. 6. Then, in Sect. 7 we present the parametric study of the energy dependence on the amount of divergence added to solenoidal magnetic fields. An analysis specific to numerical fields obtained by NLFFF extrapolations of observed vector magnetograms is presented in Sect. 8, and conclusions are finally given in Sect. 9.
2. Magnetic energy of solenoidal fields
We first consider the decomposition of the magnetic energy for perfectly solenoidal fields. By decomposing the field B as the sum of a potential, B_{p} = ∇φ, and a current carrying contribution, B_{J}, the total magnetic energy E in CGSGaussian units in a volume is given by (1)where represents the boundary of , , and is the external normal to the bounding surface.
Two conditions are classically considered:

[a]
, i.e., the potential field B_{p} is computed from the same distribution of normal field of B on the boundary of . This condition implies that and the surface integral vanishes in Eq. (1);

[b]
∇·B_{J} = 0, in which case also the rightmost volume integral in Eq. (1) vanishes.
If these two conditions hold, then (2)and the energy of a magnetic field is bounded from below by the energy of the corresponding potential field that has the same distribution of the normal component on the boundary of the considered volume. When applied to discretized fields, the above result holds under the implicit assumption that fields are numerically well resolved, yielding, in particular, continuous derivatives.
The mathematical equivalent of Eq. (2) is known as Thomson’s (or Dirichlet’s) theorem, see e.g., Lawrence (1998).
To satisfy the above requirement [a], the scalar potential φ(x,y,z) is computed as the solution of the Laplace equation (3)In practical applications, Eq. (3) can be solved numerically using standard methods. In the applications presented in this paper, the Poisson solver included in the Intel Mathematical Kernel Library was used.
3. Magnetic energy of nonsolenoidal fields
In this section we provide expressions for evaluating errors in the energy that stem from an imperfect fulfillment of the solenoidal property, as is the case for discretized magnetic fields. In deriving Eq. (1) the divergence theorem, i.e., (4)is used, which may not be fulfilled by the techniques employed in constructing the numerical representations of magnetic fields or in their analysis. Moreover, if the numerically computed potential field B_{p} and currentcarrying field B_{J} have a finite divergence, additional contributions can appear in the corresponding energy terms, E_{p} and E_{J}.
We, therefore, seek a formulation of Eq. (1) for applications to numerical, nonsolenoidal fields that includes all possible sources of errors separately, that satisfies the requirement [a], and that includes only volume integrals (thus avoiding using the divergence theorem). To obtain that, we first introduce the method of computing the potential and currentcarrying parts.
3.1. Helmholtz decomposition of the potential part of the field
The accuracy in the numerical solution of Eq. (3) is limited, which may result in a finite divergence of the potential field. To quantify its effect, we can write (5)which separates in B_{p} the solenoidal part, B_{p,s} ≡ B_{p} − B_{p,ns} from the nonsolenoidal one, B_{p,ns} ≡ ∇ζ. This is equivalent to adopting the Helmholtz decomposition for the vector B_{p}, together with the choice that all the nonsolenoidal component of B_{p} is contained in ∇ζ. Finally, the boundary condition for ζ(x,y,z) in Eq. (5) is chosen such that B_{p,s} satisfies the same boundary condition as B_{p}; i.e., they both fulfill requirement [a].
In practical applications, we first solve Eq. (3) numerically to determine φ, then we compute B_{p} = ∇φ, and finally we overwrite the values of the normal components of B_{p} on each boundary according to Eq. (3). Since the latter operation enforces the requirement [a], then any residual inaccuracy in the solution of Eq. (3), close to the boundary, implies a jump in the field, i.e., a finite divergence that adds to the divergence of the potential field discussed above. Second, we solve Eq. (5) to compute the residual nonsolenoidal component in B_{p}.
3.2. Helmholtz decomposition of the currentcarrying part of the field
Using the Helmholtz decomposition on B_{J} we define a solenoidal component, B_{J,s}, and a nonsolenoidal one, B_{J,ns}, such that (6)the nonsolenoidal part of B_{J} being: B_{J,ns} ≡ ∇ψ. The boundary condition for ψ in Eq. (6) is chosen to have the same boundary condition for B_{J,s} and B_{J}, i.e., to fulfill the requirement [a]. As for the potential field, the required values of B_{J,s} at the boundaries (i.e., zero in this case) are overwritten onto the solution of Eq. (6) which is obtained numerically, so that any error in matching these values by ψ(x,y,z) reduces to a finite jump close to the boundaries.
Finally, we notice that this method is often used to remove the divergence of a vector field (Brackbill & Barnes 1980,sometimes referred to as “projection method”), and it has the property of conserving the current, i.e., ∇ × B_{J} = ∇ × B_{J,s}.
3.3. Gaugeinvariant decomposition of the magnetic energy
We now summarize the procedure for the decomposition of the magnetic field. For a given numerical magnetic field B, we solve Eq. (3) numerically and compute the corresponding potential component B_{p} and currentcarrying component B_{J} = B−B_{p}. Next, we compute the solenoidal component B_{p,s} = B_{p}−∇ζ and the nonsolenoidal component B_{p,ns} = ∇ζ of the potential field by solving Eq. (5) numerically. Similarly, the numerical solution of Eq. (6) provides the solenoidal component B_{J,s} = B_{J}−∇ψ and the nonsolenoidal component B_{J,ns} = ∇ψ of the currentcarrying part of B. The values of the different components at the boundary are such that the condition [a] is satisfied (Sect. 2). Finally, by substituting the field decomposition in and grouping it again as in Eq. (1), we obtain (7)with
(8)All terms in Eq. (7) are positively defined, except for E_{mix}. For a perfectly solenoidal field, it is E_{p,s} = E_{p}, E_{J,s} = E_{J}, E_{p,ns} = E_{J,ns} = E_{mix} = 0, and Eq. (7) reduces to Eq. (2).
Finally, Eq. (7) is normalized such that (9)where the tilde indicates that the corresponding definition in Eq. (8) is divided by E.
Using the divergence theorem, Eq. (4), and the condition [a], several terms in the above expressions could be simplified. However, since practical test fields may be obtained with methods that do not insure that the divergence theorem holds numerically, we have kept all the terms in Eq. (8). Indeed, the simplification obtained by using the divergence theorem results in mixing other numerical issues with the issue of the finite divergence, producing cumbersome results, up to the point where Eq. (7) is not satisfied numerically. Moreover, the direct appearance in the integrals of the scalar potentials, rather then their gradients, introduces an undesired gaugedependence.
3.4. Sources of the violation of the Thomson theorem
We summarize which are the source of errors that we consider in Eq. (7). First, the energy is affected by the finite divergence of the currentcarrying part of the magnetic field, which enters the E_{J,ns} and E_{mix} terms. Additionally, the potential field may have a finite divergence, owing to the limited numerical accuracy of the solution of Eq. (3), both in the volume and close to its boundary. These effects are contained in the E_{p,ns} and E_{mix} terms.
As long as these are the only source of errors, then the sum of the terms on the righthand side of Eq. (7) must be equal to the total energy E computed using B directly, and Eq. (7) must hold numerically even for nondivergencefree fields. Equivalently but using normalized quantities, the sum on the righthand side of Eq. (9) must be equal to one. We show in Sect. 5 that the total energy is indeed retrieved by the decomposition we adopted, allowing us to identify the source and extent of the eventual violation of Thomson’s theorem, Eq. (2).
3.5. Accuracy of the decomposition of the energy equation
A further step is the assessment of the accuracy of the decomposition, Eq. (7). First, we address how effective the decomposition in the solenoidal and nonsolenoidal parts is in concrete numerical applications in Sect. 6.1.
Second, the continuity condition, implicit in the derivation of Eq. (7), implies that numerical derivatives can be computed precisely enough in the employed discretization. This may not be the case in some numerical applications, e.g., when observed values are used as boundary conditions for computing magnetic fields. The continuity of the fields in relation to small scales is discussed in Sect. 8.
Finally, our decomposition employs the numerical solution of Laplace and Poisson equations. We briefly recall the conditions for uniqueness of the general Poisson equation (10)where f(x,y,z) is a source term in , and g is the boundary value on . The use of Neumann boundary conditions implies that the solution u(x,y,z) is only unique up to an additive constant. For Eqs. (3), (5) and (6), the freedom in the additive constant is equivalent to a gauge freedom for the scalar potentials φ, ζ, and ψ, respectively. This gauge dependence is, however, irrelevant for Eq. (7), since the energy decomposition is intentionally derived in a way such that the scalar potentials only appear in conjunction with the gradient operator.
Integrating Eq. (10) in and using the divergence theorem, Eq. (4), we find that source and boundary values must satisfy (11)which is a necessary condition for the uniqueness of the solution u. This implies that, for Eq. (3) where f = 0 and , the flux of B through must vanish. For Eq. (5), where f = ∇·B_{p} and g = 0, it implies that the volume integral of ∇·B_{p} must vanish. Similarly, for Eq. (6), where f = ∇·B_{J} and g = 0, the volume integral of ∇·B_{J} must vanish. When such conditions cannot be insured, the uniqueness of the solution is not guaranteed. The effect of the violation of Eq. (11) is studied in Sect. 6.2.
4. Test fields
Fig. 1
Selected field lines of the six test cases: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. The vertical component of the magnetic field at the bottom boundary is shown on a gray scale, with the positive (respectively, negative) polarity in white (respectively, black). The different line colors outline different types of connectivities. 
To explore the effects of a finite divergence in a representative sample of practical situations, we consider six test fields B_{test} obtained from analytical models, numerical simulations, and NLFFF extrapolations. Their magnetic configuration is outlined in the fieldline plots in Fig. 1. Furthermore, we consider six additional test cases B_{test,s}, which are obtained from each of the B_{test} by removing the nonsolenoidal part of the field.
4.1. Discretized analytical test fields
The first test field that we consider is the potential field B_{test} = B_{DD} generated by a pair of vertical magnetic dipoles, located at (0, ± y_{DD},z_{DD}), see, e.g., Eq. (7) in Török & Kliem (2003) for the analytical expression of the field. We set y_{DD} = 2 and z_{DD} = −1.5, and the field is normalized such that the zcomponent has a maximum value equal to unity at the bottom boundary (z = 0). The only currents and finite divergence errors present in B_{DD} are generated by truncation errors in its discretization.
The second employed test field, B_{test} = B_{TD}, is the model of the magnetic field of an active region derived in Titov & Démoulin (1999), given by a section of a current ring surrounded by a stabilizing potential field. The employed configuration is the same as in Valori et al. (2012), to which we refer the reader for further details. In this case, the test field has an explicit currentcarrying component sustained by a flux rope. The analytical formulae defining the test field are approximate, which together with the rather coarse resolution employed here, yield relatively large finitedivergence errors.
For both test fields B_{DD} and B_{TD} the discretized volume is , with uniform resolution Δ = 0.12 in all directions.
4.2. Numerical tests fields
The next test field that we consider, B_{test} = B_{MHD}, is a snapshot of a magnetohydrodynamic numerical simulation of magnetic reconnection in a nullpoint topology (Masson et al. 2012). To use our presentstage diagnostic, we interpolated the original snapshot onto a uniform and homogeneous grid, whereas the original simulation was performed using a nonuniform one. Because the divergence values are slightly increased by the interpolation, they are not representative of the quality of the simulations presented in Masson et al. (2012). However, they still serve our purpose of providing a typical situation arising from the numerical evolution of magnetohydrodynamic equations. The considered volume is with uniform resolution Δ = 0.05 in all directions, and the field is normalized such that the vertical component is unity at its maximum.
Next, we consider three NLFFF extrapolations of Hinode/SOT vector magnetograms, obtained with the magnetofrictional method in Valori et al. (2010). The original resolution of the vector magnetograms is 0.3′′, and they can be preprocessed (Fuhrmann et al. 2011) to improve their compatibility with the forcefree assumption on which the extrapolation code is based.
Our fourth test field, B_{test} = B_{Ex1}, is the nonlinear extrapolation of a vector magnetogram of AR 11158, measured on 14 February 2011. The vector magnetogram was binned to the resolution Δ = 1.1′′ prior to extrapolation, and no preprocessing was applied in this case. The analyzed coronal model volume in arcsec is . The Hinode/SOT field of view of the measurements employed for this extrapolation cuts through the external sunspots of a quadrupolar field distribution, resulting in high field values at the lateral edges of the magnetogram. Even computing the potential field is problematic in this case, therefore we limited the considered volume to the bipolar core of the extrapolated field.
The fifth test field, B_{test} = B_{Ex2PP}, is the extrapolated coronal field model above AR 11024 on 4 July 2009. In this case, the full resolution of Hinode/SOT is used, and the vector magnetogram is preprocessed before extrapolation. The extrapolation covers a volume of arcsec, with uniform resolution Δ = 0.3′′. This model of the coronal field of AR 11024 is discussed in detail in Valori et al. (2011), where more details about extrapolation of vector magnetograms can be found.
Finally, the sixth test field, B_{test} = B_{Ex2}, is the same case as B_{Ex2PP} except that the vector magnetogram is not preprocessed prior to extrapolation. More details on the numerical implementation are given in Appendix A.
4.3. Cleaned test fields
Since a small divergence of B is one major condition for the Thomson theorem, Eq. (2), for each test field B_{test} we consider a corresponding solenoidal version of it, B_{test,s}, which is computed from B_{test} employing the divergence cleaner described in Appendix B. In Cartesian coordinates, such a solenoidal field has the same x and ycomponents as B_{test}, whereas the zcomponent is changed everywhere in the volume, except for the top boundary. Therefore, B_{test} and B_{test,s} have the same distribution of normal field on all boundaries except for the bottom one, where B_{test,s} differs from B_{test} by an amount that is related to the combined effect of ∇·B_{test} in the whole volume. Since the divergence cleaner changes the value of the normal field component on one boundary, the potential fields computed from the boundary values of B_{test} and of the corresponding solenoidal B_{test,s} are not the same. Additionally, the divergence cleaner alters the current of the field, as prescribed by Eq. (B.5), of an amount that is proportional to the divergence of B_{test}. Therefore, the field that is obtained by applying the cleaner may have drastically different properties than the original field. Finally, let us notice that different solenoidal fields can be derived from B_{test} using different methods. The divergencecleaned versions of the test fields B_{test,s} are used here as illustrative examples.
5. Numerical tests of Thomson’s theorem
In this section we apply Eq. (9) to the test cases described in Sect. 4. Table 1 summarizes the values of the divergence metric ⟨ f_{i} ⟩ defined in Appendix C and the contribution of each term to Eq. (9), for all test fields. The divergence metric spans values from 10^{21} to 10^{3}. In all cases, the rightmost column, corresponding to the sum of the righthand side of Eq. (9), is equal to unity, despite the large difference in the divergence values. Therefore, we conclude that Eq. (9) completely accounts for all relevant contributions to the energy, in all test cases. We then consider the different contributions to Eq. (9) case by case.
5.1. Results with the test fields
The top part of Table 1 refers to the test fields B_{test}. In general, the energies E of the different test fields go from the purely potential case of B_{DD}, where , to highfreeenergy cases (B_{TD} and B_{Ex1}, with ), where the field is strongly nonpotential. The main source of violation of Thomson’s theorem, Eq. (2), in all cases is the mixed currentpotential term , except for the B_{Ex2PP} case where is slightly higher in absolute value than .
More precisely, B_{DD} is nearly perfectly potential, with nonsolenoidal spurious fluctuations contributing to the total energy for few parts per thousand at most (in ). B_{TD} has a 16%energy contribution from the currentcarrying part of the field , with a 2% contribution from the nonsolenoidal field related to the currentcarrying structure (in but not in ). This is the effect of the approximate nature in the matching between currentcarrying and external potential fields in the equilibrium that defines the B_{TD} field. B_{MHD}, which has 6% free energy , has an even lower nonsolenoidal contribution (− 0.1%). In all three cases, there is very small (B_{TD}) or no significant (B_{DD}, B_{MHD}) violation of Thomson’s theorem.
We now move to the NLFFF extrapolations. These show values of ⟨ f_{i} ⟩, which are twotothree orders of magnitude greater than in the first three cases. The contribution of the nonsolenoidal part of the potential field to the total energy, , is always negligible with respect to the other terms. In the B_{Ex2PP} case, the free energy associated with the solenoidal part of the currentcarrying field is about 11%, and the potential field energy is 88% of the total energy. The sum of the potential and currentcarrying solenoidal parts accounts for 99% of the total energy, apparently verifying Thomson’s theorem accurately. However, is 14% and is − 12%; i.e., the errors related to the divergence of the currentcarrying part of the field have comparable magnitudes and compensate for each other. These are the dominant sources of error, almost three orders of magnitude more than .
The test case with the highest value of ⟨ f_{i} ⟩ is B_{Ex1}. With respect to the B_{Ex2PP} case, B_{Ex1} is characterized by three times higher free energy , twice the error on current , and almost a four times larger error on . Again, the last two are largely compensating each other. We conclude that the interpolation to one third of the resolution used for B_{Ex1} is less efficient than preprocessing (used for B_{Ex2PP}) in eliminating the source of violation of Thomson’s theorem.
This situation is even more extreme in the case of the extrapolation of the nonpreprocessed, noninterpolated magnetogram B_{Ex2}. Although this case has a value of the mean divergence ⟨ f_{i} ⟩ that is only a factor two higher than for B_{Ex2PP}, and not even the highest one, it shows the most pathological behavior: The potential field has an energy 2.29 times the energy of the test field, which is downright unphysical according to Eq. (2). Such a high value is compensated for by an equally high value of (− 2.38). On the other hand, the currentcarrying part of the field accounts for 14% of the energy, but the associated error is more than six times larger. Such large errors are related to the high values of the divergence – in particular at the bottom boundary – and their actual values are very sensitive to the numerical details of the computation. Additional analysis of B_{Ex2} and B_{Ex2PP} is discussed in Sect. 8.1.
We finally notice that in the preprocessed case B_{Ex2PP}, the error from or might be considered as still tolerable if compared with the total energy (errors on vector magnetograms are similar, after all), but it seriously compromises the reliability of the free energy estimation, each one being as high as itself.
5.2. Results with the cleaned test fields
We now consider the bottom part of Table 1 for the solenoidal fields. The values of the divergence are drastically reduced in all cases to 10^{17} or less, which shows that the cleaner in Appendix B is an effective – and fast – way of removing the nonsolenoidal component of a discretized magnetic field. For the purpose of this article, we can then consider all B_{test,s} to numerically be perfectly solenoidal. All error terms, i.e., , , and , are smaller than 1%, and we recover Eq. (2) in a numerical sense.
More precisely, the B_{DD,s} and B_{MHD,s} cases are practically identical to their corresponding test fields, as far as the energy metrics E, , and are concerned. On the other hand, B_{TD,s} shows an increase of about 1.3% of the total energy, E, as a result of the removal of the error in of B_{TD}. The error removal affects the potential field energy E_{p} more, which raises about 5% with respect to the energy of the potential field in B_{TD} (in nonnormalized values), as a consequence of the cleaner’s modification of the bottom boundary. In contrast, the relative contribution of the currentcarrying part is unaffected by the cleaner. It is true that ⟨ f_{i} ⟩ differs by 15 orders of magnitude between B_{TD} and B_{TD,s}, but it is significant anyway that the removal of a 2%error in changes the nonnormalized values of the total energy E and potential field energy E_{p,s} of 1% and 5%, respectively. We conclude that, even in relatively divergencefree fields, residual nonsolenoidal effects can be energetically significant.
In the extrapolated cases, the removal of the larger divergence has far stronger consequences. In the first place, the nonnormalized field energy E of the cleaned fields B_{Ex1,s}, B_{Ex2PP,s}, and B_{Ex2,s} is increased of 42%, 109%, and 38%, respectively, with respect to those of the corresponding test fields. As a consequence of the higher values of E, the importance of potential fields relative to the total energy is decreased (to 78%, 95%, and 40% of their testfield values, respectively). In contrast, the energy contribution related to the currentcarrying part of the field is strongly increased, as expected, since the cleaner introduces currents that are related to the cumulated divergence that is removed (see Appendix B).
We conclude that the cleaned fields that are obtained from the test ones using the method in Appendix B all comply with Thomson’s theorem accurately. However, three of them, namely B_{Ex1,s}, B_{Ex2PP,s}, and B_{Ex2,s}, are energetically very different from the original fields B_{Ex1}, B_{Ex2PP}, and B_{Ex2}, respectively. Incidentally, we notice that the removal of the finite divergence does not conserve the approximate forcefreeness of the extrapolated fields.
5.3. Contributions to for the test fields
Contributions to in Eqs. (8), (9).
In many of the test fields in Table 1, is the largest source of error. Table 2 shows the six contributions to in the order in which they appear in Eq. (8) and their sum for the six test cases B_{test}. We do not consider the solenoidal fields B_{test,s} since all terms are mostly zero and never bigger than 0.7% The following conclusions can be drawn. First, the main contribution to is in all cases. The main source of violation of Thomson’s theorem is then the divergence of the currentcarrying part of the field. More often than not, this term has a similar magnitude and opposite sign of , which is positivedefinite. However, there is no obvious reason for to be always – or predominantly – negative, and we regard this as a coincidence.
Second, the terms with residual divergence of the potential field (i.e., any term containing ∇ζ in Eq. (8)) are always negligible. Therefore, also in view of the always low values in Table 1, we can conclude that the divergence of the potential field always gives a negligible contribution to the energy.
Third, the integral , and the integral have finite values for the extrapolations in Table 2. Analytically, they should be vanishing. Using the divergence theorem, Eq. (4), the surface integral vanishes because , and the volume integral vanishes because ∇·B_{J,s} = 0. The first condition is enforced at the boundary, but the second is only approximately true numerically (see also Sect. 6.1). This is not enough to insure that and vanish numerically. This is why we adopted the decomposition of the energy of Sect. 3.3 that only contains volume integrals.
6. Source of errors in the decomposition
6.1. Values of ⟨ f_{i} ⟩ for the field decomposition in Eqs. (3)–(6)
In this section we quantify how accurate the decomposition in solenoidal and nonsolenoidal contributions is. Table 3 reports the values of the logarithm of ⟨ f_{i} ⟩, defined by Eq. (C.2), for the field decomposition used in Eq. (7). Since ⟨ f_{i} ⟩ is not additive in the field, its value for, say, B is different from the sum of its values for the potential B_{p} and currentcarrying B_{J} components.
We next consider the decomposition of the potential field given by Eq. (5) for the test fields B_{test} (upper half of Table 3). Values of ⟨ f_{i} ⟩ for the solenoidal part of the potential field B_{p,s} are better (i.e., more negative) than those for B_{p,ns}, so that the B_{p,s} is indeed more solenoidal than B_{p,ns}. However, it is only in the first three cases, B_{DD}, B_{TD}, and B_{MHD}, that log _{10}(⟨ f_{i} ⟩) has a noticeably more negative value for B_{p,s} than for B_{p}. In the other cases, the values are relatively close to each other, and B_{p,s} is only marginally more solenoidal than B_{p}. On the other hand, B_{p,ns} is always much less solenoidal than B_{p}. This is partly the effect of the nonadditivity of the metric ⟨ f_{i} ⟩, and partly because B_{p,ns} is, on average, much smaller than B_{p,s}, as the corresponding energy metrics in Table 1 show. (In particular, , i.e., the energy associated with B_{p,ns}, is always extremely small.)
Similar conclusions can be drawn looking at the decomposition of the currentcarrying part, B_{J}, where this time the energy associated with the solenoidal error (see in Table 1) is more significant. In this case, values of ⟨ f_{i} ⟩ for all three contributions B_{J}, B_{J,s}, and B_{J,ns} are of similar magnitude. Again, the nonsolenoidal part, B_{J,ns}, has a higher divergence value than the solenoidal one, B_{J,s}, but only marginally so for B_{DD} and extrapolated fields.
We consider the solenoidal test fields B_{test,s} (bottom half of Table 3). Values of ⟨ f_{i} ⟩ for a given field component belonging to B_{test} and to the corresponding B_{test,s} are very similar. For instance, the value of log _{10}(⟨ f_{i} ⟩) for B_{J,s} in, say, the test field B_{Ex1} is − 2.36, whereas for the corresponding contribution for B_{Ex1,s} it is − 2.11. Therefore, the above discussion of the contributions to B_{test} holds for those of B_{test,s} as well. In contrast, the total divergence of the field is very different in the two cases, i.e., − 2.4 and − 17.2, respectively. This is a clear indication that the accuracy of the field decomposition is determined by the accuracy in the solution of Eqs. (3)–(6) rather than by the divergence of the total field.
In conclusion, the Poisson solver provides a decomposition of the magnetic field where the solenoidal parts have a smaller divergence than the original field, as required. The limit in the accuracy of the decomposition comes from the accuracy of the solver, and not from the level of solenoidality of the initial field. One possible source of inaccuracy for the solver is the incompatibility of the boundary conditions used in Eqs. (3)–(6), which is discussed in the next section.
6.2. Compatibility of boundary conditions in Eqs. (3)–(6)
Removal of flux imbalance.
We here consider the normalized flux of the field, , computed as the surface flux through all six boundaries, normalized to the mean flux entering and exiting from the lower boundary: (12)The values of in Table 4 show that the test fields of the extrapolation cases B_{Ex1}, B_{Ex2PP}, and B_{Ex2} are not fluxbalanced. Therefore, the decomposition of Eq. (7) based on the solutions of Eqs. (3)–(6) may be inconsistent (see Eq. (11) and related text). The purpose here is to determine whether the unbalanced flux affects the accuracy of any of the terms in Eq. (7).
A fluxbalanced field, , can be computed from a fluxunbalanced one, B_{test}, by splitting the original field as (13)and assuming B_{Φ} = ∇Θ to be generated by an uniformly distributed, constant divergence; i.e., ΔΘ = constant. We choose the simple solution Θ ∝ r^{2}, and fix the constant such that the flux of B_{Φ} equals the flux of B_{test}, yielding Table 4 shows that the above modifications to B_{test} is effective, drastically reducing the net flux of the original field to very low values (compare with ). On the other hand, the effect on the field of B_{Φ} is very small. Both energy terms related to that (i.e., and ) are negligible (with a contribution below 0.01% of the total energy E computed for ).
Repeating the same analysis of Sects. 5 and 6.1 for the fluxbalanced part of the field only, , yields no significant change: all values in Tables 1–3 are identical. Inaccuracies of the Poisson solver in solving Eq. (3) are therefore related to the solver itself, not to the incompatibility of the boundary conditions.
In a similar way, the test field can be modified to have vanishing volume divergence, which is the requirement for consistency in solving Eqs. (5), (6), using The result is likewise clear: no significant change is found in the values of Tables 1–3.
Therefore, an imperfect consistency of source and boundary conditions play no role in the accuracy of the solution of the Laplace and Poisson equations employed in the decomposition, Eq. (7), for any of the test cases. Recalling the results of Sect. 6.1, we conclude that the accuracy limitation of our analysis comes from the solver itself. In this respect, we note that, when the method used in Eqs. (5) and (6) is viewed as an algorithm for removing the divergence (Projection method), it is far less efficient than our divergence cleaner described in Appendix B. On the other hand, the projection method has other advantages; for instance, it change neither the current in the volume nor the normal component of the original field at the boundaries.
7. Parametric study
In this section we study how the relative energy of the field depends on its divergence in progressively going from a solenoidal to nonsolenoidal realizations. The purpose is to offer a practical method of fixing the level of solenoidal errors that can be tolerated in a given numerical realization, based on their consequences on the energy of the field.
Fig. 2
Magnetic energy, Eq. (16), normalized to the energy of the potential field having the same distribution of normal field on the boundaries (black lines) or to the energy of the reference field (orange lines), as a function of the amplitude δ of the nonsolenoidal term, Eq. (14). Two models for the divergence are shown according to Eq. (15). Each panel shows the results of a test field: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. An important change of scale of both axes is present between the top and the bottom rows. 
7.1. Parametric models of finitedivergence fields
For a given test magnetic field B_{test}, the corresponding solenoidal field B_{test,s} is considered. A parametric, nonsolenoidal field B_{δ} is obtained by adding a nonsolenoidal component B_{div} to B_{test,s}, using a control parameter δ, as (14)We consider here two models of B_{div}, namely (15)Adding the first divergence model for δ = 1 is the inverse operation of the cleaner in Sect. B, since B_{δ}(δ = 1) = B_{test}. For other δ values, the resulting field B_{δ} only differs from the solenoidal field B_{test,s} in the zcomponent. The second model for the divergence is more general, because it changes all three components of B_{test,s} in the volume, although not on all boundaries.
Both divergence models in Eq. (15) are based on the computed ∇·B_{test}. In this way, we relate the divergence models to the source of error that is specific to the considered test field. For instance, we expect that errors in the test case B_{DD}, which are only generated by truncation errors, have a different distribution in space than those coming from the approximate nature of the B_{TD} equilibrium, or from a numerically constructed field like B_{Ex2PP}.
The influence of a finite divergence of B_{δ} on the energy value can be written as (16)where E, E_{div}, and E_{test,s} are defined as usual as proportional to the volume integrals of , , and , respectively, and .
Below, the energy dependence on δ is studied for the test fields described in Sect. 4. Since the separation in solenoidal and nonsolenoidal components is known by construction, we simplify the presentation by analyzing the energies of the total fields according to Eq. (16), and we do not separate the error sources as in Sect. 5. For each value of δ, we consider B_{δ} as the test field to analyze and compute the corresponding potential field according to Eq. (3).
Fig. 3
Magnetic energy, Eq. (16), normalized to the energy of the corresponding potential field (black lines, zoom of Fig. 2), and energy associated to the nonsolenoidal part of B_{J}, normalized to the energy of B_{test} (red lines, equal to E_{J,ns}/E in the notation of Eq. (8)), as a function of the amplitude δ of the nonsolenoidal term, Eq. (14). Two models for the divergence are shown according to Eq. (15). Each panel shows the results of a test field: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. The black dashdotted line at E/E_{p} = 1 marks the value below which the solution is unphysical. A large change of scale of both axes is present between the top and the bottom rows. 
7.2. Parametric dependence of the energy
Figure 2 shows the energy for the two divergence models in Eq. (15), as a function of the control parameter δ in a wide range of values. Due to the large difference in ∇·B_{test} between the six models, the top and bottomrows have different scales. The orange lines show the energy normalized with the energy of B_{test}, which is not dependent on δ: They follow the expected parabolic profile of Eq. (16), only scaled by the normalization factor. Model 1 (continuous orange lines) yields a smaller variation of the energy with δ (corresponding to lower values of E_{div}) with respect to Model 2 (dashed orange lines), and is centered farther away from δ = 0 (i.e., Model 1 has higher values of E_{s,div}/E_{div}).
The orange curves in the top row of Fig. 2 show that it takes very high values of δ in order to have a variation of order one of the energy in the B_{DD}, B_{TD}, and B_{MHD} cases (e.g., for the Model 2 applied to B_{TD} at δ = 15). On the other hand, the energy of the extrapolated fields shows a much steeper increase with δ, related to the much higher value of ∇·B_{test}, and particularly so for Model 2.
The location of the minimum of each of the orange curves is at δ_{min} = −E_{mix}/E_{div}, therefore its location depends on the average orientation and amplitude of the divergence field B_{div} with respect to the solenoidal field B_{s}. The orientation and amplitude of B_{div} also determines the height of the minimum (since the energy of the test field is fixed). With both divergence models, there are no general rules; i.e., the energy can increase or decrease with δ, and the location of the minimum depends on the case.
7.3. Comparison with the potential field energy
The physically meaningful quantity is represented by the energy normalized to the energy of the corresponding potential field, represented in Fig. 2 by black lines. For different δ values, the normal component of the field B_{div} at the boundary changes according to Eq. (15), hence also the energy of corresponding potential field depends – quadratically – on δ. Due to the additional δdependence, the shape of the black lines is not always parabolic in the six cases, and the actual profiles depend on the details of the spatial distribution of divergence in the test field.
To show that, we first notice that the two divergence models behave very differently, except for B_{DD} where the range in δ is too narrow to show significant differences. For instance, E/E_{p} of Model 1 (continuous black lines) is an increasing function of δ in the range (− 15, 15) in the B_{TD} case. Model 2, on the other hand, has a parabolic energy profile with minimum at δ ≈ −4. For both models, the energy variation is relatively large (1.8 and above 2 for Models 1 and 2, respectively), whereas the variation in the same range of δ is smaller for the B_{DD} and B_{MHD} cases.
The extrapolated cases yield not only much larger variations (note again the difference in scales between the top and bottom rows of Fig. 2), but also a stronger dependence on δ. In particular, Model 2 yields a relative energy that sharply increases with δ, for instance, to one order of magnitude increase for δ going from the value 0 to 1 in the B_{Ex2PP} case. A saturation at high values of δ is clearly visible in the dashed black line (Model 2) of the B_{Ex1} case, and is hinted at in the B_{Ex2PP} case. Such saturation is actually present in all three extrapolated cases, yielding values that are higher than those shown in the corresponding plots. The saturation happens when the quadratic dependence on δ of the energy of potential field compensates the quadratic term δ^{2}E_{div}.
Fig. 4
Power spectra of the twodimensional fields B_{x} (continuous line) and B_{z} (dashed line), summed over all wave numbers k_{y}, for the three cases B_{TD} (orange), B_{Ex2PP} (red), and B_{Ex2} (black). Left: at the bottom boundary. Right: at the tenth pixel in height. Spectra are normalized to their maximum value, and the spatial resolution is taken to be 1 in both directions (i.e., the wave number k_{x} has the dimension of pixel^{1} and is normalized to the total number of modes). 
On the other hand, Model 1 shows a more complex dependence on δ, which is shown in magnified scale by the black lines in Fig. 3. Counterintuitively, the largest variation in the relative energy E/E_{p} as a function of δ is found for B_{Ex2PP}, i.e., for the extrapolation case, which satisfies Thomson’s theorem better, see Table 1. The continuous black lines in Fig. 3d,f show the presence of one maximum and one minimum in the considered range of values of δ (for the B_{Ex2PP} case, these lie outside the considered range), implying that, at high values, the potential field energy grows faster than the total energy. The location of the extrema is different in the three B_{Ex1}, B_{Ex2PP}, and B_{Ex2} cases, and in none of the cases are the extrema found for the solenoidal (δ = 0) or the test (δ = 1) configurations. In general, the maximum and minimum energy configurations depend on the spatial distribution of the divergence of the test field, through E_{s,div}.
7.4. Unphysical cases.
The black dashdotted line at E/E_{p} = 1 in Fig. 3 is the value below which unphysical fields are obtained. We find that only Model 1 can produce unphysical solutions, and only for specific range of δ values in the B_{DD}, B_{TD}, B_{Ex2} cases. The latter case is known from the value of E_{p}/E in Table 1, and is considered to be an extreme case because of the large divergence that it involves. However, the possibility of also creating unphysical solutions in the far more solenoidal field B_{DD} and B_{TD} (for values of  δ > 5) is unexpected. It confirms that not just the value of the divergence is important, but also its detailed spatial distribution with respect to the solenoidal component, as evident from E_{s,div}. It is the alignment between B_{div} and B_{s}, and not just the magnitude of B_{div}, which determines how strongly the energy depends on δ. Moreover, while E > E_{p} is always satisfied for B_{Ex1} and B_{Ex2PP}, the minimum value of E is close to E_{p} (see Fig. 3d,e), showing that unphysical fields may be found relatively easily in NLFFF extrapolations.
From Table 2 and the related discussion of Sect. 5.3 we showed that the main source of violation of Thomson’s theorem is the term E_{p,s/J,ns} in E_{mix}. The dependence on δ of this term, normalized to the energy of the test field, is shown by the red curves in Fig. 3 for both models of divergence (Eq. (15)). The contribution to the total energy is negligible in the B_{DD} and B_{MHD} cases, and can be a few percent for large δ in the B_{TD} case.
In the extrapolated cases, the dependence of E_{p,s/J,ns} on δ is linear for Model 1 and parabolic for Model 2. In Model 1, the steepness of the linear curve increases, going to B_{Ex2PP} to B_{Ex1} and B_{Ex2}, as expected (see Table 1 and related text). The amplitude of the error is two to three orders of magnitude larger than in the B_{DD}, B_{TD}, and B_{MHD} cases. In the B_{Ex2PP} case the error is smaller, but it is still about a factor 20 larger than in B_{TD} for δ = 5.
If we consider the black curves in Fig. 3 for δ = 0, we can identify the energy of the solenoidal field as a natural reference value for the free energy. Starting from this reference value, for increasing  δ , the linear contribution of E_{p,s/J,ns}, together with the quadratic change in the potential field energy, creates the maximum and minimum values of E/E_{p}. If the linear contribution is large enough, the minimum lies below the threshold E/E_{p} = 1, and there is a range of values where the solution is unphysical. For even higher values of  δ , the quadratic dependence of the potential field energy dominates E. From this point onward, E_{p,s/J,ns} is not the main source of error in Eq. (7).
More generally, a parametric study like the one in Figs. 2 and 3 can be used to identify what is the level of divergence (i.e., the level of E_{mix}or E_{p,s/J,ns}) that can be tolerated and which is the threshold above which the solution becomes entirely unphysical (i.e., with E/E_{p} < 1).
In conclusion, the parametric study shows that the energy may be severely influenced by the solenoidal property of the field. The effect depends not only on the amplitude of the nonsolenoidal component, but also on the specific average orientation of the nonsolenoidal component with respect to the solenoidal one (directly affecting E_{s,div} in Eq. (16)). As a result, a singlenumber divergence metric, such as ⟨ f_{i} ⟩, is insufficient to deduce what errors should be expected in the energy. A more proper indication is found by the numerical verification of Thomson’s theorem (Sect. 5) and by a parametric study as presented in this section.
8. Source of divergence in NLFFF extrapolations
We now investigate in more detail some of the test fields discussed in Sect. 5, with emphasis on the reason for the large divergence that leads to violating Thomson’s theorem. The main source of error comes, in almost all the cases, from the mixed term E_{mix}, and is associated with the nonsolenoidal component of the currentcarrying part of the field. Also, there are markedly larger errors in the extrapolated test fields, B_{Ex1}, B_{Ex2PP}, B_{Ex2}, than in B_{DD}, B_{TD}, and B_{MHD}. Finally, the preprocessing of the vector magnetogram before extrapolation yields more solenoidal fields, whereas a simple averaging does not seem to be enough for removing errors, and yields a more severe violation of Thomson’s theorem (Eq. (2)).
8.1. Analysis of small scales
One main difference among the B_{test} cases in the upper half of Table 1 is the length scale of the magnetic field: While the first three cases are smooth fields with a magnetic field variation spanning several times the spatial resolution, the extrapolated cases have large variations on the pixel scale, especially at the bottom boundary, i.e., on the vector magnetogram that is used as a boundary condition for extrapolations. This is true to a different degree for the three cases: For B_{Ex1} the vector magnetogram was interpolated (with a fluxconserving average) at a resolution of about one third that of B_{Ex2} and B_{Ex2PP}. Such an interpolation smooths part of the small scale away, yielding results that are closer to the B_{Ex2PP} case rather than to the B_{Ex2} one. B_{Ex2PP} is not interpolated, but it is preprocessed, an operation that includes an explicit smoothing of smaller scales, especially on the transverse components. Finally, B_{Ex2} has neither interpolation nor preprocessing, and it retains all the small scales that are present at the full resolution of the Hinode/SOT vector magnetograms.
As an example, Fig. 4 shows the power spectrum of the x and zcomponents of the fields B_{TD}, B_{Ex2PP}, and B_{Ex2}, at two different heights as a function of the normalized wave number k_{x}. The lefthand panel of the figure shows that, at the bottom boundary, B_{TD} has power spectra that decrease rapidly with k_{x}, in both components. In contrast, the power spectra of B_{Ex2PP} and B_{Ex2} have higher values on all scales, which are particularly strong in the vertical component.
Ten pixels above the bottom boundary (right panel in Fig. 4), the B_{TD} power spectrum is essentially the same as at z = 0 because both planes cut through the flux rope, so a similar magnetic structure is present. In contrast, B_{Ex2PP} on the upper plane has a much more peaked spectrum, except for the distribution tail on the smallest scales which is basically as strong as at the bottom boundary. Such a component on the shortest scales comes from the forcefree condition that is enforced by the extrapolation code, which propagates into the volume the small scales that are present at the bottom boundary.
We now consider the difference between preprocessed case B_{Ex2PP} and the nonpreprocessed one B_{Ex2}. The difference in ⟨ f_{i} ⟩ between the two is about a factor 2, and it is large in the other energy metrics in Table 1. The comparison between the normalized spectra of B_{Ex2PP} and B_{Ex2} in Fig. 4 shows that there are comparable (relative) energies on small scales in both cases. Actually, by locally changing the magnetic field at the bottom boundary to enforce forcefree compatibility, preprocessing increases the small scales. The smoothing term that is present in preprocessing only has a limiting effect on such an increase. Therefore, the two cases B_{Ex2PP} and B_{Ex2} do not differ strongly as far as the presence of small scales is concerned, while Thomson’s theorem is much better satisfied for B_{Ex2PP} than for B_{Ex2} (see Sect. 5.1).
The cleaned test fields B_{test,s} are numerically solenoidal, and there is no violation of Thomson’s theorem. However, in these cases, too, small scales are increased (not shown), since the divergence cleaner introduces extra electric currents that are related to derivatives of the divergence of the original field, see Eq. (B.5). This is an additional confirmation that the presence of small scales as such is not directly at the origin of the violation of Thomson’s theorem.
8.2. Role of small scales and preprocessing
Valori et al. (2010) show that the NLFFF extrapolation of the B_{TD} vector magnetogram yields a very accurate reconstruction of the whole test field, which is also solenoidal to a very high degree. On the other hand, there is a large difference in the scale distribution between smooth fields like the B_{TD} and the extrapolated fields.
The presence of small scales inside the volume, which are induced by the small scales at the boundary, may not be correctly approximated by the discretization employed in extrapolation code, yielding local violation of the solenoidal constrain. However, when the extrapolation from a preprocessed magnetogram is considered, the extent of the violation of Thomson’s theorem is greatly reduced, even though small scales are actually increased. By partially enforcing forcefree compatibility on the bottom boundary, the preprocessing provides the extrapolation code with a boundary condition that is more compatible with the forcefree equations. Since extrapolation codes attempt to construct a solution of the forcefree equations that is simultaneously force and divergencefree, the more compatible the boundary, the more consistent (i.e., force and divergencefree) the obtained solution. Conversely, when the boundary condition is incompatible with the forcefree equation, the reduction of the Lorentz forces is at the expense of the solenoidal condition. In such cases, the divergence of the solution is higher, and Thomson’s theorem is more severely violated. We thus conclude that the incompatibility of the boundary condition with the forcefree condition is at the origin of the difference in the errors E_{J,ns} and E_{mix} between B_{Ex2PP} and B_{Ex2}.
We notice that preprocessing is a parametric method that can produce progressively more forcefreecompatible vector magnetograms for higher values of the employed parameters, at the price of larger modifications of observed values. The energy values and their relative errors therefore vary continuously as a function of the preprocessing parameters, quite independently of the particular extrapolation method that is employed (see, e.g., Schrijver et al. 2008; Metcalf et al. 2008). No unequivocal method is available in order for determining the best parameters to use (see, e.g., Wiegelmann et al. 2006, 2012; Fuhrmann et al. 2011), which leaves energy estimations subjected to uncomfortable arbitrariness.
9. Conclusions
Thomson’s theorem states that the energy of a magnetic field is given by the sum of the energy of the currentcarrying part of the field plus the energy of the potential field that has the same distribution of the normal component on the boundary of the considered volume. The field must be perfectly solenoidal for the theorem to be valid. Such a condition is often only approximately satisfied in numerical simulations, such as in MHD simulations and NLFFF extrapolations. However, it is a nontrivial task to identify a quantitative estimation of solenoidal errors that can be applied to different discretizations of magnetic fields, essentially due to the nonlocal consequences that such errors produce. Our goal has been to develop physically meaningful metrics and practical methods that can be used to judge whether the solenoidal property is fulfilled with sufficient accuracy.
To this aim, we introduced a decomposition of the energy of a discretized field into solenoidal and nonsolenoidal contributions that allowed an unambiguous and numerically welldefined estimation of the effect of the divergence in terms of associated energies. Moreover, we introduced a method of parametrizing the divergence that allows for an exploration of the nonsolenoidal effects.
In this way, the numerical verification of Thomson’s theorem offers an operational and quantitative way of checking the reliability of energy estimations in numerical computations. Since the violation of Thomson’s theorem is solely determined by the presence of magnetic charges, it is at the same time a quantitative estimation of the importance of solenoidal errors.
We applied our method to six different test cases, covering a representative sample of numerical realizations. Of the six test cases considered here, two of them (the dipolar field B_{DD} and a snapshot of an MHD simulation of nullpoint reconnection B_{MHD}) presented negligible violations, and one (a forcefree current ring B_{TD}) offered only a moderate one that, however, has finite effects on the energy. In the case of an NLFFF extrapolation of a preprocessed vector magnetogram (B_{Ex2PP}), the sum of the potential energy E_{p} and free energy E_{J,s} is very close to the total energy E, and one could draw the conclusion that almost no violation of Thomson’s theorem occurs. However, by separating all contributions in Eq. (7), our analysis showed compensating energy contributions (E_{J,ns} and E_{mix}) that are close to E_{J,s}. If the most conservative view is adopted by considering errors in absolute values, then the opposite conclusion must be drawn: The violation is large enough to compromise the estimation of the free energy, since both E_{J,ns} and E_{mix} are on the order of the free energy value E_{J,s}. The last two cases we studied, also NLFFF extrapolations but of nonpreprocessed magnetograms (B_{Ex1} and B_{Ex2}), represent cases with very large errors.
The energy of the potential field E_{p,s} is the reference value for the free energy. In our applications, the inaccuracy in its determination, E_{p,ns}, is practically never significant. The currentcarrying part of the field is responsible for the largest errors instead.
The parametric study shows that the amplitude of the nonsolenoidal component is not the only factor that generates errors in the energy. The average orientation of the nonsolenoidal component with respect to the solenoidal one (affecting directly E_{s,div} in Eq. (16)) plays an even more important role. Indeed, even using a relatively solenoidal discretized magnetic field (like B_{TD}), it is possible to create configurations where the energy of the field is lower than that of the corresponding potential field. Such unphysical solutions have also been found in some cases of NLFFF extrapolations.
More generally, in NLFFF extrapolations the energy of the reconstructed field was found to vary according to the extent of the modification that was enforced on the vector magnetogram that is used as boundary condition (by preprocessing, i.e., by smoothing and/or by enforcing forcefreeness compatibility). Our study shows quantitatively the effect of such practices on the energy, and makes it clear that the origin of the variability (and errors) in energy estimations based on NLFFF extrapolations is the presence of a large divergence, which is eventually caused by the lack of compatibility between the equations solved (solenoidal forcefree field) and the photospheric boundary conditions, rather than by noise or the small scales present in the vector magnetogram.
Finally, the parametric study is based on a numerically solenoidal field that is obtained from a given, nonsolenoidal one. We introduced a method for the complete removal of the nonsolenoidal component of a discretized field. At the price of changing boundary values and the current density, this method provides a field that is solenoidal to numerical precision. When the solenoidal versions of the test fields are considered, the Thomson theorem is found to be fulfilled with more than 99% accuracy.
We concluded that testing Thomson’s theorem in numerical realizations of magnetic fields is a powerful method quantifying the amount of nonsolenoidal contributions to a numerical magnetic field. In particular, it allows assessing the reliability of free magnetic energy estimations, a crucial quantity in phenomena such as flares and coronal mass ejections. To this purpose, we proposed a set of analytical and numerical tools that allowed us to fully test the reliability of numerical magnetic fields. Such a set includes a method for removing the divergence from a given discretized field, to numerical precision. The effect of larger and larger divergence contributions is studied by parametrically adding a known divergence to the numerically solenoidal field. In this way, it is possible to monitor the effect of the nonsolenoidal part of the magnetic field and to quantify its effect in terms of magnetic energy. Our method can be applied to any discretization of magnetic fields, e.g., in MHD simulations and in
NLFFF extrapolations, to constrain quantitatively errors due to violation of the solenoidal property.
Acknowledgments
The authors are pleased to thank Guillaume Aulanier for fruitful discussions, Bernhard Kliem and Tibor Török for making the numerical solution of the TD equilibrium available, and the referee for helpful comments that improved the clarity of the paper. G.V. is indebted to the NLFFF Consortium for stimulating discussions and collaborations. The research leading to these results has received funding from the European Commission’s Seventh Framework Program (FP7/20072013) under the grant agreement eHeroes (project n° 284461, www.eheroes.eu). S.M. gratefully acknowledges support from the NASA Postdoctoral Program, administrated by Oak Ridge Associated University through a contract with NASA, during her stay at NASA Goddard Space Flight Center. Calculations were done on the quadricore biXeon computers of the Cluster of the Division Informatique de l’Observatoire de Paris (DIO). SDO data are courtesy of NASA/SDO and the HMI science team. Hinode is a Japanese mission developed and launched by ISAS/JAXA, collaborating with NAOJ as a domestic partner, NASA and STFC (UK) as international partners. Scientific operation of the Hinode mission is conducted by the Hinode science team organized at ISAS/JAXA. This team mainly consists of scientists from institutes in the partner countries. Support for the postlaunch operation is provided by JAXA and NAOJ (Japan), STFC (UK), NASA (USA), ESA, and NSC (Norway).
References
 Amari, T., Luciani, J. F., Aly, J. J., Mikic, Z., & Linker, J. 2003, ApJ, 585, 1073 [Google Scholar]
 Asai, A., Ishii, T. T., Kurokawa, H., Yokoyama, T., & Shimojo, M. 2003, ApJ, 586, 624 [NASA ADS] [CrossRef] [Google Scholar]
 Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314 [Google Scholar]
 Aulanier, G., Janvier, M., & Schmieder, B. 2012, A&A, 543, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aulanier, G., Démoulin, P., Schrijver, C., et al. 2013, A&A, 549, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Berger, M. A. 1985, ApJS, 59, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Brackbill, J. U., & Barnes, D. C. 1980, J. Comput. Phys., 35, 426 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Chen, P. F. 2011, Liv. Rev. Sol. Phys., 8, 1 [Google Scholar]
 Emslie, A. G., Dennis, B. R., Shih, A. Y., et al. 2012, ApJ, 759, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y. 2010, ApJ, 719, 728 [NASA ADS] [CrossRef] [Google Scholar]
 Forbes, T. G. 2000, J. Geophys. Res., 105, 23153 [NASA ADS] [CrossRef] [Google Scholar]
 Forbes, T. G., Linker, J. A., Chen, J., et al. 2006, Space Sci. Rev., 123, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Fuhrmann, M., Seehafer, N., & Valori, G. 2007, A&A, 476, 349 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuhrmann, M., Seehafer, N., Valori, G., & Wiegelmann, T. 2011, A&A, 526, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Golub, L., & Pasachoff, J. M. 2009, The Solar Corona [Google Scholar]
 Karpen, J. T., Antiochos, S. K., & DeVore, C. R. 2012, ApJ, 760, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Lawrence, C. E. 1998, Partial Differential Equations (American Mathematical Society) [Google Scholar]
 Liu, R., Liu, C., Wang, S., Deng, N., & Wang, H. 2010, ApJ, 725, L84 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C. 2001, J. Geophys. Res., 106, 25141 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, S., Aulanier, G., Pariat, E., & Klein, K.L. 2012, Sol. Phys., 276, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalf, T. R., Derosa, M. L., Schrijver, C. J., et al. 2008, Sol. Phys., 247, 269 [Google Scholar]
 Priest, E. R. 2003, in Solar magnetohydrodynamics, ed. B. N. Dwivedi, 217 [Google Scholar]
 Savage, S. L., McKenzie, D. E., & Reeves, K. K. 2012, ApJ, 747, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., & Zwaan, C. 2008, Solar and Stellar Magnetic Activity (Cambridge Astrophysics) [Google Scholar]
 Schrijver, C. J., DeRosa, M. L., Metcalf, T., et al. 2008, ApJ, 675, 1637 [Google Scholar]
 Shibata, K., & Magara, T. 2011, Liv. Rev. Sol. Phys., 8, 6 [Google Scholar]
 Török, T., & Kliem, B. 2003, A&A, 406, 1043 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tajima, T., & Shibata, K. 2002, Plasma astrophysics (Boulder, CO, USA: Perseus Books) [Google Scholar]
 Taylor, J. B. 1986, Rev. Mod. Phys., 58, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707 [NASA ADS] [Google Scholar]
 Valori, G., Kliem, B., Török, T., & Titov, V. S. 2010, A&A, 519, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valori, G., Green, L. M., Démoulin, P., et al. 2011, Sol. Phys., 374 [Google Scholar]
 Valori, G., Démoulin, P., & Pariat, E. 2012, Sol. Phys., 278, 347 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wheatland, M. S., Sturrock, P. A., & Roumeliotis, G. 2000, ApJ, 540, 1150 [Google Scholar]
 Wiegelmann, T., Inhester, B., & Sakurai, T. 2006, Sol. Phys., 233, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., Thalmann, J. K., Inhester, B., et al. 2012, Sol. Phys., 281, 37 [NASA ADS] [Google Scholar]
Online material
Appendix A: Details of the numerical implementation
In the applications presented in this paper we have considered uniform Cartesian grids of resolution Δ in all directions, discretizing a rectangular volume (see Sect. 4 for the actual values of Δ and in each case). We compute derivatives using the standard secondorder, centraldifference operator, and we employ the relevant onesided (i.e., forward or backward), secondorder differences at the boundaries of . The only exception is the computation of the divergence of B_{test}, since all test fields are known in a volume that is larger than the selected (on lateral and top boundaries). In this case, ∇·B_{test} is computed using the central differences also at the location of the lateral and top boundaries of .
In the computation of volume integrals, the cell volume Δ^{3} is assigned to each internal node of the grid, whereas the cell volume is reduce to half, one fourth, and one eighth for nodes on the lateral surfaces, edges, and corners of , respectively. Similarly, in the computation of surface integrals, the cell surface Δ^{2} is assigned to each node inside each side of , whereas the cell surface is reduced to half and one fourth on edges and corners of each side, respectively. Despite the accurate computation of integrals, the divergence theorem, Eq. (4), is not insured to hold numerically, a property that requires special techniques, like finitevolume discretizations, to be fulfilled.
Appendix B: Divergence cleaner
To construct a numerically solenoidal field [B_{s}] from a field [B] let us define (B.1)where A is the vector potential computed from B in the volume . The vector potential A can be derived as in Valori et al. (2012) using the gauge ẑ·A = 0, yielding the expression (B.2)where b ≡ (A_{x}(x,y,z = z_{2}),A_{y}(x,y,z = z_{2}),0) is any solution of (B.3)A direct substitution of Eq. (B.2) into Eq. (B.1) shows that (B.4)with the property that ∇·B_{s} = 0. In other words, Eq. (B.4) naturally separates B into a solenoidal part B_{s} and a nonsolenoidal one, thus defining a divergence cleaner for B. The zcomponent of B is changed throughout the volume except on the top boundary, whereas the x and ycomponents are unchanged. The amplitude of the modification to B at a given height z is given by the cumulative effect of “magnetic charges” above that altitude.
Since only the zcomponent of the field is changed, the divergence cleaner changes the x and ycomponents of the current, but not the zcomponent, (B.5)
therefore the cleaner changes the injected magnetic flux but not the injected electric current through the bottom layer. On the other hand, since most of the test fields considered in this article have the highest values of divergence close to the bottom boundary, only the lower part of the field is changed significantly by the cleaner.
Computation B_{s} requires numerical computation of an integral of the type , as in Eq. (B.4) for f = ∇·B. To achieve numerical accuracy in the solenoidal property of B_{s}, G(z) must satisfy ∂_{z}G(z) = −f(z) numerically, i.e., must satisfy the numerical formulation of the fundamental theorem of integral calculus in the employed discretization. For the secondorder central differences that are used in the analysis, this can be obtained by the recurrence formulae (B.6)where G(z) = G(z_{1} + kΔ) ≡ G(k) with k = 0,1,2,··· ,(n_{z} − 1), and Δ is the uniform spatial resolution in z.
The constraint ∂_{z}G(z) = −f(z) in the secondorder, centraldifference discretization does not fix the value of G(n_{z} − 2). To do that, we require that the divergence of Eq. (B.4) also vanishes at the bottom boundary, i.e., (∇·B_{s})_{z = z1} = 0. Here the secondorder divergence operator is computed by using a secondorder, forward derivative in the zdirection, i.e., defining the operator , where and . By using the recurrence formula Eq. (B.6), the condition on the bottom boundary is transformed into the condition for G(n_{z} − 2), yielding where f_{os} = ∇_{os}·B. Such a numerical trick is only possible if the volume is discretized by an even number of points in the zdirection, therefore the analysis volumes employed in the article were chosen to satisfy such a requirement.
Appendix C: Measures of ∇ · B
The total divergence of a field B can be conveniently expressed by a single number using the average ⟨ f_{i} ⟩ over the grid nodes of the fractional flux (C.1)through the surface ∂v of a small volume v including the node i (Wheatland et al. 2000). Taking a cubic voxel of side equal to Δ as the small volume v centered on each node, the divergence in the discretized volume of uniform and homogeneous resolution Δ is then given by (C.2)where i runs over all N nodes in . This metric depends on the considered volume, so that values are strictly comparable only if computed on equal volumes.
All Tables
All Figures
Fig. 1
Selected field lines of the six test cases: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. The vertical component of the magnetic field at the bottom boundary is shown on a gray scale, with the positive (respectively, negative) polarity in white (respectively, black). The different line colors outline different types of connectivities. 

In the text 
Fig. 2
Magnetic energy, Eq. (16), normalized to the energy of the potential field having the same distribution of normal field on the boundaries (black lines) or to the energy of the reference field (orange lines), as a function of the amplitude δ of the nonsolenoidal term, Eq. (14). Two models for the divergence are shown according to Eq. (15). Each panel shows the results of a test field: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. An important change of scale of both axes is present between the top and the bottom rows. 

In the text 
Fig. 3
Magnetic energy, Eq. (16), normalized to the energy of the corresponding potential field (black lines, zoom of Fig. 2), and energy associated to the nonsolenoidal part of B_{J}, normalized to the energy of B_{test} (red lines, equal to E_{J,ns}/E in the notation of Eq. (8)), as a function of the amplitude δ of the nonsolenoidal term, Eq. (14). Two models for the divergence are shown according to Eq. (15). Each panel shows the results of a test field: a) the potential field of a double dipole, B_{DD}; b) the TD model, B_{TD}; c) the MHD model, B_{MHD}; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, B_{Ex1}; e) the NLFFF model of the preprocessed magnetogram of AR 11024, B_{Ex2PP}; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, B_{Ex2}. The black dashdotted line at E/E_{p} = 1 marks the value below which the solution is unphysical. A large change of scale of both axes is present between the top and the bottom rows. 

In the text 
Fig. 4
Power spectra of the twodimensional fields B_{x} (continuous line) and B_{z} (dashed line), summed over all wave numbers k_{y}, for the three cases B_{TD} (orange), B_{Ex2PP} (red), and B_{Ex2} (black). Left: at the bottom boundary. Right: at the tenth pixel in height. Spectra are normalized to their maximum value, and the spatial resolution is taken to be 1 in both directions (i.e., the wave number k_{x} has the dimension of pixel^{1} and is normalized to the total number of modes). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.