© ESO, 2013

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 long-range, magnetic interaction dominates other forces, e.g., plasma pressure and gravitational forces. A typical example is the low-corona (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 current-free (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 current-carrying 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 low-coronal 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 inertia-dominated 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 current-free 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 post-event 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 post-event 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 post-event 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 force-free 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 magneto-hydrodynamic (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 CGS-Gaussian 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 current-carrying 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 current-carrying 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 current-carrying 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. Gauge-invariant 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 current-carrying 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 current-carrying 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 gauge-dependence.

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 current-carrying 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 nondivergence-free 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, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. 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 field-line 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 z-component 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 current-carrying 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 finite-divergence 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 magneto-hydrodynamic numerical simulation of magnetic reconnection in a null-point topology (Masson et al. 2012). To use our present-stage 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 magneto-hydrodynamic 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 magneto-frictional 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 force-free 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 y-components as B_test, whereas the z-component 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 divergence-cleaned 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 high-free-energy 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 current-potential 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 current-carrying part of the field , with a 2% contribution from the nonsolenoidal field related to the current-carrying structure (in but not in ). This is the effect of the approximate nature in the matching between current-carrying 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 two-to-three 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 current-carrying field is about 11%, and the potential field energy is 88% of the total energy. The sum of the potential and current-carrying 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 current-carrying 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 current-carrying 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 non-normalized values), as a consequence of the cleaner’s modification of the bottom boundary. In contrast, the relative contribution of the current-carrying 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 non-normalized values of the total energy E and potential field energy E_p,s of 1% and 5%, respectively. We conclude that, even in relatively divergence-free 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 non-normalized 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 test-field values, respectively). In contrast, the energy contribution related to the current-carrying 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 force-freeness of the extrapolated fields.

5.3. Contributions to for the test fields

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 current-carrying part of the field. More often than not, this term has a similar magnitude and opposite sign of , which is positive-definite. 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 current-carrying 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 ₁₀(⟨    |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 current-carrying 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 ₁₀(⟨    |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)

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 flux-balanced. 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 flux-balanced field, , can be computed from a flux-unbalanced 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², 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 flux-balanced 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, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. An important change of scale of both axes is present between the top and the bottom rows.

7.1. Parametric models of finite-divergence 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 z-component. 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 BJ, normalized to the energy of Btest (red lines, equal to EJ,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, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. The black dash-dotted line at E/Ep = 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 bottom-rows 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 δ²E_div.

Fig. 4 Power spectra of the two-dimensional fields Bx (continuous line) and Bz (dashed line), summed over all wave numbers ky, for the three cases BTD (orange), BEx2PP (red), and BEx2 (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 kx 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 dash-dotted 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_mixor 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 single-number 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 current-carrying 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 flux-conserving 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 z-components 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 force-free 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 non-preprocessed 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 force-free 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 force-free compatibility on the bottom boundary, the preprocessing provides the extrapolation code with a boundary condition that is more compatible with the force-free equations. Since extrapolation codes attempt to construct a solution of the force-free equations that is simultaneously force- and divergence-free, the more compatible the boundary, the more consistent (i.e., force- and divergence-free) the obtained solution. Conversely, when the boundary condition is incompatible with the force-free 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 force-free 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 force-free-compatible 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 current-carrying 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 non-trivial task to identify a quantitative estimation of solenoidal errors that can be applied to different discretizations of magnetic fields, essentially due to the non-local 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 well-defined 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 null-point reconnection B_MHD) presented negligible violations, and one (a force-free 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 current-carrying 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 force-freeness 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 force-free 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.

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 second-order, central-difference operator, and we employ the relevant one-sided (i.e., forward or backward), second-order 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 Δ³ 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 Δ² 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 finite-volume 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₂),A_y(x,y,z = z₂),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 z-component of B is changed throughout the volume except on the top boundary, whereas the x- and y-components 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 z-component of the field is changed, the divergence cleaner changes the x- and y-components of the current, but not the z-component, (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 ∂_zG(z) = −f(z) numerically, i.e., must satisfy the numerical formulation of the fundamental theorem of integral calculus in the employed discretization. For the second-order central differences that are used in the analysis, this can be obtained by the recurrence formulae (B.6)where G(z) = G(z₁ + kΔ) ≡ G(k) with k = 0,1,2,···   ,(n_z − 1), and Δ is the uniform spatial resolution in z.

The constraint ∂_zG(z) = −f(z) in the second-order, central-difference 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 second-order divergence operator is computed by using a second-order, forward derivative in the z-direction, 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 z-direction, 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.

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/2007-2013) 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 quadri-core bi-Xeon 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 post-launch operation is provided by JAXA and NAOJ (Japan), STFC (UK), NASA (USA), ESA, and NSC (Norway).

References

All Tables

All Figures

Fig. 1 Selected field lines of the six test cases: a) the potential field of a double dipole, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. 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, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. 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 BJ, normalized to the energy of Btest (red lines, equal to EJ,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, BDD; b) the TD model, BTD; c) the MHD model, BMHD; d) the NLFFF model of the nonpreprocessed magnetogram of AR 11158, BEx1; e) the NLFFF model of the preprocessed magnetogram of AR 11024, BEx2PP; f) the NLFFF model of the nonpreprocessed magnetogram of AR 11024, BEx2. The black dash-dotted line at E/Ep = 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 two-dimensional fields Bx (continuous line) and Bz (dashed line), summed over all wave numbers ky, for the three cases BTD (orange), BEx2PP (red), and BEx2 (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 kx has the dimension of pixel-1 and is normalized to the total number of modes).

In the text