On the Application of Differential Evolution to the Analysis of X-Ray Spectra^*

The analysis of X-ray spectra is a powerful tool for plasma diagnostics. Physical parameters of the emitting coronal plasma, such as the values of temperature and emission measure, are not directly observed but can be determined from the modeling of the X-ray spectrum. For many sources the X-ray spectrum is not a unique function of a single temperature but is described using a distribution of plasma with temperature, the so-called differential emission measure distribution (DEM, φ(T)). DEM connects the observed fluxes with atomic data (characterizing the line/continuum emissivity) and the elemental abundance that is usually assumed from existing models. However, the DEM thus has to be determined independently of the elemental abundances. In this contribution we present a new approach in which DEM and elemental abundances are determined simultaneously.

The DEM distribution describes the amount of the optically thin emitting plasma present in a given temperature interval (dT) and volume (dV) and is defined as $\varphi (T)={N}_{e}^{2}\tfrac{{dV}}{{dT}}$, where N_e is the plasma electron density. Knowledge of the DEM distributions is essential for the study of plasma heating, calculation of synthetic spectra, and prediction of the intensity of UV, EUV, and X-ray radiation of the Sun, as well as plasma radiation losses in the solar atmosphere.

The form of DEM can be determined by solving a set of Fredholm's integral equations of the first case, in which the line/spectral band fluxes are defined as

$$\begin{eqnarray}&&{F}_{i}={A}_{i}{\int }_{T}{f}_{i}(T,{\lambda }_{i},{N}_{e})\varphi (T){dT},\end{eqnarray} \tag{ 1 }$$

where F_i is flux in appropriate spectral band/line i, A_i represents the abundance of an element contributing to this flux, and f_i (T, λ_i , N_e ) is the known emission function (atomic physics—CHIANTI in our case) corresponding to the observed flux F_i . The abundance values in Equation (1) are used in linear scale. Further in the paper, the values of elemental abundances are given in the "astronomical" logarithmic scale (where log A_H = 12) for easier comparison.

Equation (1) is an example of an ill-posed problem. The results of the reconstruction of DEM distributions are nonunique and subject to systematic errors (Craig & Brown 1976). However, the methods available in the literature described recently, for example, by Aschwanden et al. (2015), Plowman et al. (2013), Hannah & Kontar (2012), Su et al. (2018), and Kepa et al. (2020), which have previously been tested on synthetic data, provide DEM solutions closely reproducing the observed input spectra as tested on synthetic models. The features present on DEM shape obtained from inversion of real data were interpreted in terms of separate magnetic loop volumes making up solar flare and active region plasmas.

Since spectral line and continuum intensities depend directly on the absolute elemental abundance (see Equation (1)), the results of any DEM inversion strongly depend on the assumed chemical composition of the plasma. Usually, the elemental composition assumed in the inversion relies on earlier spectral observations or plasma composition determinations performed from in situ measurements (solar wind, solar energetic particles (SEPs)). The accuracy and applicability of these abundance determinations count on the calibration of the instruments that provided observations, physical models used in interpretation, and also the data statistics. The uncertainties in elemental abundances strongly influence the results of the inversion process to an extent similar to the uncertainties inherent to spectral fluxes input to inversion. The assumptions used in theoretical derivation of line and continuum fluxes like "Maxwellian plasma in ionization equilibrium" may not be appropriate for some (transient) sources. Spectrally determined abundances depend, in addition, on the accuracy of values for all the atomic parameters involved in spectra modeling like excitation rate coefficients, transition probabilities, etc. (Doschek & Warren 2016).

The coronal elemental abundances (derived from solar wind particle measurements and UV and X-ray emission-line observations) are known to differ from photospheric abundances according to the first ionization potential (FIP) of the element: the so-called "FIP effect." Elements with FIP values below 10 eV (e.g., Si, Mg, Ca, Fe) are enriched by a factor of 2–5 as compared with photospheric abundances (Feldman et al. 1992). Elements with FIP values of more than 10 eV (e.g., C, N, O, Ne, Ar) have abundances that are equal to photospheric values or proxies (such as interplanetary plasma abundances). Earlier work by Meyer (1985) and Feldman et al. (1992) summarizes the FIP effect as it was known then for solar wind, SEPs, cosmic rays, and the solar corona itself. Observations of magnetically active stars having hot coronae (Brinkman et al. 2001; Güdel et al. 2001) performed by the Chandra and XMM-Newton Reflection Grating Spectrometer gave information about the so-called inverse FIP effect, in which the coronal abundances of low-FIP elements are decreased rather than enhanced relative to photospherical ones. The abundance ratios of Ne (highest FIP) to Fe (low FIP) can be order of 10 compared to the photospheric abundances of the Sun (Güdel 2007). An inverse FIP effect was also found for the Sun, for example, by Doschek et al. (2015) and Doschek & Warren (2016) when analyzing the spectra of calcium and argon lines, and recently by To et al. (2021) from Hinode/Extreme-ultraviolet Imaging Spectrometer and by Katsuda et al. (2020) for silicon, calcium, sulfur, and argon abundances determined for four solar flares of GOES X class. In all cases, the inverse FIP effect was only found during flares in small areas of very complex active regions.

RESIK (REntgenovsky Spektrometr s Izognutymi Kristalami; Sylwester et al. 2005a) recorded solar spectra from 2001 to 2003 aboard the Russian Coronas-F spacecraft (Kuznetsov 2014) in four spectral bands covering the nominal range 3.3–6.1 Å. In this spectral range the emission lines of hydrogen- and helium-like ions of such elements as Si, S, Cl, Ar, and K are observed. RESIK measurements constitute a unique database for analysis of the DEM distributions of various solar structures (flares, active areas, quiet corona) and determination of the elemental abundances. The uniqueness of the RESIK data is related to very good intensity calibration of this instrument performed at Rutherford Appleton Laboratory and Mullard Space Science Laboratory, a higher accuracy than for previous spectrometers. The accurate estimation of the level of emission arising from instrumental fluorescence, which is often a problem for X-ray spectrometers, has been performed for RESIK.

DEM distributions calculated based on RESIK spectra are presented in a number of papers available in the literature (Sylwester et al. 2005b; Kepa et al. 2006; Chifor et al. 2007; Dzifčáková et al. 2008; Sylwester et al. 2015a, 2015b). The abundances of the main elements contributing to these spectra (potassium, argon, chlorine, sulfur, and silicon) were determined using two assumptions: the isothermal and multithermal X-ray emission. The results obtained in the first approach are given in Sylwester et al. (2010a, 2010b, 2011, 2012, 2013), while those for the multithermal plasma assumption are given in Sylwester et al. (2014, 2015a). In this second case the abundances of Si, S, Ar, and K were optimized using the maximum likelihood, Bayesian Withbroe–Sylwester (WS) inversion technique (Sylwester et al. 1980).

In Table 1 we summarize values of abundances determined using RESIK spectra and compare them with the values available in the literature.

Values of calculated abundances for argon and potassium values in the case of isothermal and multithermal assumptions based on RESIK observations are in agreement. For isothermal assumption the estimated sulfur abundance depends on the lines used in analysis and lies in the range from 7.13 to 7.24, with the most reliable value being 7.16 (estimated based on the S xv w4 line). These values are very close to photospheric abundance estimates (Asplund et al. 2009; Caffau et al. 2011) and to quiet-Sun solar wind (Reames 1999) and meteoritic abundances (Lodders 2003). The value of silicon abundance depends on the spectral lines used for calculations. It differs by a factor of 2.6 and 2.4 from the photospheric abundance (Asplund et al. 2009), respectively. The determined silicon and sulfur abundances are higher than those obtained with the multithermal assumption. However, abundances obtained with the multitemperature approach should be preferred owing to their smaller uncertainties. In the case of sulfur, this lower abundance indicates that, in contrast to the view expressed by Feldman et al. (1992), sulfur behaves like a high-FIP element (the flare S abundance being even lower than in the photosphere). Chlorine abundance was determined with high error, which reflects a factor of two uncertainty in measured line fluxes.

In this contribution we propose a new approach for determination of elemental abundances from EUV and X-ray spectra. The new method solves the inversion simultaneously searching for the optimum combination of DEM distribution and the set of abundance values for elements important in the RESIK spectra formation. In the described approach a differential evolution (DE) method (Storn & Price 1997; Price et al. 2005) is used. DE is a powerful technique for difficult optimization problems and has been used in dozens of studies, e.g., for the identification of stellar flares (Lawson et al. 2019), stellar classification (Contreras Sanchez et al. 2018), and the optimization of neutrino oscillation parameters (Mustafa et al. 2013). We have used this method for the analysis of RESIK spectra. The DE method has already been tested and used for determinations of DEM distributions of flaring plasma (Kepa et al. 2016, 2020). Here, this approach is applied to calculations of the elemental abundances for the first time. This novel analysis was performed for the example, well-observed RESIK flare that occurred on 2003 January 21 (SOL2003-01-21T15:26, M1.9 GOES class).

This contribution is organized as follows: A description of the investigated flare and a short introduction to the DE method are given in Section 2. The algorithm test and the results of analysis are shown in Section 3. The conclusions and discussion are provided in Section 4.

We studied the M1.9 GOES class flare that occurred on the east limb of the solar disk (S07E90) on 2003 January 21 and reached maximum in X-ray at ∼ 15:26 UT (SOL2003-01-21T15:26). RESIK observations for this event are available from 14:56:35 UT to 18:56:34 UT with a few gaps when the Coronas-F spacecraft crossed through the Van Allen belts or the South Atlantic Anomaly. The RESIK spectrometer recorded 289 high-quality spectra with exposition from 4 to 905 s, and the total time for spectra collection was about 1.5 hr.

The average RESIK spectrum for the chosen flare is shown in Figure 1. The gray strips indicate 25 wavelength bands used for analysis. The information about emission lines observed by RESIK and about wavelength ranges selected is presented in Table 2. The different colors of the spectra represent four spectral channels (1: 3.3–3.8 Å, orange; 2: 3.8–4.3 Å, red; 3: 4.3–4.9 Å, green; 4: 5.0–6.1 Å, blue) in which RESIK spectra were recorded. RESIK observations include the 1s²–1s n p and 1s–n p lines of He-like and H-like ions and such elements as silicon, sulfur, argon, and potassium. For some flares the chlorine lines were visible. The lines of argon and sulfur were observed in the first and second channels. The potassium triplet (3.53 Å) was observed in the first channel, and the chlorine line was noticed in the second channel (Sylwester et al. 2011). In the third and fourth channels the lines of sulfur and silicon were the most pronounced.

Zoom In Zoom Out Reset image size

For the present analysis the fluxes containing lines and the continuum were used, and the corresponding emission functions have been determined using CHIANTI v8.0.7 (Dere et al. 1997; Young et al. 2019) and are described in more detail in Section 3. The spectral line Si xiv 1s ² S_1/2−3p² P_1/2,3/2 observed around 5.2 Å (in channel 4) was excluded from the present analysis, as the intensity of this line is observed to be systematically too strong—for unknown reasons, probably related to partly nonthermal excitation of this line.

Here we introduce how the DE approach is adapted in order to be used for simultaneous determination of DEM distributions and elemental abundances.

DE belongs to the class of evolutionary algorithms designed to solve optimization problems with real-valued parameters. It is a very powerful method for the black-box optimization, where the function to be optimized is very complex and it is difficult to find the appropriate derivatives necessary to compute its extreme. DE is based on the biological evolutionary processes: mutation, crossover, and selection.

The values of three parameters are important to achieve the solution: a size of the population, and mutation (F) and crossover (CR) factors.

The most important characteristics of DE are the following:

• 1.  
  The mutation is provided by arithmetic combinations of individuals, not as the result of small perturbations to the genes as in genetic algorithms (Feoktistov & Janaqi 2004).
• 2.  
  DE works with three equally large populations of individuals: the population of parents S, trial, and descendants.
• 3.  
  The size of populations (a set of potential solutions)—usually does not change during the evolution process.
• 4.  
  The DE works with floating point numbers, so the crossover/mutation is easier.

A typical working scheme for DE is presented in Figure 2. At the beginning the population (S) has to be initialized. In our case, each solution (individual) s_i , where i = 1, 2, 3, ..., S, represents the values of elemental abundances or DEM. The individuals consist of the genes. This means that a few genes are related to a chemical composition (one gene corresponds to the abundance for one element) and a number of genes correspond to the values of DEM distribution. We use two models of the individuals with different representation of DEM distribution. In the first one (Figure 3(a)) genes represent the abundances, temperatures, and DEM values. In the second model (Figure 3(b)) the genes represent abundances and DEM values only. In this approach we set (define) values of temperatures at which DEM distribution will be determined.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each parameter (gene; one box in Figure 3) the limiting values have to be present and the initial population should cover the entire search space as much as possible. In our specific case, each (s_i ) solution (corresponding to one DEM distribution and one set of abundances of elements that dominate the spectra formation) is a randomly generated vector.

Next, the population is subjected to mutations. In this process, for each individual in the population S (S is often referred to as the parent population) the trial vector is created. It is important that for each individual exactly one individual of target population U is created. In the basic version of mutation the target population is created according to the following formula:

$$\begin{eqnarray}&&{\forall }_{i,j}{U}_{i,j}={S}_{i,j}+M\times ({S}_{{r}_{1},j}-{S}_{{r}_{2},j}),\end{eqnarray} \tag{ 2 }$$

where S_i,j is the jth element of the ith individual, M is the mutation factor, and ${S}_{{r}_{1},j},{S}_{{r}_{2},j}$ are the jth elements of the other vectors r₁ and r₂. The r₁, r₂ indexes are randomly chosen and different from i (r₁ ≠ r₂). The value of the M factor affects the exploration speed of the space solution; the greater value corresponds to the higher speed (Georgioudakis & Plevris 2020).

Besides the mutation scheme according to Equation (2), there are many other combinations available in the literature (e.g., Das et al. 2016). However, in this paper we restricted ourselves to using the basic version of mutation.

In order to increase the diversity of the perturbed parameter vectors, the population of target vectors (U) is subject to crossing. The crossover builds trial vectors V_i through the recombination process from the elements of the parent vector, S_i , and the elements of the mutant vector, U_i in the following way:

$$\begin{eqnarray}{\forall }_{i,j}{V}_{i,j}=\left\{\begin{array}{ll}{U}_{i,j} & \mathrm{if}\,{\mathrm{rand}}_{j}(0,1)\leqslant \mathrm{CR},\\ {S}_{i,j} & \mathrm{otherwise.}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

Here the crossover probability, CR ∈ [0, 1], is user defined and rand_j (0, 1) is a random number from the specified range.

To decide whether or not a given individual should become a member of the next generation, the trial vector V_i is compared to its direct parent S_i . If V_i has a lower or equal objective function value than that of its target vector, S_i , it replaces the target vector in the next generation:

$$\begin{eqnarray}{\forall }_{i}{S}_{i}^{k+1}=\left\{\begin{array}{ll}{V}_{i}^{k} & {\rm{if}}\,\,\,f({V}_{i}^{k})\geqslant f({S}_{i}^{k}),\\ {S}_{i}^{k} & {\rm{otherwise.}}\end{array}\right.\end{eqnarray} \tag{ 4 }$$

Here ${S}_{i}^{k+1}$ is the ith individual in the k + 1 generation and f is the objective function.

Our objective function is given as the χ² value, defined as

$$\begin{eqnarray}&&f={\chi }^{2}=\sum _{i=1}^{L}\displaystyle \frac{{({F}_{i}-{F}_{{ic}})}^{2}}{{\delta }_{i}^{2}},\end{eqnarray} \tag{ 5 }$$

where F_i is the observed flux in the ith spectral band from Table 2, F_ic is the flux in the ith line/spectral band calculated using Equation (1) (based on φ(T) distribution and abundance values taken from the individual), δ_i is interpreted as an error consisting of the uncertainties in the measurements, and L is the number of lines (or bands) used in the analysis (L = 25 in our case). Our aim is to minimize the χ² value. The evolutionary process is repeated until the convergence criterion is fulfilled.

The choice of the optimal value of parameters for the DE method is a challenging task (Storn 2008). In the literature there exist many proposals of self-adaptive methods that select the values for M and CR (Milani 2020). The parameters used in our study were adopted based on the tests performed on synthetic models. In our calculations we set mutation and crossover factors to be 0.1. The examples of the tests using DE for simultaneous calculation of DEM distributions and plasma composition are presented in the next section.

DEM distribution is derived by deconvolving φ(T) from the fluxes observed in selected spectral ranges (Equation (1)) and known from theoretical emission functions f_i (T, λ_i , N_e ). In order to be able to interpret the results of DE inversion, we performed numerous tests of spectra reconstruction for spectra generated for model sources. Testing provides insight into the reliability of the DE approach for simultaneous determination of element abundances and DEM distribution.

We have used two types of source models, the isothermal and the two-temperature, and several sets of elemental abundances. In this respect we assumed the "coronal extended" (as in CHIANTI, after Feldman et al. 1992) for most of the elements. The abundances for silicon, sulfur, argon, and potassium were drawn from the values within specified limits in various exercises. The generated X-ray spectra were perturbed within limits given by errors coming from expected count statistics. Next, these spectra have been treated as observations, which we should reconstruct using the DE method. In the analysis we used the fluxes and corresponding (integral kernel) emission functions determined for 25 spectral bands presented in Table 2 that contain the contribution of emission in line and continuum. The emission functions have been calculated using CHIANTI package ver. 8.0.7 available in SolarSoft for ionization equilibrium based on Bryans et al. (2009) and unit elemental abundances. Thanks to this, the change in the elemental abundances during each step of iteration of the DE method was associated with the multiplication of emission functions by the appropriate values only. The ranges of values in which parameters could change were assumed to be T = 1–30 MK, log EM = 40–50, A_Si = 7–9, A_S = 6–8, A_K = 5–7, and A_Ar = 5–8. Abundances were expressed logarithmically on a scale with A_H = 12. We started each test by generating a "population" of 500 trial solutions (S = 500), where each of the individuals has been defined according to the scheme presented in Figure 3(a). The values for genes were random and placed within assumed ranges.

In Figure 4 we illustrate the process of convergence between reconstructed and generated synthetic spectra at six indicated iteration steps during DE progress. For each iteration step we present the "best-fit" calculated spectrum and plot it over the synthetic "observed" one. The spectrum assumed as "observed" has been calculated for the isothermal plasma at temperature T = 9 MK and emission measure EM = 2 × 10⁴⁸ cm⁻³. We assumed the abundances of silicon, sulfur, argon, and potassium to be A_Si = 8.3, A_S = 7.12, A_Ar = 6.82, and A_K = 5.41, respectively. The generated input synthetic spectrum was randomly perturbed by ∼10% fluctuation in spectral irradiance.

Zoom In Zoom Out Reset image size

The analysis of the presented test and many other tests made with the isothermal plasmas reveals that the DE method satisfactorily reconstructs the assumed plasma parameters after 500 steps of iteration. The obtained results are stable and do not change if the iterations continue (see the results presented for the 700th and 900th steps). A similar pattern of convergence is observed for different sets of temperature and emission measure pairs, as well as for different sets of abundances.

In Figure 5 we present the sequence of spectra reconstruction obtained when the input synthetic spectrum was calculated for the two-temperature model of the source. In the model we assumed a very small amount of high-temperature plasma EM₂ = 2 × 10⁴⁵ cm⁻³ at T₂ = 12 MK "contaminating" the main cooler component of EM₁ = 1 × 10⁴⁹ cm⁻³ at T₁ =5 MK. The elemental abundances of silicon, sulfur, argon, and potassium were taken to be A_Si = 7.70, A_S = 7.10, A_K = 5.91, and A_Ar = 6.50. In this case our population consists of 500 individuals (possible solutions) and the calculation process was interrupted after 5000 steps.

Zoom In Zoom Out Reset image size

The results of the tests are similar to those obtained from the reconstruction of spectra calculated for the isothermal model of plasma. The spectrum has been reconstructed properly, and the obtained solution is in agreement with the assumed input parameters of the model. However, achievement of a good solution needs more iteration steps. In the presented case it was over 1300. The calculations take a few minutes on a personal computer.

In Figures 4 and 5 we show the best solutions obtained in selected iteration steps. However, each individual in the population is changed by mutations, crossover, and selection during the iteration. In Figures 6 and 7 we show how the solutions differentiate within the whole population at four selected iteration steps for isothermal and two-temperature models of plasma, respectively. Each parameter (temperature, emission measure, and silicon, sulfur, argon and potassium abundances values) is presented in a separate panel. Different colors represent different iteration steps.

Zoom In Zoom Out Reset image size

The differentiation of the whole population during reconstruction of the spectrum obtained with the isothermal model of plasma (from Figure 4) is presented (in Figure 6) for steps 5 (in black), 200 (in blue), 500 (in yellow), and 1000 (in red). As can be seen after 1000 steps, the entire population becomes almost homogeneous (all 500 solutions are similar), but 200 iterations are enough to find a proper solution. It is observed that the temperature structure (i.e., T and EM values) converges the fastest. The potassium abundance value is established as the last parameter owing to a very weak A_K and a weak contribution of potassium lines to the single-band flux (the H-like K xviii triplet at 3.53–3.57 Å contributing to bands 3 and 4). Potassium lines extend only slightly above the continuum, contrary to the other bands, where strong lines of Ar, S, and Si dominate the fluxes in bands.

The similar plots that present the results of spectrum reconstruction for a two-temperature model from Figure 5 are shown in Figure 7. Here, the state of "population differentiation" of solutions is presented as obtained in iteration steps 5, 500, 3000, and 5000.

Zoom In Zoom Out Reset image size

One can see that for the two-temperature source the population becomes homogeneous after a longer evolutionary process than for the isothermal model. It is clearly seen that even in the last iteration shown (5000) a few percent of solutions still are different from assumed input parameters. This reflects the more complicated two-temperature "nature" of the source. The population differentiation is slightly dependent, however, on the random draw of the initial population for the DE evolution. For various populations a satisfactory solution is obtained in different iteration steps. The process of evolution was stopped after 5000 generations because the convergence became very slow. The minimum χ² in each case was below 1.

We checked also the possibility of simultaneous calculation of elemental abundances and DEM distribution for the multitemperature approach. We have calculated DEM for 30 temperatures (30T) with ΔT = 1 [MK] in the range 1–30 [MK] (using the DE method). In this case, the population consisted of individuals created according to the scheme from Figure 3(b). The temperatures in the DE method were predefined, and the sought variables were the abundances and emission measures in these 30 temperatures. The obtained results of reconstruction for three different models of plasma distribution are given in Figure 8.

Zoom In Zoom Out Reset image size

Based on the results of the assumed model reconstructions using the multithermal approach, we can conclude that DE constitutes a viable alternative method not only for determination of DEM distributions but also for inferring values of elemental abundances from well-calibrated X-ray spectra. All unknowns (DEM and composition) are being reconstructed simultaneously, but at various paces, reproducing satisfactorily not only the shape of the DEM but also the set of input elemental abundances. The method correctly reproduces the saturation values of silicon, sulfur, and argon. However, the value of the abundance of potassium is burdened with a very large error owing to the sensitivity limit of this method (in the case of multitemperature plasma).

After testing on the synthetic models, the DE method was used for the study of the flare that occurred on 2003 January 21 with maximum in X-ray at 15:26 UT. In Figure 9, we show the temporal evolution of the soft X-ray fluxes in GOES bands 1–8 Å and 0.5–4 Å. All time intervals when RESIK operated and recorded spectra are marked by gray strips. The temporal evolutions of isothermal temperature and emission measure as calculated based on GOES data are presented too. The values of temperature and emission measure were derived using the filter-ratio technique (with the flare background level subtracted) using the standard routine goes.pro from SolarSoft package. The determined maximum values of temperature and emission measure were 14.9 MK (at 15:04 UT) and 1.37 × 10⁴⁹ cm⁻³ (about 15:29 UT), respectively.

Zoom In Zoom Out Reset image size

For further analysis we selected spectra in four time intervals. They are marked as numbers at the top of Figure 9 (the detailed timing is given in the successive Figures 10 and 11). For each of the selected time intervals, we calculated the mean RESIK spectrum and fluxes in the set of 25 wavelength ranges (see Table 2). These fluxes were the input data used for the reconstruction of the DEM distribution and elemental abundance determination. The DEM distribution and elemental abundances were calculated simultaneously using the DE method.

At the beginning we used the isothermal model of plasma (values of temperatures, emission measures, and abundances for silicon, sulfur, potassium, and argon) that reproduces the best RESIK fluxes in 25 spectral bands (Table 2) in each of four chosen time intervals. The results are plotted in the left panel of Figure 10 (in red). Based on these values, the synthetic spectrum was calculated. The comparison of our best-fit calculated spectrum with observations (in black) is given in the right panel. The same procedure has been adopted assuming the two-temperature model. The blue color is related to the two-temperature model of plasma. The results of calculations made based on unperturbed observed fluxes are given as dark-red and blue dots. We have also checked how the data errors influence the solutions. We have randomly perturbed fluxes (10% error) and determined plasma parameters. The light red and blue represent obtained solutions for 100 independent generations.

Zoom In Zoom Out Reset image size

One can perceive that for the rise phase of the flare (time intervals 1 and 2) the temperature that describes the observation in the isothermal model the best is between the two temperatures derived for the two-temperature model. For the observations during the decay phase (time intervals 3 and 4) the values for the isothermal model of plasma are closer to the values of the cooler component in the two-temperature model of the emitting plasma. A scatter of solutions obtained based on perturbed fluxes for isothermal models is smaller than for the two-temperature model. The values of total emission measure determined for two different models are similar, and the small discrepancy is related to the DE method.

Next, based on the obtained plasma models from the left panels (temperature, emission measures, and abundances), we calculated the synthetic spectra and compared them to observed ones (right panel of Figure 10). The observations are plotted in black, and the red and blue colors correspond to the spectra determined from isothermal and two-temperature models of emitting plasma, respectively. One can see that the spectra calculated based on the two types of models are almost identical. However, two differences are clearly visible. The first is associated with the continuum level in the shortest wavelengths (3.3–3.6 Å), the second one with a fit of S xvi observed at 4.729 Å in wavelength band No. 18 (see Table 2). In purple (thin line) the synthetic spectra obtained based on the analysis of GOES fluxes are plotted. The determined values of temperature and emission measure are presented in Table 3. One can see that spectra calculated based on the isothermal model of plasma from GOES data and coronal abundances by Feldman et al. (1992) (Figure 10, in gray) do not agree with RESIK-observed spectra. They are plotted for comparison only. We focused on the similarities and differences of the spectra calculated for the isothermal and two-temperature models of plasma from RESIK fluxes.

The results presented for the 30-temperature (30T) model of plasma (individuals created according to the scheme from Figure 3(b)) are shown in Figure 11. For each time interval 10 DEM distributions from 10 independent generations of the population are shown (left panel, in light green). The average DEM distributions are presented in dark green. The results of the analysis using the two-temperature model of plasma are shown in blue for comparison. One can see that these results (2 and 30 temperatures) are not in contradiction. Two DEM components obtained using 30T models are in agreement with temperature and emission measure values obtained from the two-temperature models.

Zoom In Zoom Out Reset image size

Based on each of 10 independent DEM distributions revealed by the 30T DE exercise, we calculated the respective set of 10 synthetic spectra. These spectra are shown in the right panel of Figure 11 (thin light-green lines). The shapes are so similar that they converge into the slightly wider light-green line falling atop the observed RESIK spectra (in black). It is clearly seen that spectra calculated in the 30T exercise do not reproduce RESIK observations perfectly. The discrepancy between observed and calculated spectra that is seen in the wavelength range from 3.3 to 3.6 Å may be due to some unknown instrumental effects. The dissimilarity in spectral range 5–6 Å may be associated with the unaccounted-for overlapping contribution from higher orders of crystal reflection. Work on this problem (discrepancies) is in progress.

As noted many times, in the DE analysis presented here the DEM distributions and elemental abundances were fitted simultaneously. The best-fit characteristics (temperature and emission measure) for single- and two-temperature approaches (unperturbed fluxes) are presented in Table 3, together with values obtained from the analysis of GOES fluxes. The values presented in parentheses represent limits obtained from perturbed data sets. In the case of the 30T temperature model we calculate the values of average temperature (T_av) and total emission measure (EM_tot) for each of the 10 independent DE runs using the following formula: ${T}_{\mathrm{av}}=\tfrac{{\int }_{T}T\varphi (T){dT}}{{\int }_{T}\varphi (T){dT}}$, EM_tot = ∫_T φ(T)dT. The indications "1av, 1tot" and "2av, 2tot" are related to temperature ranges 1–10 MK and 10–30 MK, respectively.

One can see that the temperature and emission measure values obtained from RESIK and GOES in the isothermal approximation differ substantially. The most probable reason is that GOES determinations (goes.pro, SSW package) use the standard "coronal extended" set of abundances while we allow for abundances to vary, which directly influences the results, and the values of emission measure the most.

The values of the average temperatures and the total emission measure (30T model) are very similar to those obtained from 2T models. It is not a surprising conclusion and reflects the inner consistency of the DE approach.

The best-fit values of abundance for the four elements considered (Si, S, Ar, and K) calculated for isothermal (1T), two-temperature (2T), and 30-temperature (30T) models of plasma are given in Table 4. For comparison the values of RESIK abundances obtained in previous analysis by Sylwester et al. (2015a) and Feldman et al. (1992) are given.

The graphical representation results given in Tables 3 and 4 are shown in Figure 12. It is clearly seen that obtained values of abundances are similar and lie within the bounds of errors. They also do not change during flare evolution. The exception is the solution obtained for time interval 2 (15:13:45 UT–15:20:29 UT). Such behavior may be related to nonthermal effects observed during the rise phase of the flare. Unfortunately, there are no RHESSI observation at these times. The obtained values of elemental abundances for silicon, sulfur, and argon are in agreement with those obtained by Sylwester et al. (2015a) using the WS approach for DEM calculations and the AbuOpt method for determinations of elemental abundances. However, potassium abundance could not be determined in our method. This is probably the consequence of the limit of sensitivity of the method and of the before-mentioned problem about disagreement between observed and calculated (based on observed continuum in shortest wavelengths RESIK channel, 3.3–3.6 Å). In addition, the fact that the potassium lines are barely observed only in two spectral bands may cause the problem.

Zoom In Zoom Out Reset image size

One of the fundamental assumptions when analyzing X-ray spectra with the aim to calculate DEM distribution is related to the elemental abundances. The method discussed in the literature relies on the adoption of the elemental abundance set before calculations of DEM distributions. In such an approach the elemental abundance errors translate into inaccuracies in the determined DEM distribution.

In this research, we introduce the DE approach that resolves this troublesome problem by simultaneously iteratively improving elemental composition and the DEM shape based on the evolutionary selection process. The DE algorithm was extensively tested, and the results of the tests confirmed the stability of the solutions and the capability of the algorithm to simultaneously reconstruct the DEM distributions and elemental abundances.

The DE method was applied for the analysis of RESIK spectra obtained for the 2003 January 21 flare. Our aim was to find the plasma model (abundances and distribution of plasma with temperature) that fits RESIK observations the best. We calculated the isothermal model, two-temperature model, and DEM distribution described in 30 temperatures (from 1 to 30 MK). Our analysis can be summarized as follows:

• 1.  
  The values of temperature and emission measure calculated based on GOES fluxes using the isothermal model of plasma do not allow us to reproduce RESIK spectra. This (most probably) indicates that the coronal set of abundances used by the GOES SSW analysis package is not generally applicable for detailed analysis of every coronal source.
• 2.  
  The two-component model of plasma better describes RESIK observations than the isothermal model (as expected).
• 3.  
  The values of the average temperatures and the total emission measure for multithermal analysis (30T model) are very similar to those obtained for the 2T model.
• 4.  
  The values of abundance for silicon, sulfur, and argon do not change during flare evolution.
• 5.  
  Best-fit DE abundance values are, within uncertainties, the same as those obtained by Sylwester et al. (2015a) for the same flare using a different approach.
• 6.  
  It was not possible to determine the potassium abundance based on the RESIK spectra used in this study (too small potassium line contribution).

We have shown how the DE algorithm can be used to make a breakthrough advance in analysis of spectra from optically thin plasmas on the example of spectra recorded by RESIK. This approach can be adopted to analysis of other data sets available from the many instruments observing the Sun (e.g., the Bent Crystal Spectrometer aboard the Solar Maximum Mission spacecraft) and other stellar objects in the EUV and X-rays in cases in which the source plasma composition is unknown.

When we finished the calculations a new version of CHIANTI (10.0.1) became available. This latest version of CHIANTI includes substantial changes as compared with the previous ones in the line intensities, as well as in the X-ray range. We found these changes for some Si line fluxes to be as large as 10%. In our future work, we will make use of the most recent version of CHIANTI.

We acknowledge financial support from the Polish National Science Centre grant No. 2017/25/B/ST9/01821. We are thankful to an anonymous referee for providing us with very useful feedback.

Software: CHIANTI (v8.0.7; Dere et al. 1997; Young et al. 2019), SolarSoft (Freeland & Handy 2012).