Smoothed Particle Inference Analysis of SNR DEM L71

TL;DR

This paper introduces a novel Smoothed Particle Inference method to analyze the complex three-dimensional structure of supernova remnant DEM L71 using X-ray data, revealing detailed properties and substructures.

Contribution

The paper presents the first application of Smoothed Particle Inference to supernova remnants, enabling comprehensive characterization of their structure, composition, and dynamics in a single analysis.

Findings

Confirmed high Fe abundance in the ejecta, supporting Type Ia classification.

Results are consistent with numerical simulations and previous studies.

Revealed irregular substructures in temperature and other parameters not seen with traditional methods.

Abstract

Supernova remnants (SNRs) are complex, three-dimensional objects; properly accounting for this complexity when modeling the resulting X-ray emission presents quite a challenge and makes it difficult to accurately characterize the properties of the full SNR volume. We apply for the first time a novel analysis method, Smoothed Particle Inference, that can be used to study and characterize the structure, dynamics, morphology, and abundances of the entire remnant with a single analysis. We apply the method to the Type Ia supernova remnant DEM L71. We present histograms and maps showing global properties of the remnant, including temperature, abundances of various elements, abundance ratios, and ionization age. Our analysis confirms the high abundance of Fe within the ejecta of the supernova, which has led to it being typed as a Ia. We demonstrate that the results obtained via this method are consistent with results derived from numerical simulations carried out by us, as well as with previous analyses in the literature. At the same time, we show that despite its regular appearance, the temperature and other parameter maps exhibit highly irregular substructure which is not captured with typical X-ray analysis methods.

Figures

Tables

Table 1. Table 1: XMM-EPIC Observation 0201840101

Detector	Events	Net Exposure (ks)
MOS1	175100	58.1
MOS2	357621	58.3
pn	1034574	58.3

Table 2. Table 2: Model Priors

Parameter	Minimum	Maximum	Global
log kT (keV)	0.08	7.0	No
nH,LMC (1020${\text{cm}}^{-2}$)	1.0	10.0	Yes
log net (cm-3s)	1010	5$\times$1012	No
O/O⊙	0.01	1.0	No
Ne/Ne⊙	0.1	1.5	No
Mg/Mg⊙	0.1	1.5	No
Si/Si⊙	0.01	3.0	No
S/S⊙	0.01	3.0	No
Fe/Fe⊙	0.01	3.0	No

Table 3. Table 3: EM-weighted summary statistics.

	Median	Mode	Standard Deviation
nH,LMC (${10}^{20}$ cm-2)	6.7$\pm$0.9	9.0$\pm$1.7	2.4$\pm$0.3
kT (keV)	0.24$\pm$0.04	0.15$\pm$0.02	0.43$\pm$0.07
log $n_{e}t$ (cm${}_{}{}^{-3}s$)	12.4$\pm$0.1	12.7$\pm$0.2	0.6$\pm$0.1
O/O⊙	0.19$\pm$0.04	0.02$\pm$0.08	0.27$\pm$0.02
Ne/Ne⊙	0.45$\pm$0.12	0.11$\pm$0.30	0.35$\pm$0.05
Mg/Mg⊙	0.75$\pm$0.14	0.95$\pm$0.36	0.37$\pm$0.05
Si/Si⊙	1.19$\pm$0.30	0.40$\pm$0.77	0.79$\pm$ 0.10
S/S⊙	1.43$\pm$0.27	1.54$\pm$0.74	0.77$\pm$0.09
Fe/Fe⊙	0.74$\pm$0.31	0.16$\pm$0.85	0.83$\pm$0.12

Equations

$${\delta }_F=\frac{1}{N_i}\sum _i^{N_i}[F_i-\overline{F}{]}^2,$$

Full text

circumstellar matter — ISM: supernova remnants — methods: data analysis — X-rays: individual (DEM L71) — X-rays: ISM

††journal: ApJ††facilities: XMM(EPIC)

1 Introduction

Supernova remnants (SNRs) are invaluable tools for studying stellar evolution, supernovae, and the interstellar medium. As a SNR expands into the interstellar medium (ISM), the shock probes the ambient medium. Heavy elements produced by nucleosynthesis in the progenitor star, and in the explosion, are carried out by the shock wave. The outgoing fast shock reaches velocities of thousands of kilometers per second, with the resulting post-shock temperatures exceeding millions of degrees. As a result, the best wavelength range to probe the bulk of the remnant is X-rays. However, SNRs have complex morphology, in 3 dimensions. They also have correspondingly complex spectral properties, typically with significant variation across the face of the remnant, as well as along the line-of-sight. Relative to this spatial and spectral complexity, X-ray observations of SNRs typically provide very few photons. These factors combine to make accurately measuring the physical conditions throughout the entire SNR volume a very challenging endeavor.

The SPIES (Smoothed Particle Inference Exploration of Supernova remnants) project aims to address this challenge by applying a fundamentally different approach to analyzing X-ray observations of SNRs. Smoothed Particle Inference (SPI) was developed specifically for XMM-Newton (Peterson et al., 2007) and it simultaneously takes advantage of both the spatial and spectral information available in the data. SPI is a Bayesian modeling process that fits a population of gas blobs (“smoothed particles”) such that their superposed emission reproduces the observed spatial and spectral distribution of photons. The results of this completely independent approach can be compared with results from numerical simulations as well as results derived from more traditional analysis techniques. Here we demonstrate the unique capabilities of SPI and compare it to more traditional X-ray analysis methods by applying SPI to the Type Ia supernova remnant DEM L71.

DEM L71 was first identified as an optical SNR in the Large Magellanic Cloud (LMC) more than 40 years ago by Davies et al. (1976). It was detected in X-rays in an Einstein survey by Long et al. (1981). Since then, it has been further observed in X-rays using ASCA (Hughes et al., 1998), Chandra (Hughes et al., 2003; Rakowski et al., 2003, hereafter H03 and R03), and XMM-Newton (van der Heyden et al., 2003; Maggi et al., 2016, hereafter vdH03 and M16).

The X-ray morphology is remarkably regular, with a bright outer rim and a faint, diffuse center, usually interpreted as the forward shock and reverse-shocked ejecta, respectively. Most previous analyses have therefore characterized the physical conditions in DEM L71 with two sets of spectral parameters: one for the outer rim, and one for the center emission. ASCA observations had indicated that the remnant ejecta were enhanced in Fe (Hughes et al., 1998), which was further confirmed by subsequent high-resolution observations with Chandra (H03, R03) and XMM-Newton (vdH03). These observations, along with associated Fe ejecta mass estimates of $\sim$1.4M*⊙*(H03, vdH03), suggest that the progenitor of DEM L71 was a Type Ia supernova (SN), which are expected to have high Fe abundances. Optical observations of the Balmer-dominated shock velocities by Ghavamian et al. (2003, hereafter G03) yield an age estimate of $\sim$4400 years. Throughout this work, we assume a distance to the LMC and DEM L71 of 50 kpc.

2 Numerical Hydrodynamic Simulations

In order to better interpret the results of the SPI analysis, we have carried out spherically symmetric one and two-dimensional numerical hydrodynamic simulations of the evolution of DEM L71. Analytical approximations can help to understand the overall evolution of the remnant, but are unable to capture the complex dynamics and kinematics. Numerical simulations are necessary to investigate the detailed morphology, deviations from symmetry, and hydrodynamical instabilities. They can be used to study the temperature, density, and velocity profiles, and the small-scale structure within the remnant, and to provide an adequate theoretical framework to analyze our SPI results. Results from the these simulations can be used to determine distinguishing properties of the different components of the SNR, such as the contact discontinuity and the ejecta, which can then be used to identify these components in our SPI results.

The explosion of a star to form a supernova (SN) gives rise to a forward shock that propagates outwards into the medium with a high velocity, and a reverse shock that moves back into the ejecta in a Lagrangian sense. The two are separated by a contact discontinuity (CD) that separates the shocked ejecta from the shocked circumstellar medium (CSM). The decelerating contact discontinuity is subject to Rayleigh-Taylor (R-T) instabilities that can result in some mixing of shocked ejecta into shocked CSM. The density profile of the ejecta and of the CSM are essential to deciphering the shock structure, and kinematics and dynamics of the evolution (see, for instance, Dwarkadas, 2011, and references within).

Observations around the rim of the remnant have detected a wide variation in the expansion velocity of the remnant, ranging from around 400 to about 1250 km s*-1* (Rakowski et al., 2003, 2009). These suggest an asymmetry in the expansion. Our spherically symmetric simulations are not designed to reproduce such an asymmetry; therefore we have not attempted to match all the parameters. Rather, our goal was to carry out simulations that allow us to adequately understand the histograms and maps generated from the SPI analysis, and provide a framework to interpret the output from the SPI analysis. Our simulations were run using the VH-1 code, a 1, 2, and 3-dimensional finite-difference numerical hydrodynamics code, based on the Piecewise Parabolic Method of Colella & Woodward (1984). The initial setup is very similar to that used in Dwarkadas & Chevalier (1998) and Dwarkadas (2000). Description of the code is available in these papers.

H03 identify the contact discontinuity in DEM L71 as the edge of an inner region of harder X-ray emission, and find that it has a mean radius of 4.3 pc, which is about halfway out to the forward shock. The reverse shock will be at a much lower radius. Since the ejecta are expanding homologously, their density decreases as t*-3*, which means that the density has decreased substantially, given the estimated age of the remnant. The reverse shock expanding into this low density would have moved back into the ejecta, suggesting that much of the ejecta has been recently shocked and the reverse shock has recently reached, or is close to, the center of the explosion.

Dwarkadas & Chevalier (1998) have shown that the ejecta density in a Type Ia SN explosion is best represented by an exponential profile. We used this profile to describe the SN ejecta expanding into a constant density medium with a density \sim$$0.5 cm*-3*, as suggested by G03. Our simulations used an expanding grid, such that the grid expands outwards with the flow; this is extremely useful in cases like this where the remnant size grows by 4 orders of magnitude over the simulation time. The main constraints are the radius of the remnant, which is somewhere between 8.6 – 10 pc, and the forward shock velocities, which are of order a few hundred km s*-1*.

We first carried out simulations assuming canonical parameters such as a standard energy of 1051 erg, and a mass of 1.4M*⊙*. However, at the remnant radius of even 10pc, we find that the velocities are higher than the highest observed velocities, at an age of about 3600 years. This implies either a much higher density, or a lower explosion energy. Following the results of G03, we lowered the energy to 0.4 $\times 10^{51}$ erg, and reran the simulations. These provide much better agreement with the available data.

In Figure 1 we show the pressure and density profile within the remnant at an age of 3987 years from our simulation. The age and the shock velocity are consistent with the results of H03, which is not very surprising since they used a similar model in their calculations.

The forward shock of the remnant is at about 8.5pc. In our expanding grid simulation, the reverse shock has crossed inward of the inner boundary, and presumably reached the center or is very close to it. The solution in a few zones near the inner boundary reflects that and should be ignored. At around 4.3pc lies the contact discontinuity (CD), which separates the shocked ejecta from the shocked circumstellar medium (CSM). Everything inward of the CD is ejecta, thus showing that even at this late stage the (shocked) ejecta are clearly visible. Everything outside the CD is the shocked CSM. We note that the maximum density is immediately behind the circumstellar shock. The density decreases inwards till it reaches the CD, where it increases again, before decreasing towards the center. The majority of the CSM can be seen to have a higher density than most of the ejecta. The pressure drops somewhat behind the outer shock but is then mostly constant throughout the remnant. The temperature structure will then be essentially the inverse of the density structure.

In Figure 2 we show the gas temperature within the remnant. It must be emphasized that this is the fluid temperature, and is obtained by assuming a constant value of the mean molecular weight throughout the remnant. One would expect that the ejecta would generally be higher in metals than the CSM and therefore the mean molecular weight would be different. More importantly, this does not reflect the electron temperature, which may be lower if electron-ion equilibration has not been established, as is most likely the case. However, it provides us some idea of the temperature profile throughout the remnant, and indicates that on average the ejecta have a higher temperature than the CSM. There will be a small region of shocked circumstellar material with a much higher temperature than the ejecta, but the mass there is very small compared to the total CSM mass.

It is clear that the material nearest to both the shocks forms the most recently shocked material. This material will consequently have the lowest ionization age $n_{e}t$. However, the low ionization age material behind the outer shock will have essentially the highest density, whereas the low ionization age material near the reverse shock, closer to the center of the remnant, will have the lowest density. Since the ejecta generally have a lower density, they will have on average lower ionization age than the CSM. Similarly, we expect that the ejecta, if this is truly a Type Ia remnant as emphasized by many of the above papers, should display higher metallicity, and specifically high abundances of Fe and Si as expected from Type Ia SNe. These distinguishing characteristics can be used to separate the ejecta from the CSM during the SPI analysis.

In the spherically symmetric simulations, the CD is a sharp surface that cleanly separates the shocked ejecta from the shocked CSM. In reality however the decelerating CD is unstable to the Rayleigh-Taylor (R-T) instability (Chevalier et al., 1992), leading to the growth of R-T ‘fingers’ extending out from the shocked ejecta into the shocked CSM. Shearing between the fingers and the surrounding fluid gives rise to the Kelvin-Helmholtz instability, and ‘mushroom caps’ can be seen atop the R-T fingers. The growth of R-T instability in SNRs of Type Ia was studied by Dwarkadas (2000) for an exponential ejecta profile. In order to study the multi-dimensional evolution and kinematics, we have run a two-dimensional (2D) simulation using approximately the same parameters as the 1-dimensional run. The simulation was run with 500 radial and 500 angular zones on an expanding grid. As expected, the CD was found to be unstable to the growth of R-T instability. In Figure 3 we show a frame from the 2D simulation at approximately the same time as the 1D run above. The image shows the temperature distribution across the grid. The temperature profile follows more or less the 1D profile in Figure 2, except at the contact discontinuity, where the turbulence due to the R-T and K-H instabilities results in considerable temperature variation around the CD, due to the presence of the R-T fingers that mix shocked ejecta with the shocked CSM. The effect of the temperature (and density) variations is to effectively widen the CD, which is now no longer a sharp spherical surface but has a width of around 10-15% of the remnant radius (0.9- 1.4pc), centered at about 4.5 pc.

3 SPI Analysis

3.1 Observations and Data Reduction

We use the XMM EPIC observation 0201840101 from 2003 December, the longest available observation of DEM L71. The exposure-corrected image is shown in Figure 4. Our SPI analysis requires only a filtered photon event list and exposure map for each detector (MOS1, MOS2, and pn). These were created using SAS 16.0, primarily the tasks emchain and epchain. To exclude non-X-ray events, the event lists were filtered to include only photon event patterns 0-12. They were filtered further to exclude events with energies outside of the 0.2-10 keV range. Periods affected by soft proton flares were removed. The resulting number of events and exposure for each detector are shown in Table 1. For the SPI analysis, only events within a $150^{\prime\prime}\times 150^{\prime\prime}$ box centered on DEM L71 were included. DEM L71 is approximately circular with a radius of $\sim$45*′′*.

3.2 Smoothed Particle Inference

Smoothed Particle Inference (Peterson et al., 2007) owes its flexibility to its modeling of the gas as a collection of independent ‘smoothed particles,’ or blobs  of gas. Emission from all blobs is combined to represent the net X-ray emission from the SNR. Each blob has its own spectral and spatial model, typically including the gas temperature, abundances, spectral normalization, Gaussian width, and spatial position. Multiple blobs can occupy the same line of sight, or even the same space, providing the ability to model a multiphase gas. While the blobs themselves are spherically symmetric, no particular morphology or symmetry is assumed in their arrangement. The size and position of each blob can be adjusted independently, which allows the multi-blob model to reproduce any arbitrary distribution (within limits imposed by the quality of the X-ray observation). This method is quite powerful, as it becomes a relatively straightforward task to characterize the distributions of any number of gas properties from the posterior distributions of the relevant blob parameters. SPI was specifically developed for the analysis of XMM-Newton RGS and EPIC observations (Peterson et al., 2007), and has been previously applied to observations of galaxy clusters (Andersson et al., 2007, 2009; Frank et al., 2013).

3.3 Fitting Procedure

SPI uses a Markov Chain Monte Carlo (MCMC) process to forward-fold the blob model and predict detector positions and energies for each photon. In this way SPI utilizes all available spectral and spatial information without imposing artificial restrictions on either the spatial or spectral distribution of the gas emission. A set of emitted photons is simulated for each blob, according to each blob’s spectral and spatial model. These are then propagated through the instrument’s response functions, including mirrors, gratings, and detector responses as necessary. Given the energy and location of each of these model photons, the probability of its detection is calculated for the given instrument response, as is its predicted location on the detector and the measured energy. The result is a set of model X-ray events which is compared to the observed events by binning both on the same three-dimensional grid (two spatial dimensions and energy). The two are then compared by calculating a two-sample likelihood statistic that quantifies the difference between the model and observed datasets. A modified Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970) is used to choose a new set of blob parameters and the process is repeated. For more details of the SPI implementation, see Peterson et al. (2007). After a number of iterations, the MCMC chain converges and the posterior parameter distributions become stable. After convergence, the fit typically has reduced $\chi^{2}$ values of 1.0-1.5. From this point onwards, each iteration returns a model that is statistically consistent with the data.

The number of blobs is fixed and chosen based on the spectral and spatial complexity of the emission, as well as the number of photons available in the data. Previous work with galaxy clusters has suggested that the optimal number of blobs is such that the number of model photons per blob should be at least $10^{4}$ (Frank et al., 2013); this number can be increased (and computation time decreased) by reducing the number of blobs. However, this criterion must be balanced with the need to have enough blobs to adequately fit the more complex emission of SNRs. For our DEM L71 analysis, 50 blobs are used, resulting in \sim$$3\times 10^{4} photons per blob.

The SPI fit of DEM L71 results in a reduced $\chi^{2}=1.15$. We ensure that 50 blobs provides sufficient accuracy by also testing with 100 blobs and comparing the posterior distributions. Examples for temperature, absorbing column density, and ionization age are shown in Figure 5. It is expected that the posteriors will not be identical, but they are extremely similar, indicating no improvement or significant change in the overall results when going from 50 to 100 blobs. The reduced $\chi^{2}$ is also $\sim$1.15 for both the 50 and 100 blob models, and the medians, modes, and standard deviations of the parameter posteriors are nearly identical. The model spectrum also reproduces all the main features of the observed spectrum (Figure 6).

3.4 Model

With SPI, a combination of spectral and spatial models may be chosen to suit the particular astrophysical context, in this case an absorbed thermal plasma which is likely not in ionizational equilibrium. We also include model components to simultaneously fit the X-ray background. The model components are described in more detail below. The model is applied independently to each blob; however some model parameters may be global (the same for all blobs) and some are frozen. The spectral normalizations of each model component are always free parameters. All prior distributions are flat. The end points of the priors are chosen via an iterative process; initial guesses are made based on typical physical conditions in SNRs and anything known from previous X-ray analyses. We then test the model by performing an SPI fit and investigating the posterior distributions. The priors are adjusted if necessary and the MCMC restarted, to ensure that the priors span the entire relevant parameter space. For example, if the fit places an excess of blobs at the upper limit of the prior (the largest allowed values), this upper limit is increased to allow exploration of higher values. Similarly, the limits are decreased if it is clear the parameter space spanned by the prior is unnecessarily wide, as this improves the efficiency of the MCMC.

3.4.1 Thermal Supernova Remnant Emission

The X-ray emission from DEM L71 results from shocked plasma. In the case of DEML71, the emission is known to be thermal (H03, vdH03). However the gas may not necessarily be in ionization equilibrium, meaning that the ionization state of the gas may not be consistent with the post-shock temperature. Therefore we use the thermal, non-equilibrium ionization (NEI) plasma model vpshock (Borkowski et al., 2001). The free spectral parameters for each blob, summarized in Table 2, include temperature and ionization age (net), as well as the abundances of O, Ne, Mg, Si, S, and Fe (relative to solar abundances from Anders & Grevesse, 1989). All other abundances are fixed at typical LMC values of 0.5 solar (Russell & Dopita, 1992). Free parameters may be different for every blob, unless marked as global in Table 2. Note that temperature and ionization age are chosen from a logarithmically uniform prior. Volume emission measures can be easily derived directly from the spectral normalization of each blob. Galactic absorption is included via the phabs absorption model, and is fixed to the Dickey & Lockman (1990) value of n{}_{H}=5.7\times 10^{20}$$\mbox{cm}^{-2}, as it is not expected to vary substantially across such a small object. Absorbing material within the LMC is accounted for with a zvphabs model, with abundances set to average LMC values of 0.5 solar. This column density is tied for all model components, but is allowed to vary from blob to blob, to allow for the possibility of small spatial variations. Several of the background components also include Galactic and LMC absorption. The column densities for these components are tied to those associated with the vpshock component. Spatially, the vpshock emission is modeled as coming from a spherical Gaussian distribution (the ‘blob’). The central coordinates in the plane of the sky are allowed to vary independently for each blob anywhere within the $150^{\prime\prime}\times 150^{\prime\prime}$ box centered on DEM L71  while the logarithmic prior for the Gaussian width of each blob is ln 2.$\prime\prime$5$<\sigma<$ ln $75^{\prime\prime}$ (thus the minimum diameter is $\sim 5^{\prime\prime}$).

3.4.2 X-ray Background

The X-ray background is a result of several sources of emission, including instrumental, cosmic, and Galactic. Each is accounted for with an additional model component. Background model components are fixed and the same for all blobs  except the spectral normalizations.

The EPIC instrumental background consists of soft protons, internal line emission, and electronic noise, and is accounted for with a custom model which includes all three components, as used in Frank et al. (2013). The model is described in detail in Andersson et al. (2007).

The Galactic X-ray background is soft, with energies $\lesssim 1$ keV, and consists of several thermal components with $kT<0.5$ keV (Lumb et al., 2002; Kuntz & Snowden, 2000; Snowden et al., 2008). We use two model components, one to account for the Local Hot Bubble (LHB) and one for slightly warmer emission from the Galactic Halo. Both employ thermal MEKAL spectral models, uniform across the field of view. The temperatures of the LHB and Halo are set to kT$=0.1$ keV and $0.25$ keV, respectively (Snowden et al., 2008). The LHB is unabsorbed, while Galactic absorption is applied to the Halo emission.

The final background component results from unresolved cosmic sources, mainly AGN, and contributes primarily at energies $\gtrsim 1$ keV. These are modeled as a powerlaw, with $\Gamma=1.47$, that is spatially uniform. The powerlaw is absorbed by both Galactic absorption and absorption within the LMC.

3.5 Post-SPI Analysis

After a successful SPI fit, information about the SNR can be derived directly from the posterior parameter distributions, which include all blobs from every iteration after convergence. Summary statistics, such as the medians and standard deviations, are the simplest to derive. Medians are similar to the averages that are measured with more traditional methods which extract a single spectrum from the entire SNR, or from several lines of sight through the SNR, and fit to a simple model. Standard deviations of each parameter can also be calculated, and are one of the advantages in using SPI, as they provide a straightforward measurement of how homogeneous the gas conditions are.

Taking a step further, it can be even more informative to investigate the parameter distributions directly. SPI allows the construction of continuous parameter distributions (the posteriors), the shapes of which contain useful information on the gas composition and state. Distributions can be constructed from the set of blobs by essentially creating histograms of the parameters of interest.

In addition, because each blob has a location and size, it is possible to create maps of any of the parameters. Maps are created by selecting a region containing the SNR and dividing it into spatial bins of a specified size (typically 1 to 5 arcsec). For each spatial bin, the contributions of each blob are combined to determine the value. A number of different methods can be used to combine the blobs within a spatial bin, depending on the purpose. The most common is to take the median value (e.g. the median of the blob temperatures), but it is also possible to use other quantities, including the maximum, minimum, sum, or standard deviation. The standard deviation is of particular interest, as it provides a means of mapping the level of homogeneity across the remnant. Note that these maps are all fundamentally different from X-ray images (including narrow-band images) or equivalent width images, which show only the intensity of related emission rather than measured values of the parameter. Narrow band and equivalent width images are fast methods for mapping the location of different elements or relative temperatures but do not measure the differences in abundances or other properties across the remnant.

Associated statistical uncertainties can also be derived, as the probability density for each parameter is represented by the collection of models from all individual iterations. The uncertainty of a value $F$, $\delta_{F}$, can thus be characterized as the standard deviation of the distribution of $F$ over converged iterations,

[TABLE]

where $N_{i}$ is the number of iterations after convergence, the sum is over these iterations, and $F$ is the measurement of interest, a function of the blob parameters, for example the overall median or mean temperature. Systematic uncertainties are much more difficult to estimate. The shape of the posterior distributions can provide some information. The posterior shape is influenced by a several factors in addition to the true parameter distribution. First is the shape of the prior; the prior distributions are all uniform, and thus the further the posterior deviates from uniform, the more confident we can be that the SPI fitting has successfully constrained the distribution. The converse, however, is not true; a uniform posterior may be the result of the true distribution also being uniform, and in this case it is not distinguishable from the prior. Second, there is some broadening of the posterior distribution due to the imperfect quality of the data and the inherently incomplete spectroscopic information about the gas properties that can be derived from atomic transitions. This broadening will be smaller for higher quality data and can be minimized by choosing an optimal number of blobs (Frank et al., 2013). Because of the difficulty in assessing the details of these systematic uncertainties, small scale features of the distributions should be interpreted with caution, but the basic shape of the distributions are sufficiently robust to be informative.

The raw summary statistics, distributions, and maps, for which all blobs are treated equally, are useful for exploring the results and identifying minor but important components of the SNR gas. However, it is usually more meaningful to first weight each blob by its emission measure (EM), to account for the different sizes and densities and thus better represent the true contribution of each blob to the overall SNR.

SPI also provides the capability to explore any subset of the parameter space by selecting a set of blobs based on any combination of parameters and investigating the summary statistics, distributions, and maps of only that subset. For example, it is possible to select only blobs in a particular temperature or ionization age range (or both) and then investigate the corresponding abundance distributions or maps. Filtering in this way has huge advantages for isolating and measuring the properties of specific physical components in the SNR, as it eliminates confusion from unrelated material that may be along the same line of site. This is one of the most powerful and unique capabilities of SPI.

4 Application of SPI to DEM L71

We detail here the results of the application of the SPI analysis, as outlined in Section 3, to the SNR DEM L71.

4.1 Summary Statistics

The EM-weighted summary statistics are shown in Table 3. The volume emission measure for the entire remnant is $EM=5.53(\pm 0.59)\times 10^{59}$ cm*-3*. These simple quantities provide a very limited view of the actual parameter distributions, shown in Figure 7, but they are useful for comparison with more standard spectral fitting methods.

In most cases, large standard deviations and/or substantial differences between the median and mode indicate the distribution is complex enough that a simple mean or median is a poor descriptor. For example, the median EM-weighted temperature is 0.54 keV, but a very different mode (0.22 keV) indicates the actual distribution is more complicated. Similarly, O, Ne, Mg, and Fe EM-weighted median abundances are all roughly similar to typical LMC values ($\sim$0.5 solar), while Si and S are higher. If this was the only information available, it might be concluded that O, Ne, Mg, and Fe are all from shocked CSM, while the Si and S is shocked ejecta. However, in the cases of O, Ne, and Fe, the other summary statistics suggest the situation is not so simple; the modes are all very low, implying there are regions with DEM L71 that are largely devoid of O, Ne, or Fe, and the standard deviation of Fe is large enough to indicate the presence of not only gas with very low Fe, but also very high Fe. The Si mode of 0.40 is more indicative of shocked CSM than ejecta.

4.2 Parameter Distributions

Emission-measure-weighted posterior distributions are shown in Figure 7. Typically some blobs are present throughout the entire allowed parameter space. The temperature distribution reveals most gas has temperatures \lesssim$$1 keV, but some does exist at higher temperatures with substantially lower EM. Ionization age tends to be above 1012, indicating much, but not all, of the gas is in ionization equilibrium111Ordinarily the fact that the ionization age peaks at the highest allowed values would suggest the prior should be extended. However, in this case all higher values would be in ionization equilibrium, and thus increasing the range is not justified.. O shows a clear peak at $\lesssim$0.3, in line with the existence of some O-deficient gas, while the majority of O emission comes from the shocked CSM. The Ne distribution is similar, but with a more gradual decline at higher values. Mg peaks near typical LMC values of $\sim 1$ and drops off sharply above solar. Si also tends to have abundances $\lesssim 1$, but with a substantial contribution from higher abundance material. S has a rather wide distribution, and likely is not well-constrained, but the distribution suggests that S abundances $>2.5$ are unlikely. Most gas in DEM L71 clearly has low Fe, $\lesssim 1$, but some exists with higher Fe abundances as well, up to \sim$$2.5.

4.3 Maps

EM-weighted parameter maps are shown in Figures 8 and 9. Overall, the center of DEM L71 is hotter, with lower EM, and larger Fe, Ne, and O. Apart from these general trends, the maps are all highly irregular. The outer rim is particularly clumpy, with azimuthal variations of at least a factor of two in all parameters.

4.4 Identifying Ejecta Emission

The substantially higher EM of the clumps in the outer rim implies that the EM-weighted distributions and maps (Figures 7-9) are dominated by the outer rim material. However, the parameter distributions (Figure 7) clearly show the presence of of some lower EM gas. Maps which are not weighted by EM, and therefore give equal weight to both high and low EM blobs, can provide some insight on the properties of the lower EM gas. Unweighted temperature and Fe maps (Figure 10) indicate that the lower EM core gas has both higher temperatures and high Fe. As demonstrated in the simulations (Figure 2), this is expected for emission from the shocked ejecta, which is also expected to have lower density and lower $n_{e}t$. Therefore to investigate the core gas, we selected only blobs with kT$>1$ keV and Fe$>$1, a combination which produces a clearly centralized spatial distribution (Figure 11, left). Selecting blobs with particular characteristics to measure properties of the ejecta dramatically reduces contamination in these measurements from overlying CSM material. Reducing either the temperature or Fe thresholds results in the inclusion of significant amounts of the high EM gas in the outer rim, and thus will not represent the core emission, while the maps and parameter distributions with higher thresholds look nearly identical to those with kT$>1$ keV and Fe$>$1. In addition to being warmer with higher Fe, the core emission has distinctly lower ionization ages, $\lesssim 10^{12}$, and higher O, as can be seen in the associated parameter distributions (Figure 12). The temperature and ionization age maps both show a gradient, with the hottest gas and highest ionization ages in the southeast and lowest in the west (temperature) and northeast (ionization age).

Our 2D simulations illustrate that the unstable CD, which is spread over a wide region due to R-T instabilities, has the largest temperature variations (Figure 3). Using SPI, we can measure and map the temperature variations, as seen in Figure 8 (center right). This is a map of the standard deviation of the temperature within each spatial bin, and clearly shows a wide hook-shaped structure in the outer core where the variations are largest. Both the width of this region, $\sim$1.7pc, and its distance from the center, $\sim$4.5pc, agree very well with the location and width of the CD from the simulations (Figure 3). Although the widths from the SPI map and the 2D simulations are both upper limits, and the hook shape itself cannot be reproduced in our spherically symmetric simulations, the agreement on the size and location of this region with maximum temperature variations allows it to be identified as the location of the CD and lends further support for the core emission inside this region to be arising from the SN ejecta.

It is also clear from the simulations that the temperature in this region will be higher than in the rest of the ejecta, and indeed both the overall (Figure 8) and ejecta (Figure 11) maps indicate that the highest temperatures are found there. M16 used two temperature components to fit the ejecta, while vdH03 used a shell like structure for the core region. Although the above discussion provides some understanding of the basis for these assumptions, we emphasize that neither of them is required when using the SPI technique. Any asymmetry in the shape of the CD region in SPI will be due to asymmetries in the ejecta or the medium into which the remnant was expanding, or due to projection effects, and are not captured in our simulations.

5 Discussion

Previous measurements of the physical properties of DEM L71 rely on extracting X-ray spectra from a specified regions and fitting a NEI spectral model, with one to three components to represent the CSM and ejecta contributions. Ejecta properties are measured by designating one or more of the spectral components to be ejecta (e.g. by allowing higher abundances). H03, R03, vdH03, and M16 all take this approach using XMM-Newton (vdH03 and M16) and Chandra (H03 and R03) observations. Both M16 and vdH03 extract and fit an EPIC-pn spectrum from the entire source region. To investigate CSM vs ejecta properties, vdH03 also extract and fit separate spectra for two regions, an outer annulus and a core region, while M16 instead use a three-component model to fit their single spectrum (one CSM component and two ejecta). Rather than fit spectra that encompass the entirety of DEM L71, H03 and R03 take advantage of the higher resolution of Chandra to extract spectra from several small regions. H03 choose 4 small regions in the southeast to represent the inner and outer CSM and the inner and outer ejecta, fitting the spectra from each independently. R03 extract spectra from regions in the bright outer shell with 5 different azimuthal positions, representing the CSM, and carry out a joint fit to determine abundances. Together, these works represent the most common approaches to X-ray analysis of SNRs, and thus make a useful basis for comparison to assess the efficacy of SPI.

All of these analyses result in discrete measurements of the properties of interest. The SPI equivalents to these measurements are the EM-weighted medians as in Table 3. We therefore expect that these values should be similar to those reported by other authors, and indeed we find this is generally true. Figures 13 and 14 compare our medians with the equivalents reported in other works, both for the ejecta and for the CSM. In addition to the simple medians, we show shaded gray bars representing the EM-weighted posterior distributions as in Figure 12. The overall picture is consistent: a bright shell of CSM surrounding a warmer core. The high EM shell is cool, kT $\sim 0.3-0.5$ keV, with ionization ages of at least a few times $10^{11}$cm*-3*s and subsolar O, Ne, Mg, and Fe abundances that suggest shocked CSM. The core has significantly lower EM and is warmer (kT $\sim 1-2$ keV) with lower ionization ages and enhanced abundances. However, there are a number of differences between previously reported results and the results reported here that are informative.

M16, vdH03, and this work all use XMM-Newton EPIC-pn observations. The total emission measure obtained from our SPI analysis is 5.53($\pm 0.59)\times$1059cm*-3*. This matches that from M16, 5.68${}^{+0.07}_{-2.90}\times$1059cm*-3*. For their fit to the full EPIC-pn spectrum, vdH03 find a total emission measure of 3.54($\pm 1.05)\times$1059cm*-3*, which is somewhat smaller but still reasonably close. Separating the CSM and ejecta emission, the total EM of the ejecta found in this work agrees with that of vdH03, but the EM of the CSM is substantially lower in vdH03 (see Figure 13). The combined EM from their separate CSM and ejecta spectral fits is only 1.13${}^{+0.23}_{-0.04}\times$1059cm*-3*, substantially lower than the total from their fit to the integrated spectrum. It is clear that such a model, where all the X-ray emission comes from two shells, ignores a substantial fraction of the remnant emission (everything outside the individual regions). M16 report a similar CSM EM, but much higher ejecta EM than either vdH03 or this work. These differences are most likely due to the very different assumptions made during the spectral fits, both in the number of components and the abundances. . All of these measurements of the total EM are based on spectra from the same instrument (EPIC-pn) extracted from the same region (all of DEM L71). The differences found here illustrate the inherent uncertainties involved when assumptions are made about the spectral state of the emitting gas, especially when only a small fraction of the gas is being considered, something that is largely avoided with SPI.

The majority of the X-ray emission is from the low temperature gas in the CSM, but there is clearly some contribution from gas with temperatures above 2 keV, as can be seen in the temperature distribution (Figure 7). Our maps show this hot gas is mostly located in the core (Figures 8, 10, and 11). However, simple fits miss this low emission measure gas. This is a byproduct of only being able to measure discrete values that are similar to a weighted average over a region. While the low EM gas may influence these measured averages, resulting in the higher reported temperatures of the ejecta for example, the properties of this gas cannot generally be determined by standard spectral fitting methods, since it is very difficult to separate the emission from the more dominant (high EM) gas. Temperature is the clearest example, but the same effect applies to all parameters. Extremely deep observations can help, but our results clearly demonstrate that the physical state of the gas is very complex. Standard spectral analyses require a discrete number of spectral components, but observations which are of high enough quality to explore the full complexity of the gas cannot be adequately fit with a manageable number of these components.

While the XMM-Newton analyses use large spatial regions, resulting in averaging over any variations, the Chandra studies take the opposite approach, with a few small regions chosen for spectral analysis. This has the advantage of a simpler spectrum, as a smaller region will be more physically homogeneous and is therefore easier to fit, although the spectrum is still integrated along the line of sight.

H03 provide the only other measurement of Ne in the ejecta, and find a very high abundance (Figure 14) compared to the CSM: 1.7 and 2.3 times solar for the inner and outer ejecta, respectively. We do not find this to be the case, instead finding LMC-like values that are similar to those in the CSM. The Ne distribution (Figure 12) does indicate a very small amount of gas with Ne $>$ 1 in the ejecta, and Ne is highest in the region used by H03 in the southeast (Figure 9), but in neither case do we find a substantial amount of Ne gas with such a high abundance as indicated by H03. However, note that the uncertainties on the reported high Ne are large, likely a consequence of choosing a small region; there are relatively few counts in the spectrum.

In contrast to Ne, we find higher Fe in the CSM than either H03 or R03. This can be explained by the choice of the spectral extraction regions. As with all parameters, we find that Fe varies irregularly across DEM L71 (Figure 9). In choosing only a few small regions, H03 and R03 are sampling different parts of this spatial distribution, and thus measure Fe values for the CSM that are different from both each other and this work. The region used by H03, in the southeast, is among the lowest Fe regions in DEM L71, and our Fe distribution (Figure 7) does peak at very low Fe (the asymmetric shape is what causes the median to be higher). R03, on the other hand, combines several regions from around the bright CSM shell and correspondingly measures a higher value for Fe, similar to the median Fe we report here. While this better matches our median Fe, it does not encompass the low Fe gas indicated by both our results and that of H03. These differences in abundance measurements highlight the selection effects resulting from the practice of choosing discrete regions. Because physical conditions vary widely and irregularly across the entire source, the results are highly dependent on choice of region and are unlikely to be representative of the object as a whole.

6 Conclusions

SPI analysis of DEM L71 reveals a more comprehensive and complicated picture than previously known, with irregular distributions of temperature, ionization age, and other properties. The cool outer shell of shocked CSM has a median temperature of $0.24$ keV but includes gas with temperatures that range from $\sim 0.2$ to $\sim 0.45$ keV. Ionization age of the CSM material varies widely around the outer rim but is high on average, indicating much of the CSM is in ionizational equilibrium. The CSM exhibits subsolar LMC-like abundances. The inner core is warmer with enhanced Fe (greater than solar), which supports a Type Ia origin for DEM L71. The core also reveals a somewhat higher, but still subsolar, abundance of O and Ne. We use the enhanced Fe and warmer temperatures to isolate the ejecta from the CSM and examine it in more detail. We find that in addition to typical temperatures of $\sim 1-3$ keV, the ionization age of the core material is much lower than that of the CSM. Both the warmer temperatures and lower ionization age of the ejecta are in agreement with expectations from our simulations. We also independently identify the contact discontinuity that separates the shocked CSM from the ejecta via a distinctive hook-shaped region where temperature variations are at a maximum, matching the R-T unstable region predicted BY our 2D simulations.

Overall, we find median values that generally agree well with what other studies have found; however, these are broad averages that do not reflect the distribution of properties over the whole remnant that SPI is able to realize. Even for a case which looks relatively ‘clean’ such as DEM L71, SPI reveals substantial and irregular variation in all properties across the remnant. We show that this introduces biases for analyses which choose specific regions and assume they are spectrally homogeneous. Small regions cannot be assumed to be representative of the SNR as a whole, and larger regions average over the variations and miss small features.

Relatively simple traditional approaches provide only very limited or no information about the wide range of physical conditions throughout the entire SNR. The most common approach of fitting the spectrum with a multi-component model will not reproduce the true variations (both spatial and spectral) of gas properties, unless an inordinately large number of components are used. Analyzing spectra of small regions covering the entire remnant is an alternative that is sometimes used (e.g. Yang et al., 2008; Lopez et al., 2013); however, this requires extremely deep observations and still requires simple spectral models for each region. SPI does not require exceptionally deep observations, and allows for arbitrarily complex spectral models. It also removes any need to make assumptions about the spatial or spectral distribution of physical properties (e.g. spherical symmetry or number of different temperatures) while providing a way to measure more than simple ‘average’ quantities. It can be used to make maps and distributions of any of the parameters present in the spectral models. SPI can therefore provide much more globally accurate conclusions about physical properties or the evolutionary state of a SNR, as compared to standard methods of X-ray analysis.

The nature of SPI makes it possible, and comparatively easy, to separate and study the properties of low EM gas, or to identify and isolate any subset of gas with given properties. Simple 1D simulations can enhance this capability by providing guidance on which parameters are most informative for a given SNR and by aiding the identification of different gas components, e.g. shocked ejecta, by considering the evolution of the SNR. In the case of DEM L71, we are able to separate emission from the ejecta based not on visual identification of surface brightness features, but rather on measured temperature and Fe abundances. This allows us to map and construct distributions of the ejecta temperature, ionization age, and other properties. We find, as expected from the 1D simulations, that this gas has lower ionization ages than the bulk of the SNR gas, as well as flatter O, Ne, and Fe distributions. For other SNRs, a similar procedure can be used to isolate and study components of interest.

SPI also offers a unique method for identifying the CD without relying on surface brightness features. 2D simulations clearly predict a shell-like region of R-T instabilities at the CD and provide an estimate of its location and width. The largest variations in temperature are expected in this region. With SPI, it is possible to visualize such variations by mapping the standard deviation of the temperature distribution within each spatial bin. We find the largest temperature variations in a distinct hook-shaped feature that matches the expected size and location of the R-T unstable CD region. This demonstrates the power of SPI and its potential for extracting new physical insights from existing X-ray observations of SNRs.

Standard approaches are very valuable because the methods and uncertainties are generally very well understood, which is not yet the case for SPI. They are also very fast compared to SPI, which is resource-intensive and slow. However, SPI can provide unique and valuable insight, and is a good complement to traditional methods of X-ray analysis of SNRs. In future papers we will apply the SPI technique to several more SNRs for which observations with XMM-Newton exist in the database.

This work was supported by NASA ADAP grant NNX15AH70G to Pennsylvania State University, with subcontracts to the University of Chicago and Northwestern University. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.