1 Introduction
The RHESSI instrument (Lin et al., 2002) was designed to investigate particle acceleration and energy release in solar flares through imaging spectroscopy of hard Xray and gammaray emissions from keV to MeV. Over the course of its 15plusyear lifetime, RHESSI has provided thousands of hard Xray images of solar flares. Modeling of the often complex structures thus revealed (Emslie et al., 2003), and relating them to theoretical models of the spatial distributions of hard Xrays in solar flares (e.g., Brown, 1971, 1973; Brown & McClymont, 1975; Emslie & Machado, 1987) 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 (e.g., Krucker et al., 2007, 2008; Krucker & Lin, 2008), “thicktarget” nonthermal sources confined to the corona (Xu et al., 2008; Guo et al., 2012a, b) and chromospheric footpoints (e.g., Brown et al., 2002; Aschwanden et al., 2002a; Sato, 2004; Kontar et al., 2008).
A particularly interesting example of the importance of such subsource 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 Sui et al. (2002) show images of the event produced using the Clean algorithm (Dennis & Pernak, 2009) and the MEMSato maximum entropy method (Sato et al., 1999). 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 (Masuda et al., 1995).

Figure 2 of Aschwanden et al. (2004) shows further evidence for a twofootpoint structure, using five different imaging methods, including forwardfit (Aschwanden et al., 2002b), Clean (Dennis & Pernak, 2009), Pixon (Metcalf et al., 1996), and two different maximumentropy methods, one operating in space and the other in visibility (i.e., Fourier) space.

Figure 4 of Aschwanden et al. (2002a) shows the source centroids at different hard Xray energies all concentrated in these two footpoints.

Adding to the complexity, Figure 1 of Veronig et al. (2005) shows “two distinct footpoints and a faint looptop source” located some 20 arcseconds toward the limb compared to the footpoints.

In their analysis of this event, Guo et al. (2011) discuss (cf. their Figure 3) “a complicated pattern of low and highenergy sources.”

Krucker & Lin (2002) 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 (their Figure 4; see also Figure 2 of Vilmer et al., 2002). 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 Xray kernel” (Wang et al., 2002), providing evidence for energy release to due to magnetic reconnection/reconfiguration and possible evidence for an acceleration site in the chromosphere (Brown et al., 2009).

Piana et al. (2007) 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 countbased images. It is therefore especially interesting that while the countbased images in Piana et al. (2007) do suggest the possibility of a third footpoint source (see especially the 3438 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, Krucker et al. (2003) 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 Xray 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 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, rotatingcollimatorbased imaging techniques (Hurford et al., 2002). 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” (twodimensional spatial Fourier components of the image), with each detector sampling Fourier components tracing a circle in 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 space, which is then improved using methods like Maximum Entropy (Bong et al., 2006), Clean (Högbom, 1974; Dennis & Pernak, 2009), or Pixon (Metcalf et al., 1996);

Inversion of count modulation profiles by means of optimally stopped iterative schemes (Benvenuto et al., 2013);

Interpolating the (relatively sparse) visibility data to obtain a smooth surface in space, which can then be Fouriertransformed to produce an image in the domain (Massone et al., 2009);

Forwardfitting the visibilities with a parametric model in order to obtain the bestfit parameters
(Aschwanden et al., 2002b).
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 Xray 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 simplyshaped geometrical objects. The parameters of these objects are then estimated through nonlinear optimization techniques. To the best of our knowledge, to date this visibility forward fitting (VFF) algorithm
(Xu et al., 2008) 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 (Dennis & Pernak, 2009) 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 subregions in a hard Xray 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 2 we describe the method. In Section 3 we describe the RHESSI observations used for the analysis, i.e. the relatively wellstudied 2002 February 20 11:06 UT event, and we compare the obtained results with the ones described in (Sui et al., 2002; Krucker & Lin, 2002). Section 4 discusses the impact of the prior distribution on the results. Our conclusions are offered in Section 5.
2 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).
2.1 Parametrization of flare sources
Hard Xray images in solar flares (e.g., Emslie et al., 2003; Xu et al., 2008) 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 twodimensional 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 (Aschwanden et al., 2002b), for the sake of completeness we provide here a brief description of it.
The various model sources used are characterized as follows (see Figure 1):

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

Ellipse: an ellipse is defined by the location of its center, the semimajor axis corresponding to the 50% intensity level that extends furthest from the center, the rotation angle of that direction, and the eccentricity
where and correspond to the length of the semimajor and semiminor axes, respectively.

Loop: a “loop” is defined by the above ellipse parameters plus an additional parameter which indicates the degree of curvature of the loop shape (see Figure 1). 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
(1) 
where the last three parameters can possibly (for the case of a circle or an ellipse) be zero or indeterminate. Figure 1 provides a visual description of the parameters.
The parametrization in Equation (1) has already been presented in Aschwanden et al. (2002b), where a forwardfitting method to reconstruct hard Xray 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 Aschwanden et al. (2002b) 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 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:
(2) 
where represents the number of sources, the source types and the source parameters for each source.
2.2 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 lowercase 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 (
2), conditioned on the data, namely a set of observed visibilities. From Bayes’ theorem we have
(3) 
where
is the posterior probability distribution representing the conditional probability density of
given the data , is the likelihood function that expresses the likelihood that the data is indeed due to the image ,is the prior probability density which describes the properties of the source configuration that we want to recover, and
is a normalization factor, independent of . 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 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 variabledimension model, i.e., the state space of the sources is defined as follows:
(4) 
where is the state space for source number of the sources present.
2.2.1 The prior distribution
As in all Bayesian problems, the choice of the prior distribution , 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 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 can be separated as follows:
(5) 
where the () 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
:(6) 
and we will limit . In addition, we shall limit , the maximum number of sources, to , because past observations (e.g., Emslie et al., 2003) 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
(7) 
where . Now, in general the morphology of solar flare hard Xray sources is dependent on the photon energy observed: images at lower energies keV generally are dominated by extended looplike (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 (see, however, Xu et al., 2008; Guo et al., 2012a, b), while the latter are generally considered as chromospheric footpoint sources, which are compact because of the much higher density and so electron stopping distance (Brown, 1972; Emslie, 1978). In our analysis we will use if the energy is high and otherwise. We stress that these are only assumed prior distributions; they do not represent the final (datainformed) result.
The individual source parameters are a priori
assumed to be uniformly distributed in their respective ranges:
(8)  
where is the map center, indicates the field of view of the map in each direction and represents the (complex) visibilities. Combining Equations (6), (7), and (8), the prior distribution is given by
2.2.2 Likelihood function
Since visibilities are linear combinations of independent photon count numbers (Hurford et al., 2002), we can make the standard approximation that noise on the visibilities is Gaussian additive; the likelihood function is thus given by
(9) 
where
indicates the forward operator, i.e. the projector from the imaging parameter space onto the visibility space. The standard deviations
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.2.3 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 Del Moral et al. (2006) for the original article in which the SMC method is introduced, and to Sorrentino et al. (2014) 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
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;
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(10) 
3 Application to RHESSI data  the 2002 February 20 event
3.1 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 1, 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 (Xu et al., 2008) 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 (Dennis & Pernak, 2009) 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 Sui et al. (2002) and Krucker & Lin (2002), in our experiments, we use RHESSI detectors 3 through 9, a field of view arcseconds; further, we choose mapcenters at and for Sections 3.2 and 3.3, which compare our results with those of the respective previous authors. The number of particles is set to , the maximum number of sources admitted in a map is and, coherently to what done in VFF, the noise level is given by adding a systematic error to the error computed by the RHESSI SSW routines, i.e.,
where is the absolute value of the visibilities. Therefore, the only parameter that needs to be tuned is , and we assumed values in the range . In most cases slightly changing will not influence the reconstruction, but only the convergence speed; for more details, we refer the reader to Section 4.
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 hsi_vis_fwdfit_func.
3.2 Comparison with the results in Sui et al
In this section, we focus on the time interval 11:06:1011:06:24 UT and several energy intervals spanning the range from 6 keV to 70 keV. We compare the results found in Sui et al. (2002) with those using our Bayesian approach.
In Sui et al. (2002), the Clean (see Hurford et al., 2002) and the maximum entropy (see Sato et al., 1999) reconstruction techniques have been used to study the number of sources with different energy bands at a fixed time interval. In particular, Sui et al. (2002) found (their Figures 3 and 4) that the map in the lowenergy band (610 keV) shows a single structure which is located between the two footpoints observed at higher energies. In the interval 1014 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 looptop source, while in the last energy band 5070 keV only the two footpoints are present and the loopsource is not visible.
For the prior distribution used in our approach (see equation (5)), we set the number of sources according to the results found in Sui et al. (2002), i.e., for the energy band 610 keV, for 1014 keV and 5070 keV, and for 1450 keV.
In Figure 2 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 3 and 4, together with the reconstructions obtained by Clean and MEMNJIT, respectively; the values in the color bands represent the weights of each particle, with , with . 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 2 and 3 we can draw some quantitative conclusions about the number and locations of sources in the map. For the lowest energy band (610 keV) the probability of having two sources is almost (Figure 2). At first glance, this result seems to run counter to the results of Sui et al. (2002), 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 3, top left panel and Figure 5). 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 (7), the algorithm recovers a single loopsource with a significantly higher probability and with parameters and uncertainties close to the ones obtained by VFF (see Table 1). Thus, we can conclude that only a single source is found in the energy band 610 keV, which we interpret as a thermal coronal source situated at the apex of a magnetic structure connecting the footpoints.
r  

Bayesian approach  Loop  908.17  263.79  7784.65  14.05  0.74  182.0  106.7 
Standard deviation  
Visibility forward fitting  Loop  908.78  263.82  7905.22  14.87  0.72  177.6  124.8 
Standard deviation  0.85  0.42  278.12  1.01  0.11  5.5  41.9 
In the next energy band (1014 keV), the map (Figure 3(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 2, we see a higher probability of three sources (the two footpoints and the looptop source); the probability of having four sources present (red band in Figure 2) is around . By comparing the weights of the particles in Figure 3, 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 1420 keV energy band shows a predominance of three sources (Figure 2): Figure 3(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 2030 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 Sui et al. (2002) and with a significant probability of four sources present, including a third footpoint, as found by Krucker & Lin (2002) (see next Section). The positions of the two main footpoints are completely identified, while the location of the looptop source is a bit uncertain. In addition, our results suggest that the looptop source is relatively weak, which is coherent with the figures in Sui et al. (2002), and relatively large. To show this, in Figure 6 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 (3050 keV) the number of sources in the next case decreases, with a high probability, to two; only of the particles now pertain to the looptop 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, 5070 keV, the looptop source has totally disappeared. Only less than of the particles constitute a third source, which is not even identified as the looptop source but probably represents only noise.
3.3 Comparison with the results of Krucker & Lin
A similar analysis to the one in Section 3.2 was done for the results found in Krucker & Lin (2002), which revealed a source structure with evolving complexity, including nonsimultaneous brightening of different footpoints. Figure 4 of Krucker & Lin (2002) shows that the Southernmost footpoint (“Source 2”) brightened some 8 s after the Northernmost one (“Source 1”); see also Vilmer et al. (2002). 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 (Krucker & Lin, 2002): 11:06:0411:06:12, 11:06:1111:06:19, and 11:06:1611: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 7, 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 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 2 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.
Time  Bayesian method  VFF  

r  r  
11:06:08 UT  Ellipse  901.95  271.80  7.32  3.41  0.12  92.0  0.4  902.58  271.30  11.82  5.53 
Ellipse  905.63  252.94  17.25  19.72  0.18  27.5  0.2  906.81  252.15  12.11  8.67  
11:06:15 UT  Circle  903.00  271.96  24.20  4.03  0.10  22.8  0.1  903.19  271.99  25.13  3.99 
Ellipse  906.99  249.38  25.18  8.39  0.17  40.6  0.2  907.14  249.53  24.94  7.80  
11:06:24 UT  Ellipse  902.95  273.22  14.39  4.68  0.19  46.3  0.2  902.89  273.16  14.96  4.67 
Ellipse  908.15  247.60  24.27  6.20  0.30  66.2  0.0  908.04  247.60  24.96  6.18 
To explore this more fully, it is desirable to examine the source structure at finer time resolution. The shortest practical time resolution for visibilitybased 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 Krucker & Lin (2002) 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 8 shows the distributions of the particles at the last iteration, and Figure 9 shows the estimated coordinate (arcseconds West), coordinate (arcseconds North), and flux (photons s cm) of the Southern footpoint source as functions of time.
The coordinate of this Southern footpoint, identified as “Source 2” by Krucker & Lin (2002), has a relatively sudden Southerly (i.e., towards more negative values) jump around 11:06:13 UT. This is the time at which Krucker & Lin (2002) suggest that two different Southern footpoint sources may be present. However, at this time our method (top right panel of Figure 8) shows only a small probability of an additional source a few arcseconds to the North of the main source of emission in this region. Krucker & Lin (2002) 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 Krucker & Lin (2002) is very shortlived (only one time interval; Figure 8) 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 (, , 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 5.
4 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.
4.1 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 and with steps of .
The first case was selected as an example of a simple case: it is the 3080 keV interval, with a 8 second interval centered at 11:06:15 (see Section 3.3). The results are reported in Figure 10: 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 2030 keV interval (see Section 3.2). 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 when – still lower than the of the three–source configuration, but close. On the other hand, the MAP estimate is always a three–source model, even with the lowest .
4.2 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 2030 keV interval (see Section 3.2) using three different sets of values for in the prior distribution (see eq. (7)). These values were set to favour, each time, a different source type: the preferred source type was given a prior probability of , while the other two were given a prior probability of . The results are shown in Figure 11. 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 3, panel (d)).
5 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 (Krucker & Lin, 2002; Sui et al., 2002; Aschwanden et al., 2002a; Wang et al., 2002; Aschwanden et al., 2004; Veronig et al., 2005; Piana et al., 2007; Guo et al., 2011; Falewicz, 2014). Most of these studies have centered around the morphology of the source, and in particular that of the hard Xray emission. It is generally agreed that there are several distinct sources of hard Xray emission in the event, including a looptop 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 (Sui et al., 2002) or three (Krucker & Lin, 2002; Vilmer et al., 2002) such footpoints. The presence of three separate locations of footpoint emission, especially if each brightens at a different time (Krucker & Lin, 2002; Vilmer et al., 2002) would imply a much more complicated magnetic field geometry than a simple twochromosphericfootpointplusthermalcoronalsource geometry. More specifically, Figure 4 of Krucker & Lin (2002) (see also Figure 2 of Vilmer et al. (2002)) 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 Krucker & Lin (2002) 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 (see also Vilmer et al., 2002) 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 loweccentricity (
) 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 Krucker & Lin (2002) and Vilmer et al. (2002), 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 (Figure 4 of Krucker & Lin, 2002), 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 8 and 9). While Krucker & Lin (2002) 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 Krucker & Lin (2002) as the continuation of “Source 2,” is, however, very weak and shortlived; the other, identified by Krucker & Lin (2002) 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 9). 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 9), drifting about 7 arc seconds in a SouthSouthwest direction over a period of about 15 seconds, corresponding to a velocity 350 km s. Krucker et al. (2003) 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 (comparable to those inferred by Sakao et al., 1998). 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 interfootpoint 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 wellknown CSHKP model (Carmichael, 1964; Sturrock, 1966; Hirayama, 1974; Kopp & Pneuman, 1976).
Krucker et al. (2003) also showed that the speed of the footpoint motion in the 2003 July 23 event correlates well with the hard Xray 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: . Since the 2002 February 20 event was fairly small (a GOES Cclass) flare, we take B = 100 G as a representative value; taking cm, the associated energy content in the magnetic field is then of order ergs, consistent with the GOES classification.
The electric field associated with the inferred 350 km s velocity is thus statvolt cm V cm. This is several orders of magnitude larger than the Dreicer (1959) field V cm (using cm and 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 superDreicer acceleration in current sheets (see, e.g., Somov, 1992; Litvinenko, 1996).
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 , such that the first distribution is the prior distribution , the last distribution is the target posterior distribution , and the sequence is “smooth.” We draw a random sample from , make the corresponding sample points evolve according to properly selected Markov kernels, and reweigh them; eventually, the sample points at the final iteration will be distributed according to the target distribution, which in our case is the posterior distribution.
At the th iteration we use
(11) 
where , and , with . In other words, we start from the prior distribution and we go progressively towards the posterior distribution by systematically increasing the exponent .
The choice of the is very important for the success of the algorithm; here we employ an adaptive approach in which the sequence of the distributions is determined online, with an empirical measure of the distance between two subsequent distributions (i.e., it depends on how much the distribution at the iteration differs from the one at the iteration ). For more details we refer the reader to Sorrentino et al. (2014).
Since the Bayesian model considered here has a variable dimension (see Equation (4)), 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 . The type of the source is chosen according to in Equation (7) and the parameters of the newly born source are uniformly distributed following Equation (8). 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 MetropolisHastings 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.
References
 Aschwanden et al. (2002a) Aschwanden, M. J., Brown, J. C., & Kontar, E. P. 2002a, Sol. Phys., 210, 383
 Aschwanden et al. (2004) Aschwanden, M. J., Metcalf, T. R., Krucker, S., et al. 2004, Sol. Phys., 219, 149
 Aschwanden et al. (2002b) Aschwanden, M. J., Schmahl, E., & RHESSI Team. 2002b, Sol. Phys., 210, 193
 Benvenuto et al. (2013) Benvenuto, F., Schwartz, R., Piana, M., & Massone, A. M. 2013, A&A, 555, A61
 Bong et al. (2006) Bong, S.C., Lee, J., Gary, D. E., & Yun, H. S. 2006, ApJ, 636, 1159
 Brown (1971) Brown, J. C. 1971, Sol. Phys., 18, 489
 Brown (1972) —. 1972, Sol. Phys., 26, 441
 Brown (1973) —. 1973, Sol. Phys., 31, 143
 Brown et al. (2002) Brown, J. C., Aschwanden, M. J., & Kontar, E. P. 2002, Sol. Phys., 210, 373
 Brown & McClymont (1975) Brown, J. C., & McClymont, A. N. 1975, Sol. Phys., 41, 135
 Brown et al. (2009) Brown, J. C., Turkmani, R., Kontar, E. P., MacKinnon, A. L., & Vlahos, L. 2009, A&A, 508, 993
 Carmichael (1964) Carmichael, H. 1964, NASA Special Publication, 50, 451
 Del Moral et al. (2006) Del Moral, P., Doucet, A., & Jasra, A. 2006, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68, 411
 Delaboudinière et al. (1995) Delaboudinière, J.P., Artzner, G. E., Brunaud, J., et al. 1995, Sol. Phys., 162, 291
 Dennis & Pernak (2009) Dennis, B. R., & Pernak, R. L. 2009, ApJ, 698, 2131
 Dreicer (1959) Dreicer, H. 1959, Physical Review, 115, 238
 DuvalPoo et al. (2017) DuvalPoo, M. A., Massone, A. M., & Piana, M. 2017, in Sampling Theory and Applications (SampTA), 2017 International Conference on, IEEE, 677–681
 DuvalPoo et al. (2017) DuvalPoo, M. A., Piana, M., & Massone, A. M. 2017, ArXiv eprints, arXiv:1708.03877
 Emslie (1978) Emslie, A. G. 1978, ApJ, 224, 241
 Emslie et al. (2003) Emslie, A. G., Kontar, E. P., Krucker, S., & Lin, R. P. 2003, ApJ, 595, L107
 Emslie & Machado (1987) Emslie, A. G., & Machado, M. E. 1987, Sol. Phys., 107, 263
 Falewicz (2014) Falewicz, R. 2014, ApJ, 789, 71
 Felix et al. (2017) Felix, S., Bolzern, R., & Battaglia, M. 2017, ApJ, 849, 10
 Guo et al. (2012a) Guo, J., Emslie, A. G., Kontar, E. P., et al. 2012a, A&A, 543, A53
 Guo et al. (2012b) Guo, J., Emslie, A. G., Massone, A. M., & Piana, M. 2012b, ApJ, 755, 32
 Guo et al. (2011) Guo, J., Liu, S., Fletcher, L., & Kontar, E. P. 2011, ApJ, 728, 4
 Hirayama (1974) Hirayama, T. 1974, Sol. Phys., 34, 323
 Högbom (1974) Högbom, J. A. 1974, A&AS, 15, 417
 Hurford et al. (2002) Hurford, G. J., Schmahl, E. J., Schwartz, R. A., et al. 2002, Sol. Phys., 210, 61
 Kontar et al. (2008) Kontar, E. P., Hannah, I. G., & MacKinnon, A. L. 2008, A&A, 489, L57
 Kopp & Pneuman (1976) Kopp, R. A., & Pneuman, G. W. 1976, Sol. Phys., 50, 85
 Krucker et al. (2003) Krucker, S., Hurford, G. J., & Lin, R. P. 2003, ApJ, 595, L103
 Krucker & Lin (2002) Krucker, S., & Lin, R. P. 2002, Sol. Phys., 210, 229
 Krucker & Lin (2008) —. 2008, ApJ, 673, 1181
 Krucker et al. (2007) Krucker, S., White, S. M., & Lin, R. P. 2007, ApJ, 669, L49
 Krucker et al. (2008) Krucker, S., Battaglia, M., Cargill, P. J., et al. 2008, A&A Rev., 16, 155
 Lin et al. (2002) Lin, R. P., Dennis, B. R., Hurford, G. J., et al. 2002, Sol. Phys., 210, 3
 Litvinenko (1996) Litvinenko, Y. E. 1996, ApJ, 462, 997
 Massone et al. (2009) Massone, A. M., Emslie, A. G., Hurford, G. J., et al. 2009, ApJ, 703, 2004
 Masuda et al. (1995) Masuda, S., Kosugi, T., Hara, H., et al. 1995, PASJ, 47, 677
 Metcalf et al. (1996) Metcalf, T. R., Hudson, H. S., Kosugi, T., Puetter, R. C., & Pina, R. K. 1996, ApJ, 466, 585
 Piana et al. (2007) Piana, M., Massone, A. M., Hurford, G. J., et al. 2007, ApJ, 665, 846
 Sakao et al. (1998) Sakao, T., Kosugi, T., & Masuda, S. 1998, in Astrophysics and Space Science Library, Vol. 229, Observational Plasma Astrophysics : Five Years of YOHKOH and Beyond, ed. T. Watanabe & T. Kosugi, 273
 Sato (2004) Sato, J. 2004, in COSPAR Meeting, Vol. 35, 35th COSPAR Scientific Assembly, ed. J.P. Paillé, 3787
 Sato et al. (1999) Sato, J., Kosugi, T., & Makishima, K. 1999, PASJ, 51, 127
 Somov (1992) Somov, B. V., ed. 1992, Astrophysics and Space Science Library, Vol. 172, Physical processes in solar flares.
 Sorrentino et al. (2014) Sorrentino, A., Luria, G., & Aramini, R. 2014, Inverse Problems, 30, 045010
 Sturrock (1966) Sturrock, P. A. 1966, Nature, 211, 695
 Sui et al. (2002) Sui, L., Holman, G. D., Dennis, B. R., et al. 2002, Sol. Phys., 210, 245
 Veronig et al. (2005) Veronig, A. M., Brown, J. C., Dennis, B. R., et al. 2005, ApJ, 621, 482
 Vilmer et al. (2002) Vilmer, N., Krucker, S., Lin, R. P., & Rhessi Team. 2002, Sol. Phys., 210, 261
 Wang et al. (2002) Wang, H., Ji, H., Schmahl, E. J., et al. 2002, ApJ, 580, L177
 Xu et al. (2008) Xu, Y., Emslie, A. G., & Hurford, G. J. 2008, ApJ, 673, 576
Comments
There are no comments yet.