Spectral Evolution of the X-Ray Remnant of SN 1987A: A High-resolution Chandra HETG Study

Based on the previous Chandra data, it has been amply shown that spectral model fits with two characteristic components are required to adequately describe the observed Chandra gratings spectrum of SNR 1987A (e.g., Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009).

Adopting such an approach, we performed simultaneous spectral fitting for all the deep HETG (MEG, HEG) and LETG (LEG) first-order dispersed spectra taken at five different epochs including our new data (Table 1).

Significant non-equilibrium effects have been recognized in the X-ray-emitting plasma of SNR 1987A based on previous Chandra grating spectral analyses (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009). We model the X-ray emission spectrum with a two-component plane-parallel shock model (vpshock in XSPEC) that takes into account the non-equilibrium ionization (NEI) effects in a hot optically thin plasma (Borkowski et al. 2001). We adopt the NEI version 3.0 in XSPEC, based on ATOMDB v3.0.9. ⁷ In our latest HETG 2018 broadband spectra (Figure 1), we do not find any significant emission features corresponding to trace elements. This is consistent with previous analyses of the X-ray emission from SNR 1987A between 2004 and 2007 (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009), where X-ray emission from these trace elements were negligible. This is also noted with XMM-Newton data from 2012–2017 (Alp et al. 2021). Thus, we set the abundances of all trace elements with vpshock (version ATOMDB v3.0.9) to zero, using the XSPEC command NEI_TRACE_ABUND. In our spectral model fits, we varied the electron temperature (kT) and the ionization age (τ = n_e t, where n_e is the post-shock electron density and t is the time since the gas first entered the shock). We also keep the normalization parameter free, which corresponds to the volume-emission measures (EM). We also convolved the physical shock model with a Gaussian profile (gsmooth) to represent the broadening of the spectral lines due to the bulk and turbulent motion of the X-ray-emitting plasma (see Dewey et al. 2012 for a detailed discussion). Within gsmooth, both the magnitude of the smoothing measured by the Gaussian kernel and the α-index (exponent for gsmooth energy dependence) are kept free.

We fit the total foreground absorption column by applying two-component multiplicative absorption models, tbabs and tbnew. ⁸ The absorption component tbabs accounts for the Galactic absorption column N_H,Gal, and it is fixed at N_H,Gal = 6 × 10²⁰ cm⁻² (Dickey & Lockman 1990). The second absorption component, tbnew, accounts for the LMC absorption column, N_H,LMC and is varied in the fits. We use tbnew as the absorption column for LMC as it allows us to set individual elemental abundances associated with N_H,LMC at their respective LMC values. Recent measurements of the LMC abundances based on the X-ray spectral analysis of SNRs (Maggi et al. 2016; Schenck et al. 2016) suggest ∼50% lower LMC abundance values for O, Ne, Mg, and Fe than the previously estimated values (Hughes et al. 1998). Our derived shock parameters are consistent (within statistical uncertainties), assuming either set of these LMC abundances in tbnew. We tie N_H,LMC among all the fitted spectra, assuming no temporal variation in the LMC column for SNR 1987A. We obtain a best-fit total absorption column density, ${N}_{{\rm{H}}}={N}_{{\rm{H}},\mathrm{Gal}}\,+{N}_{{\rm{H}},\mathrm{LMC}}={2.17}_{-0.22}^{+0.22}\times {10}^{21}\,{\mathrm{cm}}^{-2}$. This value is comparable to estimates of N_H obtained by other X-ray analyses: ${2.35}_{-0.08}^{+0.09}\,\times {10}^{21}\,{\mathrm{cm}}^{-2}$ (Park et al. 2006), ${1.44}_{-0.12}^{+0.16}\times {10}^{21}\,{\mathrm{cm}}^{-2}$(Zhekov et al. 2009), and ${2.60}_{-0.05}^{+0.05}\times {10}^{21}\,{\mathrm{cm}}^{-2}$ (Alp et al. 2021). We note that more recent H i surveys have suggested much higher N_H,Gal values toward the LMC, i.e., ∼(2.5–4) × 10²¹ cm⁻² (Kalberla et al. 2005; Willingale et al. 2013; HI4PI Collaboration et al. 2016). By adopting these high Galactic columns, the overall fits are equally good, but this results in a negligible LMC column. The implied negligible LMC column does not appear to be reasonable because optical extinction estimates show that the LMC contribution is greater than the Milky-Way contribution (Fitzpatrick & Walborn 1990; France et al. 2011). A detailed analysis of the contribution of the Galactic N_H toward total absorption column density is beyond the scope of our work. While this issue was similarly outlined by Alp et al. (2021) in relation to their analysis of the XMM-Newton data of SNR 1987A, we note that the total columns (Galactic + LMC) are generally consistent between either values of the Galactic column, and thus that the best-fit parameters in our spectral model fits (i.e., electron temperatures, ionization age, abundances, and volume-emission measures) are not significantly affected (within statistical uncertainties).

We note that the low-energy emission (E ≤ 2 keV) is likely to be dominated by radiative shocks (Pun et al. 2002; Gröningsson et al. 2008). For a given temperature, these radiative shocks will be softer than adiabatic shocks (Nymark et al. 2006). This, along with uncertainties in the measured absorption column and the decreasing low-energy sensitivity of the ACIS detectors, adds to the systematic uncertainties in our analysis. Additional systematic uncertainties arise from the two-temperature fit to our broadband spectra, which is a crude approximation of a complex multi-temperature shocked medium. In Section 4.4, we compare our inferred results from the broadband spectral fits with a continuous distribution of plasma characteristics from self-consistent MHD simulations. We note that our data are somewhat insensitive to high-energy/temperature components (compared to XMM-Newton or NuSTAR observations) due to the relatively low sensitivity of Chandra grating spectra above ∼6 keV, which also contributes to systematic uncertainties, particularly regarding any high-energy components.

Previous works have shown no evidence for significant temporal changes of the elemental abundances in the observed X-ray spectrum of SNR 1987A (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009; Bray et al. 2020). Based on our spectral model fits to each individual spectrum taken at every epoch between 2004 and 2018, we confirm that the estimated abundances are consistent with each other within statistical uncertainties. Thus, we assumed constant abundances for Ne, Mg, Si, S, and Fe, for which various emission lines from the K- and L-shell electron transitions are evident in the observed spectra of HETG 2018 (Figure 1) and tied the abundance of each element together for joint fits to all five data sets (LETG 2004, HETG 2007, LETG 2007, HETG 2011, and HETG 2018) in order to obtain the best-constrained abundance values for all epochs between 2004 and 2018. With XMM-Newton data, Sun et al. (2021) reported decreasing average abundances of N, O, Ne, and Mg over a similar period as our analysis. They have suggested that this may reflect the different chemical compositions between the two plasma components. However, based on our deep Chandra HETG data, we find no evidence for a statistically significant evolution of elemental abundances in either of our soft or hard component spectral model fits.

The solar abundance table angr (Anders & Grevesse 1989) was adopted in all our previous Chandra spectral analyses of SNR 1987A (Michael et al. 2002; Park et al. 2002, 2004, 2006; Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009). In this work we use an updated solar abundance table, aspl (Asplund et al. 2009). We found that this change in the assumed solar abundances does not significantly affect our results. The best-fit values of the spectral model parameters are consistent within statistical uncertainties between model fits, adopting these two different solar abundance tables. Thus, hereafter, we quote all elemental abundance values with respect to the aspl solar abundances. We fix other elemental abundances for SNR 1987A at: He = 1.98 (Mattila et al. 2010), C = 0.12 (Fransson et al. 1996), N = 0.92, O = 0.14 (Zhekov et al. 2009). We also fix Ar = 0.776, Ca = 0.354, and Ni = 0.662 at their respective LMC values (Russell & Dopita 1992). We tie the elemental abundances between the soft and hard component shock models for all epochs.

In Table 2, we show the results from our simultaneous two-component shock model fits to the MEG and HEG (±1) spectra from HETG 2007, HETG 2011, HETG 2018, and LEG (±1) spectra from LETG 2004, LETG 2007. In Figure 1, we present our broadband best-fit model overlaid over HETG 2018 (MEG ± 1). We also present our best-fit spectral model in several sub-bands around strong emission lines from Si, Mg, and Ne for these data sets (Figure 2). Between 2004 and 2011, the observed electron temperature (kT) for the soft component showed no significant changes (within uncertainties), while for the hard component it gradually decreased from ∼2.5 to ∼1.8 keV. In contrast, between 2011 and 2018, both the soft- and hard-component kT have increased from ∼0.62 to ∼0.86 keV (∼38%) and from ∼1.96 to ∼2.41 keV (∼23%), respectively. Electron temperatures and ionization ages measured in Table 2 for epochs of LETG 2004, HETG 2007, LETG 2007 are consistent (within statistical uncertainties) with previous measurements (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009). From such a broadband spectral model fit, our best-fit elemental abundances are: Ne = 0.34 ± 0.01, Mg = 0.25 ± 0.01, Si = 0.36 ± 0.01, S = 0.40 ± 0.04, and Fe = 0.19 ± 0.01. These best-fit abundances are consistent (within statistical uncertainties) with those for the ER as measured in the previous Chandra grating spectral analyses of SNR 1987A (Zhekov et al. 2006; Dewey et al. 2008, 2012; Zhekov et al. 2009).

Zoom In Zoom Out Reset image size

Table 2. Simultaneous Two-shock Fit Results

Parameters	LETG 2004	HETG 2007	LETG 2007	HETG 2011	HETG 2018
χ2/dof	⋯	⋯	⋯	⋯	3454/5181
NH(1021 cm−2) a	2.17	2.17	2.17	2.17	${2.17}_{-0.22}^{+0.22}$
kTsoft (keV)	${0.54}_{-0.05}^{+0.04}$	${0.61}_{-0.03}^{+0.02}$	${0.57}_{-0.01}^{+0.02}$	${0.62}_{-0.03}^{+0.03}$	${0.86}_{-0.08}^{+0.07}$
kThard (keV)	${2.49}_{-0.28}^{+0.46}$	${2.25}_{-0.21}^{+0.30}$	${1.88}_{-0.13}^{+0.14}$	${1.96}_{-0.15}^{+0.18}$	${2.41}_{-0.16}^{+0.20}$
τsoft b	${3.29}_{-0.77}^{+1.25}$	${3.36}_{-0.41}^{+0.58}$	${3.80}_{-0.47}^{+0.59}$	${7.51}_{-1.53}^{+2.17}$	${6.20}_{-1.42}^{+1.84}$
τhard b	${1.10}_{-0.25}^{+0.50}$	${2.43}_{-0.51}^{+0.63}$	${2.27}_{-0.36}^{+0.46}$	${3.20}_{-0.34}^{+0.72}$	${5.99}_{-0.50}^{+1.15}$
EMsoft c	${4.46}_{-0.46}^{+0.34}$	${12.58}_{-0.43}^{+0.49}$	${14.41}_{-0.93}^{+0.77}$	${23.59}_{-1.92}^{+0.96}$	${18.26}_{-1.61}^{+1.58}$
EMhard c	${2.26}_{-0.43}^{+0.12}$	${5.02}_{-0.40}^{+0.58}$	${7.13}_{-0.45}^{+0.58}$	${13.51}_{-1.20}^{+1.05}$	${17.85}_{-1.27}^{+0.96}$
He (fixed)	1.98	1.98	1.98	1.98	1.98
C (fixed)	0.12	0.12	0.12	0.12	0.12
N (fixed)	0.92	0.92	0.92	0.92	0.92
O (fixed)	0.14	0.14	0.14	0.14	0.14
Ne	0.34	0.34	0.34	0.34	${0.34}_{-0.01}^{+0.01}$
Mg	0.25	0.25	0.25	0.25	${0.25}_{-0.01}^{+0.01}$
Si	0.36	0.36	0.36	0.36	${0.36}_{-0.01}^{+0.01}$
S	0.40	0.40	0.40	0.40	${0.40}_{-0.04}^{+0.04}$
Ar (fixed)	0.776	0.776	0.776	0.776	0.776
Ca (fixed)	0.354	0.354	0.354	0.354	0.354
Fe	0.194	0.194	0.194	0.194	${0.19}_{-0.01}^{+0.01}$
Ni (fixed)	0.662	0.662	0.662	0.662	0.662

Notes. 90% confidence intervals are provided as error bars and all abundances are expressed as ratios to the aspl solar abundance table. The two shock components are represented by soft and hard labels.

^a N_H = N_H,Gal + N_H,LMC. ^b Ionization age (τ = n_e t) in units of 10¹¹ cm⁻³ s. ^c EM = ∫n_e n_H dV in units of 10⁵⁸ cm⁻³ for an adopted distance of 51 kpc.

Download table as:  ASCIITypeset image

High-resolution HETG/LETG spectroscopy provides an excellent opportunity for line diagnostics to investigate the thermal condition of the X-ray-emitting hot gas, complementary to the broadband spectral model fits. In the broadband spectrum (Figure 1), we identify several emission lines from He-like and/or H-like ions of S, Si, Mg, and Ne. We fit the Heα triplet lines (from He-like ions) with three Gaussian profiles and an underlying constant continuum, in the "narrow" bands around the emission lines. Similarly, for the Lyα and Lyβ doublet lines (from H-like ions), we use two Gaussian profiles and an underlying constant continuum. For strong Heα lines of Si and Mg we clearly resolve resonance (r), and forbidden (f ) lines, while their intercombination (i) lines are faint and are not clearly distinguishable (Figure 3). Thus, we allow the line ratios (f /r and i/r) to vary in the fits. For these emission lines from H- and He-like ions, we fixed the line centers at their lab-measured values from ATOMDB v3.0.9. We include a redshift parameter for each emission line to account for the line-of-sight velocity effects on the line-center wavelengths. This redshift parameter is free for each individual Lyα, Lyβ, and Heα multiplet, to account for any possible differential line shifts between them. Thus, free parameters in our Heα line model are the FWHM, normalizations for Gaussians, f /r and i/r ratios, continuum normalizations, and redshifts. For our models of the Lyα and Lyβ lines, we assumed an equal flux from each of the sub-component lines based on ATOMDB v3.0.9 because our HETG data cannot resolve them individually. Free parameters are the FWHM, Gaussian and continuum normalizations, and redshifts.

Zoom In Zoom Out Reset image size

Based on the best-fit redshift parameter, we find that the centroid shifts for all the strong emission lines are consistent (within statistical uncertainties) with the optical estimates of SNR 1987A redshift of ∼286 km s⁻¹ (Gröningsson et al. 2008). We also find that the fluxes between positive and negative arms are consistent within statistical uncertainties. Thus, we fit all four of the HETG ± 1 (MEG ± 1 and HEG ± 1) spectra simultaneously, sharing normalizations, line centers, line-flux ratios, and redshifts. For the LETG spectra, we adopt a similar method for the LEG ± 1 spectral fits. We untied their FWHMs to account for differential widths that are caused by the tilted orientation of the ER with respect to the plane of the sky (inclination angle i = 44°–45°; Plait et al. 1995) with the northern rim toward Earth. In Figure 3 we show example best-fit line profiles from our Gaussian model fitting for Heα and Lyα lines from Si, Mg, S, and Ne ions. In Table 3, based on these line-profile fits, we present our measured line-emission fluxes. Our estimated fluxes for LETG 2004, HETG 2007, and LETG 2007 observations are consistent with the published values (Zhekov et al. 2005; Dewey et al. 2008; Zhekov et al. 2009) within statistical uncertainties. In Figure 4, we show these temporal variations of Heα and Lyα line fluxes for Si, Mg, and Ne in 2004–2018. Between 2011 and 2018, Si Lyα line flux has increased ∼51% while Si Heα line flux has stayed roughly constant. On the other hand, Mg Heα line flux has decreased by ∼37% while Mg Lyα line flux has stayed roughly constant. These changing trends of the Heα and Lyα line fluxes of the strongest line profiles (Si and Mg) can also be quantified with the Heα/Lyα line-flux ratios. For both Si and Mg, Heα/Lyα line-flux ratios have significantly decreased (by ∼40%) between 2011 and 2018 (Figure 5). Thus, between 2004 and 2018, line fluxes and line-flux ratios (Si, Mg: Heα/Lyα) have undergone statistically significant changes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Combining information from all five deep Chandra grating spectroscopic data of SNR 1987A, including our new HETG 2018 data, we explore the X-ray spectral evolution of SNR 1987A and changing thermal conditions of the X-ray-emitting plasma over the 14 yr period from 2004 to 2018. Volume-Emission Measure (EM $\sim {n}_{e}^{2}V$, where n_e is the electron density and V is the X-ray-emitting volume) is considered to be a tracer of density of the shocked gas. Based on our broadband spectral model fits (Table 2), we show temporal changes in EM for the best-fit soft- and hard-component spectral models (Figure 6). While the EM increase for both the soft and hard components between 2004 and 2007 is consistent (within statistical uncertainties) with results from previous works (Zhekov et al. 2006, 2009), we note significant deviations in this trend since 2011. The soft component EM has declined (by ∼23%) between 2011 and 2018 (Figure 6(a)), while the hard component EM continues to increase (by ∼32%), but at a slower rate than before 2011 (Figure 6(b)).

Zoom In Zoom Out Reset image size

While the changes in EM after 2011 are intriguing, we note that it is based on a minimal sample of our measurements, covering only two epochs (HETG 2011, HETG 2018) with a relatively long time separation of 7 yr between them. To complement this small number of measurements, we add the published EM measurements by Bray et al. (2020) in our sample. Based on similar two-component spectral model fits to our annual Chandra HETG monitoring observation data of SNR 1987A, Bray et al. (2020) measured EMs for the soft- and hard-component emission models. Their data were taken with shorter exposure times (∼50–70 ks), but covered a similar time period (from 2011 to 2019) to that of our interest with a significantly higher cadence of observations at 11 epochs. In Figure 6, we overlay EM values measured by Bray et al. (2020). We find that the soft component EMs measured by Bray et al. (2020) between 2011 and 2019 are well aligned with our measured decline, filling the gap between 2011 and 2018. While the statistical uncertainties are large due to poor count statistics in Bray et al.'s (2020) measurements, they provide us confidence for the declining soft-component EM since 2011. The hard-component EM measurements by Bray et al. (2020) also show a trend similar to our measurements (Figure 6(b)). Their hard-component EM values are systematically higher (by ∼30%) than our corresponding measurements of 2011 and 2018. We find that this discrepancy is the model-dependent systematics between Bray et al. (2020) and this work (e.g., differences among plasma models adopted for the soft component, fitting the foreground absorption column, adopted versions of CALDBs, etc.). Applying exactly the same spectral model fits (as those we use in this work) for the Bray et al. (2020) data, we confirmed that such a discrepancy in the hard-component EMs (between this work and Bray et al. 2020) are removed.

Based on these EM variations, we estimate the density profiles of the X-ray-emitting gas. We construct a toy model in the form of a simple power law (PL), EM ∼ t^β (where β is the PL index). We fit the temporal changes of the soft- and hard-component EMs with this simple PL model (Figure 6). The best-fit PL indices for the soft-component EM variation are β = 4.8 ± 0.8 until 2011 and β = −2.4 ± 0.8 after 2011, indicating the clear turnover since 2011. For the hard-component EM variation, β = 5.5 ± 0.7 until 2011 and β = 1.1 ± 0.1 afterwards. While the change in the PL slope for the hard component is less pronounced than for the soft component, the decrease in the PL index is evident for the hard-component EM after 2011. Based on our current physical picture with previous Chandra spectral analyses (e.g., Park et al. 2002, 2004, 2006; Zhekov et al. 2006, 2009; Dewey et al. 2008; Helder et al. 2013; Frank et al. 2016; Bray et al. 2020), the soft component is dominated by emission from the shocks interacting with high-density clumps in the ER, and the hard component may represent a combination of shock interactions with lower-density inter-clump regions in the ER and less dense CSM (which may include regions outside of the ER). Thus, the X-ray-emitting volumes represented by the soft and hard components could be different, and we estimate their associated density profiles independently.

As the X-ray emission characterized by the soft component is dominated by emission from the dense ER, we assume the soft component X-ray-emitting volume approximately to be the volume of the ER. Early HST images of SN 1987A show that the ER has a thin disk-like ring geometry (Plait et al. 1995). The volume for such a disk-like ring can be expressed as $V\sim h\times ({r}_{\mathrm{out}}^{2}-{r}_{\mathrm{in}}^{2})$, where h represents the vertical thickness (height) of the disk, and r_in and r_out are the inner and outer radii of the disk, respectively. The height h of the ER was estimated to be confined within ±30° of the equatorial plane based on modeling of the Lyα line in the HST spectra (Michael et al. 1998, 2003). As no significant evolution of this thickness is observed, we can assume that h has been relatively constant. The soft X-ray spectrum of SNR 1987A is dominated by emission from a disk-like volume of the ER (e.g., Zhekov et al. 2005, 2006). Thus, we consider that the X-ray emission originates from a disk-like ring for which the emission volume would be proportional to ${r}_{\mathrm{out}}^{2}-{r}_{\mathrm{in}}^{2}$ ∼ r². Assuming n_e ∼ r^ν (where r is the distance from the SN site, and ν is the density-profile index) and V_soft ∼ r², EM for the soft component can be expressed as EM_soft ∼ r^2ν+2. From morphological studies of ACIS images of SNR 1987A, the expansion rate of the X-ray ring is linear in time (r ∼ t) since ∼2004 (Racusin et al. 2009; Frank et al. 2016). Since all our epochs in this work (2004–2018) fall within this linear regime, the volume-emission measure is EM ∼ t^β ∼ t^2ν+2, and thus, β = 2ν + 2. Our best-fit values of β for the soft component thus indicate the associated density profiles of n_e ∼ r^1.4±0.4 in 2004–2011 and n_e ∼ r^−2.2±0.4 in 2011–2018. While our PL analysis with a disk-like emission volume might be an oversimplification of a physically more complicated geometry, a significant change in the density profile of the shocked gas responsible for the soft-component X-ray emission since 2011 is evident.

In Figure 7 we present the time evolution of the the soft-to-hard component EM ratio. It shows a significant decrease (by a factor of ∼2) in 2018, compared to those in 2004–2011. Such a drop indicates that the soft X-ray emission from the shocked dense clumpy gas in the ER is less significant in 2018 compared to 2004–2011. We note that our derived density profiles associated with the soft component post 2011 (n_e ∼ r^−2.2±0.4) are in plausible agreement with the standard density profile of red supergiant (RSG) winds (ρ ∼ r⁻²). Thus, the X-ray emission associated with the soft shock component appears to have increasing contributions from the low-density CSM (beyond the dense ER) since 2011, probably created by the massive progenitor's stellar winds in its RSG stage before the SN explosion.

Zoom In Zoom Out Reset image size

The observed X-ray spectrum fitted by the hard-component shock model has been likely dominated by the emission from the low-density inter-clump regions in the ER. In our recent Chandra observations, this hard-component X-ray emission may also include some contribution from the shocked gas beyond the ER (Frank et al. 2016; Alp et al. 2021; Sun et al. 2021). As a zeroth-order approximation, we consider a spherical shell-like volume with a radius of the blast wave (r_b ) for the hard-component X-ray-emitting gas. In such a case, the emission volume for the hard component is approximated to be ${V}_{\mathrm{hard}}\sim {r}_{b}^{3}$. While changes in the temporal variation of the hard-component EM (after 2011) are less prominent than those in the soft component EM, a single PL model fit significantly underestimates the measured hard component EM (by ∼60%) in 2018. We find that a two-component PL fit is required to adequately fit the entire EM variation including the 2018 data (e.g., ${\chi }_{\nu }^{2}\sim 17$ for the single PL model fit, and ${\chi }_{\nu }^{2}\sim 3$ for the two-component PL fit, Figure 6(b)). Thus, we consider the two-component PL model fit to adequately describe the temporal variation of the hard component EM between 2004 and 2018, with the break around 2011, similar to its softer counterpart. As discussed in Section 4.1, we note that there are systematic differences between hard-component EM values between deep Chandra grating (our work) and monitoring grating spectra (Bray et al. 2020). Based on the hard-component EM measurements by Bray et al. (2020), we obtain the best-fit PL index β = 1.4 ± 0.2, which is in general agreement (within statistical uncertainties) with that fitted to our deep Chandra HETG spectra taken in 2011 and 2018. Thus, best-fit values of β for the hard-component EM evolution based on measurements with the deep HETG/LETG observations (β = 5.5 ± 0.7 in 2004–2011 and β = 1.1 ± 0.1 in 2011–2018) and the annual monitoring grating observations (β = 1.4 ± 0.2 in 2011–2019) both indicate that the shocks responsible for the hard component have been interacting with less dense CSM since 2011.

For our assumed simple spherical geometry (${V}_{\mathrm{hard}}\sim {r}_{b}^{3}$) for the X-ray-emitting gas outside the ER, we consider self-similar solutions for the SNR dynamics (Chevalier 1982). The radius of the blast wave in such a solution can be expressed as r_b ∼ t^{(n−3)/(n−s)}, where n is the PL index for the SN ejecta density profile (i.e., ${\rho }_{\mathrm{ejecta}}\sim {r}_{b}^{-n}$) and s is the PL index for the density profile of the CSM (i.e., ${\rho }_{\mathrm{CSM}}\sim {r}_{b}^{-s}$). For SN 1987A, we assume n = 9 (Eastman & Kirshner 1989; Suzuki et al. 1993; Borkowski et al. 1997a, 1997b). For the two representative cases s = 0 (constant density) and s = 2 (RSG wind density), the radius of the blast wave would be r_b ∼ t^2/3 (constant density) and r_b ∼ t^6/7 (RSG wind density). With our adopted spherical X-ray-emitting geometry, the hard component EM is EM_hard $\sim {\rho }_{\mathrm{CSM}}^{2}{V}_{\mathrm{hard}}$ $\sim {\rho }_{\mathrm{CSM}}^{2}{r}_{b}^{3}$. Thus, EM_hard ∼ t² and EM_hard ∼ t^−6/7 are inferred for the two representative cases of constant and RSG wind densities for the interacting CSM, respectively. Our derived PL index post 2011 (β = 1.1 ± 0.1) for the hard component is in between what we expect for the constant density and the RSG wind-density cases, suggesting a transitionary phase for the blast wave interacting with the ER densities and lower ambient CSM densities outside it. A new phase of SNR 1987A's evolution, in which the shock starts interacting with the low-density CSM beyond the dense ER, has been proposed to have started at around ∼2016 (∼10,400 days) based on our previous Chandra monitoring data (Frank et al. (2016). HST observations show that the optical emission from the ER peaked around 2009 (∼8000 days) and has been fading since then, signaling that the blast wave has started to leave the ER (Fransson et al. 2015). While the physical processes leading to X-ray and optical emission might be different, our estimated "switching" period at around 2011 (∼8500 days) for both soft and hard components (Figure 6) is in plausible agreement with those previously suggested phase transitions.

Based on the broadband spectral model fits (Table 2), we show the evolution of electron temperatures (kT) in 2004–2018 (Figure 8). The soft-component electron temperature has clearly increased between 2011 and 2018 (by ∼38%). While it is statistically marginal due to large uncertainties, the hard-component electron temperature shows a similar increase (∼23%) during this period. In contrast, such a large increase in kT was not evident in 2004–2011 for either the soft or hard components. Similar temporal changes of kT between 2011 and 2018 have been reported in the literature based on the annual monitoring data of SNR 1987A with Chandra HETG (Bray et al. 2020). Our results in 2004–2007 for both the soft and hard components are consistent (within statistical uncertainties) with those previously reported by Zhekov et al. (2009) for the same epochs. Such increases in kT for both the soft and hard components in the new HETG 2018 data are in plausible agreement with the overall decreasing density of the shocked gas since 2011.

Zoom In Zoom Out Reset image size

We present the evolution of the ionization age over the same period in Figure 9. The soft-component ionization age increases between 2004 and 2011 as the dense gas is shocked within the ER. While statistical uncertainties are large, the soft component ionization age appears to have leveled off or is even marginally decreasing since 2011 (Figure 9(a)). These results generally support the decreasing density for the shocked gas responsible for the soft component. Such a leveling off or a marginal decrease of the soft component ionization age since 2011 may also indicate that some of the densest material in the ER (representing the highest τ) becomes cold (possibly through radiative cooling), that its emission falls out of the X-ray band, making only the material close to the shock front dominate the observed X-ray spectrum. The hard-component ionization age increases linearly (by ∼87%) from 2011 to 2018 (Figure 9(b)). Thus for both kT and τ, the soft and hard components have become less distinguishable since 2011. Such a "merging" trend could be a combined effect of increasing contribution from the low-density shocked gas around the ER (beyond and/or above and below the ER) and of the gradual thermalization of the shocked plasma of the ER due to Coulomb collisions. Future high-resolution X-ray spectroscopic observations would be required to adequately trace and constrain the changing physical conditions in the evolution of SNR 1987A.

Zoom In Zoom Out Reset image size

We compare the results from our Chandra HETG spectral analysis of SNR 1987A with the 3D MHD simulations by Orlando et al. (2019). It is worth noting that these MHD simulations are not the result of fitting to the X-ray data; rather, they have been constrained by X-ray data collected before 2016. In particular, the model parameters have been explored within ranges of values inferred from optical and X-ray observations, performing an iterative process of trial and error to converge on model parameters that best reproduce the X-ray light curves (soft and hard bands) and moderate-resolution spectra of SN 1987A (see Orlando et al. 2015, 2019, 2020). Then the X-ray emission has been synthesized from the best model by adopting the approach described in a number of papers (e.g., Orlando et al. 2015; Miceli et al. 2019) and based on the NEI emission model vpshock, available in the XSPEC package. Thus the model is not fully independent of the X-ray data (having been constrained with them) but it provides a single 3D physical (not phenomenological) model that reproduces self-consistently the X-ray light curves, spectra (even high-resolution dispersed spectra not used to constrain the model), and morphology of the remnant at all epochs (even after 2016, the date of the last X-ray observations used to constrain the model) since the SN event.

In the left panels of Figure 10, we show the simulated EM distribution relative to the electron temperature and ionization age of the shocked gas in SNR 1987A from 2004 (t = 17 yr) to 2018 (t = 31 yr). The contributions from each of the three assumed components (H ii region, ER, and SN ejecta) are shown in the right panels of Figure 10. The model predicts a peak of EM for the shocked clumps of the ER at higher ionization age and lower temperature than the isothermal components fitting the observed spectra. This discrepancy might have been caused in part by our simple approach in the broadband spectral model fits, where we fit the observed spectra with only two characteristic components rather than the physically more realistic "continuous" distributions of electron temperatures. A continuous distribution of electron temperatures has been considered to describe the X-ray spectrum of SNR 1987A in previous works (Chandra: Zhekov et al. 2006, 2009; XMM-Newton: Alp et al. 2021). These models, even if they allow continuous temperature distributions, show peaks around the characteristic temperatures of the discrete models. The peak of the modeled EM gradually shifts to higher kT and τ as the shock evolves. This is in general agreement with the observed increases in kT and τ for both of the soft- and hard-component X-ray emission spectra, indicating the progressive equilibration of the shocked gas. Increasing kT from the MHD models and observed decreasing EM ratios between the soft and hard components since 2011 (Figure 9) suggests that the CSM responsible for the soft component is less dense in 2018 than before. This also supports the scenario of the CSM associated with the soft-component emission being less clumpy than before.

Zoom In Zoom Out Reset image size

In 2004 (17 yr after the SN) the hard-component kT and τ as observed with Chandra data (the upper white cross in the upper right panel of Figure 10) are in between two local peaks of the MHD-modeled EM distribution, indicating contributions from both the shocked H ii region (blue in the upper right panel of Figure 10) and the low-density component of the shocked ER region (red in the upper right panel of Figure 10). There is also some contribution expected from the shocked-ejecta (in green: upper right panel) component, though at this stage of SNR 1987A evolution (t = 17 yr), this component may include only the ejected mantle, which is not metal-rich. As the blast wave evolves, the second local peak (low-density component of the shocked ER) gradually becomes dominant over the first peak (shocked H ii region), due to the increase in the amount of the shocked ER gas. In fact, from 2004 (t = 17 yr) to 2011 (t = 24 yr), the hard-component kT in the MHD models slightly decreases because it is most affected by the second peak (the kT component migrates to fit the second peak). At later times, once the ER is fully shocked, no new freshly shocked ER material contributes to the EM. Hence the peak of the EM shifts to higher values of kT due to the plasma thermalization through Coulomb collisions. Thus, EM maps from MHD simulations in 2018 (t = 31 yr: lower left and right panels) suggest that the blast wave has started propagating beyond the ER. These maps are consistent with previous findings of Frank et al. (2016) in the X-ray band, of Cendes et al. (2018) in the radio band, and of Fransson et al. (2015) and Larsson et al. (2019) in the optical band, showing that the blast wave has propagated beyond the ER. These effects in the X-ray band have also been corroborated by recent results from Sun et al. (2021) and Alp et al. (2021). The ionization age gradually increases during this time period due to the progressive equilibration of the shocked plasma.

These EM distribution maps based on our MHD calculations suggest that our characteristic soft-component kT represents emission from the shocked high-density clumps of the ER, whereas emission from the low-density component of the shocked ER gas is mainly responsible for the hard component (although a contribution from the shocked H ii region is also evident in 2004; t = 17 yr). Shifting of the EM peak toward higher kT is consistent with the scenario of the blast wave propagating into a low-density region beyond the dense ER. Thus, these results from the 3D MHD calculations support our conclusions on the changing density profiles of the X-ray-emitting gas as we discussed in Section 4.1. Unveiling matter beyond the ER is important to characterize the structure of the CSM sculpted by the winds of the progenitor star and thus to trace the final phases of its evolution.

The ionic Heα/Lyα line-flux ratios are functions of both the electron temperatures (kT) and ionization ages (τ). It has been shown that a single shock component cannot explain the observed line ratios from SNR 1987A spectra (Zhekov et al. 2005) as these are products of a distribution of shocks with a range of ionization ages and post-shock electron temperatures (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009). Utilizing a simple NEI plane-parallel shock model with a foreground absorbing column (N_H = 2.17 × 10²¹ cm⁻², Table 2), we show the ranges of electron temperature and ionization age which are consistent with observed Heα/Lyα line-flux ratios for Si and Mg (Figure 11). For these comparisons with observed line-flux ratios, we have considered HETG 2007, HETG 2011, and HETG 2018.

Zoom In Zoom Out Reset image size

Based on broadband spectral model fits (Table 2), the observed range of τ is ∼1 × 10¹¹ – 8 × 10¹¹ s cm⁻³. In this range, we note that the kTs corresponding to the observed line-flux ratios for any given τ have steadily increased from 2007 to 2018 for both Si and Mg. In Figure 11, except for a very narrow range of τ (around ∼1 × 10¹¹ s cm⁻³), the allowed kTs for X-ray-emitting plasma are likely to be in the range ∼0.5–2.0 keV (HETG 2007), ∼0.6–2.0 keV (HETG 2011), and ∼0.8–2.5 keV (HETG 2011). While there are significant overlaps in the allowed kTs for HETG 2007, HETG 2011, and HETG 2018, these ranges have consistently shifted to higher values from 2007 to 2018. This is consistent with the broadband spectral model fits showing a continuous increase of kT from 2007 to 2018 (Figure 8). For a direct comparison, we overlay the best-fit soft- and hard-component kT and τ from Table 2 over Si and Mg isolines in HETG 2007, HETG 2011, and HETG 2018 (Figure 11(a), (b)). We note that for all three epochs, the locus of temperatures and ionization ages suggested by the single NEI component is in between the values of kT and τ measured with the more detailed two-component shock model fits. Allowed electron temperatures derived from line-flux ratios for HETG 2007, HETG 2011, and HETG 2018 are consistent with what we find from broadband spectral model fits (Table 2). Thus, our observed Heα/Lyα line-flux ratios are in plausible agreement with the overall increase in the post-shock electron temperatures.

Adopting a similar approach to that used for previous Chandra grating data of SNR 1987A (Zhekov et al. 2005), we simultaneously fit the widths of individual lines detected in both grating dispersion arms (positive and negative). We assume the same physical characteristics but different line-broadening parameters (due to spatial–spectral effects of the extended source along the dispersion direction) between the positive and negative arm spectra. As a first approximation, we assume that the line profiles are Gaussian and have two sources of broadening: (1) the angular extent of SNR 1987A and (2) the bulk motion of the shocked gas. In general, lines for a given element in higher ionization states form at shorter wavelengths. Also, the higher ionization states are more populated at higher plasma temperatures. The higher derived plasma temperatures may be considered to be associated with faster shock velocities. Thus, we adopt a stratification in the line-formation zone. We include this broadening from bulk motion as $2{z}_{0}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }\lambda$, where z₀, a redshift-like parameter, fits the line broadening at some fiducial wavelength, λ₀, and α measures the degree of stratification. The total line broadening can be expressed as

$$\begin{eqnarray}&&{\rm{\Delta }}{\lambda }_{\mathrm{total}}=2{\rm{\Delta }}{\lambda }_{0}\pm 2{z}_{0}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }\lambda .\end{eqnarray} \tag{ 1 }$$

Here, the first term (2Δλ₀) accounts for the broadening caused by the angular extent of SNR 1987A, which is independent of wavelength, while the second term ($\pm 2{z}_{0}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }\lambda$) accounts for broadening (or compression) by bulk motion of the gas (broadening in the MEG -1 arm and compression in the MEG +1 arm due to the 45° tilted geometry of the ER along the line of sight). We fit this line-broadening model (Equation (1)) to the measured FWHM (Section 3.2) as a function of wavelength for both the positive and negative grating spectra. Free parameters of this fit are Δλ₀, z₀, and α. An additional source of broadening due to the shock heating of heavy ions was reported in Miceli et al. (2019). As the effects of such a thermal broadening are significantly smaller than those of the bulk motion of the gas (see Miceli et al. 2019), as a first approximation, we do not include it in this model. We discuss the effects of thermal broadening in Section 4.5.

Our best-fit parameters (with 1σ errors) are: Δλ₀ = 0.0172 ± 0.0006 Å, z₀ = 0.00049 ± 0.00007, α = −0.7 ± 0.5 (MEG 2011); Δλ₀ = 0.0204 ± 0.0004 Å, z₀ = 0.00047 ± 0.00006, α = −1.3 ± 0.3 (MEG 2018). Using the grating dispersion of 0.0226 Å arcsec⁻¹ for MEG, we derive the source "half-size" θ_R to be 076 ± 001 (MEG 2011) and 090 ± 001 (MEG 2018). These angular extents are in good agreement (within statistical uncertainties) with those measured in the X-ray imaging analysis of the ACIS data taken in 2011 (∼079; Frank et al. 2016) and in 2018 (∼084; A. P. Ravi et al. 2021 in preparation). We also apply this method for the MEG 2007 data and obtain consistent (within statistical uncertainties) values of Δλ₀, z₀, and α when compared with the previously published values (Zhekov et al. 2009). For comparing shock velocities measured with emission-line widths, we consider the three available sets of MEG (±1) spectra as representative data sets from 2007 until 2018, respectively. Best-fit model plots to line widths of MEG 2007, MEG 2011, and MEG 2018 are shown in Figure 12.

Zoom In Zoom Out Reset image size

Based on these best-fit values for z₀ and α, we estimate the radial shock velocities encountered by the X-ray-emitting gas in 2007–2018. For the simple case of a strong shock with an adiabatic index γ = 5/3, the shock velocity (v_s ) is related to the bulk velocity (v_b ) of the shock-heated gas as v_s = (4/3)v_b . Assuming a radially expanding circular ring with velocity v_b , we adopt an average value for the azimuthal ϕ ($\sin {\phi }_{\mathrm{mean}}$ = 2/π, where 0 ≤ ϕ ≤ π/2) and inclination angle, i = 45°, to account for the geometric effects from viewing angles. Then, the shock velocity as a function of wavelength is

$$\begin{eqnarray}&&{v}_{s}=\displaystyle \frac{4}{3}\times 2{z}_{0}\times 1.10c\times {\left(\lambda /{\lambda }_{0}\right)}^{\alpha },\end{eqnarray} \tag{ 2 }$$

where c is the velocity of light. λ₀ is a fiducial wavelength in our broadband range, for which we choose λ₀ = 10 Å (roughly a mean between 4.5 and 15 Å). The derived average shock velocities are: ∼431 ${\left(\lambda /10\right)}^{-0.7}$ km s⁻¹ (MEG 2011) and ∼413 ${\left(\lambda /10\right)}^{-1.3}$ km s⁻¹ (MEG 2018). In Figure 13, we show the estimated shock velocity profile as a function of wavelength, plotted for MEG 2007, MEG 2011, and MEG 2018. Despite relatively large errors, the average shock velocities decreased between 2007 and 2011, which is a continuation of the trend detected in 2004–2007 (Zhekov et al. 2009). The shock velocities increase between 2011 and 2018, particularly at short wavelengths (λ < 9 Å) for Mg, Si, and S ions (Figure 13). As an example, shock velocities encountering Si ions (6–7 Å) are in the range: ∼553–616 km s⁻¹ (MEG 2011) and ∼656–802 km s⁻¹ (MEG 2018). These faster shocks since 2011 are consistent with an increase in the post-shock electron temperature (see Section 4.1). A caveat of this method is that, while we are assuming a single bulk velocity from each line, these lines more likely form from multi-velocity gas flows.

Zoom In Zoom Out Reset image size

For a fast-moving shock (with the velocity v_s ) entering a stationary packet of gas, the immediate post-shock ion temperature (kT_i ) can be expressed as ${{kT}}_{i}\approx 1.4{\left({v}_{s}/1000\,\mathrm{km}\,{{\rm{s}}}^{-1}\right)}^{2}\,\,\mathrm{keV}$, where the assumed mean molecular mass of the gas (μ) is 0.72 for the abundances of SNR 1987A. Based on the shock velocity ranges encountered by Si ions, their post-shock ion temperatures are: 0.42–0.53 keV (MEG 2011) and 0.60–0.90 keV (MEG 2018). In the non-equilibrium plasma condition, we find that these ion temperatures from our simple line-profile modeling are still significantly lower than our estimated post-shock electron temperature ranges (from broadband spectral model fits): ∼0.6–1.9 keV (MEG 2011) and ∼0.8–2.4 keV (MEG 2018). This discrepancy suggests that emission from reflected shocks in the ER (heating the inter-clump regions multiple times) may significantly contribute in the observed X-ray spectrum as proposed by Zhekov et al. (2009). The increase in the ion temperatures from 2011 to 2018 may indicate that emission from the shocked low-density gas outside of the ER contributes in the observed X-ray spectrum more in 2018 than in 2011.

Based on the high-resolution HETG (MEG) spectrum (taken in 2011) and hydrodynamic (HD) models of SNR 1987A (Orlando et al. 2015), evidence for thermal broadening of X-ray emission lines was reported (Miceli et al. 2019). Thermal broadening can be modeled by a Gaussian, whose width is directly proportional to the ion temperature and inversely proportional to the ion mass. We assume that the ion temperature is equal to the product of atomic mass of ions and proton temperature (Miceli et al. 2019). This is because at the shock front, the heavy ions have been heated up to a temperature proportional to their mass. In such a case, the thermal motion of all ions should be relatively the same, regardless of their mass. Thus, for simplicity, we assume a similar functional form to that of our Doppler-broadening term, for which the thermal-broadening parameter (z₁) is included as an additional broadening term. We modify our model (Equation (1)) to include thermal broadening, and thus, the overall broadening is expressed as:

$$\begin{eqnarray}&&{\rm{\Delta }}{\lambda }_{{\rm{t}}{\rm{o}}{\rm{t}}{\rm{a}}{\rm{l}}}\,=\,\sqrt{{\left(2{\rm{\Delta }}{\lambda }_{0}\pm 2{z}_{0}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }\lambda \right)}^{2}+(2{z}_{1}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }){\left.\lambda \right)}^{2}},\end{eqnarray} \tag{ 3 }$$

where the plus (minus) sign refers to the negative (positive) order of the MEG spectrum. Free parameters in the fit are Δλ₀, z₀, z₁, and α. We define the broadening due to the third term in Equation (3) per Å (i.e., $2{z}_{1}{\left(\lambda /{\lambda }_{0}\right)}^{\alpha }$) as the thermal broadening contribution. Observed line widths are significantly broader in the negative dispersion arm than in their positive counterpart (see Section 3.2). Therefore, the relative contribution of thermal broadening (which is the same in the two orders) is higher in the positive arm dispersion spectrum. Miceli et al. (2019) showed that, because of this effect, thermal broadening can be detected with high statistical significance in the positive order (MEG +1), while it is barely needed in the negative order (MEG -1).

In agreement with these findings, we find that thermal broadening is not statistically needed when including both the positive and negative first-order dispersed MEG 2018 spectra in our fits with the model described by Equation (3). Thus, we place an upper limit of z₁ < 0.0018 (90% C.L.) on the thermal-broadening parameter with α ∼ −0.57. The wavelength-dependent upper limits on the thermal broadening contribution in 2018 are expressed as 0.0036 ${\left(\lambda /10\right)}^{-0.57}$. In Figure 14(a), we show the upper limits on thermal broadening in the observed individual lines based on our MEG 2018 data. The assumed thermal broadening in our synthetic spectrum (as described below) is lower than the observed upper limits in 2018, which is consistent with the non-detection of thermal broadening in 2018.

Zoom In Zoom Out Reset image size

Based on a similar approach to that adopted in Miceli et al. (2019), we synthesized MEG +1 X-ray spectra of SNR 1987A as of 2018 from the HD simulation presented in Orlando et al. (2015). In this HD simulation, the thermal condition of the shocked gas is calculated self-consistently in each cell of the spatial domain during the whole evolution from the SN event to the remnant interacting with the inhomogeneous CSM. For comparisons, we considered two different cases of including and excluding the thermal broadening effects in the synthesis of the X-ray spectra. We note that these HD models predicted a significant contribution from the reverse-shocked high-velocity ejecta in the synthetic X-ray emission spectrum of SNR 1987A by 2018. This ejecta contribution resulted in previously unseen broad wings in our synthetic spectral lines calculated for 2018, even without the thermal-broadening effects. However, because our deep HETG spectrum taken in 2018 shows no evidence for enhanced elemental abundances (see Section 3.1), SNR 1987A has unlikely arrived at such an ejecta-dominated phase yet. Thus, we fit these extra broad-wing features in the synthetic spectra of 2018 with a separate Gaussian component. We exclude this extra Gaussian component from our discussion.

Figure 14(b) shows the comparisons between the observed line widths and the synthetic spectral line widths (with and without thermal broadening), for MEG +1 in 2018. We note that including thermal broadening in the HD simulated spectrum better describes these observed line widths (in agreement with Miceli et al. 2019; see Figure 14(b), with an exception at ∼8.5 Å). Detailed studies of thermal broadening contributions will be important for the forthcoming generation of X-ray spectrometers based on microcalorimeters.

HD/MHD models of SN 1987A suggest that the contribution of the reverse-shocked outer layers of ejecta would soon (∼35 yr after the SN) become dominant in the soft (0.5–2 keV) X-ray band (Orlando et al. 2015, 2019, 2020). These models predict that such an ejecta-dominated phase will be marked by a significant change in the observed X-ray spectrum, evidently showing strong enhancements in the individual line fluxes and thus in the metal abundances (compared to those of the ER emission). Our estimated metal abundances based on the broadband spectral model fits are consistent (within statistical uncertainties) with those for the ER as measured in previous works (Zhekov et al. 2006; Dewey et al. 2008; Zhekov et al. 2009). Thus, we conclude that the onset of such an ejecta-dominated era in the evolution of SNR 1987A has not been ensued as of 2018. We note that the detection of Fe K emission lines possibly from the reverse-shocked overabundant ejecta has been suggested based on earlier XMM-Newton observations (Sun et al. 2021). We cannot confirm the detection of such Fe K lines based on our deep Chandra HETG observations due to the low responses of the HETG spectrometer in the ∼6 keV energy band.

Recently, NuSTAR observations of SNR 1987A showed evidence for hard X-ray emission up to E ∼ 20 keV (Alp et al. 2021; Greco et al. 2021). The NuSTAR-detected hard X-ray spectrum may be attributed to either nonthermal (PL with photon index Γ ∼ 2.5, Greco et al. 2021) or thermal (kT ∼4 keV, Alp et al. 2021) origins. Our Chandra grating data do not have significant count statistics at E ≥ 5 keV and a third hard-component X-ray emission is not statistically required to model our Chandra HETG/LETG spectra of SNR 1987A. In fact, the contribution from the third component (either the nonthermal PL component with Γ ∼ 2.5 in Greco et al. 2021 or the thermal component with kT ∼ 4 keV in Alp et al. 2021) in our observed Chandra HETG/LETG bandpass is negligible (≤4% of the observed flux in the 0.5-5 keV band), and thus it does not affect our conclusions. Discussing the true nature of the NuSTAR-detected hard X-ray emission at E ≥ 10 keV from SNR 1987A is beyond the scope of this work. Nonetheless, we note that electron temperatures estimated for the soft X-ray spectrum (E ∼ 0.5–8 keV) in Greco et al. (2021; kT ∼ 0.6 and 2.6 keV) are in better agreement with our results than those in Alp et al. (2021; kT ∼ 0.5 and 0.9 keV), although these comparisons need to be taken with caution due to different spectral modeling among them (i.e., two-component model fits in this work versus three-component models in Greco et al. (2021) and Alp et al. (2021).