Identification of multiple hard X-ray sources in solar flares: A Bayesian analysis of the February 20 2002 event
Federica Sciacchitano Alberto Sorrentino A Gordon Emslie Anna Maria Massone Michele Piana
astro-ph.SR astro-ph.IM stat.CO 62F15, 65R32, 68U10




Michele Piana
piana@dima.unige.it



Dipartimento di Matematica, Università di Genova, via Dodecaneso 35 16146 Genova, Italy


Dipartimento di Matematica, Università di Genova, via Dodecaneso 35 16146 Genova, Italy and CNR - SPIN Genova, via Dodecaneso 33 16146 Genova, Italy


Department of Physics and Astronomy, Western Kentucky University, Bowling Green, KY 42101, USA


CNR - SPIN Genova, via Dodecaneso 33 16146 Genova, Italy


Dipartimento di Matematica Università di Genova via Dodecaneso 35 16146 Genova Italy and CNR - SPIN Genova, via Dodecaneso 33 16146 Genova, Italy





The hard X-ray emission in a solar flare is typically characterized by a number of discrete sources, each with its own spectral, temporal, and spatial variability.  Establishing the relationship amongst these sources is critical to determine the role of each in the energy release and transport processes that occur within the flare.  In this paper we present a novel method to identify and characterize each source of hard X-ray emission.  The method permits a quantitative determination of the most likely number of subsources present, and of the relative probabilities that the hard X-ray emission in a given subregion of the flare is represented by a complicated multiple source structure or by a simpler single source. We apply the method to a well-studied flare on 2002 February 20 in order to assess competing claims as to the number of chromospheric footpoint sources present, and hence to the complexity of the underlying magnetic geometry/toplogy.  Contrary to previous claims of the need for multiple sources to account for the chromospheric hard X-ray emission at different locations and times, we find that a simple two-footpoint-plus-coronal-source model is the most probable explanation for the data.  We also find that one of the footpoint sources moves quite rapidly throughout the event, a factor that presumably complicated previous analyses.  The inferred velocity of the footpoint corresponds to a very high induced electric field, compatible with those in thin reconnecting current sheets.







§ INTRODUCTION
 

The RHESSI instrument <cit.> was designed to investigate particle acceleration and energy release in solar flares through imaging spectroscopy of hard X-ray and gamma-ray emissions from ≃ 3 keV to ≃ 20 MeV.  Over the course of its 15-plus-year lifetime, RHESSI has provided thousands of hard X-ray images of solar flares. Modeling of the often complex structures thus revealed <cit.>, and relating them to theoretical models of the spatial distributions of hard X-rays in solar flares <cit.> typically involves the identification of subregions of the image and the classification of each subregion into a particular “type” of source, such as extended coronal sources <cit.>, “thick-target” nonthermal sources confined to the corona <cit.>  and chromospheric footpoints <cit.>.

A particularly interesting example of the importance of such sub-source identification is the C7.5 flare observed by RHESSI on 2002 February 20, at ∼11:06 UT. This was one of the first solar flare events observed by RHESSI and its characteristics, particularly its morphology, have been extensively studied, and we now summarize the results of these studies.




  * Figure 3 and 4 of <cit.> show images of the event produced using the Clean algorithm <cit.> and the MEM-Sato maximum entropy method <cit.>.  Both image reconstruction methods reveal the presence of two chromospheric footpoints (which become more prominent as the observing energy increases), in addition to a remote looptop source <cit.>.


  * Figure 2 of <cit.> shows further evidence for a two-footpoint structure, using five different imaging methods, including forward-fit <cit.>, Clean <cit.>, Pixon <cit.>, and two different maximum-entropy methods, one operating in (x,y) space and the other in visibility (i.e., Fourier) (u,v) space.


  * Figure 2 of <cit.> again shows evidence for two chromospheric footpoints, a morphology also supported by SOHO EIT <cit.> observations.


  * Figure 4 of <cit.> shows the source centroids at different hard X-ray energies all concentrated in these two footpoints.


  * Adding to the complexity, Figure 1 of <cit.> shows “two distinct footpoints and a faint loop-top source” located some 20 arcseconds toward the limb compared to the footpoints.


  * In their analysis of this event, <cit.> discuss (cf. their Figure 3) “a complicated pattern of low- and high-energy sources.”


  * <cit.> discuss the possibility of a more complicated morphology with three footpoint sources, plus the looptop source.  Their Figure 3 shows that one of the footpoints (“Source 2”) disappears between 11:06:11 UT and 11:06:17 UT, and is promptly replaced by a “Source 3” located some 10 arcseconds to the South <cit.>. Both Source 2 and Source 3 are close to a location involving the disappearance of a small sunspot “located near the edge of the dominant hard X-ray kernel” <cit.>, providing evidence for energy release to due to magnetic reconnection/reconfiguration and possible evidence for an acceleration site in the chromosphere <cit.>.


  * <cit.> show images of the event based both on observed counts (their Figure 1); they also performed a spectral inversion of the observed count visibilities to obtain “electron visibilities” and so produced images of the source in the electron domain (their Figure 5).  Because of the regularization in the spectral inversion step, the electron images vary much more smoothly with energy than the original count-based images.  It is therefore especially interesting that while the count-based images in <cit.> do suggest the possibility of a third footpoint source (see especially the 34-38 keV image in their Figure 1), there is little evidence for this additional footpoint in the electron images.



In light of these various results, assessing the quantitative likelihood (and significance) of a third footpoint source is an intriguing problem.  Further, <cit.> have shown that one of the footpoints in another event on 2003 July 23 moves across the solar disk with a speed that correlates well with the hard X-ray intensity of that footpoint source, suggesting that the rate of electron acceleration is correlated with the speed of the footpoint motion.  Indeed, one can easily hypothesize that higher footpoint speeds result in higher induced electric fields 𝐄 = (𝐯/c) ×𝐁 and hence greater acceleration rates.
There is therefore a clear need for a robust algorithm which can, amongst other things, distinguish between the alternative hypotheses of (a) the disappearance of one source and appearance of another, and (b) the drift of a single footpoint source. The development of just such an algorithm is the focus of the present work.

RHESSI uses indirect, rotating-collimator-based imaging techniques <cit.>. The “native” output of the imaging data is in the form of count modulation profiles from each of RHESSI's nine detectors, which can be temporally stacked to produce a set of “visibilities” (two-dimensional spatial Fourier components of the image), with each detector sampling Fourier components tracing a circle in (u,v) space over the course of the ≃4 s spacecraft rotation period.  Reconstruction of images then proceeds in one of several ways:



  * Fourier transformation (“backprojection”) of the data to produce a “dirty” image in (x,y) space, which is then improved using methods like Maximum Entropy <cit.>, Clean <cit.>, or Pixon <cit.>;

  * Inversion of count modulation profiles by means of optimally stopped iterative schemes <cit.>;

  * Interpolating the (relatively sparse) visibility data to obtain a smooth surface in (u,v) space, which can then be Fourier-transformed to produce an image in the (x,y) domain <cit.>;

  * Reconstruction methods from visibilities using compressed sensing approaches <cit.>

  * Forward-fitting the visibilities with a parametric model in order to obtain the best-fit parameters <cit.>.



The first four items essentially refer to image reconstruction methods that involve the reconstruction from data of an entire image, i.e., a set of pixel values that correspond to the intensity of the hard X-ray emission, at the photon energy of observation, in an area element on the solar disk.  By contrast, the last method represents an alternative approach, i.e. it is a parameter estimation method, where it is assumed that the data (again, count modulation profiles or visibilities) have been generated by a very small number of simply-shaped geometrical objects. The parameters of these objects are then estimated through non-linear optimization techniques. To the best of our knowledge, to date this visibility forward fitting (VFF) algorithm <cit.> is the only algorithm that provides statistical uncertainties of fitted source parameters; however, the VFF method is limited by the need for a priori knowledge of the number and the shapes of sources.  Furthermore, in the current default version of VFF provided in the SSW software, only a single circular or elliptical Gaussian, or single loop, or two circular Gaussians can be reconstructed, although more recently the VVF algorithm has been modified <cit.> such that it can address two elliptical sources simultaneously.

We here present a Bayesian subsource identification and characterization method that is ideally suited to addressing the issue of determining the number, and structure, of various sub-regions in a hard X-ray image.  The method quantitatively assesses the likely number of discrete sources in an image and the parameters that characterize each source, together with their quantitative uncertainties. This approach works with an arbitrary number of sources, the shapes of which are automatically determined by the method, together with the source parameters and their uncertainties. Importantly, the method also provides a quantitative way of determining the relative probability of configurations with different numbers of sources. It is therefore particularly useful in assessing the complexity (or lack thereof) in the source image, such as the number of footpoint sources in the 2002 February 20 event.

In Section <ref> we describe the method.  In Section <ref> we describe the RHESSI observations used for the analysis, i.e. the relatively well-studied 2002 February 20 11:06 UT event, and we compare the obtained results with the ones described in <cit.>. Section <ref> discusses the impact of the prior distribution on the results.  Our conclusions are offered in Section <ref>.

 
The RHESSI instrument <cit.> was designed to investigate particle acceleration and energy release in solar flares through imaging spectroscopy of hard X-ray and gamma-ray emissions from ≃ 3 keV to ≃ 20 MeV.  Over the course of its 15-plus-year lifetime, RHESSI has provided thousands of hard X-ray images of solar flares.  Modeling of the often complex structures thus revealed <cit.>, and relating them to theoretical models of the spatial distributions of hard X-rays in solar flares <cit.> typically involves the identification of subregions of the image and the classification of each subregion into a particular “type” of source, such as extended coronal sources <cit.>  and chromospheric footpoints <cit.>.

RHESSI uses indirect, rotating-collimator-based imaging techniques <cit.>. The “native” output of the imaging data is in the form of count modulation profiles from each of RHESSI's nine detectors, which can be temporally stacked to produce a set of “visibilities” (two-dimensional spatial Fourier components of the image), with each detector sampling Fourier components tracing a circle in (u,v) space over the course of the ≃4 s spacecraft rotation period.  Reconstruction of images then proceeds in one of several ways:



  * Fourier transformation (“backprojection”) of the data to produce a “dirty” image in (x,y) space, which is then improved using methods like Maximum Entropy <cit.>, Clean <cit.>, or Pixon <cit.>;

  * Inversion of count modulation profiles by means of optimally stopped iterative schemes <cit.>;

  * Interpolating the (relatively sparse) visibility data to obtain a smooth surface in (u,v) space, which can then be Fourier-transformed to produce an image in the (x,y) domain <cit.>;

  * Reconstruction methods from visibilities using compressed sensing approaches <cit.>

  * Forward-fitting the visibilities with a parametric model in order to obtain the best-fit parameters <cit.>.



The first four items essentially refer to image reconstruction methods that involve the reconstruction from data of an entire image, i.e. a set of pixel values, each corresponding to the intensity of the hard X-ray emission at the photon energy of observation in the corresponding area element on the solar disk.  The last method represents an alternative approach, i.e. it is a parameter estimation method, where it is assumed that the data (again, count modulation profiles or visibilities) have been generated by a very small number of simply-shaped geometrical objects. The parameters of these objects are then estimated through non-linear optimization techniques. However, the currently available version of this algorithm is limited by two constraints: first, the pre-defined source configurations that can be used for the fitting are just one circular source, one elliptical source, one loop source or two circular sources; second, the employed configuration must be selected in advance. We here present a Bayesian method that determines the optimum way of constituting a given data set from an assembly of such simple source structures. This approach overtakes the previous limitations, inasmuch it works with an arbitrary number of sources, whose type is automatically determined by the method, together with the individual parameter uncertainties. Importantly, the method also provides a quantitative way of determining the relative probability of configurations with different number of sources, and therefore is particularly useful in assessing the complexity (or lack thereof) in the source image.  In Section <ref> we describe the method.  In Section <ref> we describe the RHESSI observations used for the analysis, i.e. the relatively well-studied 2002 February 20 11:06 UT event. Section <ref> and Section <ref> discuss the obtained results and compare them with the ones described in <cit.>. Our conclusions are offered in Section <ref>.

 




§ BAYESIAN PARAMETRIC ESTIMATION OF SOLAR FLARE STRUCTURES


In this section we describe the source parametrization, the Bayesian model and the Monte Carlo algorithm that we use to estimate the flare parameters from the visibilities measured by RHESSI (this same approach can be extended to the Bayesian fitting of RHESSI count modulation profiles).



 §.§ Parametrization of flare sources


Hard X-ray images in solar flares <cit.> can involve either a single source (compact or extended), compact sources, extended sources, or combinations thereof. Very often, such maps can be effectively represented as the superposition of different two-dimensional distributions, which we assume to be each described by a simple geometric object, such as a circle, an ellipse, or a “loop”.  Even though the parametrization we shall use to describe these objects is not entirely new <cit.>, for the sake of completeness we provide here a brief description of it.

The various model sources used are characterized as follows (see Figure <ref>):




  * Circle: a circular source is defined by the x and y positions of its center, measured in arcseconds from the center of the image, the flux ϕ>0 (which represents the integrated intensity over the source), and the full width at half maximum (FWHM) intensity 2r>0 (where r measures the radius of the source).


  * Ellipse: an ellipse is defined by the location (x,y) of its center, the semi-major axis a corresponding to the 50% intensity level that extends furthest from the center, the rotation angle α∈ [0,180] of that direction, and the eccentricity


    ε= √(1- b^2/a^2)    ,


where a>0 and b>0 correspond to the length of the semi-major and semi-minor axes, respectively.


  * Loop: a “loop” is defined by the above ellipse parameters plus an additional parameter β∈ [-180,180] which indicates the degree of curvature of the loop shape (see Figure <ref>). The angle β is inversely proportional to the radius of curvature; we choose to use the reciprocal curvature parameter β to avoid the problem of dealing with infinite values of curvature for the elliptical and circular sources.



All these sources are thus defined by the following vector


    θ= (x,y, ϕ, r, α, ε, β),

where the last three parameters can possibly (for the case of a circle or an ellipse) be zero or indeterminate.  Figure <ref> provides a visual description of the parameters.



The parametrization in Equation (<ref>) has already been presented in <cit.>, where a forward-fitting method to reconstruct hard X-ray images of the solar flares was described. In their study, the authors recovered the best estimates of the various source parameters together with error bars on the estimated parameters. However, a major limitation of the <cit.> method is that the number and shape of the sources have to be assumed a priori. Furthermore, the algorithm is able to recover only the following limited configurations: at the time of writing, a single circle, a single ellipse, a single loop, or two circles. Here we overcome such limitations, by explicitly embodying the number N of sources and the different source types among the unknowns to be determined.  Specifically, in the approach to be described the goal is to probabilistically infer the following object:


    X=(N,T_1:N, θ_1:N)

where N=0, …, N_max represents the number of sources, T_1:N=(T_1, …, T_N) the source types and θ_1:N=(θ_1, …, θ_N) the source parameters for each source.



 §.§ Bayesian model


In a Bayesian setting, all the variables of the problem are considered as random variables, with information on their values encoded as their probability distributions. In the following we will for simplicity use lower-case letters to indicate both the random variables and their specific values in a given instance. We are interested in the posterior distribution for flare source parameters, as described by Equation (<ref>), conditioned on the data, namely a set of observed visibilities y.  From Bayes' theorem we have


    π_po(x|y)=ℒ(y|x)π_pr(x)/p(y),

where  π_po(x|y) is the posterior probability distribution representing the conditional probability density of x given the data y, ℒ(y|x) is the likelihood function that expresses the likelihood that the data y is indeed due to the image x, π_pr(x) is the prior probability density which describes the properties of the source configuration that we want to recover, and p(y) is a normalization factor, independent of x. In the following we will describe our choice of the prior density and the algorithm we use to then approximate the posterior density. We notice that the normalization factor p(y) is generally unknown, but its knowledge is not necessary to apply most Monte Carlo algorithms and, specifically, the Monte Carlo algorithm described below.

Since the number of sources is not known a priori, we need to consider a variable-dimension model, i.e., the state space χ of the sources is defined as follows:


    χ∪_k=1^N {k}×Θ_k     ,

where Θ_k is the state space for source number k of the N sources present.



  §.§.§ The prior distribution


As in all Bayesian problems, the choice of the prior distribution π_pr(x), which is the starting point for the calculation of subsequent improved distributions using the available data, is critical.  Our choice of suitable prior probability densities π_pr(x) reflects basic expectations that are based on available knowledge of solar flare geometries. By assuming independence between the number of sources, the source types and the source parameter values, the prior distribution π_pr can be separated as follows:


    π_pr(X)= π_pr(N, T_1:N, θ_1:N)= ρ_1(N)  ∏_k=1^N ρ_2(T_k)  ρ_3(θ_k)     ,

where the ρ_i (i ∈{ 1,2,3 }) indicate the prior distributions of the number of sources, the source types, and the source parameters, respectively. We remark that the assumption of independence of these quantities is not simply useful; rather, it corresponds to the minimization of the information content of the prior.

We now further assume that the number of sources in a map is Poisson distributed with mean λ>0:


    ρ_1(N)=e^-λλ^N/N!    ,

and we will limit λ∈ [0,3]. In addition, we shall limit N, the maximum number of sources, to N ≤ 5, because past observations <cit.> suggest that, given the spatial resolution of RHESSI, a higher number of distinct sources is unlikely.

For the source type we adopt the following prior distribution


    ρ_2(T_k)=
    
    p_C      if  T_k  is a circle
    
    p_E      if  T_k  is an ellipse
    
    p_L      if  T_k  is a loop,

where p_C,p_E,p_L ∈ [0,1].  Now, in general the morphology of solar flare hard X-ray sources is dependent on the photon energy observed: images at lower energies 15 keV generally are dominated by extended loop-like (or elliptical) structures, while images at higher photon energies are more likely to include compact sources (possibly exclusively so).  The former are typically interpreted as extended thermal sources <cit.>, while the latter are generally considered as chromospheric footpoint sources, which are compact because of the much higher density and so electron stopping distance <cit.>.  In our analysis we will use (p_C,p_E,p_L)=(1/2, 1/4, 1/4) if the energy is high and (p_C,p_E,p_L)=(1/4, 1/4, 1/2) otherwise.  We stress that these are only assumed prior distributions; they do not represent the final (data-informed) result.

The individual source parameters are a priori assumed to be uniformly distributed in their respective ranges:


    (x,y)    ∼𝒰 (  [ x_c-F/2,x_c+F/2 ],  [ y_c-F/2,y_c+F/2 ]  ) 
    
    r    ∼𝒰 (  [ 0,F/3 ]  ) 
    ϕ   ∼𝒰 (  [ 0,max Re(V)  ]  ) 
    ε   ∼𝒰 (  [ 0,1  ]  ) 
    α   ∼𝒰 (  [ 0,180 ]  ) 
    β   ∼𝒰 (  [ -180 ,180 ]  )     ,

where (x_c,y_c) is the map center, F indicates the field of view of the map in each direction and V represents the (complex) visibilities.  Combining Equations (<ref>), (<ref>), and (<ref>), the prior distribution π_pr is given by


    π_pr   = ρ_1(N) ( ∏_k=1^N ρ_2(T_k))( ∏_k=1^N ρ_3(θ_k))
       = e^-λλ^N/N! p_C^N_C p_E^N_E p_L^N_L( 1/F)^2N( 3/F)^N
               ( 1/max Re(V))^N( 1/180)^N_E+N_L( 1/360)^N_L.




  §.§.§ Likelihood function

Since visibilities are linear combinations of independent photon count numbers <cit.>, we can make the standard approximation that noise on the visibilities is Gaussian additive; the likelihood function is thus given by

    ℒ(y|x) ∝exp(- ∑_i |y-A(x)|_i^2/2 σ_i^2)

where A indicates the forward operator, i.e. the projector from the imaging parameter space onto the visibility space.  The standard deviations σ_i^2 are assumed to be known; in the data analysis below, we will be using the uncertainty values associated to the measured visibilities as provided by the SolarSoftWare (SSW) IDL routines.



 §.§ Approximation of the posterior distribution


The posterior probability density is a complicated function on a relatively high–dimensional space, so that dealing with it analytically is impractical.  We therefore use a Sequential Monte Carlo (SMC) algorithm that produces a sample set which is approximately distributed according to the posterior, and can thus be used to make inference on the values of the various parameters. A brief description of the SMC algorithm is given in the Appendix below; for more details, we refer to <cit.> for the original article in which the SMC method is introduced, and to <cit.> for an example of its application to a mathematically similar inverse problem. 

Once the algorithm converges, the parameters of the reconstructed map are computed by using the following point estimates. For the number of sources N̂_i we use the most probable value, given by the MAP (Maximum a Posteriori) estimate. Similarly, the estimated locations of the sources are given by the local modes of the posterior distribution for the source locations, conditioned on the estimated number of sources; ϕ, r, and ε are estimated by the mean values of the conditional distribution conditioned on the source location and number of parameters; the angles α, β are estimated as the modes of the corresponding conditional distributions. Finally, the types of the sources are estimated a posteriori by using the estimates of the parameters. In particular, for our experiments we use the following classification


    circle;    if ε≤ 0.1 
    ellipse;    if  0.1 ≤ε < 0.4, |β| ≤ 45
        and ε≥ 0.4   and  |β| ≤ 10
    loop;   otherwise.





§ APPLICATION TO RHESSI DATA - THE 2002 FEBRUARY 20 EVENT




 §.§ General


We now apply our analysis to the C7.5 flare observed by RHESSI on 2002 February 20, at ∼11:06 UT. This was one of the first solar flare events observed by RHESSI and, as discussed in Section <ref>, its characteristics, particularly its morphology, have been extensively studied, particularly with regard to the number, and nature, of different sources present as time progresses.

Our Bayesian subsource identification and characterization method is ideally suited to addressing the issue of the number of discrete subsources in the region. Using the method as described above, we can quantitatively assess both the number of subsources in an image and the parameters that characterize each source, together with their uncertainties. To the best of our knowledge, to date only the visibility forward fitting (VFF) algorithm <cit.> is able to provide statistical uncertainties of fitted source parameters; however, as we have already pointed out, the VFF method is limited by the need for a priori knowledge of the number and the shapes of sources. Furthermore, in the current default version of VFF provided in the SSW software, only a single circular or elliptical Gaussian, or single loop, or two circular Gaussians can be reconstructed, although more recently the VVF algorithm has been modified <cit.> such that it can address two elliptical sources simultaneously.

In this section, we will compare the values obtained through both the new Bayesian method and the VFF. As we will see, the results are very similar; however the Bayesian method has two very important advantages: it does not require any a priori knowledge of the shapes and number of sources. Furthermore, it can recover every plausible configuration with the corresponding probability, and not only the simplest one as in the VFF algorithm.

In order to be consistent with the results in <cit.> and <cit.>, in our experiments, we use RHESSI detectors 3 through 9, a field of view F=64 arcseconds; further, we choose map-centers at [918,270] and [912,263] for Sections <ref> and <ref>, which compare our results with those of the respective previous authors.  The number of particles is set to N=10000, the maximum number of sources admitted in a map is  N_max=5 and, coherently to what done in VFF, the noise level is given by adding a systematic error  σ_sys=0.05 to the error e computed by the RHESSI SSW routines, i.e.,


    σ_noise = √(e^2  + σ_sys^2  V_amp^2)    ,

where V_amp is the absolute value of the visibilities. Therefore, the only parameter that needs to be tuned is λ, and we assumed values in the range [1,3].  In most cases slightly changing λ will not influence the reconstruction, but only the convergence speed; 
for more details, we refer the reader to Section <ref>.

In order to calculate the visibilities generated by a given source, which is needed to compute the likelihood of the Monte Carlo samples, we employ the SSW function .


























 §.§ Comparison with the results in Sui et al


In this section, we focus on the time interval 11:06:10-11:06:24 UT and several energy intervals spanning the range from 6 keV to 70 keV. We compare the results found in <cit.> with those using our Bayesian approach.

In <cit.>, the Clean  <cit.> and the maximum entropy <cit.> reconstruction techniques have been used to study the number of sources with different energy bands at a fixed time interval.  In particular, <cit.> found (their Figures 3 and  4) that the map in the low-energy band (6-10 keV) shows a single structure which is located between the two footpoints observed at higher energies. In the interval 10-14 keV, a thermal source and a possible loop source are present. In the energy bands between 14 and 50 keV there are two footpoints and the loop-top source, while in the last energy band 50-70 keV only the two footpoints are present and the loop-source is not visible.



For the prior distribution used in our approach (see equation (<ref>)), we set the number of sources λ according to the results found in <cit.>, i.e., λ=1 for the energy band 6-10 keV,  λ=2 for 10-14 keV and 50-70 keV, and λ=3 for 14-50 keV.
































In Figure <ref> we show the posterior probability distribution of the number of sources in the six different energy bands. 

The distribution of the particles during the last iteration is shown in Figures <ref> and <ref>, together with the reconstructions obtained by Clean and MEM-NJIT, respectively; the values in the color bands represent the weights ω_i of each particle, with ω_i ∈ [0,1], with i ∈ [1,N]. Please notice that this plot of the Monte Carlo samples and weights only summarizes information on the posterior probability of the source locations, and not on the source size or strength; specifically, a source whose location is very certain will show up as a very small cloud of points, irrespective of the source size, while a source with an uncertain position will appear as a rather widespread cloud. Combining the information from Figures <ref> and <ref> we can draw some quantitative conclusions about the number and locations of sources in the map. For the lowest energy band (6-10 keV) the probability of having two sources is almost 90% (Figure <ref>). At first glance, this result seems to 
run counter to the results of <cit.>, where only a single source, which lies between the two footpoints, was described. However, the reconstruction provided by Clean shows a rather complex shape for that source, which cannot be well explained by a combination of simple geometric forms. The Bayesian approach addresses this problem by approximating the source with the combination of a circle and an ellipse characterized by very close (almost superimposed) locations (see Figure <ref>, top left panel and Figure <ref>). As a confirmation of this, we note that if we increase the probability of having a loop instead of a circle or an ellipse in (<ref>), the algorithm recovers a single loop-source with a significantly higher probability and with parameters and uncertainties close to the ones obtained by VFF (see Table <ref>). 
Thus, we can conclude that only a single source is found in the energy band 6-10 keV, which we interpret as a thermal coronal source situated at the apex of a magnetic structure connecting the footpoints.






In the next energy band (10-14 keV), the map (Figure <ref>(b)) appears  to have possibly three chromospheric sources  (the one in the middle having a much smaller weight than the others), but by looking at the distribution on the number of sources in Figure <ref>, we see a higher probability of  three sources (the two footpoints and the loop-top source); the probability of having four sources present (red band in Figure <ref>) is around 3%. By comparing the weights of the particles in Figure <ref>, we see that most particles are at one of two footpoints, while only a few others constitute  the putative source in the middle of the two footpoints.

The map in the 14-20 keV energy band shows a predominance of three sources (Figure <ref>): Figure <ref>(c) shows these to be two footpoint sources plus the looptop source. The position of the southern footpoint is somewhat uncertain due to its relatively low intensity.

The energy band from 20-30 keV is probably the one which describes best the flexibility of the Bayesian approach, with the highest probability of having three sources as found by <cit.> and with a significant probability of four sources present, including a third footpoint, as found by <cit.> (see next Section). The positions of the two main footpoints are completely identified, while the location of the loop-top source is a bit uncertain. In addition, our results suggest that the looptop source is relatively weak, which is coherent with the figures in <cit.>, and relatively large. To show this, in Figure <ref> we plot the histograms that represent the distribution of the positions, FWHM and flux for each source. In order to provide a good estimate of the FWHM of the smaller sources, we included in the analysis detectors 1 and 2, that take the spatial resolution to 2.26 arcsec.


For the next energy band (30-50 keV) the number of sources in the next case decreases, with a high probability, to two; only 24% of the particles now pertain to the loop-top source.  The position of the (two) footpoints is clearly determined, while the third footpoint source is no longer clearly identified.  

Finally, for the last interval, 50-70 keV, the loop-top source has totally disappeared.  Only  less than 10% of the particles constitute a third source, which is not even identified as the loop-top source but probably represents only noise.








 §.§ Comparison with the results of Krucker & Lin


A similar analysis to the one in Section <ref> was done for the results found in <cit.>, which revealed a source structure with evolving complexity, including non-simultaneous brightening of different footpoints.  Figure 4 of <cit.> shows that the Southernmost footpoint (“Source 2”) brightened some 8 s after the Northernmost one (“Source 1”); see also <cit.>. Source 2 then slowly disappeared and a third source (“Source 3”), located more Southward than Source 2, reaches maximum brightness very shortly after the time of peak brightness of Source 2. To further explore these results, we first computed the source structure for the same three time intervals considered in Figure 4 of <cit.>: 11:06:04-11:06:12, 11:06:11-11:06:19, and 11:06:16-11:06:32.  The first two intervals correspond to two full rotation periods of the RHESSI spacecraft, while the third one corresponds to four rotations. In Figure <ref>, the distributions of the particles at the last iteration and the histograms of the positions of the footpoints are shown. The regularization parameter has been set to λ=3 for all the three time intervals. Note that the particles in the first time interval (11:06:08 UT) are more spread than in the others two most probably due to the noise level in the data.  We can clearly see that the southern source has changed its coordinates, while the northern one is more stable. Table <ref> gives the values of the parameters  obtained by employing the proposed method and the VFF at different time intervals. It is worth noting that the publicly available version of the VFF algorithm works only with circular shapes in case of multiple sources. Therefore, with the VFF algorithm we cannot give an estimate of the eccentricity, rotation angle and loop angle. The parameter values are very similar for both recovering methods, the only difference regards the FWHM, due to the fact that the sources reconstructed by the VVF algorithm  are  circular. Clearly, we can conclude that for both methods the southern source is given by very close sources which appear at different locations and time intervals. Because these sources do not appear simultaneously, it seems likely to conclude that the southern source is a single footpoint moving towards south–southwest.





To explore this more fully, it is desirable to examine the source structure at finer time resolution. The shortest practical time resolution for visibility-based imaging using RHESSI is ∼2 s, corresponding to half the satellite rotation period and so the shortest time interval for which a complete set of visibilities can be found.  Figure 3 of <cit.> shows results every 2 s, but it should be noted that in the interests of improving the statistics, each image corresponds to a (full rotation) 4 s integration time, so that there is some overlap among the images.  To further explore this evolving structure with our Bayesian method, we computed the source structure with the same 2 s time resolution and 4 s integration time.  Figure <ref> shows the distributions of the particles at the last iteration, and Figure <ref> shows the estimated x coordinate (arcseconds West), y coordinate (arcseconds North), and flux (photons s^-1 cm^-2) of the Southern footpoint source as functions of time.

The y-coordinate of this Southern footpoint, identified as “Source 2” by <cit.>, has a relatively sudden Southerly (i.e., towards more negative y-values) jump around 11:06:13 UT. This is the time at which <cit.> suggest that two different Southern footpoint sources may be present.  However, at this time our method (top right panel of Figure <ref>) shows only a small probability of an additional source a few arcseconds to the North of the main source of emission in this region. <cit.> identify this weak source as the continuation of “Source 2”  (which then promptly fades in intensity), and they identify the Southernmost bright patch as the appearance of a new “Source 3,” which rises in intensity very quickly (their Figure 4). However, the weak source identified as the continuation of “Source 2” by <cit.> is very short-lived (only one time interval; Figure <ref>) and therefore most likely a small artifact of noise. We instead contend that the Southernmost source, and not the one just to the North of it, is the continuation of “Source 2,” appearing at a slightly new location 3 arcseconds to the South. (This relatively sudden change in location of this source might also be an artifact of noise, considering that we are using a relatively short integration interval.)  The fact that all the parameters (x, y, and flux) associated with the brightest Southerly source all vary smoothly throughout the event confirms with a high probability that they are all manifestations of a single source in this region. We shall return to the significance of this result in Section <ref>.



















§ THE IMPACT OF THE PRIOR DISTRIBUTION



When using a Bayesian approach to solve an inverse problem, the results may depend on the choice of the prior distribution.
In our tests, we observed that the impact of the prior distribution is negligible whenever the data are highly informative, while, as expected, becomes more relevant when the data carry only little information on the estimated quantity. In order to show this, we present here results obtained by repeated analysis with different values of the prior probabilities for the number of sources and for the source types. 




 §.§ The prior on the number of sources


In order to assess the sensitivity of the results to the parameter λ of the prior, we select two data sets among those already analyzed. For both cases, we let λ vary between 0.5 and 4 with steps of 0.5. 

The first case was selected as an example of a simple case: it is the 30-80 keV interval, with a 8 second interval centered at 11:06:15 (see Section <ref>). The results are reported in Figure <ref>: the posterior probabilities change very little, compared to the prior probabilities, in this case; this is likely due to the good information content of the data.

As an example of a more challenging case, we report the posterior probabilities for the 20-30 keV interval (see Section <ref>). Here the number of sources is higher, and the results are more affected by a change of λ; specifically, the posterior probability for a four–source configuration reaches 43% when λ = 3.5 – still lower than the 53% of the three–source configuration, but close.
On the other hand, the MAP estimate is always a three–source model, even with the lowest λ.




 §.§ The prior on the source type


In order to assess the sensitivity of the results to the prior on the source type, we analyzed the 20-30 keV interval (see Section <ref>) using three different sets of values for p_C, p_E, p_L in the prior distribution (see eq. (<ref>)). These values were set to favour, each time, a different source type: the preferred source type was given a prior probability of 1/2, while the other two were given a prior probability of 1/4. The results are shown in Figure <ref>. In each box, we plot the posterior probability for the source type of each estimated source. The results show some expected dependence of the posterior value on the prior, when the source shape is more ambiguous: for instance, the North source is estimated to be a circular source in the first and in the third case, while it is estimated to be an elliptical source in the second case, when the elliptical sources are favoured. On the other hand, the South source is always estimated as an ellipse, though with different probabilities, for all the tested cases. Finally, the loop top is always estimated to be a circular source, though again with different probabilities; however, this result should be taken with care, because the looptop is the weaker source and, as such, the most difficult to estimate; this is also confirmed by the higher uncertainty on its location (cfr. Figure <ref>, panel (d)).





§ DISCUSSION AND CONCLUSIONS
 


The 2002 February 20 event at around 11:06 UT has been the source of intense study by a number of authors  <cit.>. Most of these studies have centered around the morphology of the source, and in particular that of the hard X-ray emission. It is generally agreed that there are several distinct sources of hard X-ray emission in the event, including a loop-top source that dominates at lower energies (presumably indicative of the thermal source), and a number of chromospheric footpoints that dominate at higher energies (presumably associated with the impact of accelerated electrons onto the dense chromosphere).

More interesting is the question as to the number and location of these chromospheric footpoints, in particular whether there are two <cit.> or three <cit.> such footpoints.  The presence of three separate locations of footpoint emission, especially if each brightens at a different time <cit.> would imply a much more complicated magnetic field geometry than a simple two-chromospheric-footpoint-plus-thermal-coronal-source geometry. More specifically, Figure 4 of <cit.> (see also Figure 2 of <cit.>) suggests a relatively steady Northern footpoint, followed by the appearance of two more southerly footpoints that brighten at different times and appear at different locations.  Although <cit.> claim (their Figure 4) that “Source 3” (the most Southerly footpoint source) replaces “Source 2” (the South footpoint source closest to the Northern footpoint), Figure 4 of that paper <cit.> shows that “Source 2” only gradually diminishes in intensity, and is still present even at the time of peak intensity of “Source 3.”  It is therefore still a possibility  that “Source 2” and “Source 3” are simply different manifestations of a single source that moves in a Southerly direction as the event progresses.

So is this event characterized by (1) a Northern footpoint and two separate Southern footpoint sources located at different locations and that brighten at different times (making three footpoint sources in all), or (2) a Northern footpoint and a single Southern source that moves (making two footpoint sources in all)?  The Bayesian method used in this paper is ideally suited to distinguishing between these two possibilities. The results (Sections 3 and 4) show that the most likely scenario is that there are only two principal sources at any given time, each of which has the highest probability of being a low-eccentricity (ϵ 0.3) elliptical source (in one case – for the Northern source at 11:06:15 UT – the inferred eccentricity is even sufficiently low that the source is best classified as “circular.”)

Consistent with the analyses of <cit.> and <cit.>, we find that the Northern footpoint source is fairly stationary, moving less than 1 arcsecond over the 16 second duration between the first and last images in Figure 5. However, even for the 11:06:15 UT image (middle panel in top row of Figure 5) that is the best case for showing three footpoint sources visible at the same time <cit.>, the Bayesian method reveals that the most likely geometry involves only two sources, with a single elliptical source (rather than two distinct footpoints) at the Southern location.  Analysis at higher (2 s) time resolution is fully consistent with this interpretation (Figures <ref> and <ref>).  While <cit.> suggest that two different Southern foopoint sources (“Source 2” and “Source 3”) may be present near the peak of the event, the time (∼ 11:06:13 UT) at which this interpretation is most convincing corresponds to a time at which the Southern source is split into two components.  One of these, previously identified by <cit.> as the continuation of “Source 2,” is, however, very weak and short-lived; the other, identified by <cit.> as the sudden appearance of a new “Source 3,” has a location, size, and intensity that, given the statistical uncertainties, is more probably consistent with a continuation of “Source 2” (Figure <ref>).  Our analysis therefore favors an interpretation in which the weak source is to be disregarded, so that the emission in the Southern part of the flare arises from only one footpoint, with smoothly varying size, intensity, and location.

This Southern source does move quite rapidly (left and middle panels of Figure <ref>), drifting about 7 arc seconds in a South-Southwest direction over a period of about 15 seconds, corresponding to a velocity ∼350 km s^-1. <cit.> studied footpoint motions in the 2003 July 23 (∼00:30 UT) event; Figure 3 of that paper shows that the fastest motions (in the Northern footpoint of that event) had a speed of order 50 km s^-1 <cit.>.  The relative motion of the various footpoints in the 2002 July 23 event is qualitatively similar to what we have found for the 2002 February 20 event; the motion in both events is such that the inter-footpoint separation increases with time.  Such a behavior is consistent with the successive activation of higher and higher field lines in an arcade, as in the well-known CSHKP model <cit.>.

<cit.> also showed that the speed of the footpoint motion in the 2003 July 23 event correlates well with the hard X-ray intensity of the footpoint source, suggesting that the rate of electron acceleration is correlated with the speed of the footpoint motion. Indeed higher footpoint speeds result in higher induced electric fields perpendicular to both the magnetic field direction and the direction of motion: 𝐄 = (𝐯/c) ×𝐁.  Since the 2002 February 20 event was fairly small (a GOES C-class) flare, we take B = 100 G as a representative value; taking V = 10^27 cm^-3, the associated energy content (B^2/8π) V in the magnetic field is then of order 10^29 ergs, consistent with the GOES classification.

The electric field associated with the inferred 350 km s^-1 velocity is thus E ≃ 0.12 statvolt cm^-1 ≃ 35 V cm^-1.  This is several orders of magnitude larger than the <cit.> field E_D ≃ 10^-8 n/T ≃ 10^-5 V cm^-1 (using n ≃ 10^10 cm^-3 and T ≃ 10^7 K) that is necessary to accelerate electrons out of a background plasma. It is also sufficiently large that it can accelerate electrons to the required ≃30 keV energies in a distance less than 10 meters, which is orders of magnitude smaller than the characteristic size of a flare loop.  Such electric field magnitudes and acceleration lengths are, however, nicely compatible with those associated with super-Dreicer acceleration in current sheets <cit.>.


















 
<cit.> also showed that the speed of the footpoint motion in the 2003 July 23 event correlates well with the hard X-ray intensity of the footpoint source, suggesting that the rate of electron acceleration is correlated with the speed of the footpoint motion. Indeed higher footpoint speeds result in higher induced electric fields perpendicular to both the magnetic field direction and the direction of motion: 𝐄 = (𝐯/c) ×𝐁.  Since the 2002 February 20 event was fairly small (a GOES C-class) flare, we take B = 100 G. The electric field associated with the inferred 200 km s^-1 velocity is thus E ≃ 0.07 statvolt cm^-1 ≃ 20 V cm^-1.  This is several orders of magnitude larger than the <cit.> field E_D ≃ 10^-8 n/T ≃ 10^-5 V cm^-1 that is necessary to accelerate electrons out of a background plasma of density n (cm^-3) and temperature T (K) and indeed is sufficiently large that it can accelerate electrons to the required ≃ 20 keV energies in a distance of only 2 meters, which is orders of magnitude smaller than the characteristic size of a flare loop.  Such electric field magnitudes and acceleration lengths are, however, nicely compatible with those associated with super-Dreicer acceleration in current sheets <cit.>.
 



§ ACKNOWLEDGEMENTS

AGE was supported by grant NNX17AI16G from NASA's Heliophysics Supporting Research program. FS, AS, AMM and MP have been supported by the H2020 grant Flare Likelihood And Region Eruption foreCASTing (FLARECAST), project number 640216. The authors kindly thank Dr. Kim Tolbert for providing access to results from previous studies.


§ APPENDIX: THE ADAPTIVE SEQUENTIAL MONTE CARLO ALGORITHM



The main reason to introduce SMC samplers is that, due to the complexity of the posterior distribution, drawing random samples from it is not trivial, so that basic Monte Carlo approaches would either fail or be extremely inefficient.
The main idea behind SMC samplers is to construct an auxiliary sequence of distributions {π_i(x)}_i=1,…,I, such that the first distribution is the prior distribution π_1(X)=π_pr(X), the last distribution is the target posterior distribution π_I(X)=π_po(X|y), and the sequence is “smooth.”  We draw a random sample from π_1(x), make the corresponding sample points evolve according to properly selected Markov kernels, and re-weigh them; eventually, the sample points at the final iteration I will be distributed according to the target distribution, which in our case is the posterior distribution.

At the i-th iteration we use


    π_i(X)= ℒ(y|X)^γ_iπ_pr(X),

where γ_1=0, γ_I=1 and γ_i ∈ [0,1], with γ_i+1 > γ_i, ∀ i. In other words, we start from the prior distribution and we go progressively towards the posterior distribution by systematically increasing the exponent γ_i.

The choice of the γ_i is very important for the success of the algorithm; here we employ an adaptive approach in which the sequence of the distributions is determined on-line, with an empirical measure of the distance between two subsequent distributions (i.e., it depends on how much the distribution at the iteration (i+1) differs from the one at the iteration i). For more details we refer the reader to <cit.>.

Since the Bayesian model considered here has a variable dimension (see Equation (<ref>)), we need to introduce the concept of reversible jump moves, which allow us to explore the space χ of the parameters with which flares are modeled. In our analysis we consider the following moves:



  * birth and death;

  * change;

  * split and merge; and

  * update.

Birth and death moves are attempted with probability 1/3. The type of the source is chosen according to ρ_2 in Equation (<ref>) and the parameters of the newly born source are uniformly distributed following Equation (<ref>). A change move is done during each iteration and transforms a circle into an ellipse (and vice versa), or an ellipse into a loop (and vice versa). (We do not consider the transition from a circle to a loop (and vice versa) since we want to have smooth transformations.) Split and merge moves are applied only to the ellipses and allow one ellipse to split in two (or two ellipses to merge in one) and are attempted at each iteration. The update move is a classical move of the Metropolis-Hastings algorithm and it is done at the end of each iteration. This move does not change the type of the source, but it does determine how each parameter of each source evolves independently from the others


In summary, this process allows the description of physical situations in which flaring sources may appear, disappear, and change their shape and physical properties. From the mathematical viewpoint, the split, merge and change moves are introduced here for the first time.  


In Table <ref>, we give a sketch of the ASMC algorithm. For more details, we refer the reader to  <cit.>.








aasjournal



