© The Authors 2024

1. Introduction

The X-ray telescope capability of focusing protons with energies lower than a few hundred kiloelectron volts, so-called soft protons, was discovered just after the launch of the Chandra X-ray observatory (Weisskopf et al. 2002), when a rapid degradation of the front illuminated CCDs at the focal plane of the Wolter type-1 telescope occurred (Lo & Srour 2003). The damages were in fact caused by soft protons encountered along the orbit that, hitting the X-Ray optics surface at grazing angles, were transmitted to the focal plane.

XMM-Newton (Jansen et al. 2001) has been exposed to the same risk, carrying three Wolter type-1 mirrors to focus X-rays into the focal plane instrument EPIC (European Photon Imaging Cameras), consisting of three cameras, covering the 0.2–15 keV energy range: two MOS (Turner et al. 2001) and one PN CCDs (Strüder et al. 2001). The focal plane detectors are protected from the background induced by IR, visible, and UV light with three optical blocking filters with different thicknesses of both the aluminum film and its support: the thin, medium, and thick filters. They are mounted on a filter wheel and can be changed depending on the requirement on the observations. In addition, the two MOS telescopes are equipped with reflection grating spectrometer (RGS, den Herder et al. 2004) nested behind the mirrors. These instruments, operated simultaneously with EPIC, allow for high resolution (E/ΔE ranging from 100 to 500) measurements in the soft X-ray range of 0.35–2.5 keV.

The soft protons transmitted by the optics, after crossing the optical filters, release a part of their energy in the EPIC CCDs, generating a local charge indistinguishable from the signal produced by photons. In the case of a sudden increase in the input flux, the final effect is the generation of flaring events, lasting from hundreds of seconds to hours, which degrade the quality of science observations because of the increased background, resulting in an average loss of 30–40% of exposure (Marelli et al. 2017).

However, what is a disadvantage for the observation of astronomical sources would be a good opportunity to study the environment crossed by the satellite along its orbit, provided that a proton energy-transfer function is available. As is explained in detail in Fioretti et al. (2024, hereafter Paper I), this function in the format of a response matrix has been built from detailed Geant4 simulations modeling all the proton interactions with the telescope elements: the transmission from optics, the energy degradation from the filters, and the total release of the soft proton energy in the detectors. The RGS was not included in the simulator. However, considering that it shades ∼50% of the photon flux in the EPIC/MOS (den Herder et al. 2004), the same attenuation factor was assumed for protons and the total effective area of this instrument was halved independently of the energy (Paper I). This assumption is confirmed by the simulations performed by Nartallo (2002) to quantify the level of soft protons focused by the optics to XMM-Newton and Chandra focal planes.

In summary, three matrices for the MOS and three for the PN, one for each filter, have been generated. They are relative to a single event pattern with no spatial selection at the focal plane. The working energy range is 2–11.5 keV for the MOS and 2–19.5 for the PN, with a systematic uncertainty on the absolute flux of ∼20%. As is explained in Paper I, a first validation of the matrices has been performed using representative MOS spectra obtained from the EXTraS¹ (Exploring the X-ray Transient and variable Sky) archive of focused background spectra (Salvetti et al. 2017), accumulated in the first 13 years of XMM-Newton life.

The proton response matrices² have been produced to be compliant with XSPEC (Dorman & Arnaud 2001), an interactive X-ray spectral-fitting program, and are given as two separate files: i) the ancillary response file that stores the summed contributions of all efficiencies; and ii) the redistribution matrix that contains the probability that an incoming proton of energy, E, is detected in the output detector channel, PHA. The effective area relative to MOS and PN for the thin and medium filter are shown in Fig. 1. The two peaks in the MOS curves are generated by the complex structure of the open electrode with two regions of different thickness (Paper I).

Fig. 1. Proton response matrix effective area for MOS (top panel) and PN (bottom panel). The application of the thin (medium) filter is represented by a continuous (dash-dotted) curve.

Fig. 2. Light curves relative to observation ID 0844860601 (2020-01-06) for MOS1, MOS2, and PN. The dashed vertical lines delimit the interval in which the background rate was evaluated, the dash-dotted lines indicate the highest background rate, and the dotted lines indicate the lowest flare rate.

It is important to note that our matrices are a first example of response matrices for soft protons and have never been applied to observational data. This paper, then, represents the first example of spectral analysis of flares detected by XMM-Newton. Here, presenting our work and results in fitting each spectrum independently, we discuss with particular attention the intercalibration between the MOS and PN matrices under the consideration that the input spectra measured by the three detectors must be statistically equivalent. The data reduction and a model-independent analysis of the considered observations are presented in Sect. 2 and Sect. 3, respectively. The spectral analysis is described in Sect. 4; the results and systematics are shown in Sect. 5 and discussed in Sect. 6. Conclusions are drawn in Sect. 7.

2. Observations and data reduction

Our analysis considered as the signal the events collected in the intervals with flaring activities and extracted the background from quiescent emission intervals. XMM-Newton data were obtained from the archive³, selecting observations with no strong sources in the field of view and full frame mode available either for MOS and PN, where the light curves and the focal plane images showed an evident presence of flares. These criteria were applied for data selection in four epochs: December 2001–January 2002 (named epoch A and identified by the red color throughout the paper), July–August 2002 (epoch B, blue color), December 2019–January 2020 (epoch C, green color), and July–August 2020 (epoch D, gray color). They involve different environmental conditions due to seasons or to solar cycle phases, the years 2001–2002 being at the maximum of solar cycle 23 and 2019–2020 at the minimum of solar cycle 24. This choice would allow us to investigate any eventual effect on flares, as the dependence of its density on the angle between the satellite orbit and the magnetosphere axis (Fioretti et al. 2016). Table 1 and Table 2 present the log of the selected observations at the maximum and minimum solar activity, respectively.

All data were reprocessed using the standard tools, epproc and emproc (SAS ver. 21.0.0). Light curves and spectra (both from flares and quiescent emission) were accumulated by selecting only single events (PATTERN = 0) to be compliant with the available proton matrices. Bad pixels (XMMEA_EM and XMMEA_EP set to on and FLAG = 0) were excluded and no spatial selection was applied, again in agreement with the matrices that are relative to the whole field of view.

Flare spectra were accumulated, filtering EPIC events for periods of high background flaring activity⁴, while intervals with steady emission constitute our background. This includes the sources in the field of view, the diffuse X-ray sky, and the instrumental dark components. In particular, in the 50 s binned light curves of each observation, we first identified by eye periods with low steady rates (see the interval between the two vertical dashed lines in Fig. 2 as an example), where we computed the average rate (ave) and the root mean square (rms) relative to it. These were used to produce good time intervals (GTIs) to accumulate spectra. Background GTIs were created by selecting intervals with rates lower than two times the rms above the average [ave + 2 × rms], where up to ∼95% of the background fluctuations are included. Flares GTIs instead are relative to periods where background fluctuations have a probability lower than ∼0.3% of contaminating the signal [rate > ave + 3 × rms]. The fraction of events above or below the two thresholds in units of rms assumes a Gaussian distribution of the count rates. Figure 2 shows the background interval for observation ID 0844860601 with two dashed vertical lines, and the rate thresholds for the background and for the flares with dash-dotted and dotted lines, respectively.

The average spectrum per observation was accumulated without taking into account eventual spectral variations in the observations. The spectral analysis was performed on background-subtracted spectra with channels grouped in order to have a minimum of 50 counts per bin, allowing for the application of the χ² statistic.

3. Model-independent analysis

A model-independent analysis on the selected observational data is required before applying the response matrices, considering that our sample spans 20 years of satellite life. This was performed using the 2–11.5 keV flare rates measured by the three instruments in the four epochs after subtracting the contemporary background (see Table 1 and Table 2).

It is important to note that data relative to MOS1 CCD3 and CCD6 are not available after 2005, due to damage from micrometeorite impacts. The main effects expected from losing the two detectors should be an average reduction of MOS1 rates on the order of ∼30%. However, looking at the rates listed in the two tables, a more complex behavior between different epochs appears.

We investigated, at first, the variation in the count rate between the three cameras in each epoch. To this purpose, the averages plus errors and the dispersion relative to them evaluated from the standard deviation (rms) were computed. From the resulting values, presented in Table 3, it is evident that there are no statistically significant variations in each instrument rate during epochs because of the large rms.

However, the comparison of the average values shows some inconsistency between instruments: PN average rates in the interval C and D are significantly higher with respect to those in A and B, while the MOS2 ones are lower, even if no variation in this detector response is expected. In addition, the average rates observed by MOS1, where a decrease of ∼30% is expected, are almost constant, even if values in A and B are higher than those in C and D. A more precise evaluation of these differences was computed from the ratios between rates in the three instruments. The corespondent average values, listed in Table 4, confirm the inconsistency.

As a further analysis, the correlations between the rates detected by the three instruments were investigated, plotting the values relative to MOS1 as a function of the corresponding MOS2 ones and the rates from MOS1 and MOS2 against the PN ones (see Fig. 3). Each of the curves was then fit with a line fixing the intercept to zero (y = k * x). The fit used the least squares method for the minimization because the statistical errors are negligible with respect to the systematic ones. The values of the resulting slopes, k, are shown in Table 5; they are in agreement with the values in Table 4, confirming again the inconsistencies detected between the rates in epochs C and D and those in A and B.

Fig. 3. Correlations between rates detected by the three instruments in the four epochs. Colors indicate epochs (see text). The superimposed dashed black line is relative to the best fit of epoch B. The top panel is relative to MOS1 vs. MOS2 rates, the middle panel to MOS1 vs. PN, and the bottom panel to MOS2 vs. PN.

As a check, correlations between rates detected by instruments were also investigated for photons using the background intervals identified in our analysis. We computed these rates in the range of 2–5 keV to avoid contamination from the internal background component, which is dominant above this range⁵ and which is not affected by changes in the effective area. The results are shown in Table 5. In this case, the correlation factors are coherent with what was expected: the slope of the linear fit of MOS1 vs. MOS2 rates in epochs C and D decreases by about 30% with respect to the one in epochs A and B, while no significant variation is observed for MOS2 vs. PN rates.

To complete the model-independent analysis, an average hardness ratio was computed to investigate spectral differences among observations. The hardness was computed using counts in the ranges 5–11.5 keV and 2–5 keV. Values as a function of the 2–11.5 keV rates for the three instruments are shown in the three panels of Fig. 4. Fitting all the points with a line, no correlations between the hardness values and rates were obtained for both MOS1 and PN being the best fit value of the slope compatible with zero. The MOS2 rate instead shows a weak positive correlation with hardness (slope = 0.010 ± 0.008).

Fig. 4. Hardness ratio between counts in the ranges 5–11.5 keV and 2–5 keV vs. 2–11.5 keV rate for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel). Colors indicate the four epochs (see text).

4. Flare spectral analysis

A large collection of XMM-Newton EPIC blank sky observations for the years 2000–2012 was provided by EXTraS, allowing one to characterize the focused background and to evaluate the fractions of time that the satellite spent in the regions of the Earth magnetosphere crossed by its orbit (Ghizzardi et al. 2017). However, information on the input spectral model expected for the soft protons detected by XMM-Newton is still not accurate, despite the detailed studies performed on the subject (Kronberg et al. 2020; Gastaldello et al. 2017; Kuntz & Snowden 2008).

The ESA AREMBES project (Macculi et al. 2021) predicted the soft proton environment for the Athena mission at the second Lagrange point. Lotti et al. (2018) and Fioretti et al. (2018) derived power law distributions above 50 keV for different regions of the magnetosphere and interplanetary space. These models can also be considered valid for the XMM-Newton environment. Paper I presents the spectral analysis, using the proton response matrices, of EPIC MOS2 averaged spectra obtained from the EXTraS archive, normalized for the maximum intensity expected in 90% of background observations (Molendi et al. 2017). The best-fitting model is in agreement with the predicted spectrum of the proton environment encountered by XMM-Newton’s orbit, predicted for the same 90% exposure fraction.

From our analysis, performed on single flares, we first started using power law models that describe well the proton environment above 50 keV. However, because XMM-Newton EPIC detectors are also sensitive to lower-energy protons (see Fig. 1), alternative models must also be considered to investigate changes in the intrinsic slope. In this section, we show that a single power law is not able to fit the considered spectra and that, among all of the models tried, a power law plus a black body is the only one that better describes all of the spectra and gives coherent results between the three detectors.

4.1. Fitting with a power law

The analysis was performed within XSPEC (version: 12.12.1) and fitting all spectra individually in the common energy range of 2–11.5 keV. The χ² statistic was applied to define the goodness of the fits. However, the reduced χ² is not a good estimator because the selected spectra have different number of grouped channels (from 200 to 2000) depending on the statistics. We instead used the null-hypothesis probability (NHP) provided by XSPEC, which is the probability of obtaining the observed value of χ² if the model is a correct representation of the data.

We defined as acceptable fits those with NHPs higher than 2.7 × 10⁻³, which corresponds to 3σ of a Gaussian distribution. Errors are reported at a 90% confidence level.

The power law model is not a valid model for most of the spectra, given that 72% of the 165 spectra have an NHP lower than the adopted threshold. Among those consistent with a single power law, about 60% are detected by the PN, and only in 6 observations of the 55 selected do the three detectors agree in observing a single power law. Their observation IDs are: 0037981201, 0037981601, 0037982701, 0111281501, and 0111281601 in interval B and 0844860801 in interval C. These spectra are generally characterized by low count rates so that the larger errors allow us to fit with a single power law.

Reducing the fitting energy range to 5–11.5 keV significantly improved the fits. In this case, 96% of the spectra are consistent with a single power law model, having an NHP within the 3σ Gaussian threshold. The remaining 4% spectra have in any case a probability within 4.5σ. The measured spectral indices range in the intervals of 2–7 for MOS and 1–6 for PN and, on average, the values in epochs A and B are lower that those in C and D. We observe no correlation between indices and rates.

In order to investigate possible alternative models for the whole energy range, we computed, for each spectrum and each instrument, the differences between the observed rates in the 2–5 keV band and the rates derived from the extrapolations of the 5–11.5 keV best-fitting power law. We found that for all MOS spectra and for most of the PN ones, the power law significantly underestimates the detected rates with statistically significant differences (> 3σ). As an example, Fig. 5 shows in the upper plot the extrapolation of the 5–11.5 keV power law up to 2 keV for the MOS1 spectrum of observation ID 0108060601, where a low energy excess is evident.

Fig. 5. Spectral fit with a power law. The upper plot shows the MOS1 spectrum relative to observation ID 0108060601. The top panel presents with a continuous line the 5–11 keV best-fit power law extrapolated down to 2 keV, while residuals over the whole range are presented in the lower panel. The lower plot is relative to the same figure for the 0852190101 PN spectrum.

For the PN spectra where the extrapolation does not underestimate the rates, some (0112551101 in the epoch A, 0037981601, 0037981701, and 0056022201 in B, 0844860301 and 0844860801 in C, and 0863810301 in D), present an excess consistent with zero and, in these cases, PN spectra are those fit by a single power law in the whole 2–11.5 keV range. In a few other observations (0052140201 and 0108061901 in epoch A, 0827241201, 0844210101, 0852190101, and 0854590401 in C, and 0862400101 in D), the extrapolation of the power law led to higher rates than the detected ones. In the lower plot of Fig. 5, an example one of these cases is presented. For these spectra, a different spectral model should be considered, and consequently we excluded them from the following analysis.

The excess rates were then correlated with the best-fit spectral indices and with the observed 2–5 keV rates. While no significant correlation was observed with the spectral indices for any of the three detectors, the correlation with the observed rates is instead well evident, with no significant differences between epochs. Figure 6 shows the 2–5 keV excess rate with respect to the 5–11.5 keV power law plotted as functions of the 2–5 keV rates for the three instruments. MOS1 and MOS2 show a strong correlation, while the correlation is weaker for PN data. From fitting all points with a line, we found that the 2–5 keV excess rate is 21% of the detected rate for MOS1 and MOS2, while it is only 5% for the PN.

Fig. 6. Excess rates with respect of the 5–11.5 keV best-fit power law in the range of 2–5 keV as a function of the total 2–5 rates for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel).

The power law obtained from fitting in the 5–11.5 keV range extrapolated to lower energies constitutes the most of the detection for almost all spectra, the residuals comprising a small fraction of the total rate. Taking into account this, and considering that this is a first check of the proton matrices with real data, to verify that the matrices correctly represent the spectral variations measured by the detectors, the best-fit values of the spectral indices are correlated with the hardness ratios. Figure 7 shows the correlation plots for the three detectors. It is evident that the spectral indices derived from the matrices closely correlate with the hardness ratios, with higher spectral indices corresponding to lower hardness ratios, and again no systematic differences are present between the four epochs.

Fig. 7. Power law spectral index computed from fitting the 5–11.5 keV ranges vs. the hardness ratios for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel).

4.2. Fitting with a power law plus a second component

According to the results presented in the previous subsection, a correct spectral modeling of the whole range must take into account the small excess at low energy. The 2–11.5 keV spectra were then fit, adding a second component: a second power law and black body were investigated. The fitting was performed in two steps. The statistics of the data, in fact, is not sufficient to limit the second component to the small fraction of the excess below 5 keV. Leaving all the parameters free in the 2–11.5 range, a strong reduction in the 5–11.5 keV power law contribution and a loss of correlation between the MOS and PN spectral indices were caused. The power law was then fit, at first, limiting the range to 5–11.5 keV, and then the second component was added to the model after fixing the power law parameters.

4.2.1. Power law plus black body

Using a black body as a second component, 83% of the spectra have NHPs > 2.7 × 10⁻³ and, in any case, the goodness of the fit improves for all of the spectra, increasing the values of the NHP, even for those compatible with a single power law model. The MOS spectra not compatible with this model have positive residuals below 3 keV, showing a more complex energy-dependent structure of the excess with respect to the simple black body model that we adopted. Most of these spectra have high count rates (0093160201 in epoch B, 0844860201 and 0844860401 in C, and 0827060801, 0862400101, and 0862730601 in D), but this effect is also present in the MOS2 spectrum 0844860601 with a moderate count rate.

For spectra with an acceptable fit, the kT values are always around 1 keV. Specifically, the two MOSs present an average value of 1.0 keV with an rms of 10%, while for the PN a slightly lower average value (0.9 keV) is obtained, coupled with a larger standard deviation of the best-fitting values (20%). Fitting the best-fit values with a line, no correlation either with the excess rate or with the spectral indices is observed for this parameter in MOS data. A weak correlation with the excess rate is instead detected for the PN (slope = 0.038 ± 0.011).

4.2.2. Fitting with two power laws

The composite model of two independent power laws is representative of only 37% of the spectra.

4.3. Fitting with a power law with an energy-dependent spectral index

In order to investigate the possibility that the low energy excess with respect to the power law is the effect of a spectral index steepening, we fit all spectra with a broken power law and with a power law whose index varies with energy as a log parabola (Massaro et al. 2004). Both models give input proton spectra that cannot be considered physically reliable.

4.3.1. Broken power law

One of the simplest models with a variable spectral index is the broken power law that foresees two constant spectral indices, respectively, above and below the break energy, E_break. In our dataset, this model gives reduced χ² lower than 2.7 × 10⁻³ for 86% of the spectra. The energy indices detected above E_break are consistent within errors with those computed in the 5–11.5 keV energy range with a simple power law in all epochs except in D, where MOS detectors shows a discrepancy of about 10%. However, all E_break values are close to the lower boundaries of the soft proton matrices. Their best-fit values are in the range of 22–27 keV for MOS and 10–19 keV for PN, with the consequence that the fitting procedure is not able to constrain the energy indices below E_break. Applying this model to the more statistically significant MOS2 spectra relative to 90% of the observing time of the flaring background used to validate the matrices gives unacceptable χ² (14 855 for 1896 d.o.f.) and underestimates the soft excess.

4.3.2. Log-parabola power law

The other model with an energy-dependent spectral index that we investigated was proposed by Massaro et al. (2004) for X-ray nonthermal emission. It follows the spectral distribution

$$\begin{array}{c}\hfill F(E)=K\times E^{-(\alpha +\beta \times log(E))},\end{array}$$(1)

where the parameter α is the photon index at 1 keV and β measures the curvature of the parabola. The energy-dependent photon index, Γ(E), formula, given by the log-derivative of Eq. (1), is

$$\begin{array}{c}\hfill \mathrm{\Gamma }(E)=\alpha +2\beta \times log(E).\end{array}$$(2)

This model fits 68% of the 165 spectra well, but only for 51% of the observation IDs is this true in all three detectors at the same time. In these 51% of cases, the parameters α and β measured by the two MOSs are always consistent with each other, while no coherence with the PN ones is obtained. Furthermore, computing the spectral indices over the matrix energy working range shows that above 50 keV the spectral indices cannot be considered in any way constant. Figure 8 shows the results for the spectral index measured by MOS1 (blue line) and by PN (red line) for observation ID 0110980601 (epoch D) computed with Eq. (2). From the plot, it is evident that the two instruments measure completely different α and β values and that the spectral index is not constant above 50 keV.

Fig. 8. Spectral index computed with Eq. (2) for MOS1 (continuous blue line) and PN (continuous red line) spectra of observation ID 0110980601. The dash-dotted lines indicate the 90% confidence level. For comparison, the spectral index relative to the 90% flux in the active magnetosphere given in Lotti et al. (2018) is indicated with the black line.

5. Spectral analysis systematics

To test the validity of the results and to quantify the systematic between the matrices, we compared the best-fit parameters obtained by the three detectors by analysing their ratios in pairs. This approach allows us to take into account errors detected by both detectors, while a simple fitting of one parameter with respect to the other would only account for errors on the fit variable.

5.1. Power law

Figure 9 displays, in the three panels, the ratios between the spectral indices measured by MOS1 and MOS2, MOS1 and PN, and MOS2 and PN as a function of the MOS2 and PN rates in the 5–11.5 keV range. The averages of the ratios have also been computed for each epoch and the resulting values are listed in Table 6, where it is evident that no significant variations are detected between epochs. Considering this, we assume the average values for the entire dataset to be systematic in the spectral indices. They are as follows: 0.97 ± 0.05 for the spectral index ratio MOS1–MOS2, and 1.22 ± 0.09 and 1.25 ± 0.11 for the MOS1 and MOS2 ratio with PN, respectively. These constants are plotted with a black line in the three panels of Fig. 9, while its error computed as rms is indicated by the dashed lines.

Fig. 9. Ratio of the spectral indices measured by MOS1 and MOS2 (top panel), MOS1 and PN (middle panel), and MOS2 and PN (bottom panel) as a function of the rates in the 5–11.5 keV band. The dashed black line is the best-fit constant, while the dotted lines indicate the relative error range. Colors indicate the four epochs: epoch A in red, epoch B in blue, epoch C in green, and epoch D in gray.

The same procedure was applied to quantify the systematic in the fluxes measured by the three instruments. To this purpose, the number of input protons at the telescope pupil in the energy range of 30–300 keV and its errors were computed with the XSPEC command flux. The plots of the ratios as a function of the rates are shown in the three panels of Fig. 10. Computing the average values in each epoch, we find that, in this case, consistent fluxes within the statistical uncertainties are measured in epochs A and B, and C and D. For this reason, two constants have been obtained, one for the year 2001–2002 and one for the 2019–2020. In the first case, the ratio MOS1–MOS2 has an average of 1.08 ± 0.08, which becomes 1.48 ± 0.22 for epochs C–D. The average values for the ratios MOS1–PN and MOS2–PN are: 2.17 ± 0.57 and 1.99 ± 0.49 for A and B, while they are 1.06 ± 0.27 and 0.74 ± 0.23 for C and D, respectively. The two values are shown in Fig. 10 with dashed and dash-dotted lines, respectively.

Fig. 10. Ratio of the 30–300 proton input flux measured by MOS2 and by PN as a function of the PN rates in the same band. The dashed blue line is the average value for the interval A and B, while the dash-dotted green line is the average for C and D. The dotted lines indicate the error ranges. Colors on the points are relative to the four epochs: epoch A in red, epoch B in blue, epoch C in green, and epoch D in gray.

5.2. Black body

As was already stated, the black body component is a phenomenological model introduced to take into account the small excess in fitting the whole 2–11.5 keV range with a simple power law. For this model, the best-fit spectral parameters are evaluated with large errors; however, the ratios of the kT and flux values are consistent with those in Table 6 relative to the spectral index and flux of the power law.

6. Discussion

From the analysis performed on the selected 55 observations of flares, we detected variation in the MOS response to protons with time, identified two components for the spectral model, and evaluated the systematic errors in the best-fit values of both components. The MOS data collected from flares in the epoch 2001–2002 appear to be inconsistent with those from 2019-2020. MOS1 and MOS2 exhibit statistically equivalent counts a couple of years after launch, while a noticeable reduction occurs after two decades, causing a different count ratio (MOS2 registers lower counts than MOS1). This discrepancy is not reconcilable with the loss of two MOS1-CCDs. The reduction is proportional across all energy bins, as is evidenced by the lack of dependency of the spectral index ratio as a function of time or count rate (see Fig. 9). Considering that the matrices are relative to only single events, we checked for some observations if this reduction could be due to a different distribution of the event grades in 2019-2020 flare photon lists, but no significant differences were obtained.

The effect might be related to the variation in the electrode layer structure observed primarily for MOS2 due to an increase in the surface layers of nitrogen and oxygen and due to the addition of a carbon layer deposited on the electrodes over the course of the mission (Plucinsky et al. 2017). From the average rate ratios of the two MOS cameras with respect to PN, as is shown in Table 4, we assessed the average reduction in the proton effective area between the C and D epochs compared to the A and B epochs, accounting for the loss in MOS1 due to damage in the two CCDs. We found that the MOS1 effective area decreases to 68 ± 5% of its value in 2001, while MOS2 has decreased to 38 ± 6%. Our results are, then, in agreement with studies of the contamination layers that have have formed on the MOS cameras since launch, which indicate higher contamination in the MOS2 camera compared to MOS1, while no contamination was observed on PN (Plucinsky et al. 2017).

In any case, the spectral analysis of MOS flares necessitates a new version of the matrices, with variations mainly involving the computation of the total flux at the telescope input pupil. However, such corrections can only be implemented if the real causes of this phenomenon are understood. If this is not possible, the discrepancy of the matrices with respect to the real transmission can be considered a systematic in the input flux measurements.

Concerning the spectral analysis, the most important information is that a single power law is not able to adequately fit the entire energy range of 2–11.5 keV for most of the spectra, while it is a suitable model in the range of 5–11.5 keV, and its extrapolation to lower energies constitutes the most of the emission. In fact, the differences in the power law rate with respect to the total detection are ∼20% for the MOS and ∼5% for the PN, and the 5–11.5 keV power law spectral indices are well correlated with the hardness ratios (see Fig. 7). To take into account the soft excess, a second spectral component was included. The most suitable model was a phenomenological black body.

However, we also investigated the possibility of a physical steepening of the power law using two different models: a broken power and a power law with a log parabola spectral index. Neither model gives acceptable results. For the broken power law, the unreliability comes from the impossibility of confining the spectral indices below the energy breaks, probably because of the limited validity range of the matrices. The log-parabola power law produces acceptable fits for about half of the spectra, but no coherence between the MOS and PN best values of α and β was obtained. In summary, with the present version of the proton matrices and from the analysis presented in this paper, it is not possible to state that there is any steepening of the spectrum. The new version of the matrices or the analysis of a larger sample of spectra could add more items to this point.

To assess the spectral variability of the black body as a second physical component, we computed a color-color diagram from a few observations (0011830201, 0037980401, and 0052140201). With this aim, we selected light curves in different energy ranges – soft (0.5–2 keV), intermediate (2–5 keV), and hard (5–11.5 keV) – and computed the intermediate-soft and hard-intermediate ratios with a time-bin of 50 seconds. Considering that the black body component contributes mostly below 5 keV, we expect that 0.5–2 keV rates could give a strong indication of its variability that we cannot obtain from the spectral analysis being the matrices not valid in this range. As a further check, we also correlated the 0.5–2 keV rates with the 5–11.5 keV ones. The top panel of Fig. 11 shows as an example the color-color diagram of observation ID 0037980401, with the correlation between rates shown in the bottom panel. Throughout the color-color diagram, points are concentrated around an average with variation always lower than 30% for both the intermediate-soft ratio and the hard-intermediate one. Moreover, a very strong correlation between rates in the ranges of 0.5–2 keV and 5–11.5 keV is always observed. We conclude, then, that there is no significant difference in the spectral variability of the two components, and that if this second component is real, it must vary coherently with the power law. However, because we are limited by the validity range of the matrices, a more physical origin – that is, a component coming from heavier ions – cannot be entirely excluded.

Fig. 11. Analysis of spectral variability in the observation ID 0037980401: the top panel shows the color-color diagram and the correlation between rates in the ranges 0.5–2 keV and 5–11.5 keV is presented in the bottom panel.

We strongly think that all these results suggest that the power law is the main model for the proton input spectra, and that the second component, necessary to fit the data well at lower energies, is an artifact of the matrices. As is shown in Paper I, the relative flux of the second component in the MOS analysis increases with the passive material of the optical filters. This, together with the lower excess found in the PN, which is back-illuminated and does not include the electrode structure in front of the depletion region, brings us to the conclusion that it can be attributed to imprecise modeling of the proton transmission at the focal plane in the simulator used for generating the response matrices. In this case, considering its small contribution to the total count, it can be accounted for by a systematic factor in the input flux derived from the power law.

It is worth pointing out that in some of the PN spectra, the power law overestimates the low energy count rates, and the corresponding MOS spectra have a slightly lower excess with respect to spectra with the same rates. For these, the input proton spectrum is probably flattening below 2 keV, and the simple power law is not applicable.

This work has also allowed us to quantify the systematic uncertainties in the spectral analysis of flares with the proton matrices. These results are quite simple: the three detectors always measure coherent spectral indices independently from the epoch. The differences are of ∼3% between MOS1 and MOS2 and about 20% between MOS and PN. This higher discrepancy could be an effect of incorrect modelling of transmission at lower energies, which is more severe in the MOS response matrix.

The measurement of the total flux input at the telescope pupil gives less clear results because of two different behaviors observed in the data: one in the years 2001–2002 (epochs A and B) and the other in 2019–2020 (C and D). This makes it difficult to disentangle effects arising from the response matrix from those attributed to variations in the instrument response. The discussion of the total 30–300 keV input flux is then addressed for the 2001–2002 dataset and the 2019–2020 dataset separately.

In epochs A and B, MOS1 and MOS2 exhibit statistically equivalent proton fluxes at the telescope pupil (see Fig. 10). However, discrepancies arise in the fluxes measured by the MOS and by the PN, revealing a systematic difference of approximately a factor of two (MOS measures twice the PN flux). The systematic increase of 20% in the MOS spectra indicates that the different excess at lower energies is likely affecting the fitting of the primary power law, with an impact on the input proton fluxes.

In the C and D datasets, the two MOS detectors do not give coherent fluxes, even when accounting for the loss of the two CCDs in MOS1. The systematics of MOS with respect to PN has a more complex behavior. Apart from the second factor identified from 2001–2002 data, an additional factor need to be included. In our context, this factor can be computed from the ratio between the average rates in this period and the average ones in 2001–2002. However, without understanding the observed rate decrease of the two MOS detectors, it is difficult to correct the response matrix.

7. Conclusions

This paper presents an analysis of 55 flare spectra detected during December-January and July-August of 2001–2002 and 2019–2020, with the main task being to test the matrices on spectra of proton flares and to study the differences in the measured best-fit parameters between the MOS and PN proton response matrices. The spectra have been deconvolved within XSPEC using as a model a power law plus a black body. The main results of this analysis can be summarized in the following points:

• MOS data show some inconsistency in the detected rates across the 20 years of the satellite’s operation with respect to the corresponding PN ones;

• The physical model representative of the proton spectra at the input of the telescope is a power law, while the black body is a phenomenological component introduced to take into account the low energy excess coming from the extrapolation to lower energies. In any case, this black body component contributes 21% for the MOS and 5% for the PN to the total flux in the 2–5 keV energy range;

• The systematic error in the spectral indices measured by the two MOS detectors is 3%, while that of the MOS with respect to the PN is ∼24%;

• A systematic difference within a factor of two between MOS and PN input fluxes was determined across the four epochs, as a result of inaccuracies in the proton response matrices and changes in the instrument proton detection across 20 years of operations. This factor was computed directly from rates, although a more detailed analysis than that presented in this paper is necessary to clarify the causes.

A future version of the matrices would include all the elements pointed out by this work and some more specific analysis would be required for us to better understand and model the 2019–2020 loss of transmission efficiency. A systematic analysis of flare events across all the operation life of the satellite could give some clarification. Moreover, some laboratory measurements on the proton transmission from the filters would allow us to validate Geant4 for this kind of interaction. To conclude, this work is also necessary for a future X-ray observatory carrying on-board grazing incidence optics, such as Athena, to produce with the same physical model matrices necessary for predicting the proton flux at the focal plane.

Acknowledgments

Authors thank the anonymous referee for the comments and suggestions on the paper that really improved it. TM thanks G. Cusumano and A. D’Aí for their useful comments on the paper. The research leading to these results has received funding from European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158), the Accordo Attuativo ASI-INAF n. 2019-27-HH.0 and the EXACRAD – CTP ESA contract n. 4000121062/17/NL/LF.

References

All Tables

All Figures

	Fig. 1. Proton response matrix effective area for MOS (top panel) and PN (bottom panel). The application of the thin (medium) filter is represented by a continuous (dash-dotted) curve.
In the text

	Fig. 2. Light curves relative to observation ID 0844860601 (2020-01-06) for MOS1, MOS2, and PN. The dashed vertical lines delimit the interval in which the background rate was evaluated, the dash-dotted lines indicate the highest background rate, and the dotted lines indicate the lowest flare rate.
In the text

	Fig. 3. Correlations between rates detected by the three instruments in the four epochs. Colors indicate epochs (see text). The superimposed dashed black line is relative to the best fit of epoch B. The top panel is relative to MOS1 vs. MOS2 rates, the middle panel to MOS1 vs. PN, and the bottom panel to MOS2 vs. PN.
In the text

	Fig. 4. Hardness ratio between counts in the ranges 5–11.5 keV and 2–5 keV vs. 2–11.5 keV rate for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel). Colors indicate the four epochs (see text).
In the text

	Fig. 5. Spectral fit with a power law. The upper plot shows the MOS1 spectrum relative to observation ID 0108060601. The top panel presents with a continuous line the 5–11 keV best-fit power law extrapolated down to 2 keV, while residuals over the whole range are presented in the lower panel. The lower plot is relative to the same figure for the 0852190101 PN spectrum.
In the text

	Fig. 6. Excess rates with respect of the 5–11.5 keV best-fit power law in the range of 2–5 keV as a function of the total 2–5 rates for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel).
In the text

	Fig. 7. Power law spectral index computed from fitting the 5–11.5 keV ranges vs. the hardness ratios for MOS1 (top panel), MOS2 (middle panel), and PN (bottom panel).
In the text

	Fig. 8. Spectral index computed with Eq. (2) for MOS1 (continuous blue line) and PN (continuous red line) spectra of observation ID 0110980601. The dash-dotted lines indicate the 90% confidence level. For comparison, the spectral index relative to the 90% flux in the active magnetosphere given in Lotti et al. (2018) is indicated with the black line.
In the text

	Fig. 9. Ratio of the spectral indices measured by MOS1 and MOS2 (top panel), MOS1 and PN (middle panel), and MOS2 and PN (bottom panel) as a function of the rates in the 5–11.5 keV band. The dashed black line is the best-fit constant, while the dotted lines indicate the relative error range. Colors indicate the four epochs: epoch A in red, epoch B in blue, epoch C in green, and epoch D in gray.
In the text

Fig. 10. Ratio of the 30–300 proton input flux measured by MOS2 and by PN as a function of the PN rates in the same band. The dashed blue line is the average value for the interval A and B, while the dash-dotted green line is the average for C and D. The dotted lines indicate the error ranges. Colors on the points are relative to the four epochs: epoch A in red, epoch B in blue, epoch C in green, and epoch D in gray.

In the text

	Fig. 11. Analysis of spectral variability in the observation ID 0037980401: the top panel shows the color-color diagram and the correlation between rates in the ranges 0.5–2 keV and 5–11.5 keV is presented in the bottom panel.
In the text