RECONSTRUCTION OF NONLINEAR FORCE-FREE FIELDS AND SOLAR FLARE PREDICTION

A brief review is presented of methods for calculating nonlinear force-free fields, with emphasis on a new, fast current-field iteration procedure. The motivation is to reconstruct coronal magnetic fields using high resolution vector magnetic field boundary data from a new generation of spectro-polarimetric instruments. Methods of solar flare prediction are also reviewed, with focus on the need to reproduce observed solar flare statistics. The event statistics method is described, as well as an extension of the method to incorporate additional information, based on Bayesian predictive discrimination.


Introduction
Coronal magnetic fields around sunspots provide the energy for large scale solar activity, in particular solar flares and coronal mass ejections, so there is considerable interest in accurate descriptions of that field.Space weather effects associated with flares and CMEs motivate solar flare prediction.Flare prediction is currently in a nascent state, with existing prediction methods being imprecise and of limited practical use.Improved knowledge of the coronal magnetic field holds the promise of improved flare prediction.
It is difficult to infer the magnitude and direction of the coronal magnetic field.Measurements of the polarisation state of certain magnetically sensitive spectral lines permit inference of the vector magnetic field in the low solar atmosphere (at the photosphere or chromosphere) for regions on the solar disk, subject to a number of uncertainties.2][3] Also, the magnetic field transverse to the line of sight is determined up to a 180 degree ambiguity, which is usually resolved by an ad hoc method. 4Despite these problems, photospheric and chromospheric vector magnetic field values provide the most detailed information available on the state of the magnetic field in solar active regions.
A new generation of spectro-polarimetric instruments will soon provide a wealth of photospheric vector magnetic field determinations with high spatial resolution.The National Solar Observatory's ground-based Synoptic Long Term Solar Investigations of the Sun Vector Spectromagnetograph (SOLIS/VSM), as well as the space-based Solar Dynamics Observatory Helioseismic and Magnetic Imager (SDO/HMI) and Solar-B Solar Optical Telescope (Solar-B/SOT) feature detectors with thousands of pixels on a side.In principle measurements with these instruments may be used to better understand magnetic energy storage in solar active regions, and to improve flare prediction.However, the use of the data for these purposes hinges on the ability to accurately 'reconstruct' the magnetic field in the solar corona, based on the field values at the photosphere.This problem presents a considerable challenge.
This paper reviews the approach to the problem based on the assumption that the magnetic field in the corona (and at the boundary) is force-free, i.e. has a vanishing Lorentz force.A new, fast approach to the solution of the nonlinear force-free equations 5 is described in Section 2, as well as the prospects for application of the method to solar data.In Section 3, the state of flare prediction is discussed, with emphasis on the problem of reproducing observed flare statistics.The event statistics method 6,7 of prediction is described.This is a simple approach which uses only the past history of flare occurrence to make a prediction.A general approach for incorporating additional information into this method is also outlined.

Nonlinear force-free magnetic fields
A force-free magnetic field B satisfies as well as ∇ • B = 0. Force-free fields provide a simple model for the coronal magnetic field.The justification is that in most situations the Lorentz force density J × B = µ −1 0 (∇ × B) × B is expected to dominate over other forces so that in the static situation the Lorentz force must be close to zero.
The force-free equations may be re-written in the form: where α(r) describes the current distribution.The restrictive cases in which α is zero and α is a constant describe potential and 'linear' force-free fields respectively.The general case in which α varies with position describes 'nonlinear' force-free fields, which provide a minimal model for solar coronal magnetic fields.
A key question is whether the force-free equations may be solved for boundary conditions derived from solar spectro-polarimetric measurements.This provides an approach to the problem identified in Section 1, i.e. the reconstruction of the coronal field from lower boundary values.The appropriate boundary conditions for Equations (2) are the normal component of B in the boundary, and the value of α prescribed over one polarity of B, i.e. over one sign of the normal component of B. In principle α may be estimated from differencing of observationally inferred vector magnetic field values.For example, assuming the vector field is available at the plane z = 0, then where A basic difficulty to be met is that the field at the photosphere, the height of the measurements, is not force-free.The field is believed to become force-free in the chromosphere, typically at a height of around 500 km. 8One possible solution to this problem is 'preprocessing', i.e. alteration of the boundary conditions to make them more consistent with necessary forcefree conditions. 9 variety of methods of solution of nonlinear force-free equations have been investigated, including current-field iteration, 10 magneto-frictional relaxation, 11 and the optimization method. 12A recent test of methods 13 examined their performance on a directly calculable axially symmetric nonlinear force-free field. 14The best performing method was a specific implementation of optimization. 15lthough the different methods have been demonstrated to work on test cases, 13 they are computationally intensive, and hence slow.One simple measure of the speed of a method is the time taken, as a function of the size of the problem.Since we are concerned with three-dimensional calculations we consider the time taken as a function of N , for a calculation performed on a grid with N 3 points.As reported in Ref. 13, the time taken by specific implementations of the three methods mentioned above scaled as N 5 (optimization, magnetofrictional) and N 6 (current-field iteration).
If full resolution data from the new instruments is to be used, then the methods must cope with N ≈ 1k-2k.Based on the reported scalings, calculations of this size may be unfeasible, although it should be noted that most of the instruments provide full-disk data, and data extracted for individual active regions will then be smaller in size.In any case there is considerable interest in faster implementations of nonlinear force-free methods.Recently Ref. 16 reported implementations of both optimization and current-field iteration which scale as N 4 .

A faster current-field iteration method
Various implementations of current-field iteration 10 have been devised.The general approach consists of a Picard iteration solution of Equations (2).Specifically, at iteration k the equations are solved, subject to the boundary conditions and where B obs and α obs denote the observed values of these quantities at the boundary z = 0.For simplicity we restrict attention to a problem in the half space z > 0. We also specify the boundary conditions using α obs on the positive polarity, although either polarity may be used.
The different implementations of current-field iteration differ in their methods of solution of Equations ( 5) and ( 6), as well as in their handling of the boundary conditions both at the lower boundary, and at the side and top boundaries of the computational grid.
Recently a particularly simple and fast current-field iteration procedure has been described. 5The method involves separating the field at a given , where B 0 is the potential field satisfying The non-potential component B k+1 c may be constructed via solution of where The construction (11), which is due to Ref. 17, ensures that as required by Eqs. ( 7) and ( 9).The vector potential A k+1 c corresponding to B k+1 c may be obtained by solving the vector Poisson equation via the Fourier transform method.Since the current density J k c is specified in all space, no explicit boundary conditions are required, which makes this step particularly fast -the time taken scales as N 3 log N .For comparison, two earlier implementations of currentfield iteration using an integral solution to Ampere's law and solution by finite differences both perform this step with a N 6 scaling. 13he other equation in the current field iteration procedure, Equation (6), may be solved by field line tracing.For each point r i in the computational grid, the field line threading the point is traced in both directions.If the field line leaves the grid via the sides or top, then α(r i ) is set to zero.If the field line connects to the lower boundary at both ends, then α(r i ) is assigned based on the boundary values α obs at one end of the field line.The time taken by field line tracing for all points on the grid scales as N 4 .Hence solution of Equation ( 6) is slower than solution of Equation ( 5), and determines the overall time taken by the implementation.
The starting point for the method is the potential field B 0 , i.e.B k=0 = B 0 .The potential field is obtained from the boundary condition ẑ • B obs z=0 and needs to be calculated just once.More complete details of the method are provided in Ref. 5.
The new method has been tested 5 by application to a known nonlinear force-free field, 14 and to a simple bipole.The tests have confirmed the accuracy of the method, and the N 4 scaling of the time taken by the method.Typically the method converges in about 10 iterations.Figure 1 presents an example of two calculations performed with the new method.The boundary conditions on ẑ • B| z=0 are the same in each case, and consist of two Gaussian patches of field with positive polarity and two with negative polarity, representing two nearby bipoles.The color in the background indicates the value of the normal component of the field in the boundary, with light green showing regions with positive polarity and dark green showing regions with negative polarity.The boundary conditions on α are different in the two cases.In the upper case, there are two patches of positive α centered on the positive poles of each bipole, and α on the positive polarity is otherwise zero.In the lower case, there is a positive patch of α centered on the positive pole at lower right, and a negative patch of α centered on the positive pole at upper left, and α on the positive polarity is otherwise zero.The values of α on the negative polarity are an outcome of the calculation.In each case field lines are shown as blue curves, traced from boundary points with non-zero current density.Also shown, by transparent orange surfaces, are isocontours of the magnitude of the current density.
The two cases shown in Figure 1 represent nearby bipoles with the same, and with opposite, sign of current helicity (J • B), and the figure illustrates that the connectivity of current carrying bipoles depends on the currents which flow.The figure also illustrates a breaking of symmetry.For the upper case, the boundary conditions have a 180-degree rotational symmetry, which is reflected in the symmetry of the calculated nonlinear force-free field.In the lower case the boundary conditions no longer have the 180-degree symmetry, and consequently the field is non-symmetric.
Preliminary tests of the application of the new method to solar boundary data have also been performed.When the method is applied to solar vector magnetic field data, it is found that the current-field iteration does not converge.This may be due to the large errors in the estimated boundary values of α, which can lead to large (spurious) localised values of the current density.It also may be due to the non-force-free nature of the boundary data, as mentioned in Section 2.1.When the boundary data is preprocessed, or just smoothed, the localised currents are reduced and the method is found to converge.However, in this case the boundary conditions have been altered, so it is unclear how accurate the resulting coronal field model is.A basic problem with the application of nonlinear force-free methods to solar boundary data is that there is no unambiguous measure of the accuracy of the result.The application of the new method to solar boundary data continues to be investigated.

Existing methods of flare prediction
A variety of properties of solar active regions are known to correlate with flare occurrence.For example, active regions with certain categories of magnetic complexity 18 and certain sunspot classes 19 are more likely to produce flares, in particular large flares.Other properties associated with flaring include the length of the sheared neutral line, 20 flux emergence, 21 moments of quantities derived from vector magnetic field maps, 22 and the power-law index of the spectrum of line-of-sight magnetic field values. 23lthough these (and other quantities) correlate with flaring, there is no certain indicator that an active region will flare, and existing methods of flare prediction are probabilistic.The US National Oceanic and Atmospheric Administration uses an 'expert system' to issue flare predictions.This system incorporates sunspot classification, as well as a variety of observations and rules of thumb.Predictions are made for the probability of occurrence of at least one flare within a day with a peak 1-8 Å flux between 10 −5 W m −2 and 10 −4 W m −2 as measured by the Geostationary Observational Enviromental Satellite (an M class flare), and for the probability of occurrence of at least one flare with a within a day with a peak 1-8 Å flux greater than 10 −4 W m −2 (an X class flare).These probabilities may be labelled ǫ M−X and ǫ X respectively.
Observed solar flare statistics provide an important constraint on flare prediction.It is well known that flares follow a power-law peak flux distribution. 24In other words the number N (S) of events per unit time and per peak flux S is described by N (S) = const × S −γ , or equivalently by where dS is the total rate of events above size S i .The power-law index γ is typically found to be slightly less than two. 25ssuming Poisson occurrence in the prediction interval T = 1 day, the probabilities ǫ M−X and ǫ X may be related to corresponding rates λ M−X and λ X : Using Equation ( 13), the rates λ M−X and λ X are given by and where S M and S X are the peak fluxes associated with M and X events respectively.From Equations ( 14)-( 16) the quantity is equal to which is independent of λ i and hence constant in time.This provides a simple check of the compatibility of predictions with the power-law size distribution.
Figure 2 plots the quantity R = ln(1 − ǫ M−X )/ln(1 − ǫ X ) for the predictions made by NOAA for the year 2005.The power-law index for the GOES peak fluxes was estimated for all events in 2005 above size M using a maximum likelihood method, 26 and found to be γ = 1.88 ± 0.08.The corresponding value of R predicted by Equation ( 18) is shown in Figure 2 by the horizontal line.The NOAA predictions are rounded, i.e. are restricted to the values 0.01, 0.05, 0.10, 0.15, etc., and some of observed variation is due to this rounding.However, there is more variation than expected on this basis, and hence the NOAA predictions are inconsistent with power-law statistics.

Event statistics method
The event statistics method of flare prediction 6,7 is a particularly simple approach which is consistent with observed flare statistics.A prediction is made just on the basis of events already observed.
The method requires an estimate λ 1 of the rate of small events (events above size S 1 ) and of γ, and then predictions are made for the occurrence of big events, i.e. events above size S 2 .According to Equation (13) the rate of big events is given, in terms of the rate of small events, by This estimate may be made even if no big events have been observed.According to Poisson statistics, the probability of at least one big event in a time T is then Equations ( 19) and ( 20) provide the required prediction.The advantage of the method is that, if many small events are observed, the prediction (20) will be accurate.In particular, if M events are involved in the estimation of the rate λ 1 , then it is easy to show that the fractional error in the prediction is In Ref. 6 a Bayesian version of the event statistics method is developed.Given a sequence of events with sizes s 1 , s 2 , ..., s M above size S 1 which occur at times t 1 , t 2 , ..., t M , the Bayesian version permits calculation of a posterior distribution P (ǫ) for the quantity ǫ.If the rate of small events is determined based on a recent interval T ′ during which M ′ events occurred with a constant mean rate, then the posterior distribution is given by 6 where Λ(λ 1 ) is the prior distribution for λ 1 , and C is the normalization constant, determined by the requirement Ref. 7 describes an implementation of the method for whole-Sun prediction of GOES events, and a test of the method on the GOES record for 1976-2003.For each day in this period predictions ǫ M−X and ǫ X were made based on preceding events, and the results were compared with the historical record of whether M and X events did or did not occur.Comparison was also made with the corresponding NOAA predictions for 1987-2003.The event statistics method was found to out-perform the NOAA method for prediction of X class flares.However, on average the predictions made by the method (as well as by the NOAA method) were slightly too large, and the predictions were also conservative: for all days in 1976-2003, ǫ X was less than 0.5.

Incorporating additional information
The event statistics method uses only the past history of occurrence of small flares.It neglects all of the flare indicators discussed above, including sunspot classification and magnetic complexity.It should be possible to improve the method by incorporating additional information.
A systematic method for classifying whether an active region is flareproducing or not has recently been presented. 22The method uses discriminant analysis, an orthodox statistical technique for classification.In Ref. 22  discriminant analysis was applied to moments of quantities derived from vector magnetic field maps to classify active regions as flare-producing or non-flare producing.
The Bayesian version of discriminant analysis is predictive discrimination. 27Predictive discrimination assigns a probability to membership of a class based on observed properties, and a training sample of class members with corresponding properties.In common with discriminant analysis, it assumes continuous distributions for the properties, so it is restricted to observed properties which vary continuously (which excludes categorical properties such as sunspot classification).Predictive discrimination is more accurate than discriminant analysis, in particular for small training samples, because it takes into account variability in the training sample. 28redictive discrimination offers a Bayesian approach to incorporating additional information into the event statistics method of flare prediction.In brief this may be achieved as follows.We consider the relevant classes to be 'flare producing above size S 1 within time T ' and 'not flare producing above size S 1 within time T '.These classes may be denoted i = f 1 and i = f 1 , in an obvious notation.Following the event statistics approach, the observed rate λ 1 of events above size S 1 provides an initial guess ǫ 1 for membership of the class f 1 .The initial guess, which we denote P (i = f 1 ), plays the role of a a prior probability.We assume that the observed data x used to make the classification is a d-dimensional vector, with each element taking a continuous value.A training sample {x i,jk }, where j = 1, 2, ..., N i labels the number of vector data points and k = 1, 2, ..., d labels the components of each vector, is assumed to be available.Each j refers to observations for a particular active region.The components enumerated by k are then relevant properties of the active region, for example moments of quantities derived from vector magnetic field maps, following Ref.22. Assuming the data is multinormally distributed, predictive discrimination assigns the probability for membership of the class i as 27 where x i and C i denote the unbiased estimators of the sample means and covariance matrices: and For the general case of unequal sample means and covariance matrices for different classes the term P i (x|x i , C i ) on the right hand side of Equation (23) may be expressed as 27 where is the squared Mahalanobis distance, a measure of the distance of the point x from the mean of the sample.Equations ( 23)-( 27) provide an estimate for the probability of at least one flare above size S 1 during time T , based on both observed event statistics and the available additional information.Equations ( 19) and ( 20) may then be used to convert this into a prediction for the probability of at least one event above size S 2 during the time T .
Although the method outlined above gives a general prescription for incorporating additional information into the event statistics method, to date the approach has not been implemented and tested on data.It remains to be seen whether improved predictions result.

Summary
This paper presents a brief review of the problem of reconstructing coronal magnetic fields from vector magnetic field boundary data based on the nonlinear force-free model, and of the problem of flare prediction.
A new fast method for calculating nonlinear force-free fields is reviewed. 5he time taken by the method scales as N 4 , for calculations on grids with N 3 points.The motivation for developing fast methods of calculation of nonlinear force-free fields is the reconstruction of coronal magnetic fields based on high resolution vector magnetic field boundary data, such as that shortly to be available from new spectro-polarimetric instruments.At present there are unresolved difficulties in applying nonlinear force-free methods to solar data, and these are briefly described.
The event statistics method of flare prediction is also reviewed. 6,7This simple procedure uses only the past time history of flaring to make predictions, which have the advantage of being consistent with observed solar flare statistics.A generalization of the method to include additional information, based on Bayesian predictive discrimination, is also described.The generalization has not appeared in the literature before.The prescription is general, and should permit the incorporation of a variety of additional observational information.It remains to be implemented and tested on solar data.

Fig. 1 .
Fig.1.Examples of nonlinear force-free fields calculated with a new, fast current-field iteration procedure.This case involves two current-carrying bipoles.In the upper panel, the boundary conditions on the current are that the two bipoles have the same sign of current helicity.In the lower panel, the bipoles have the opposite sign of current helicity.