How to deal with measurement errors and lacking data in nonlinear forcefree coronal magnetic field modelling?
(Research Note)
T. Wiegelmann  B. Inhester
MaxPlanckInstitut für Sonnensystemforschung, MaxPlanckStrasse 2, 37191 KatlenburgLindau, Germany
Received 9 March 2010 / Accepted 4 May 2010
Abstract
Context. The measured solar photospheric magnetic field
vector is extrapolated into the solar corona under the assumption of a
forcefree plasma. In the generic case this problem is nonlinear.
Aims. We aim to improve an algorithm for computing the nonlinear
forcefree coronal magnetic field. We are in particular interested to
incorporate measurement errors and to handle lacking data in the
boundary conditions.
Methods. We solve the nonlinear forcefree field equations by
minimizing a functional. Within this work we extend the functional by
an additional term, which allows us to incorporate measurement errors
and treat regions with lacking observational data. We test the new code
with the help of a well known semianalytic test case. We compare
coronal magnetic field extrapolations from ideal boundary conditions
and boundary conditions where the transversal magnetic field
information is lacking or has a poor signaltonoise ratio in weak
field regions.
Results. For ideal boundary conditions the new code gives the
same result as the old code. The advantage of the new approach, which
includes an error matrix, is visible only for nonideal boundary
conditions. The forcefree and solenoidal conditions are fulfilled
significantly better and the solutions agrees somewhat better with the
exact solution. The new approach also relaxes the boundary and allows a
deviation from the boundary data in poor signaltonoise ratio areas.
Conclusions. The incorporation of measurement errors in the
updated extrapolation code significantly improves the quality of
nonlinear forcefree field extrapolation from imperfect boundary
conditions.
Key words: Sun: corona  Sun: photosphere  magnetic fields
1 Introduction
Because routine measurements of the coronal magnetic field are not available, we
have to rely on measurements of the photospheric magnetic field vector
to estimate the coronal magnetic field distribution. The
extrapolation of this surface field into the corona constitutes a
boundary value problem, which can be solved numerically if some
simplifying model assumptions for the coronal field are made.
Because of the low plasma
in the solar corona
nonmagnetic forces are often neglected to lowest order.
We have to solve the corresponding nonlinear forcefree boundary value problem
where is the magnetic field. Several methods have been developed to solve these equations, (see, e.g., Schrijver et al. 2006; Metcalf et al. 2008, for an overview and a comparison and evaluation of codes.) The comparison of different nonlinear forcefree extrapolation codes revealed that the algorithms can reliably reconstruct model fields from consistent boundary data.
However, in a recent joined study by DeRosa et al. (2009) which deals with an observed dataset (vector magnetogram taken with Hinode/SOT embedded in a lineofsight magnetogram from SOHO/MDI) the forcefree codes did not find consistent solutions. A major problem was that only for a part of the model region vector magnetograms were observed (in the Hinode fieldofview, FOV) and the transverse magnetic field component was unknown in the remaining photospheric area. In the study, the lacking field components were replaced by zeros, which was a simple way to treat the lacking data. In view of the inconsistency of the results DeRosa et al. (2009) concluded that a successful nonlinear forcefree reconstruction requires
 1.
 A large computational domain with a high resolution, which accommodates most of the connectivity within the coronal region under study.
 2.
 Taking account of measurement uncertainties, in particular for the transverse field component.
 3.
 ``Preprocessing'' of the observed vector field that approximates the physics of the photospheretochromosphere interface as it transforms the observed, forced, photospheric field to a more realistic approximation of the near forcefree field in the upper chromosphere,
2 Method
To solve the forcefree equations, we extended the optimisation approach
introduced by Wheatland et al. (2000):
The first integral contains the conditions (1) and (2) as quadratic forms and obviously for that the forcefree conditions (1, 2) are fulfilled when the functional reaches its minimum at L=0. This was the approach used in before (see Wiegelmann 2004; Wheatland et al. 2000, for details) where the boundary conditions entered directly as iterative improvements of to minimize L were constrained to at the photospheric boundary.
Here, we propose to extend the functional by adding the surface integral over the photosphere (4) and iterating to minimize L, which is otherwise unconstrained. In this integral, is a diagonal error matrix, the elements of which are chosen inversely proportional to the local measurement error of the respective photospheric field component at x,y. In principle should be specified by the instrument, and it is likely that this error distribution will be provided along with the data for SDO/HMI vector magnetograms. Because the lineofsight photospheric magnetic field is measured with much higher accuracy than the transverse field , we typically set the component to unity, while the transverse components of are typically small but positive. In regions where has not been measured or where the signaltonoise ratio is very poor, we set .
The boundary value problem (1) and (2) is nonlinear, and there is no guarantee that a solution exists for all sets of boundary values in the sense that we can obtain L=0exactly. Boulmezaoud & Amari (2000) have shown that the full vector field as boundary condition overdetermines the problem, and Wiegelmann et al. (2010) investigated the influence of photospheric measurement errors. Observed and noisy boundary fields are very probably inconsistent and belong to the set of boundary values for which a strict solution does not exist. In these cases the first integral in (4) alone could often not be iterated to values of the order of the numerical roundofferror. With the new formulation (4) we allow deviations between the model field and the observed and so the model field can be iterated closer to a forcefree field even if the observations are inconsistent. This balance is controlled by the Lagrangian multiplier .
2.1 Encoding of old code
The previous version of our optimisation code worked as follows (see also Wiegelmann 2004, for details):
 Set initial to the potential field computed from the normal component of at the photospheric boundary.
 Replace the transverse field component of the initial by the observed field components of at the photospheric boundary.
 Minimise L (Eq. (4)) iteratively, keeping unchanged at the photospheric boundary. Only the first two terms are influenced and the surface integral vanishes by construction at each iteration step.
2.2 Encoding of new code
 Set initial to the potential field computed from the normal component of at the photospheric boundary as above.
 Minimise L (Eq. (4)) iteratively without constraining at the photospheric boundary. The transverse magnetic field is gradually driven towards the observations while the field relaxes to a forcefree field. If the observed field is inconsistent, the difference remains finite depending on the control parameter . Where the observed field is automatically ignored.
 Different from the old code the magnetic field is also relaxed in the bottom boundary. Consequently the boundary values of in regions with poor signal to noise ratio ( ) are automatically relaxed towards forcefree consistent values.
 The state L=0 corresponds to a perfect forcefree and divergencefree state and exact agreement of the boundary values with observations in regions with . For inconsistent boundary data the forcefree and solenoidal conditions can still be fulfiled, but the third surface term will remain finite. This results in some deviation of the bottom boundary data from the observations in regions, where is small.
Table 1: Evaluation of the reconstruction quality.
Figure 1: Magnetic fieldline plots for the reference equilibrium (panel a)), the initial potential field b) and the other panels show nonlinear forcefree reconstructions with the old and new code for different cases, see text. 

Open with DEXTER 
Figure 2: Evolution of the entire functional L (solid line, see Eq. (4)) and its three terms (see Eq. (5)(7)). 

Open with DEXTER 
3 Code testing
We tested the new code with the well known semianalytic equilibrium from Low & Lou (1990). These are forcefree axissymmetric equilibria in spherical geometry, but the symmetry has broken for our tests by cutting a rectangular obliquely orientated box out of the spherical solution. By means of a linear shift l and a rotation angle , we positioned the centre of the equilibrium solution and tilted its axis orientation relative to the cartesian geometry of the computational domain (see Low & Lou 1990, for details). These equilibria have become a standard for testing nonlinear forcefree extrapolation codes, (see, e.g., Tadesse et al. 2009; Wiegelmann & Neukirch 2003; Amari et al. 2006; Inhester & Wiegelmann 2006; Schrijver et al. 2006; Valori et al. 2007; Wheatland et al. 2000, for earlier studies with this equilibrium). We here chose l=0.3 and for our tests and computed the equilibrium in a box with nx=ny=80, nz=72.
Table 1 shows the results of these tests, where Col. 1 names the used code. The second and third column list the values of the new parameters for the code in (4), the Lagrangian multiplier and the matrix element while was set to unity throughout^{}.
These parameters were not present in the old code.
We treated the equilibrium field on the bottom boundary as the observed field input for both our old and new extrapolation code. After the extrapolation, we evaluated the quality of our reconstruction by comparing the extrapolation result with the exact reference solution in the center 64^{3} cube by a number of figures of merit in Table 1 Cols. 48. These are the vector correlation (VC), Chauchy Schwarz (CS), normalized vector error (NE), mean vector error (ME) and the ratio of the computed magnetic energy and the energy of the Low and Lou reference field (Energy). For a perfect reconstruction NE and ME are zero and all other quantities are unity (see Schrijver et al. 2006, for a detailed explanation of these figures and the corresponding mathematical equations).
We also evaluated how well the forcefree condition
the divergencefree condition
and the boundary data term
are satisfied. The respective numbers are listed in Table 1 in Cols. 9 and 10. The last column displays the number of iteration steps needed to obtain the reconstruction. In practical computations L_{1}, L_{2} and L_{3} never vanish exactly, but the code stops when becomes more or less stationary. The final values of L_{1}, L_{2} may be high if the imposed boundary conditions are inconsistent. In this case a forcefree solution that matches the boundary fields does not exist. However, even for consistent boundary conditions L_{1} and L_{2} remain finite for numerical reasons. As a typical error we expect , where is the grid spacing and N the number of grid points along each axis. For N=64 this yields . The typical discretisation error can be estimated from the entries for the analytical solution in the first row and for the potential field in the second row. They are of the same order of magnitude as the estimated error above except for L_{1} of the potential field, which diminishes to roundoff errors because it was obtained from a numerical gradient of a potential. An iterative solution of the potential field gave residual values for both, L_{1} and L_{2}slightly below the estimated numerical error of 0.1.
3.1 Tests with ideally measured boundary conditions
The first two rows in Table 1 show for reference the figures of merit for the original reference solution and the initial potential field, respectively. Rows 39 display test results with ideal boundary conditions from the exact solution. For the new code we varied the Lagrangian multiplier and the matrix elements . For all pixels in the magnetogram are treated equally, which is reasonable for perfectly measured data. Figure 1c and 1d show the respective reconstructed field lines with the old and new code , respectively. Within a reasonable range of parameters the new code gives the same result as the old one. If the Lagrangian multiplier is chosen too small or is nonzero mainly only for strong field regions, the reconstructed solution is gradually decoupled from and the resulting field deviates from the exact solution towards a potential field. For the code should stop immediately and return the initial potential field for which L_{1}=L_{2}=0. For the bottom boundary field is forced to agree with at every iteration step and the iteration should perform as for the old code.
A main difference between old and new code is the evolution of the functional L during the iteration as shown in Fig. 2 panel a) and b). For the new code the forcepart L_{1} and divergencepart L_{2} remain small during the entire iteration. The final state is almost identical for both implementations, however.
3.2 Tests with boundary conditions with lacking data
The main reason for the new implementation of our code is that we need to deal with boundary data of different noise levels and qualities or even miss some data points completely. This occurs e.g. due to the limited field of view of some vector magnetographs like HinodeSOT. The lineofsight magnetic field is usually available for the entire solar disk, e.g. from SOHO/MDI, but the transverse components are often not available in parts of an active region. SDO/HMI will observe fulldisk vectormagnetograms, but will suffer from a poor signaltonoise ratio in weak field regions. This makes measurements less reliable in these regions.
We mimiced this effect by removing the information regarding the transverse field in certain areas, e.g. where and , respectively. For our test field, this affects and of the bottom boundary, respectively. Finally as an extreme test we assigned only in a very restricted area close to a local field concentration. This simulates measurements with a vector magnetograph of a very small FOV. For our test field is then available for only of the bottom boundary. In the new code, the boundary area with lacking data is marked by , in the old code the lacking transverse field components are filled with zeros, which implies a vertical magnetic field (see DeRosa et al. 2009, for details.)
The fieldline plots from these latter two reconstructions are shown in Fig. 1eh. The comparison with the exact solution and numerical residuals of the force and the divergence are displayed in rows 1018 of Table 1. We find that the new code is closer to a force and divergencefree field in all cases and L_{1} and L_{2} are significantly smaller, in particular for low values of the Lagrangian multiplier . The result of the new code agrees also somewhat better with the original solution than the output of the old code, but the difference is not as significant as the difference in forcefree and solenoidal conditions.
The reason for these results is that the old and new code react differently to inconsistent boundary conditions caused by noisy data and by assuming a vertical field in regions where transversal field measurements are lacking. As seen in Fig. 2 panel c) the functional L iterated in the old code soon reaches a stationary state at a finite value of L if the boundary data are inconsistent. This is different with the new code, shown in panels d) and e). Here L_{1} and L_{2} remain very small during the entire evolution and the inconsistency in the boundary data is absorbed in the L_{3} term. L_{3} can also become relatively low however, because regions with (where the transverse field has a poor signaltonoise ratio or is even unknown) do scarcely or not at all contribute to the functional. For a reasonably small Lagrangian multiplier as shown in panel e) the functionals L_{1} and L_{2} stay at a very low level (discretisation error). The new code also relaxes the bottom boundary values and changes the field there in order to find forcefree consistent boundary conditions. If information on are lacking in a significant part of the magnetogram (up to of the pixels in our example), there is no unique solution for consistently specifying , and one cannot expect the final equilibrium to agree with the model solution. Yet we have the most probable magnetic field model, because it is close to force and divergencefree and agrees with the sparse observed boundary data.
3.3 Influence of noise
Finally we investigate how noise influences the quality of the reconstructed magnetic field in the last three rows of Table 1. We added a uniform noise of of the maxim transversal magnetic field onto ^{} For the old code the inconsistency in the boundary data leads to significant deviations from the forcefree L_{1} and solenoidal condition L_{2}. Computations with the new code lead to significantly better agreement of the force and solenoidal condition, e.g., by almost two orders of magnitude for a small Lagrangian multiplier . The effect of noise in the transversal field is similar to lacking data, which is natural, as regions with a very poor signaltonoise ratio in can be as well considered as regions where we do not know . The new code still finds a forcefree solution for these cases, which are, however, more potential field like than the solution with full information. But even for this relatively high noise level the reconstructed solution still has similarities to the original Low and Lou reference field (error of in the vector correlation and an underestimation of the energy of .) and is a significantly better approximation than a potential field, which as an error of and for the vector correlation and energy estimate, respectively.
4 Conclusions
Different from previous implementations our new code allows us to deal with lacking data and regions with poor signaltonoise ratio in the extrapolation in a systematic manner because it produces a field which is closer to to a force and divergencefree field and tries to match the boundary only where it has been reliably measured. For ideal and consistent data the extrapolation result is identical with previous implementations of our code. This old result could also artificially be enforced by the choice and .
For finite , the new code also relaxes the boundary and allows us to fulfill the solenoidal and forcefree condition significantly better because it allows deviations between the extrapolated boundary field and inconsistent boundary data. This deviation can be controlled by setting the weight factors inversely proportional to the measurement error or , where the field information is lacking.
If no transverse magnetic field has been observed on a significant part of the vector magnetogram, the magnetic energy of the final relaxed state is underestimated. In an extreme case when no at all is available in the entire region, the code would produce the initial potential field.
AcknowledgementsThe work of T. Wiegelmann was supported by DLRgrant 50 OC 0501. This work was inspired by discussions with our colleagues in the NLFFFconsortium (a group of scientist sharing their experience on nonlinear forcefree field modelling) during six annual workshops in the years 20042009.
References
 Amari, T., Boulmezaoud, T. Z., & Aly, J. J. 2006, A&A, 446, 691 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Boulmezaoud, T. Z., & Amari, T. 2000, Z. Angew. Math. Phys., 51, 942 [CrossRef] [Google Scholar]
 DeRosa, M. L., Schrijver, C. J., Barnes, G., et al. 2009, ApJ, 696, 1780 [NASA ADS] [CrossRef] [Google Scholar]
 Fuhrmann, M., Seehafer, N., & Valori, G. 2007, A&A, 476, 349 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inhester, B., & Wiegelmann, T. 2006, Sol. Phys., 235, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Low, B. C., & Lou, Y. Q. 1990, ApJ, 352, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Metcalf, T. R., Derosa, M. L., Schrijver, C. J., et al. 2008, Sol. Phys., 247, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J., Derosa, M. L., Metcalf, T. R., et al. 2006, Sol. Phys., 235, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Tadesse, T., Wiegelmann, T., & Inhester, B. 2009, A&A, 508, 421 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valori, G., Kliem, B., & Fuhrmann, M. 2007, Sol. Phys., 245, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Wheatland, M. S., Sturrock, P. A., & Roumeliotis, G. 2000, ApJ, 540, 1150 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T. 2004, Sol. Phys., 219, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., & Neukirch, T. 2003, Nonlin. Processes Geophys., 10, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., Inhester, B., & Sakurai, T. 2006, Sol. Phys., 233, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., Thalmann, J. K., Schrijver, C. J., Derosa, M. L., & Metcalf, T. R. 2008, Sol. Phys., 247, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Wiegelmann, T., Yelles Chaouche, L., Solanki, S. K., & Lagg, A. 2010, A&A, 511, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Footnotes
 ... throughout^{}
 We investigated several forms of as indicated in the third column of the table. 1,0 means that we chose in the valid area and 0 elsewhere.
 ... ^{}
 This would correspond to a noise level of about 100 G for real measurements.
All Tables
Table 1: Evaluation of the reconstruction quality.
All Figures
Figure 1: Magnetic fieldline plots for the reference equilibrium (panel a)), the initial potential field b) and the other panels show nonlinear forcefree reconstructions with the old and new code for different cases, see text. 

Open with DEXTER  
In the text 
Figure 2: Evolution of the entire functional L (solid line, see Eq. (4)) and its three terms (see Eq. (5)(7)). 

Open with DEXTER  
In the text 
Copyright ESO 2010