Identification of multiple hard X-ray sources in solar flares: A Bayesian analysis of the February 20 2002 event

by   Federica Scacchiano, et al.

Hard X-ray emission in solar flares 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. In particular, 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.



There are no comments yet.


page 11


eBASCS: Disentangling Overlapping Astronomical Sources II, using Spatial, Spectral, and Temporal Information

The analysis of individual X-ray sources that appear in a crowded field ...

Automatic Detection of Occulted Hard X-ray Flares Using Deep-Learning Methods

We present a concept for a machine-learning classification of hard X-ray...

Visibility Interpolation in Solar Hard X-ray Imaging: Application to RHESSI and STIX

Space telescopes for solar hard X-ray imaging provide observations made ...

Machine learning applied to single-shot x-ray diagnostics in an XFEL

X-ray free-electron lasers (XFELs) are the only sources currently able t...

Identification of high-energy astrophysical point sources via hierarchical Bayesian nonparametric clustering

The light we receive from distant astrophysical objects carries informat...

Toward AI-enhanced online-characterization and shaping of ultrashort X-ray free-electron laser pulses

X-ray free-electron lasers (XFELs) as the world`s most brilliant light s...

A neural simulation-based inference approach for characterizing the Galactic Center γ-ray excess

The nature of the Fermi gamma-ray Galactic Center Excess (GCE) has remai...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 X-ray and gamma-ray emissions from  keV to  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 (Emslie et al., 2003), and relating them to theoretical models of the spatial distributions of hard X-rays 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), “thick-target” 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 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.

  1. Figure 3 and 4 of Sui et al. (2002) show images of the event produced using the Clean algorithm (Dennis & Pernak, 2009) and the MEM-Sato 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).

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

  3. Figure 2 of Falewicz (2014) again shows evidence for two chromospheric footpoints, a morphology also supported by SOHO EIT (Delaboudinière et al., 1995) observations.

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

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

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

  7. 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 X-ray 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).

  8. 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 count-based images. It is therefore especially interesting that while the count-based images in Piana et al. (2007) 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, 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 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 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 (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” (two-dimensional 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 Fourier-transformed to produce an image in the domain (Massone et al., 2009);

  • Reconstruction methods from visibilities using compressed sensing approaches (Duval-Poo et al., 2017, 2017; Felix et al., 2017)

  • Forward-fitting the visibilities with a parametric model in order to obtain the best-fit 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 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

(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 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 2 we describe the method. In Section 3 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 (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 X-ray 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 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 (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 semi-major 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 semi-major and semi-minor 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


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.

Figure 1: The three source types: circle, ellipse and loop

The parametrization in Equation (1) has already been presented in Aschwanden et al. (2002b), 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 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:


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 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 (

2), conditioned on the data, namely a set of observed visibilities

. From Bayes’ theorem we have



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 variable-dimension model, i.e., the state space of the sources is defined as follows:


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:


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



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


where . Now, in general the morphology of solar flare hard X-ray sources is dependent on the photon energy observed: images at lower energies  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 (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 (data-informed) result.

The individual source parameters are a priori

assumed to be uniformly distributed in their respective ranges:


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



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


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 map-centers 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:10-11: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 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.

Figure 2: Posterior probability for the number of sources, color-coded: in blue the probability of a single source, in green that of two sources, in yellow that of three sources, and so on. The parameter has been set as follows: for the energy band 6-10 keV, for 10-14 keV and 50-70 keV, and for 14-50 keV.

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 6-10 keV, for 10-14 keV and 50-70 keV, and for 14-50 keV.

(a) 6-10keV
(b) 10-14keV
(c) 14-20keV
(d) 20-30keV
(e) 30-50keV
(f) 50-70keV
Figure 3: Distribution of the particles in the last iteration and (black lines) Clean reconstruction. The parameter has been set as in Figure 2. Clean contour levels are 0.08, 0.2, 0.3, and 0.7.
(a) 6-10keV
(b) 10-14keV
(c) 14-20keV
(d) 20-30keV
(e) 30-50keV
(f) 50-70keV
Figure 4: Distribution of the particles in the last iteration and (black lines) MEM-NJIT reconstruction. The parameter has been set as in Figure 2. MEM-NJIT contour levels are 0.08, 0.2, 0.3, and 0.7.

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 MEM-NJIT, 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 (6-10 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 loop-source 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 6-10 keV, which we interpret as a thermal coronal source situated at the apex of a magnetic structure connecting the footpoints.

(a) Clean
(b) VFF
(c) Bayes
Figure 5: Images reconstructed in the 6-10 keV interval by Clean (left), forward fit (middle) and the Bayesian approach with and , (right).
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
Table 1: Values of the source parameters provided by VFF and by the Bayesian method when the probability of having a loop source is increased to and conditioned to the case when a single loop is reconstructed. The keV channel is considered.

In the next energy band (10-14 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 loop-top 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 14-20 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 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 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 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 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 (30-50 keV) the number of sources in the next case decreases, with a high probability, to two; only 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 of the particles constitute a third source, which is not even identified as the loop-top source but probably represents only noise.

Figure 6: Histograms of the Monte Carlo samples for the three sources in the energy band 20-30 keV; these histograms can be interpreted as (un-normalized) posterior distributions for the respective parameters.

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 non-simultaneous 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: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 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.

Figure 7: Distribution of the particles (first row) for the 20 February flare at 11:06:08 UT, 11:06:15 UT, 11:06:24 UT (30-80 keV, 7” resolution). The time intervals are s, s, and s, respectively. The histogram of the positions of northern and southern source are shown in the second and third row, respectively.
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
Table 2: Values of the source parameters provided by VFF and by the Bayesian method for the three time intervals considered in Figure 7

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 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 short-lived (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.

(a) 11:06:05
(b) 11:06:07
(c) 11:06:09
(d) 11:06:11
(e) 11:06:13
(f) 11:06:15
(g) 11:06:17
(h) 11:06:19
(i) 11:06:21
(j) 11:06:23
(k) 11:06:25
(l) 11:06:27
(m) 11:06:29
(n) 11:06:31
(o) 11:06:33
Figure 8: Distribution of the particles obtained by applying the Bayesian approach to visibilities taken every 2 s, with ∼ 4 s integration time. Energy interval: 30-80 keV, 7” resolution. All times in UT. The parameter has been set to 2 for all the time intervals.
Figure 9: Temporal evolution of the coordinate, the coordinate and the flux of the southern source. Error bars represent one standard deviation of the posterior distribution. The -axis represents the seconds after 11:06 UT.

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 30-80 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 20-30 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 .

(a) 11:06:15 UT
(b) 20-30 keV
Figure 10: Posterior probabilities for the number of sources, for different values of the prior parameter .

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 20-30 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)).

Figure 11: Posterior probabilities for the source types, using different a priori probability sets. Left: . Middle: . Right: .

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 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 (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 two-chromospheric-footpoint-plus-thermal-coronal-source 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 low-eccentricity (

) 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 short-lived; 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 South-Southwest 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 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 (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 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: . Since the 2002 February 20 event was fairly small (a GOES C-class) 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 super-Dreicer acceleration in current sheets (see, e.g., Somov, 1992; Litvinenko, 1996).


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 re-weigh 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


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 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 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 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 3, we give a sketch of the ASMC algorithm. For more details, we refer the reader to Del Moral et al. (2006).

  Initialization of the sample
  for  do
     draw from
  end for
  Set and
  while  do
     Apply a possible resampling
     for each particle  do
        propose birth/death/split/merge/change moves, then accept/reject
        for each parameter, update the value, then accept/reject
        compute the weights
     end for
     Adaptive determination of the next exponent
  end while
Table 3: ASMC algorithm


  • 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
  • Duval-Poo et al. (2017) Duval-Poo, M. A., Massone, A. M., & Piana, M. 2017, in Sampling Theory and Applications (SampTA), 2017 International Conference on, IEEE, 677–681
  • Duval-Poo et al. (2017) Duval-Poo, M. A., Piana, M., & Massone, A. M. 2017, ArXiv e-prints, 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