Since 1996 the blast wave driven by SN 1987A has been interacting with the dense circumstellar material, which provides us with a unique opportunity to study the early evolution of a newborn supernova remnant (SNR). Based on the XMM-Newton RGS and EPIC-pn X-ray observations from 2007 to 2019, we investigated the post-impact evolution of the X-ray-emitting gas in SNR 1987A. The hot plasma is represented by two nonequilibrium ionization components with temperatures of ∼0.6 keV and ∼2.5 keV. The low-temperature plasma has a density ∼2400 cm−3, which is likely dominated by the lower-density gas inside the equatorial ring (ER). The high-temperature plasma with a density ∼550 cm−3 could be dominated by the H ii region and the high-latitude material beyond the ring. In the last few years, the emission measure of the low-temperature plasma has been decreasing, indicating that the blast wave has left the main ER. But the blast wave is still propagating into the high-latitude gas, resulting in the steady increase of the high-temperature emission measure. Meanwhile, the average abundances of N, O, Ne, and Mg are found to be declining, which may reflect the different chemical compositions between the two plasma components. We also detected Fe K lines in most of the observations, showing increasing flux and centroid energy. We interpret the Fe K lines as originating from a third hot component, which may come from the reflected shock heated gas or originate from Fe-rich ejecta clumps shocked by the reverse shock.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
As the nearest supernova (SN) observed since Kepler's SN of 1604, SN 1987A provides us with a unique opportunity to observe in detail the onset of SN remnant (SNR) formation and subsequent evolution for an SN whose observational properties are known in detail. SN 1987A is located in the Large Magellanic Cloud (LMC) and was first detected on 1987 February 23. The progenitor star was identified as a blue supergiant (Sanduleak −69° 202). SN 1987A is surrounded by a peculiar circumstellar material (CSM) system represented by three coaxial rings, the origin of which may be explained by the merging of a binary system at about 20,000 yr before the explosion (Podsiadlowski et al. 1992; Morris & Podsiadlowski 2007) or by collisions of a fast wind with material from an earlier slow wind, during the progenitor's post-main-sequence evolution (Chita et al. 2008). Today, over 30 yr after the SN outburst, the blast wave has encountered and heavily interacted with the CSM structures, especially with the dense equatorial ring (ER). The impact of the blast wave on the ER resulted in a complex shock system, which compressed and heated the gas, giving rise to radiation in different energy bands. As a result, SN 1987A has made its transition to the remnant phase (see McCray & Fransson 2016 for a recent review).
As a newborn remnant, SNR 1987A has in its short life already gone through several dynamics phases, as monitored regularly across the electromagnetic spectrum. Observations in the X-ray band are in particular crucial for understanding the underlying physics of the shock interaction processes, as it contains most of the emitting energy. The X-ray emission of SNR 1987A was first detected by ROSAT ∼1500 days after the explosion and it steadily brightened until ∼3000 days, as the blast wave began interacting with the H ii region interior to the ER (Hasinger et al. 1996). Starting from ∼4000 days, the blast wave encountered the dense clumps protruding from the inner edge of the ER, as evidenced by the appearance of several "hot spots" in the optical band (e.g., Sonneborn et al. 1998; Lawrence et al. 2000). Corresponding to this—as observed by Chandra since ∼4600 days—the soft X-ray flux of SNR 1987A was found to exceed the linear extrapolation of the ROSAT light curve and brighten more rapidly than observed by ROSAT (Burrows et al. 2000; Park et al. 2002). At ∼6000 days the blast wave impacted the main body of the ER, as indicated by a dramatic increase in the soft X-ray flux (e.g., Park et al. 2005) and a sudden decrease in the shock velocity (e.g., Racusin et al. 2009). Thereafter, the soft X-ray flux continued increasing—first exponentially and then linearly (Helder et al. 2013)—until day ∼9500, when the light curve started to level off (Frank et al. 2016). Together with the optical/infrared light curves, this indicates that the blast wave had traversed through the entire dense ring (Fransson et al. 2015; Arendt et al. 2016).
Recently, the overall X-ray light curve has been well reproduced by 3D hydrodynamic (HD)/magnetohydrodynamic (MHD) simulations (e.g., Orlando et al. 2015, 2019, 2020). These simulations indicate that the X-ray emission is dominated first by the shocked H ii region, then by the shocked dense ring, and that SNR 1987A will enter a third phase around 32–34 yr after the explosion, which will be dominated by the SN ejecta heated by the reverse shock. According to these simulations, we are expecting to find the X-ray features from the shocked ejecta in the most recent observations, or in the near future. On the other hand, recent optical observations have revealed new features of spotlike and diffuse emission outside the ER (since day ∼9500), which may be related to high-latitude material that extends from the ER toward the outer rings (Larsson et al. 2019). Therefore, we are also expecting to observe the associated X-ray emission from this material. In addition, a recent (2016–2018) enhancement of gigaelectronvolt emission from the SNR 1987A region has been reported based on Fermi/LAT observations (Malyshev et al. 2019). X-ray observations could be crucial to clarifying the nature of the gigaelectronvolt signal (a detection of nonthermal emission in the 10–20 keV energy band has been recently reported by Greco et al. 2021, which may be indicative of a pulsar wind nebula).
In this work, we report our analysis of XMM-Newton observations of SNR 1987A taken during the last decade, aiming at the post-impact evolution of the X-ray-emitting gas. We performed detailed spectral modeling and constrained multiple physical parameters to unveil the plasma properties. The recent changes in the plasma properties clearly show that the blast wave has left the ER and is propagating into the high-latitude gas. This paper is organized as follows: Section 2 describes the observations and data reduction procedure, Section 3 presents the spectral modeling and the main results, Section 4 discusses further implications of the results, and Section 5 gives a brief summary.
2. Observations and Data Reduction
SNR 1987A was regularly monitored by XMM-Newton from 2007 to 2017 (principal investigator (PI): F. Haberl). In 2019 November, we obtained a new observation (PI: Malyshev, for the H.E.S.S. Collaboration), aiming at its latest evolution. Therefore, we got a set of data containing 12 observations in total, covering a time interval over 12 yr, which is summarized in Table 1. In this work, we utilized these data to analyze the X-ray spectra of SNR 1987A taken in the last decade. The high-resolution RGS spectra provide tight constraints on the plasma properties. We also included the EPIC-pn spectra for a complete energy coverage to 10 keV following Sturm et al. (2010).
Table 1. Observations
|texp (ks) a||(ks) b|
|0406840301||2007 Jan 17||7267||106.9||111.3||61.1||109.8|
|0506220101||2008 Jan 11||7627||110.1||114.3||70.7||102.4|
|0556350101||2009 Jan 30||8012||100.0||101.9||66.4||101.8|
|0601200101||2009 Dec 11||8327||89.9||91.8||82.4||91.7|
|0650420101||2010 Dec 12||8693||64.0||65.9||52.7||65.9|
|0671080101||2011 Dec 2||9048||80.6||82.5||64.2||80.6|
|0690510101||2012 Dec 11||9423||68.0||69.9||59.4||69.8|
|0743790101||2014 Nov 29||10,141||78.0||79.6||56.4||79.4|
|0763620101||2015 Nov 15||10,492||64.0||65.9||58.0||65.8|
|0783250201||2016 Nov 2||10,845||72.4||74.3||50.3||74.2|
|0804980201||2017 Oct 15||11,192||77.5||79.4||27.8||79.3|
|0831810101||2019 Nov 27||11,964||32.4||34.9||11.1||34.8|
Notes.a Total exposure times. b Total good time intervals after background flare removal.
We used the XMM-Newton Science Analysis Software (SAS, version 16.1.0) 7 to process the data. For the EPIC-pn data, the task epchain was used to produce calibrated photon event files, and then pn-filter was used to filter out soft-proton flares and to remove affected time intervals. The EPIC-pn spectra of SNR 1987A were extracted from a circular region centered on the source with a radius of 25'', 8 while the background spectra were extracted from a source-free region nearby (shown in Figure 1). We note that the background contributes only ∼1% to the EPIC-pn spectrum of SNR 1987A, so the background subtraction has only a minor contribution to the systematic error budget. We selected only single events (PATTERN = 0) for the source spectra in order to get rid of the photon pileup effect. The RGS data were processed using rgsproc, and then the spectra were extracted within the low-background time intervals (RATE < 0.2). For a reliable application of the χ2 statistics, both the EPIC-pn and RGS spectra were rebinned using the FTOOL 9 task grppha to obtain at least 25 counts per bin.
3. Spectral Analysis
In order to better trace the evolution of the X-ray properties of the remnant, we used a consistent method to analyze all the observations. We used XSPEC (version 12.10.1) 10 with AtomDB 3.0.9 11 for spectral analysis. Unless otherwise specified, in this paper we present metal abundances with respect to their solar values (Wilms et al. 2000, hereafter W00).
3.1. Individual Emission Lines and Their Fluxes
Taking advantage of the high spectral resolution of RGS, we were able to identify the individual emission lines and to constrain their fluxes. The 0.35–2.5 keV RGS spectra of SNR 1987A are dominated by the emission lines, which mainly come from the α-elements such as O, Ne, Mg, and Si, and from the Fe L transitions around 0.7–0.9 keV (Figure 2). In order to identify these lines, we constructed a heuristic model and fit it to the first-order and second-order RGS1/RGS2 spectra simultaneously.
The spectral model consists of a thermal continuum component (described by the nlapec model in XSPEC and including thermal bremsstrahlung, radiative recombination continua, and two-photon emission) and several Gaussian line profiles (Gauss). Both the continuum and the emission lines were subjected to foreground absorption (tbabs), for which the Galactic absorption was fixed at NH,Gal = 6 × 1020 cm−2 with abundances set to be solar, and the LMC absorption was fixed at NH,LMC = 2.2 × 1021 cm−2 with abundances set to be the average LMC abundances given by Russell & Dopita (1992). Previous studies have revealed line broadening features in the dispersed X-ray spectra of SNR 1987A (e.g., Zhekov et al. 2005; Miceli et al. 2019). Despite its complex origins, the line broadening can be approximately described with a power-law function of energy (e.g., Sturm et al. 2010; Dewey et al. 2012). Therefore, we simply adopted a Gaussian smoothing component (gsmooth) for the broadening effect. The gsmooth function convolves the spectrum with a Gaussian function, of which the Gaussian sigma is energy dependent: σ = σ6(E/6 keV)α , where σ6 is the width at 6 keV. We also introduced a redshift component (vashift) to account for the line shift which may result from the systemic motion of the remnant.
With the heuristic spectral model described above, we identified 36 emission lines in the RGS spectra of SNR 1987A. The emission lines, together with their rest-frame centroids and (absorption-corrected) fluxes, are listed in Table 2. Among these emission lines, the He-like triplets and the H-like Ly lines are especially useful for thermal X-ray plasma diagnostics. In Figure 3, we plot the light curves of various emission lines from the He- and H-like ions of O, Ne, Mg, and Si. Most of the emission lines follow a similar flux evolution pattern, showing initially a brightening, reaching a peak value at around 2011–2015, and after that exhibiting a drop-off in flux. This is consistent with the overall soft X-ray light curve of SNR 1987A (e.g., Frank et al. 2016; see also Section 3.3 below for our updates with the new data). For the oxygen lines, we calculated the G-ratios (defined as (f + i)/r, where f, i, and r stand for the fluxes of the forbidden line, the intercombination lines, and the resonance line, respectively), the Lyα/Heα-r flux ratios, and the Lyβ/Lyα flux ratios. For Ne, Mg, and Si, we calculated the Heα-f/Heα-r flux ratios and the Lyα/Heα-r flux ratios. As shown in Figure 3, we find that while the G-ratios (or Heα-f/Heα-r) show no significant variation during this 12 yr interval, the Lyα/Heα flux ratios are continuously increasing. This could result from the increase in the average plasma temperature, but may also be induced by nonequilibrium ionization (NEI) effects.
Table 2. Identified Lines in RGS Spectra
|Line||Centroid (keV)||Flux (10−5 photons cm−2 s−1)|
|2007 Jan||2008 Jan||2009 Jan||2009 Dec||2010 Dec||2011 Dec||2012 Dec||2014 Nov||2015 Nov||2016 Nov||2017 Oct||2019 Nov|
|Fe xvii||0.7271||19.06 ± 2.09|
|Fe xviii||0.7715||4.76 ± 0.49|
|O Lyβ||0.7746 a|
|Fe xix||0.9172 b|
|Fe xxi||1.0004||1.17 ± 0.30|
|Fe xxi||1.0093 c|
|Fe xxii||1.0534 d|
|Fe xxiii||1.1291 e||1.39 ± 0.18|
|(eV)||...||84 ± 19||5 ± 1|
|α||...||2.15 ± 0.12|
Notes. Errors are for the 1σ confidence range.a Could also be Fe xviii line at 0.7747 keV. b Could also be Fe xxi line at 0.9179 keV. c Could also be Fe xvii line at 1.0108 keV. d Could also be Fe xxiii line at 1.0565 keV. e Could also be Ne Heγ line at 1.1270 keV.
Our spectral analysis also provides measurements of the line broadening and the systemic velocity—see σ6 and v87A in Table 2. Except for an obviously divergent value from the 2019 observation (which has the shortest effective exposure and thus the lowest photon statistics), the line shift indicates a systemic velocity of SNR 1987A around 250 km s−1, consistent with the value derived from optical observations (286.74 ± 0.05 km s−1, e.g., Gröningsson et al. 2008). We note that the systematic error could be significantly greater than the statistical error in this case for systemic velocity measurement. The RGS wavelength scale accuracy is 5 mÅ and 4 mÅ in the first and second order, respectively. Since most of the bright emission lines lie in the range of 10–20 Å, a ∼5 mÅ uncertainty in the line position will transfer to a ∼70–150 km s−1 uncertainty in the systemic velocity. Taking this systematic error into consideration, we estimated the overall uncertainty in systemic velocity as ∼100 km s−1, which may explain the rather large dispersion of the measured v87A values. On the other hand, the line shift could also be affected by the asymmetric distribution and expansion of the X-ray-emitting gas. As revealed by Chandra observations (Frank et al. 2016, Figure 5 therein), the surface brightness distribution of the ER has been dynamically changing since the impact, which may have also led to the large dispersion of the values we obtained here for v87A. The line broadening effect is characterized by a Gaussian smooth function, by which we found the power-law index α ∼ 2, similar to the value obtained by Sturm et al. (2010), and consistent with the Chandra LETG results (Zhekov et al. 2009). The power-law index provides us further information about the origin of the line broadening (see also the discussion in Dewey et al. 2012). The source extent will result in a broadening that is constant in wavelength and thus gives α = 2 in energy. This is nearly the case for LETG/HETG, but not for RGS since SNR 1987A is spatially unresolved by XMM-Newton. On the other hand, the Doppler broadening (including the bulk Doppler broadening and the thermal broadening) gives α = 1 assuming a constant velocity over all the spectral lines. Our result of α ∼ 2 (>1) indicates that the spectral lines at higher energy (which are from the heavier elements) have higher velocities. This is consistent with the results obtained by Miceli et al. (2019), in the sense of a proportional relation between the post-shock ion temperature and the ion mass.
3.2. Global Spectral Fitting
In order to better characterize the physical properties of the X-ray-emitting gas in SNR 1987A, we fitted the RGS (0.35–2.5 keV) and EPIC-pn (0.3–10.0 keV) spectra simultaneously with detailed NEI plasma emission models.
3.2.1. Two-temperature vnei Model
While the actual distribution of temperatures, densities, and other plasma properties can be quite complex, previous studies (e.g., Sturm et al. 2010; Frank et al. 2016) have shown that the X-ray emission of SNR 1987A can be well approximated by a spectral model containing two plasma components: a low-temperature component associated with the dense ring material, and a high-temperature component associated with the lower-density gas. Therefore, we first fit the RGS and pn spectra with a two-temperature NEI plasma model, described by two vnei models in XSPEC. The vnei model characterizes the average thermal and ionization state of the plasma by a single electron temperature kTe and a single ionization parameter τ = ne t, respectively, where ne is the electron density and t is the elapsed time since the gas was shocked. Following Sturm et al. (2010) and the references therein, with abundances converted to those in W00, the abundances of N, O, Ne, Mg, Si, S, and Fe were set as free parameters, while those of He (set to 2.57), C (0.14), Ar (0.76), Ca (0.49), and Ni (0.98) were fixed, since they produce no significant features in our data. The abundances of the two components were assumed to be the same. The Galactic absorption and the LMC absorption were fixed at NH,Gal = 6 × 1020 cm−2 and NH,LMC = 2.2 × 1021 cm−2, respectively. 12 Similar to the model described in Section 3.1, the spectra were convolved with the gsmooth function and shifted by vashift to account for the line broadening and systemic velocity, respectively. The parameters of gsmooth and vashift were fixed at the values obtained above by the RGS spectral fitting (Table 2) for each observation.
This two-temperature vnei model gives acceptable fits to the spectra for all the observations, with reduced chi-square values of ∼ 1.21–1.41 (dof ∼ 1234–2433). The low-temperature component is characterized by kTe,LT ∼ 0.6 keV and the high-temperature component is kTe,HT ∼ 2.5 keV, which show no significant variations during the last decade. The detailed fitting results are summarized in Table 3, and an example of the fitted spectra is shown in Figure 4.
Table 3. Spectral Fitting Results of the Two-temperature vnei Model
|2007 Jan||2008 Jan||2009 Jan||2009 Dec||2010 Dec||2011 Dec||2012 Dec||2014 Nov||2015 Nov||2016 Nov||2017 Oct||2019 Nov|
|(keV)||0.56 ± 0.01||0.55 ± 0.01||0.57 ± 0.01||0.60 ± 0.01||0.61 ± 0.01||0.62 ± 0.01||0.67 ± 0.01||0.69 ± 0.01||0.71 ± 0.03|
|(1010 cm−3 s)|
|Norm (10−3 cm−5)||6.66 ± 0.27|
|kTHT (keV)||2.35 ± 0.10||2.42 ± 0.07|
|τHT (1010 cm−3 s)|
|Norm (10−3 cm−5)||2.20 ± 0.11||2.68 ± 0.08||3.19 ± 0.11|
|O||0.08 ± 0.01||0.08 ± 0.01||0.08 ± 0.01||0.07 ± 0.01||0.09 ± 0.01||0.08 ± 0.01||0.07 ± 0.01||0.07 ± 0.01||0.06 ± 0.01||0.07 ± 0.01||0.06 ± 0.01||0.05 ± 0.01|
|Ne||0.31 ± 0.02||0.30 ± 0.02||0.29 ± 0.01||0.27 ± 0.01||0.28 ± 0.02||0.29 ± 0.02||0.24 ± 0.02|
|Mg||0.39 ± 0.03||0.37 ± 0.02||0.37 ± 0.02||0.33 ± 0.03||0.33 ± 0.03||0.30 ± 0.04|
|Si||0.61 ± 0.04||0.66 ± 0.04||0.59 ± 0.04|
|S||0.61 ± 0.09||0.58 ± 0.09||0.55 ± 0.07||0.57 ± 0.07|
|Fe||0.23 ± 0.01||0.26 ± 0.01||0.29 ± 0.01||0.29 ± 0.01||0.28 ± 0.01||0.26 ± 0.01||0.27 ± 0.01||0.26 ± 0.02|
Note. Errors are for the 90% confidence range.
3.2.2. Two-temperature vpshock Model
In addition to the vnei model presented above, we modeled the global spectra of SNR 1987A with a plane-parallel shock model (Borkowski et al. 2001; described by the vpshock model in XSPEC). This model describes the X-ray emission from a hot gas heated by a plane-parallel shock with a constant post-shock temperature. It differs from the vnei model in its linear distribution of the ionization parameter τ versus the emission measure, instead of assuming a single ionization parameter. The NEI plasma in the vpshock model is characterized by a lower limit τl = ne tl at the shock front and an upper limit τu = ne tu for the earliest-shocked gas. We fixed τl at zero and let τu vary freely. Other parameters were set in the same way as those in the vnei model.
We find that the two-temperature vpshock model fits the data better, which is indicated by an average drop of the reduced chi-square Δχ2 r ∼ 0.07. This suggests that the single ionization state scenario in the vnei model is oversimplified, while the plane-parallel shock model better represents the actual conditions. However, the best-fit values of the parameters show no significant discrepancy from those in the vnei fitting, except for the values of τu, which are typically one order of magnitude higher than the average τ given by the vnei fitting (Table 4). We further discuss the physical implications of the difference between τ and τu in Section 4.2.
Table 4. Spectral Fitting Results of the Two-temperature vpshock Model
|2007 Jan||2008 Jan||2009 Jan||2009 Dec||2010 Dec||2011 Dec||2012 Dec||2014 Nov||2015 Nov||2016 Nov||2017 Oct||2019 Nov|
|kTLT (keV)||0.58 ± 0.01||0.60 ± 0.01||0.62 ± 0.01||0.63 ± 0.01||0.63 ± 0.01||0.68 ± 0.01||0.66 ± 0.01||0.69 ± 0.01||0.63 ± 0.02||0.63 ± 0.03|
|(1011 cm−3 s)|
|NormLT (10−3 cm−5)|
|(keV)||2.52 ± 0.21||2.44 ± 0.09||2.47 ± 0.07|
|(1011 cm−3 s)||1.64 ± 0.01||1.67 ± 0.01||1.72 ± 0.01||1.72 ± 0.01|
|Norm (10−3 cm−5)||0.97 ± 0.12||3.74 ± 0.14|
|O||0.16 ± 0.01||0.16 ± 0.01||0.21 ± 0.02||0.19 ± 0.01||0.17 ± 0.01||0.19 ± 0.01||0.16 ± 0.01||0.17 ± 0.01||0.11 ± 0.01||0.13 ± 0.02|
|Ne||0.41 ± 0.03||0.41 ± 0.02||0.52 ± 0.03||0.45 ± 0.03||0.46 ± 0.03||0.34 ± 0.03|
|Mg||0.47 ± 0.03||0.47 ± 0.03||0.59 ± 0.04||0.49 ± 0.04|
|Si||0.79 ± 0.05||0.72 ± 0.09|
|S||0.66 ± 0.10||0.62 ± 0.08||0.53 ± 0.09||0.45 ± 0.10|
|Fe||0.26 ± 0.02||0.35 ± 0.02||0.36 ± 0.02||0.33 ± 0.02|
Note. Errors are for the 90% confidence range.
3.3. X-Ray Light Curve
The X-ray light curves of SNR 1987A up to 2015 have been reported by Frank et al. (2016) based on both Chandra and XMM-Newton observations. Here we update the reader on recent changes in the X-ray fluxes with the new XMM-Newton data.
Based on our two-temperature vpshock fitting described above, we measured the 0.5–2.0 keV and 3.0–8.0 keV fluxes of the remnant from 2015 to 2019, and show them in Figure 5 together with the earlier results. As noticed by Frank et al. (2016), the soft-band flux of SNR 1987A reached a plateau at around 9000 days, and has remained approximately constant since then. The new observations show that after around 10,500 days, the soft-band flux started to decline. In the latest observation (11,964 days), the soft-band flux was ∼6.4 × 10−12 erg cm−2 s−1, a drop of ∼18% from its peak value (∼7.8 × 10−12 erg cm−2 s−1). In the meantime, the hard-band flux kept increasing over the last few years, with a latest value of ∼1.5 × 10−12 erg cm−2 s−1. As a result, the hard-to-soft flux ratio continued climbing up as before (see Figure 5).
3.4. Fe K Lines
We note that although the global spectra of SNR 1987A can be approximately characterized with a two-temperature vnei or vpshock model, both of these models leave some residuals around 6.6 keV for the EPIC-pn spectra, corresponding to Fe K line emission (see Figure 4). The detection of Fe K lines has been reported by several studies (e.g., Sturm et al. 2010; Maggi et al. 2012), but their origin remains unclear. Here, we present a comprehensive investigation of the Fe K lines based on the XMM-Newton observations.
We extracted the 5.0–8.0 keV EPIC-pn spectra, and fitted them with a combination of a thermal bremsstrahlung continuum (bremss) and a Gaussian profile (Gauss). As shown in Figure 6, the Fe K line gradually emerged from the underlying continuum. With the Gaussian component, we constrained the line fluxes and the centroids of the Fe K lines (listed in Table 5). There are several approaches to evaluating the significance of the lines. As listed in Table 5, we first provided the measured line fluxes and their 1σ uncertainties as estimates. Following Maggi et al. (2012), we also performed F-tests and calculated the p-values. Lastly, we introduced the Akaike information criterion (AIC; Akaike 1974) as a third approach. We derived the AIC value as follows:
where n is the number of measurements, p is the number of free parameters, yi denotes the observations, mi denotes the model-predicted values, and ωi denotes the weight factors. The more appropriate model is the one with the smaller AIC value. Therefore, we calculated the differences in AIC values between the models with and without a Gaussian component, and listed them in Table 5 as ΔAIC. Based on the significance analysis, we find that the first few observations taken in 2007–2009 show only modest (1σ–2σ) detections. But since 2010, we have continuously detected Fe K lines with ≳3σ significance—the less significant detection in 2019 is likely caused by the rather short effective exposure time of ∼11 ks. Despite a few less significant detections, we find that both the flux and the centroid energy of the Fe K lines have been increasing during the last decade, as shown in Figure 7. Assuming that the flux and the centroid energy are constant leads to reduced chi-square values χr 2 ≈ 4.05 and χr 2 ≈ 3.36, respectively. Both values of are clearly inconsistent with a constant flux and centroid energy. This indicates that there is an increase in the emission measure and the average ionization state of the Fe K emitting plasma as a function of time. However, the equivalent width (EW) of the line, which is more indicative of the Fe abundance, shows no significant variation during the period.
Table 5. Fitting Results of the Fe K Line
|Centroid (keV)||Flux (10−6 photons cm−2 s−1)||EW (eV)||(/dof)||F-test p-value||Significance|
|2007 Jan||235||0.97 (578/595)||−2.9||≳2σ|
|2008 Jan||207||1.04 (618/595)||−4.3||≳2σ|
|2009 Jan||93||1.07 (637/595)||0.24||1.1|
|2009 Dec||169||1.07 (637/595)||−6.3||≳2.5σ|
|2010 Dec||278||1.10 (653/595)||−13.6||≳3.5σ|
|2011 Dec||180||1.05 (625/595)||−9.4||≳3σ|
|2012 Dec||157||1.19 (709/595)||−5.7||≳2.5σ|
|2014 Nov||226||1.09 (646/595)||−17.6||≳4σ|
|2015 Nov||231||1.03 (609/595)||−19.0||≳4σ|
|2016 Nov||263||1.11 (662/595)||−23.5||≳4.5σ|
|2017 Oct||243||1.04 (619/595)||−8.8||≳3σ|
|2019 Nov||99||0.99 (587/595)||0.37||2.0||∼1σ|
Note. Errors are for the 1σ confidence range.
4.1. Temporal Evolution of the Plasma
Since the SN blast wave impacted the main body of the ER, the X-ray emission from SNR 1987A has been dominated by the shocked gas inside the ring. The XMM-Newton observations presented in this work started after the main impact and covered a time interval over 12 yr, which helped us to follow the post-impact evolution and provide further insights into the composition and structure of the ER.
The electron temperatures of the two plasma components show no significant variation during the last decade (kTe,LT ∼ 0.6 keV and kTe,HT ∼ 2.5 keV), indicating basically constant shock velocities. Given the classic relation between the shock velocity and the average post-shock temperature (where μ is the average particle mass in units of the proton mass; μ ∼ 0.72 for SNR 1987A abundance), we obtained Vs,LT ∼ 650 km s−1 and Vs,HT ∼ 1330 km s−1 for the low- and high-temperature plasma, respectively. On the other hand, Frank et al. (2016) constrained the expansion velocities of the ER as ∼1850 km s−1 for the soft band (0.5–2 keV) and as ∼3070 km s−1 for the hard band (2.0–10 keV) by fitting spatial models to Chandra ACIS images, which also provided estimates of the typical shock velocities, generally higher than what we find here. This difference might further reflect the nonequilibrium of the temperatures between different particle species behind the collisionless shock front. In that case, the electron temperature could be much lower than the proton/ion temperatures or the average temperature of the plasma (e.g., Vink 2020), and thus cause an underestimate of the shock velocities in the calculations above. Evidence for this temperature nonequilibrium effect has been found by both simulations and observations (e.g., Cargill & Papadopoulos 1988; Shimada & Hoshino 2000; Vink et al. 2003; Ghavamian et al. 2007). Recently, Miceli et al. (2019) measured the post-shock temperatures of protons and ions in SNR 1987A, and found that the ratio of ion temperature to proton temperature is significantly higher than 1, providing observational clues to the temperature nonequilibrium in this young remnant.
We find that the upper limit on the ionization parameter of the low-temperature plasma (τu,LT) increased linearly from 7300 to 10,500 days. A linear fit to τu,LT suggests an average density of ne,LT = 2450 ± 150 cm−3. By using a similar method, Sturm et al. (2010) fitted a linear function to τu,LT from 2007 to 2009 (7300–8000 days, based on the first three XMM-Newton observations listed in Table 1), which yielded a density of 3240 ± 640 cm−3. They also found that for the high-temperature plasma, τu,HT increased steadily from 2003 to 2009 (6000–8000 days). However, our results indicate that after 2009, τu,HT stayed constant at around 2 × 1011 cm−3 s, and showed no further evolution. Our fit to τu,LT also indicates that most of the low-temperature gas was first shocked at around 6000 days after the explosion, which gives us an estimate of the date of impact on the main ER. This is consistent with previous results based on X-ray light-curve analysis and expansion velocity studies (e.g., Racusin et al. 2009; Helder et al. 2013; Frank et al. 2016).
The normalization parameters of the two plasma components show different variations during the last decade. For the low-temperature plasma, the normalization NormLT increased in the first few years, leveled off at around 9000 days, and then started to decrease steadily at ∼10,500 days. On the other hand, NormHT kept increasing during the whole interval, with an approximately constant rate. Linear fits to the data gave a decreasing rate of NormLT (after 10,500 days) of (0.54 ± 0.03) × 10−3 cm−5 yr−1, and an increasing rate of NormHT of (0.29 ± 0.01) × 10−3 cm−5 yr−1.
With the spectral analysis, we have also constrained the abundances of various metal species, including N, O, Ne, Mg, Si, S, and Fe, as shown in the left panels of Figure 9. We note that the metal abundances, which are tied between two plasma components in our spectral fitting, should be taken as the (weighted) average values. We also calculated the abundance ratios with respect to Ne (Ne abundances are best constrained among all species), as shown in the right panels of Figure 9. In order to investigate their variations during the period, we derived the mean abundances and abundance ratios, and performed χ2 analysis based on these mean values (shown in Figure 9). We note that for a p-value of 0.002 (∼3σ level), the critical reduced chi-square is given as χ2 ν,c ≈ 2.67 for 11 degrees of freedom. Thus the χ2 analysis indicates that the abundances of N, O, Ne, Mg, and Fe, as well as the N/Ne and Fe/Ne ratios, have significantly changed in the last decade. The variations of the N, O, Ne, and Mg abundances follow a similar track, which stayed at their mean values until around 10,500 days, and then started to decrease recently. This is well fitted with the recent decline of the low-temperature plasma normalization after ∼10,500 days. The increasing Fe abundance conforms to the increasing high-temperature plasma normalization. As for the abundance ratios, N/Ne was decreasing while Fe/Ne was increasing during the time, consistent with the variation of the single abundances. Therefore, the overall changes of the single metal abundances and their ratios could be naturally explained with the different chemical compositions for the two plasma components. In our spectral fitting, the abundances of N, O, Ne, and Mg were determined mainly by their emission lines lying in the 0.5–2 keV energy band, which were dominated by the low-temperature component at the beginning. From the recent decline of NormLT, the high-temperature plasma started to contribute a significant part of the emission even in the soft band, and thus the abundances became more representative of the high-temperature plasma. In that case, our results indicate that the low-temperature plasma has enhanced metal abundances as compared to the high-temperature plasma, especially for N, O, Ne, and Mg, while the high-temperature plasma may be more abundant in Fe. We will further discuss the possible origin of metal differences in Section 4.3. It is also possible that the increasing Fe abundance and Fe/Ne ratio result from the recent emergence of Fe-rich ejecta. However, this needs to be further examined, and we leave the discussions to Section 4.4. Lastly, we note that in a recent paper, Bray et al. (2020) studied the metal abundances of SNR 1987A based on Chandra HETG data (taken in 2011–2018) and found no significant variations during the period. However, the spectra they adopted were limited to 0.5–3 keV; thus the hot component might not be well constrained. Specifically, they obtained a temperature of 1.2–1.6 keV for the hot plasma, while our results and many other works suggest ∼2.5 keV.
4.2. Electron Density, Ionization Age, and Filling Factor
In this section, we estimate the electron density, the ionization age, and the filling factor of the X-ray-emitting gas based on the spectral fitting results.
Assuming that the shocked ER has a toroidal geometry, its volume can be derived as , where R and r are the angular radius and the half-width of the torus, respectively, and d is the distance to SNR 1987A. We then took the total X-ray-emitting volume as V = ftot VER, where ftot is the total filling factor. With this volume given, the electron density can be derived from the normalization parameter , where nH ≈ ne/1.2 for a fully ionized plasma with near solar abundance. The ionization age t can then be obtained from τ = ne t. By further assuming that the two plasma components are in pressure balance (i.e., ne,LT kTe,LT = ne,HT kTe,HT) and share the whole X-ray-emitting volume (i.e., ftot = fLT + fHT), we were able to constrain these parameters separately for each component.
The geometry parameters R and r of the shocked ER in different epochs have been constrained by fitting a "torus + lobes" spatial model to the deconvolved Chandra images (e.g., Racusin et al. 2009; Frank et al. 2016). Here, we adopt a half-width r = 022 for all the observations according to Racusin et al. (2009), who found a nearly constant ring width after ∼6000 days. On the other hand, the average radius R increases due to the expansion of the blast wave. For each observation, we set R to the value obtained by Frank et al. (2016) using Chandra data taken in a similar epoch (for observations taken after 2015, R was estimated by an extrapolation based on the expansion velocity). We adopted a distance to SNR 1987A of d = 51.4 kpc based on Panagia (2005). We used the vpshock fitting results for calculations. As for the ionization ages, we estimated both the average t from the average τ obtained in the vnei fit, and the upper limit tu from the τu obtained in the vpshock fit. The derived electron densities, average/maximum ionization ages, and filling factors are listed in Table 6.
Table 6. Derived Parameters Based on Spectral Fitting Results
|2007 Jan||2008 Jan||2009 Jan||2009 Dec||2010 Dec||2011 Dec||2012 Dec||2014 Nov||2015 Nov||2016 Nov||2017 Oct||2019 Nov|
|ne,LT ( cm−3)|
|tLT ( days)|
|tu,LT ( days)|
|ne,HT ( cm−3)||629 ± 35||646 ± 38||650 ± 31|
|tHT ( days)||1197 ± 103||1149 ± 100|
|tu,HT ( days)|
The electron densities of the two plasma components are derived to be ne,LT ∼ (1600– cm−3 (with an average value of cm−3) and ne,HT ∼ cm−3 (with an average value ∼ cm−3). The unknown filling factor ftot leads to some uncertainties in our results. Ideally, it should be around 1 if the total X-ray-emitting volume is comparable to the shocked ER volume. However, it will be significantly lower than 1 if the X-ray-emitting gas is clumpy and concentrates in a small part of the shocked ring. It could also go beyond 1 if the hot, lower-density gas extends to a larger region outside the ring. We find that the derived gas densities appear to be increasing during the period, which may be indicative of a density gradient throughout the ring. But this is still debatable due to uncertainties in the time-variable filling factor and the assumed geometry. As mentioned above in Section 4.1, we get a density of ne,LT = 2450 ± 150 cm−3 for the low-temperature plasma by fitting a linear function to τu,LT. For comparison, Sturm et al. (2010) obtained average densities of ne,LT = 3240 ± 640 cm−3 and ne,HT = 480 ± 140 cm−3 for the two components from a similar analysis of the 2003–2009 XMM-Newton observations. On the other hand, expansion velocity studies suggest a shock velocity of ∼6800 km s−1 before the impact, while it dropped to ∼1850 km s−1 as viewed in the 0.5–2 keV band and to ∼3070 km s−1 in the 2–10 keV band after 6000 days (e.g., Frank et al. 2016). This requires that the density increase by a factor of ∼13 and ∼5 (given the relation between the shock velocity and the gas density of ). Given that the H ii region has a density ∼102 cm−3, this indicates an average density ≳103 cm−3 for the soft-band-emitting gas and ∼500 cm−3 for the hard-band-emitting gas. Taking all these results together, we suggest ne,LT ∼ 2400 cm−3 and ne,HT ∼ 550 cm−3 are reasonable estimations of the low- and high-temperature plasma densities.
The ionization age characterizes the elapsed time since the gas was shocked. The average ionization age t and its upper limit tu were constrained based on the vnei and vpshock spectral fitting results, respectively. For the low-temperature plasma, we find tu,LT gradually increased from days to days during the period of 7300–10,500 days after explosion, consistent with the physical picture that the ER started to be shocked by the blast wave at ∼6000 days and dominated the soft-band X-ray emission since then. Under the plane-parallel shock assumption, we would expect the average ionization age to be tLT = tu,LT/2 ∼ (1000– days. However, it is found to be significantly lower, as indicated by the vnei fitting, with tLT ∼ (300– days, and does not show evidence of an increase similar to that of tu,LT. This may suggest that the low-temperature component has always been dominated by the "youngest" hot gas, which has just been shock-heated within the most recent 1–2 yr. This very dynamic shock-heating process is also indicated by Chandra imaging observations. As shown by Frank et al. (2016, Figure 5 therein), the surface X-ray brightness distribution of the ER can significantly change within a few hundred days, which means the brightening and the fading of the gas take place on a rather short timescale. This short average ionization age may be related to the shock transmission timescale and/or the destruction timescale of the shocked high-density clumps. However, it should be recalled that the vnei model fitting optimizes for the line emission, and the line emissivity has a nonlinear relation with the ionization state. As a result, the line emission is suppressed when the gas is close to equilibrium at these temperatures, so the vnei parameters can be skewed to lower τ values. We have also considered whether the short timescale could reflect the radiative cooling timescale, which can be estimated as τcool ≈ 5.7k T/(neΛ(T)), where Λ is the volumetric cooling rate (e.g., Vink 2020). For the low-temperature plasma, taking , ne,LT ∼ 2400 cm−3, and Λ(T) ∼ 10−22 erg s−1 cm−3 (e.g., Schure et al. 2009), we obtain a cooling timescale of τcool,LT ∼ 720 yr, which makes radiation cooling unlikely to be the cause of the short average ionization age. The optical emission is likely to come from even denser clumps (≫104 cm−3), for which the post-shock temperature is in a regime for which τcool is an order of magnitude smaller, with timescales of months to years.
For the high-temperature plasma, we obtain tu,HT ∼ (2800– days (with an average value ∼ days), which shows no significant variation during this period. The average ionization age is derived to be tHT ∼ (1000– days—with an average value ∼ days—consistent with the expected value under the plane-parallel shock assumption.
The filling factors have also been constrained for the two plasma components. The ratio fLT/fHT, which no longer depends on the uncertain total filling factor ftot, is found to have been ∼0.2 during 7300–9000 days, and it has been continuously decreasing since then. By the latest observation at ∼12,000 days, it has dropped to ∼0.04.
With the derived electron densities and filling factors, we also estimated the gas masses for each of the components, listed in Table 6.
4.3. Origin of the X-Ray-emitting Plasma
Similar to many previous studies, our work reveals that the X-ray spectra of SNR 1987A can be approximately described by the thermal emission from a combination of two major plasma components, one with a low temperature, ∼0.6 keV, and the other one with a high temperature, ∼2.5 keV. However, the physical origins of these plasmas could be rather complicated due to the complex CSM structure and the shock system. The main structure of the shocked CSM includes the ER, an ionized H ii region surrounding the ring, and the high-latitude material beyond the ring. The ER consists of several dense clumps ("fingers") distributed mainly in its inner portion and smooth interclump gas. It has been suggested that the low-temperature plasma (or the soft X-ray emission) comes from the slowly transmitted shock into the dense clumps, while the high-temperature plasma (or the hard X-ray emission) mainly comes from less decelerated shocks (which may also include reflected shocks) moving through the lower-density ring material. This physical picture was first indicated by observations taken in the early stage of the shock–ring interaction (e.g., Park et al. 2006; Zhekov et al. 2010), and later supported by 1D to 3D HD simulations (e.g., Dewey et al. 2012; Orlando et al. 2015). However, our analysis is based on more recent observations (after ∼7500 days), and thus represents the later-stage evolution. Keeping this in mind, we revisit here the origin of the X-ray-emitting plasma.
Based on early optical analysis of ionized gas in the ER, the dense clumps have a density ∼(1–3) × 104 cm−3, and the interclump gas has a density ∼(1–6) × 103 cm−3 (e.g., Lundqvist & Fransson 1996; Mattila et al. 2010). Similar densities are also indicated by HD simulations (e.g., Dewey et al. 2012; Orlando et al. 2015). We find that the low-temperature plasma has a density of ne,LT ∼ 2400 cm−3, and thus it is more indicative of lower-density gas inside the ring. It may consist of both shocked interclump gas and evaporated clump material (or fragments with lower density). On the other hand, the high-temperature plasma has a density of ne,HT ∼ 550 cm−3, which is lower than the expected density of the ring material. However, taking into account the shock compression ratio of ∼4, it is consistent with the interpreted density of the H ii region (∼102 cm−3, e.g., Lundqvist & Fransson 1996; Mattila et al. 2010; Orlando et al. 2015).
The shocked dense clumps, which have a typical density of ≳104 cm−3, are expected to contribute a rather low temperature (≲0.1 keV, given the density ratio between the dense clump and the interclump material) plasma component. However, we have not found solid evidence for such a component in our spectral analysis. We note that at the early stage of the shock–ring interaction, the soft (low-temperature) component is found to be close to collisional ionization equilibrium with a much lower electron temperature ∼0.22 keV, the density of which was estimated as ∼7500 cm−3 (Park et al. 2004). As the blast wave penetrated further into the ER, the temperature of this soft component gradually increased to ≳0.5 keV (e.g., Park et al. 2006; Sturm et al. 2010). This indicates that the soft emission was dominated by the shocked dense clumps at the beginning, but was later overtaken by the lower-density gas inside the ring. The dense clumps themselves could have faded out due to evaporation/destruction, or they contribute little to the overall emission, which can hardly be detected.
The total mass of the ER is estimated to be ∼0.06 M⊙ in optical observations and HD simulations (e.g., Mattila et al. 2010; Orlando et al. 2015), and most of the mass is distributed in the lower-density component (∼60%–80%). Despite the large uncertainty of our mass estimate, it provides hints with regard to the plasma origin. We find a peak mass value for the low-temperature plasma ∼ M⊙, which is similar to the total mass of the ER. On the other hand, the mass of the high-temperature plasma has been continuously increasing, and has been found to be ≳ M⊙ in the latest observations, which exceeds the expected ER mass by a factor of 2. Therefore, the low-temperature plasma could be dominated by shocked ring material, while the high-temperature plasma should also contain the material around/beyond the ring, which may extend to a larger radius and has been continuously shocked by the blast wave.
As described in Section 4.1, the metal abundances of N, O, Ne, and Mg have significantly declined in the past few years, which may have resulted from the different chemical compositions between the two plasma components. In that case, it is more reasonable to interpret the two plasma components as representing the different parts of the CSM structure. The ER is found to have enhanced He and N abundances according to optical observations (e.g., Lundqvist & Fransson 1996; Mattila et al. 2010). As the enhancement in He and N implies a processing of hydrogen, it affects the abundances of alpha- and Fe-group elements with respect to hydrogen. The H ii region, as indicated by early X-ray observations taken before the impact, has relatively lower abundances of N, O, Ne, and Mg (Park et al. 2002). In a hydrodynamics-based analysis of X-ray spectra, Dewey et al. (2012) found different metal abundances in the ER and the H ii region, where the ER had generally higher abundances for most of the species. Our results suggest that the low-temperature plasma dominated the soft-band emission until ∼10,500 days, and resulted in enhanced average abundances of N, O, Ne, and Mg during this period. Given the decline of the low-temperature emission measure, the average abundances of N, O, Ne, and Mg started to decrease in the last few years, which in turn indicated lower metal abundances for the high-temperature plasma. Therefore, the low-temperature plasma could be more representative of the ER, while the high-temperature plasma could be more representative of the H ii region and the high-latitude material with low metal abundances.
Given all the discussions above, our results indicate that the low-temperature plasma is dominated by the lower-density gas inside the ER, which may include the smooth interclump gas and evaporated/destroyed clump material, while the high-temperature plasma is dominated by the H ii region and the low-density gas around/beyond the ring. The recent decrease in the low-temperature emission measure indicates that the blast wave has left the main ER, while it is still propagating into the high-latitude material, resulting in the continuously increasing high-temperature emission measure. Such a physical picture is also reflected by the drop of the filling factor ratio fLT/fHT. However, the real CSM structure must be rather complicated, and the actual density and temperature distributions cannot be fully represented by a simple two-component model (below we discuss the need for at least one additional component to fit the Fe K line emission). Alternative approaches could be differential emission measure analysis based on high-resolution spectra (e.g., Zhekov et al. 2006, 2009), and detailed full-scale HD/MHD simulations (e.g., Orlando et al. 2015, 2019, 2020).
4.4. Possible Origins of the Fe K Lines
The EPIC-pn camera reveals Fe K features in SNR 1987A's spectra, the origin of which remains unclear. Recent observations indicate an Fe K centroid energy ≳6.65 keV. Comparison with the theoretically predicted line centroid shown in Figure 10 (based on SPEX code 13 ; Kaastra et al. 1996) reveals that the high centroid energy requires either a high plasma temperature ≳2.0 keV or a large ionization parameter ≳1011 cm−3 s. Since the Fe K line emission is not well reproduced with the standard two-temperature-component scenario, we invoked a third, hot plasma component as the origin of the Fe K lines. Taking the 2015 observation as an example, we refit the RGS and pn spectra with a three-temperature vnei model. It provides an overall improved fit ( versus for a two-component fit, or Δχ2 = 320 for three additional degrees of freedom) to the data and reproduces Fe K lines successfully. The Fe K features are characterized by a high-temperature plasma component with keV and τ = (1.4 ± 0.1) × 1011 cm−3 s. The other two components have lower temperatures and similar ionization parameters: keV and τ = (3.2 ± 0.1) × 1011 cm−3 s for the low-temperature one, and keV and cm−3 s for the mid-temperature one. On the other hand, the temporal variation of the Fe K centroid energy provides further information about the plasma properties. Given a certain temperature and density of the plasma, the relation shown in Figure 10 results in a modeled curve of the Fe K centroid as a function of time, which can be compared with observations. As shown in Figure 7, the observed Fe K centroid energies basically match the predicted evolution of a plasma with kTe ∼ 3.2 keV and ne ∼ 500 cm−3, shocked at 7000 days after the explosion. This Fe K emitting component may therefore also come from the low-density gas around/beyond the ER, albeit with an even higher temperature. It may be related to the reflected shock, which further compresses and heats the already shocked gas (e.g., Zhekov et al. 2009, 2010). It is also possible that the Fe K emission comes from one or several reverse-shocked Fe-rich ejecta clumps, as a reverse shock with high ejecta-frame velocity could result in high post-shock temperature. Although most of the Fe-rich ejecta may concentrate in the inner part of the ejecta material, there could still be some high-speed debris that has already been approached by the reverse shock. We note that the gamma-ray data taken a few months after the SN explosion provided surprising evidence for fast core-ejecta material mixing into the outer envelopes (see Leising & Share 1990 and references therein).
Since entering the remnant phase, SNR 1987A has been dynamically evolving on timescales of months to years. In this work, we utilized the XMM-Newton observations of SNR 1987A from 2007 to 2019 to investigate its post-impact evolution. The high energy resolution of RGS helped us identify and measure individual emission lines, and the EPIC-pn camera provided us with an overall energy coverage from 0.3 to 10 keV and thus with a better constraint of the hard-band emission. We performed detailed modeling to the spectra; constrained the electron temperatures, the ionization parameters, and the emission measures of the plasma; and tracked their recent evolution. We summarize the main results and conclusions below.
- 1.The post-impact X-ray emission of SNR 1987A can be approximately described with a model containing two NEI plasma components with different electron temperatures (kTe,LT ∼ 0.6 keV and kTe,HT ∼ 2.5 keV). The plasma temperatures showed no significant variations during the last decade.
- 2.The emission measure of the low-temperature plasma started to decrease around 10,500 days after the explosion, resulting in a recent decline of the soft-band (0.5–2.0 keV) flux. On the other hand, the high-temperature emission measure, as well as the hard-band (3.0–8.0 keV) flux, kept increasing with a steady rate.
- 3.The average abundances of N, O, Ne, and Mg have decreased in the last few years, which could have resulted from the different chemical compositions between the two plasma components.
- 4.The densities of the X-ray-emitting gases were estimated as ne,LT ∼ 2400 cm−3 and ne,HT ∼ 550 cm−3 for the low-temperature and the high-temperature components, respectively.
- 5.The low-temperature plasma is indicated to be dominated by the lower-density gas inside the ER, which may include the interclump gas and evaporated/destroyed clump material. The high-temperature plasma could be dominated by the H ii region and the low-density gas around/beyond the ring. The blast wave has now left the main ER, but is still propagating into the high-latitude material.
- 6.Fe K line emission has been detected in most of the observations. Its centroid energy has been increasing during the last decade, with a latest value of ≳6.65 keV. The Fe K lines may originate from a third, high-temperature component, which could be the reflected shock-heated CSM, or be related to the reverse-shocked Fe-rich ejecta.
SNR 1987A has now entered a new evolutionary stage in which the blast wave is starting to shock the high-latitude material beyond the ER. On the other hand, as the reverse shock penetrates deeper into the inner ejecta, the shocked ejecta will start to contribute to the X-ray emission as well. Further observations of the remnant in the next few years will be crucial to monitoring these new changes and unveiling the underlying physics.
This work is partly supported by the National Key R&D Program of China (No. 2017YFA0402600) and the NSFC under grants 11773014, 11633007, and 11851305. This project has also received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 101004131. L.S. acknowledges the financial support of the China Scholar Council (No. 201906190108).
During the refereeing process of the present paper, two new analyses regarding the X-ray emission from SN 1987A were presented (Alp et al. 2021; Greco et al. 2021), incorporating NuSTAR hard X-ray spectra. Both papers agree that there is a need for a third hard X-ray component, but the Alp et al. (2021) paper (based on XMM-Newton RGS + NuSTAR data) attributes this third component to thermal emission, and the Greco et al. (2021) paper (based on Chandra + NuSTAR data) attributes it to nonthermal emission from a heavily absorbed emerging pulsar wind nebula. The results presented by us are indicative of a third, high-temperature thermal component, and thus qualitatively agree with Alp et al. (2021). However, they found a temperature as high as ∼4 keV is needed in order to fit the NuSTAR data, which is significantly higher than what we get in Section 4.4. In fact, the temperature and ionization parameter we get for the third component are closer to those obtained by Greco et al. (2021) (see also the analysis by Tsuji et al. 2021) for their high-temperature thermal component. We note that our decision to include a third thermal component was based exclusively on the soft X-ray data (below 10 keV), specifically the detection of Fe K emission, which could not be reproduced by the conventional two-component model. On the other hand, our analysis, based solely on XMM-Newton observations, cannot rule out or confirm the presence of a highly absorbed nonthermal component. We verified that adding an absorbed power-law component as described in Greco et al. (2021) does not affect much the soft X-ray spectrum (producing no significant X-ray emission below 7 keV), and thus makes no significant improvement to the fit. It remains then to be investigated whether a fourth component representing nonthermal emission is still required if we would like to augment our analysis by incorporating hard X-ray data. But it is beyond the scope of this paper to make strong conclusions regarding this point.
SNR 1987A is a pointlike source as viewed by EPIC-pn. However, in view of the point-spread function (FWHM ∼ 125), a circular region with a radius of 25'' can get most of the source photons involved while avoiding large background contamination.
As previously pointed out by Dewey et al. (2012), there could be a measurable absorption component local to SNR 1987A: a density of 1 × 104 H cm−3 and a path length of 0.01 pc would result in a local absorption NH,local ∼ 0.3 × 1021 cm−2. This local absorption is also affected by clumping, the ionization state, and the geometry of the remnant, and thus may change as the system evolves. By setting the LMC absorption as a free parameter among different epochs, we obtained NH,LMC varying in the range of 1.7–2.5 × 1021 cm−2, with a rather monotonically increasing trend. This may reflect the local absorption, but will introduce additional parameter degeneracy, especially between NH, kTe, and the normalization. On the other hand, the changes and calibration uncertainties in the RGS effective area may also affect the measurement of NH. As a result, we decided to keep NH,LMC fixed to its average value (2.2 × 1021 cm−2) in our spectral fitting, and actually this will not significantly affect our conclusions.