The Post-impact Evolution of the X-ray Emitting Gas in SNR 1987A Viewed by XMM-Newton

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 non-equilibrium ionization components with temperature of $\sim0.6$ keV and $\sim2.5$ keV. The low-temperature plasma has a density $\sim2400$ cm$^{-3}$, which is likely dominated by the lower density gas inside the equatorial ring (ER). The high-temperature plasma with a density $\sim550$ 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 steadily increase of the high-temperature emission measure. In the meantime, the average abundances of N, O, Ne, and Mg are found to be declining, which may reflect the different chemical compositions between two plasma components. We also detected the Fe K lines in most of the observations, showing increasing flux and centroid energy. We interpret the Fe K lines as 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.


INTRODUCTION
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 supernova remnant (SNR) formation and subsequent evolution for an SN whose observational properties are known in detail. SN1987A is located in the Large Magellanic Cloud (LMC) and detected on February 23, 1987. 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 colliding with material from an earlier slow wind, during the progenitor's post main sequence evolution (Chita et al. 2008). Today, over 30 years after the SN outburst, the blast wave has encountered and heavily interacted with the CSM structures, especially 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 during 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 first 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 three-dimensional hydrodynamic/magnetohydrodynamic simulations (e.g., Orlando et al. 2015Orlando et al. , 2019Orlando et al. , 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 years 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 revealed new features of spot-like and diffuse emission outside the ER (since ∼ 9500 days), 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)(2017)(2018) enhancement of the GeV emission from the SNR 1987A region has been reported based on Fermi/LAT observations (Malyshev et al. 2019). The X-ray observations could be crucial to clarify the nature of the GeV signal (a detection of the non-thermal 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 the further implications of the results, and Section 5 gives a brief summary. Recently, we obtained a new observation in November 2019 (PI: Malyshev, for the H.E.S.S. collaboration), aiming at its latest evolution. Thereby, we got a set of data containing 12 observations in total, covering a time interval over 12 years, 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 to the plasma properties. We also included the EPIC-pn spectra for a complete energy coverage to 10 keV following Sturm et al. (2010).
We used XMM-Newton Science Analysis Software (SAS, version 16.1.0) 1 to process the data. For 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 (SP) 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 2 , 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 is only minor contribution to the systematic error budget. We selected only the single events (PATTERN = 0) for the source spectra in order to get rid of the photon pile-up 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 of the EPIC-pn and RGS spectra were rebinned using FTOOL 3 task grppha to obtain at least 25 counts per bin.

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) 4 with AtomDB 3.0.9 5 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).

Individual emission lines and their fluxes
Taking advantage of the high spectral resolution of RGS, we are 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 an heuristic model and fit it to the first-order and the second-order RGS1/RGS2 spectra simultaneously.
The spectral model consists of a thermal continuum component (described by the nlapec model in XSPEC, 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 N H, Gal = 6 × 10 20 cm −2 with abundances set to be solar, and the LMC absorption was fixed at N H, LMC = 2.2 × 10 21 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 this 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 (absorptioncorrected) 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 plotted 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][2012][2013][2014][2015], and after that 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-year interval, the Lyα/Heα flux ratios have been continuously increasing. This could result from the increase in the average plasma temperature, but may also be affected by non-equilibrium ionization (NEI) effects.
Our spectral analysis also provides measurements of the line broadening and the systemic velocity, see σ 6 and v 87A 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 ∼ 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 v 87A 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 also lead to the large dispersion of the values we obtained here for v 87A . 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 ). 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 which 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 larger 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.

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.

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 kT e and a single ionization parameter τ = n e t, respectively, where n e 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 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 N H, Gal = 6 × 10 20 cm −2 and N H, LMC = 2.2 × 10 21 cm −2 , respectively 6 . Similar to the model described in Section 3.1, the spectra were convolved with 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 observation, with reduced chi-square values of χ 2 r ∼ 1.21-1.41 (dof ∼ 1234-2433). The low-temperature component is characterised by kT e,LT ∼ 0.6 keV and the high-temperature component has kT e,HT ∼ 2.5 keV, which show no significant variations during the last decade. The detailed fitting results are summarized in Table 3.

Two-temperature vpshock model
In addition to the vnei model presented above, we also modeled the global spectra of SNR 1987A spectra 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 the 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 emission measure, instead of assuming a single ionisation parameter. The NEI plasma in the vpshock model is characterized by a lower limit τ l = n e t l at the shock front and an upper limit τ u = n e t u 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 suggested that the single ionization state scenario in vnei model is over-simplified, 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

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 recent changes in 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 10500 days, the soft band flux started to decline. In the latest observation (11964 days), the soft band flux was ∼ 6.4 × 10 −12 erg cm −2 s −1 , a drop of ∼ 18% from its peak 6 As previously pointed out by Dewey et al. (2012), there could be a measurable absorption component local to SNR 1987A: a density of 1 × 10 4 H cm −3 and a path length of 0.01 pc would result in a local absorption N H, local ∼ 0.3 × 10 21 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 obtain N H, LMC varying in the range of 1.7-2.5 × 10 21 cm −2 , with a rather monotonic increasing trend. This may reflect the local absorption, but will introduce additional parameter degeneracy, especially between N H , kT e , and the normalization. On the other hand, the changes and the calibration uncertainties in the RGS effective area may also affect the measure of N H . As a result, we decide to keep N H, LMC fixed to its average value (2.2 × 10 21 cm −2 ) in our spectral fitting, and actually this will not significantly affect our conclusions.
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).

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 two 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 the Fe K lines has been reported by several studies (e.g., Sturm et al. 2010;Maggi et al. 2012), but its 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 out 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 evaluate 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. At last, we introduced the Akaike information criterion (AIC, Akaike 1974) as a third approach. We derived the AIC value as: where n is the number of measurements, p the number of free parameters, y i the observations, m i the model predicted values, and ω i the weight factors. The more appropriate model is the one with the smaller AIC value. Thereby, we calculated the differences in AIC values between the models with and without the Gaussian component, and listed them in Table 5 as ∆AIC. Based on the significance analysis, we found that the first few observations taken in 2007-2009 show only modest (1-2σ) detections. But since 2010, we have continuously detected the 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 of 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 χ 2 r ≈ 4.05 and χ 2 r ≈ 3.36, respectively. Both values of χ 2 r 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 function of time. However, the equivalent width of the line, which is more indicative of the Fe abundance, shows no significant variation during the period.

Temporal evolution of the plasma
Since the SN blast wave impacted on 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 years, which helps us to follow the post impact evolution and provide further insights into the composition and structure of the ER.
With the vpshock fitting (Section 3.2.2), we have constrained the electron temperatures, the upper limits of the ionization parameters, and the normalizations of the gas, as shown in Figure 8.
The electron temperatures of the two plasma components show no significant variation during the last decade (kT e,LT ∼ 0.6 keV and kT e,HT ∼ 2.5 keV), indicating basically constant shock velocities. Given the classic relation between the shock velocity and the average post shock temperature as V 2 s = (16/3)(kT )/(µm p ) (where µ is the average particle mass in units of the proton mass, µ ∼ 0.72 for SNR 1987A abundance), we obtained V s,LT ∼ 650 km s −1 and V s,HT ∼ 1330 km s −1 for the lowand 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 ∼ 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 and are generally higher than what we found here. This difference might further reflect the non-equilibrium 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 non-equilibrium effect have 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 one, providing observational clues to the temperature non-equilibrium in this young remnant.
We found that the upper limit on the ionization parameter of the low-temperature plasma (τ u,LT ) has been increasing linearly from 7300 to 10500 days. A linear fit to τ u,LT suggests an average density of n e,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 to 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 hightemperature plasma, τ u,HT has been increasing steadily from 2003 to 2009 (6000 to 8000 days). However, our results indicate that after 2009, τ u,HT stayed constant at around 2 × 10 11 cm −3 s, and show 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 the 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 Norm LT increased at the first few years, leveled off at around 9000 days, and then started to decrease steadily since ∼ 10500 days. On the other hand, Norm HT kept increasing during the whole interval, with an approximately constant rate. Linear fits to the data gave a decreasing rate of Norm LT (after 10500 days) as (0.54 ± 0.03) × 10 −3 cm −5 yr −1 , and an increasing rate of Norm HT as (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 ratio, have significantly changed in the last decade. The variations of N, O, Ne, and Mg abundances follow a similar track, which stayed at their mean values till around 10500 days, and then started to decrease recently. This is well fitted with the recent decline of the low-temperature plasma normalization after ∼ 10500 days. The increasing Fe abundance conforms to the increasing high-temperature plasma normalization. As for the abundance ratios, N/Ne has been decreasing while Fe/Ne been increasing during the time, consistent with the variation of single abundances. Therefore, the overall changes of the single metal abundances and their ratios could be naturally explained with the different chemical compositions for two plasma components. In our spectral fitting, the abundances of N, O, Ne, and Mg are determined mainly by their emission lines lying in 0.5-2 keV energy band, which were dominated by the low-temperature component at the beginning. Since the recent decline of Norm LT , 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 lowtemperature 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 the Fe/Ne ratio result from the recent emergence of Fe-rich ejecta. However, this needs to be further examined, and we leave the discussions in Section 4.4. At last, 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][2012][2013][2014][2015][2016][2017][2018] and found no significant variations during the period. However, the spectra they adopted were limited in 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 suggested ∼ 2.5 keV.

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 spectral fitting results.
Assuming that the shocked ER has a toroidal geometry, its volume can be derived as V ER = 2π 2 Rr 2 d 3 , 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 = f tot V ER , where f tot is the total filling factor. With this volume given, the electron density can be derived from the normalization parameter Norm = 10 −14 /(4πd 2 )n e n H V , where n H ≈ n e /1.2 for a fully ionized plasma with near solar abundance. The ionization age t can then be obtained from τ = n e t. By further assuming that the two plasma components are in pressure balance (i.e., n e,LT kT e,LT = n e,HT kT e,HT ) and share the whole X-ray emitting volume (i.e. f tot = f LT + f HT ), 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 = 0.22 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 is increasing due to the expansion of the blast wave. For each observation, we set R to be the value obtained by Frank et al. (2016) using the Chandra data taken in a similar epoch (for the observations taken after 2015, R were estimated by an extrapolation based on the expansion velocity). We adopt a distance to SNR 1987A of d = 51.4 kpc based on Panagia (2005). We use the vpshock fitting results for calculations. As for ionization ages, we estimated both the average t from the average τ obtained in the vnei fit, and the upper limit t u from τ u obtained in the vpshock fit. The derived electron densities, average/maximum ionization ages, and filling factors are listed in Table 6.
The electron densities of the two plasma components are derived to be n e,LT ∼ (1600- The unknown filling factor f tot led 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 got a density of n e,LT = 2450 ± 150 cm −3 for the low-temperature plasma by fitting a linear function to τ u,LT . As a comparison, Sturm et al. (2010) obtained the average densities as n e,LT = 3240 ± 640 cm −3 and n e,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 0.5-2 keV band and ∼ 3070 km s −1 in 2-10 keV band after 6000 days (e.g., Frank et al. 2016). This requires that the density increased by a factor of ∼ 13 and ∼ 5 (given the relation between the shock velocity and the gas density as v s ∝ n 0.5 e ). Given that the H II region has a density ∼ 10 2 cm −3 , this indicates an average density 10 3 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 n e,LT ∼ 2400 cm −3 and n e,HT ∼ 550 cm −3 as reasonable estimations of the lowand 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 t u were constrained based on the vnei and vpshock spectral fitting results, respectively. For the low-temperature plasma, we found t u,LT gradually increased from ∼ 2000f 0.5 tot days to ∼ 4800f 0.5 tot days during the period of 7300-10500 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 t LT = t u,LT /2 ∼ (1000-2400)f 0.5 tot days. However, it is found to be significantly lower, as indicated by the vnei fitting, with t LT ∼ (300-500)f 0.5 tot days, and does not show evidence for a similar increase as t u,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 remembered that the vnei model fitting optimizes for the line emission, and the line emissivity has a non-linear 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.7kT /(n e Λ(T )), where Λ is the volumetric cooling rate (e.g., Vink 2020). For the low-temperature plasma, taking kT e,LT ∼ 0.6 keV, n e,LT ∼ 2400 cm −3 and Λ(T ) ∼ 10 −22 erg s −1 cm −3 (e.g., Schure et al. 2009), we obtained 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 ( 10 4 cm −3 ), for which the post-shock temperature is in a regime for which τ cool is an order of magnitude smaller, with time scales of months to years.
For the high-temperature plasma, we obtained t u,HT ∼ (2800-4600)f 0.5 tot days (with an average value ∼ 3400f 0.5 tot days), which shows no significant variation during this period. The average ionization age is derived to be t HT ∼ (1000-1500)f 0.5 tot days -with an average value ∼ 1300f 0.5 tot days-consistent with the expected value under the plane-parallel shock assumption.
The filling factors have also been constrained for the two plasma components, respectively. The ratio f LT /f HT , which is no longer depending on the uncertain total filling factor f tot , is found to have been ∼ 0.2 during 7300-9000 days, and it has been continuously decreasing since then. By the latest observation at ∼ 12000 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, as listed in Table 6.

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 the smooth inter-clump gas. It has been suggested that the low-temperature plasma (or the soft X-ray emission) comes from the slow transmitted shock into the dense clumps, while the high-temperature plasma (or the hard X-ray emission) mainly comes from the less-decelerated shocks (which may also include reflected shocks) moving through the lower density ring material. This physical picture was first indicated by the 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 the 1D to 3D hydrodynamic simulations (e.g., Dewey et al. 2012;Orlando et al. 2015). However, our analysis is based on the more recent observations (after ∼ 7500 days), thus represents the later stage evolution. Keeping this in mind, we revisit here the origin of the X-ray emitting plasma.
Based on the early optical analysis of the ionized gas in the ER, the dense clumps have a density ∼(1-3)×10 4 cm −3 , and the inter-clump gas has a density ∼(1-6)×10 3 cm −3 (e.g., Lundqvist & Fransson 1996;Mattila et al. 2010). Similar densities are also indicated by hydrodynamic simula-tions (e.g., Dewey et al. 2012;Orlando et al. 2015). We find that the low-temperature plasma has a density of n e,LT ∼ 2400 cm −3 , and thus it is more indicative of the lower density gas inside the ring. It may consist of both the shocked inter-clump gas and the evaporated clump material (or the fragments with lower density). On the other hand, the high-temperature plasma has a density of n e,HT ∼ 550 cm −3 , which is lower than the expected density of the ring material. However, taking into account a shock-compression ratio of ∼ 4, it is consistent with the interpreted density of the H II region (∼ 10 2 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 10 4 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 (lowtemperature) component was 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 either 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 was estimated to be ∼ 0.06 M from optical observations and hydrodynamic 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 to the plasma origin. We found a peak mass value for the low-temperature plasma ∼ 0.08f 0.5 tot M , which is similar to the total mass of the ER. On the other hand, the mass of the hightemperature plasma has been continuously increasing, and was found to be 0.15f 0.5 tot M in latest observations, which exceeds the expected ER mass by a factor of 2. Thereby, the low-temperature plasma could be dominated by the 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 recent few years, which may result from the different chemical compositions between 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 was found to have enhanced He and N abundances according to the 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 respects to hydrogen. The H II region, as indicated by the early X-ray observations taken before the impact, has relatively lower abundances of N, O, Ne, and Mg (Park et al. 2002). In the hydrodynamics-based analysis of the 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 suggested that the low-temperature plasma dominated the soft band emission until ∼ 10500 days, and resulted in enhanced average abundances of N, O, Ne, and Mg during this period. As 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 the lower metal abundances of 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 inter-clump gas and the evaporated/destructed 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 f LT /f HT . 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 the differential emission measure analysis based on the high resolution spectra (e.g., Zhekov et al. 2006Zhekov et al. , 2009, and detailed full-scale HD/MHD simulations (e.g., Orlando et al. 2015Orlando et al. , 2019Orlando et al. , 2020.

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 a Fe K centroid energy 6.65 keV. By comparing with the theoretically predicted line centroid shown in Figure 10 (based on SPEX code 7 , Kaastra et al. 1996), the high centroid energy requires either a high plasma temperature 2.0 keV or a large ionization parameter 10 11 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 (χ 2 r ∼ 1.19, versus χ 2 r ∼ 1.33 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 kT e = 2.9 +0.2 −0.1 keV and τ = (1.4 ± 0.1) × 10 11 cm −3 s. The other two components have lower temperatures and similar ionization parameters: kT e = 0.35 +0.04 −0.01 keV, τ = (3.2 ± 0.1) × 10 11 cm −3 s for the low-temperature one, and kT e = 0.87 +0.03 −0.02 keV, τ = 1.9 +0.3 −0.2 × 10 11 cm −3 s for the mid-temperature one. On the other hand, the temporal variation of the Fe K centroid energy provided 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 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 with the predicted evolution of a plasma with kT e ∼ 3.2 keV and n e ∼ 500 cm −3 , shocked at 7000 days after the explosion. This Fe-K emitting component may thereby 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. 2009Zhekov et al. , 2010. It is also possible that the Fe K emission comes from one or several reverse shocked Fe-rich ejecta clumps, as the 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 which 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).

CONCLUSION
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, as well as the EPIC-pn camera provided us with an overall energy coverage from 0.3 to 10 keV and thus 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 (kT e,LT ∼ 0.6 keV and kT e,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 since around 10500 days after the explosion, resulting in the 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 result from the different chemical compositions between two plasma components.
4. The densities of the X-ray emitting gases were estimated as n e,LT ∼ 2400 cm −3 and n e,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 inter-clump gas and the evaporated/destructed 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 penetrating 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 monitor these new changes and unveil the underlying physics.

NOTE ADDED IN PROOF
During the refereeing process of the present paper, two new analyses regarding the X-ray emission from SN 1987A were presented (Greco et al. 2021;Alp et al. 2021), incorporating NuSTAR hard Xray 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 non-thermal 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 got in Section 4.4. In fact, the temperature and ionization parameter we got 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 non-thermal 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 the non-thermal emission, is still required if we would 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. Software: XSPEC (Arnaud 1996), SPEX (Kaastra et al. 1996), SAS (Gabriel et al. 2004 Dec. 2010Dec. 2011Dec. 2012Dec. 2014Nov. 2015Nov. 2016Nov. 2017 (2015), with bright emission lines labeled. Observed spectra are denoted by data points, with different colors standing for different instruments and grating orders. The best-fit model (nlapec + gauss) is denoted by the black line, and the residuals are showed in the bottom panel.     (2015), with different spectral models and residuals. Observed spectra are denoted by data points, with different colors standing for different instruments and grating orders. The best-fit models are denoted by black lines, and the high-and low-temperature components by orange and blue dashed lines, respectively. Top: the two-temperature vnei model fitting. Bottom: the two-temperature vpshock model fitting.  Frank et al. (2016), while the blank ones are newly published data obtained in this work. Error bars represent the 90% uncertainties, note that for most of the data points the error bars are too small to be visible.  2019 Nov.