Proper Motion of the High-velocity Pulsar in SNR MSH 15-56

Composite supernova remnants (SNRs) are characterized by an SNR shell produced by a supernova (SN) blast wave and a pulsar wind nebula (PWN) generated by a rapidly rotating, highly magnetized neutron star. The evolution of these two components is closely coupled. The PWN initially drives a shock into the innermost SN ejecta and eventually encounters the SN reverse shock that propagates toward the explosion center (e.g., Blondin et al. 2001). The signatures of the interaction that follows provide insight into the evolution of the system, as well as constraints on the SN progenitor, density structure of the ambient ISM, and the evolution of the spectrum of the relativistic particle population injected by the pulsar that ultimately escapes into the ISM.

Composite SNRs in which the reverse shock has already begun interacting with the PWN exhibit complex morphologies consisting of irregularly shaped PWNe that are displaced from the geometric centers of their SNRs (e.g., Temim et al. 2009; Slane et al. 2012). Recent multi-wavelength studies and hydrodynamical (HD) modeling of these evolved systems have significantly increased our understanding of how these morphologies form (Temim et al. 2015; Kolb et al. 2017). The structure of the composite SNR G327.1-1.1 was successfully reproduced by an HD simulation, invoking a density gradient in the surrounding ISM and a rapidly moving pulsar that effectively cause the reverse shock to interact with the PWN asymmetrically. Kolb et al. (2017) used the same HD model to explore what parameters govern the evolution of composite SNRs that contain PWNe generated by pulsars with high kick velocities. They found that the structure of these SNRs is most strongly influenced by the structure of the ambient medium and the total spin-down energy of the pulsar. Secondary factors include the pulsar's velocity, direction of motion, and its initial spin-down luminosity. The simulations also show that the passage of the reverse shock leaves a trail of PWN material behind the pulsar and that the direction of this trail reflects the flow of the SN ejecta at that location. Since a number of factors drive the morphology of composite SNRs, observational constraints on some of these input parameters are necessary to produce a more realistic model and learn about other SNR properties. In the case of SNR G327.1-1.1 (Temim et al. 2015), the radio and X-ray observations provided constraints on the pulsar's spin-down luminosity and the orientation of the ISM density gradient, which in turn allowed us to infer the direction of the pulsar's motion from the HD simulation. Here, we apply the HD model to another evolved composite SNR, MSH 15-56.

The radio and X-ray morphologies of MSH 15-56 are shown in Figure 1. The 843 MHz radio image from the Molonglo Observatory Synthesis Telescope (MOST; Whiteoak & Green 1996) is shown in the left panel of the figure. It consists of a filamentary, circular SNR shell, approximately 18' in radius. The PWN is outlined by white contours and is offset to the southwest of the geometric center of the SNR shell. Its radio luminosity of ${L}_{{10}^{7}-{10}^{11}\mathrm{Hz}}=5\times {10}^{34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ is the third highest PWN luminosity after the Crab Nebula and SNR G328.4+0.2 (Dickel et al. 2000), assuming a distance of 4.1 kpc obtained from Hα velocity measurements (Rosado et al. 1996). The presumed pulsar is located at the tip of the radio PWN, close to the southwestern edge of the SNR shell. Its position is marked by the yellow cross in the left panel of Figure 1. The pulsations from the source have not been detected, but the presence of the pulsar is inferred from the surrounding PWN. The location and structure of the displaced PWN itself suggest that the SNR reverse shock has already interacted with the nebula. Based on the spectral properties of the thermal X-ray emission in the SNR shell attributed to swept-up interstellar medium (ISM) material, Temim et al. (2013) estimated an age of 16,500 years, consistent with the SNR being an evolved system. In this paper, we measure the proper motion of the presumed pulsar in MSH 15-56 and use it to constrain an HD model of this system in order to understand the SNR's evolution, ambient environment, and the progenitor star.

Zoom In Zoom Out Reset image size

Our analysis includes Chandra X-ray observations of MSH 15-56 carried out at two different epochs. The first observation was taken on 2001 August 18 for a total exposure time of 56.8 ks (observations ID 1965, PI: Plucinsky). The second set of five observations was carried out in 2016 May 30–June 4, for a total exposure time of 175 ks (ObsIDs 17897, 18862–18865, PI: Temim). All observations were acquired using the Advanced CCD imaging Spectrometer (ACIS-S) in the very faint (VFAINT) mode. We performed the standard cleaning and data reduction on both sets of observations using CIAO version 4.9 and calibration database version CALDB 4.7.3.

In order to search for the pulsar proper motion, we compared the 2001 Chandra image with the longest of the five observations from 2016 (ObsID 17897), carried out on May 31 for 78 ks. This provided a baseline of 14.78 years. The other four observations from 2016 were not used in the analysis, since they all had exposure times shorter than the 2001 observation. Since the detection sensitivity for reference sources to be used in aligning the two images is set by the shorter 2001 observation, it was not necessary to merge the 2016 images and use all 175 ks, especially since merging images taken at slightly different roll angles may introduce additional astrometric uncertainties.

While the pointing accuracy of Chandra is ∼04, the relative astrometry between images can be improved by cross-matching sources in the two observations. The CIAO tool reproject_aspect, a script that runs the tools wcs_match and wcs_update, can be used to calculate astrometric shifts between two source lists and then apply the offsets to the observations. We first used the CIAO command wavdetect to search for point sources in both images. Eighteen sources were identified in the 2016 image. We removed the pulsar from this source list, as well as six additional sources that we identified in VizieR to have proper motions of more than 6 mas yr⁻¹. The remaining 11 sources were cross-matched with the source list from the 2001 observation using reproject_aspect. The script deleted bad matches that primarily consisted of far off-axis and blended sources, resulting in four matches that were used to calculate the correction to the aspect solution. These sources are shown in Figure 1. We then ran wavdetect on the corrected images and calculated the difference in R.A. and decl. for the four cross-matched sources and the pulsar. The final coordinates of the sources and their residuals are listed in Table 1. The shifts in R.A. and decl. measured for the pulsar are 0134 and 0532, respectively. The average absolute values of ΔR.A. and Δdecl. for the cross-matched sources are 0182 and 0134, while the average uncertainty on the source positions measured by wavdetect is 0115 in R.A. and 0077 in decl. These two sets of uncertainties were combined to derive the final uncertainty estimate in the residuals measured for the pulsar. Assuming a distance of 4.1 kpc, this translates to a pulsar velocity of ${720}_{-215}^{+290}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and a direction of motion of 14° ± 22° west of south (yellow arrow in the right panel of Figure 1). This velocity is on the high end of the distribution of kick velocities estimated for other neutron stars inside SNRs (e.g., Holland-Ashford et al. 2017).

Figure 1 (right) shows the 2016, 175 ks Chandra X-ray image of the region surrounding the pulsar, with the 0.3–1.5 keV emission in red and 1.5–9.0 keV emission in blue. The pulsar is surrounded by an extended, but compact X-ray PWN. Temim et al. (2013) analyzed the XMM-Newton and 2001 Chandra data of the PWN region, and found the total X-ray luminosity of the nonthermal emission to be L_X (0.3–10 keV) = 4 × 10³⁴ erg s⁻¹ and that the compact PWN has a photon index of ${1.84}_{-0.11}^{+0.13}$. The soft X-ray emission to the west that is outside of the radio PWN contours has a significant thermal component that may arise from SN ejecta that have been encountered by the reverse shock. Thermal emission with enhanced abundances was also present throughout the region surrounding the PWN. An arc of hard X-ray emission traces the southern edge of the radio PWN. While the spectrum of the arc showed that both thermal emission with enhanced abundances and nonthermal emission are present in the region (Temim et al. 2013), the hard nonthermal component is more prominent along the arc than in the western region that appears red in Figure 1(right). The nature of this arc is still unknown, but one possibility is that its emission originates from a pulsar's jet that has been swept back by the SN reverse shock. In this case, the thermal ejecta emission arises from projection along the line of sight.

We modeled the evolution of MSH 15-56 using the VH-1 hydrodynamics code (Blondin et al. 2001) that has been modified to simulate a composite SNR evolving into a non-uniform ISM and hosting a PWN generated by a moving pulsar. The modified code is discussed in detail in Kolb et al. (2017). The PWN is modeled as an expanding bubble of relativistic gas (adiabatic index γ = 4/3) driven by the pulsar's time-dependent spin-down luminosity $\dot{E}\,=\dot{{E}_{0}}{(1+t/{\tau }_{0})}^{-(b+1)/(b-1)}$, where $\dot{{E}_{0}}$ is the pulsar's initial spin-down luminosity, τ₀ is the spin-down timescale, and b is the braking index. The braking index and the magnetic field are both assumed to be constant (see Blondin et al. 2001). The SNR shell is modeled as a self-similar wave that expands into an SN ejecta density profile that falls off exponentially with a power-law index of n = 12 (e.g., Chevalier 1982), a value that is typical for a red supergiant progenitor (Matzner & McKee 1999). The surrounding ISM density profile includes a gradient of the form $n(z)={n}_{0}[1-(2-x)\tanh (z/H)/x]$, where z is the distance from the center for the explosion, n₀ is the average (central) density, H is the characteristic length scale over which the density varies, and x is the parameter that defines the density contrast (Blondin et al. 2001). The simulation is designed in the rest frame of the pulsar and the pulsar's motion through the SNR is simulated by giving a velocity to the SNR and the surrounding ISM.

5.1. Model Parameters

The HD model input parameters are summarized in Table 2. The primary observational constraint included in the model is the pulsar's speed, measured in Section 3 to be 720 km s⁻¹. The SN explosion energy was set to 10⁵¹ erg, and the other input parameters varied such that the morphology of the PWN and the SNR resembles the observed morphology at the time when the SNR radius reaches a size of ∼21 pc (the observed radius of 18' at a distance of 4.1 kpc). The pulsar's initial spin-down luminosity ${\dot{E}}_{0}$ and timescale τ₀ were chosen to obtain a current spin-down luminosity $\dot{E}$ that is consistent with the total X-ray luminosity of the PWN measured in Temim et al. (2013). In this case, the ratio ${L}_{X}/\dot{E}=0.00235$ is within the range of luminosities observed for other pulsars (Possenti et al. 2002). The mass of the SN ejecta was set to 10 M_⊙, the low end value predicted by the models of Sukhbold et al. (2016). The minimum ISM density was set to n₀ = 0.1 cm⁻³ (Temim et al. 2013), increasing by a factor of eight in a direction perpendicular to the direction of the pulsar's motion.

Table 2.  HD Model Parameters

Parameter	Description	Value
Input Parameters
*vp (km s−1)	Pulsar velocity	720
*vp direction	West of south	10°
E51 (1051 erg)	Explosion energy	1.0
Mej (M⊙)	SN ejecta mass	10
n	SN ejecta density profile index	12
$\dot{{E}_{0}}$ (erg s−1)	Initial spin-down luminosity	7.0 × 1038
b	Pulsar braking index	3.0
τ0 (year)	Spin-down timescale	2000
*n0 (cm−3)	Minimum ambient density	0.1
x	n0 contrast of 8.1	1.125
H (pc)	n0 characteristic length scale	3.24
n0 gradient angle	North of west	10°
Simulation Output
*RSNR (pc)	SNR radius	21
*${L}_{X(2-10)}$ (erg s−1)	PWN X-ray luminosity	4.0 × 1034
*$\dot{E}$ (erg s−1)	Current spin-down luminosity	1.7 × 1037
t (year)	SNR age	11,000

Note. The pulsar velocity v_p is calculated from the pulsar proper motion measured in this work. The ejecta mass M_ej, braking index of the pulsar, initial spin-down luminosity, the spin-down timescale, and the SNR age were varied until the simulation produced the observed PWN morphology and an SNR radius of 21 pc, for an assumed distance of 4.1 kpc. The input and output properties marked by an asterisk are those that have been constrained by observations. The equation for the density gradient of the medium into which the SN expands is given in Blondin et al. (2001).

Download table as:  ASCIITypeset image

5.2. Results

The simulation produced from the input parameters described in Section 5.1 is shown in Figure 2. The left panel is a density map at an SNR age of 11,000 years, which is when the SNR reaches a radius of 21 pc and the SNR morphology most closely resembles the observed radio emission from the PWN shown in the left panel of Figure 1. An increase in the ambient density would increase the SNR age, but would also cause the reverse shock to interact with the PWN sooner, resulting in a fully disrupted nebula at a time when the SNR shell reaches the observed radius. The early reverse shock interaction can be mitigated by increasing the SN ejecta mass. A combination of higher ambient density and higher SN ejecta mass can therefore produce a similar morphology as seen in Figure 2 at an age greater than 11,000 years.

Zoom In Zoom Out Reset image size

Using the SNR age of 11,000 years, a distance of 4.1 kpc, and the measured velocity of the pulsar, we can estimate the location of the SN explosion center, marked by the black cross in Figure 1 (left). The error box for the explosion center, reflecting the uncertainties in the pulsar's velocity and direction of motion, is also shown. As expected, the site of the explosion is displaced to the west of the geometric center of the SNR shell, toward the direction of the higher ambient density. The radio morphology of the SNR shell also shows evidence that the radius of the shell is smaller in the western half of the SNR. The higher ambient density in this direction and the displacement of the PWN to the south that results from the pulsar's motion cause the SN reverse shock to interact with the PWN asymmetrically. This can be seen in Figure 3, which shows four different time snapshots of the SNRs evolution. The PWN becomes displaced from the explosion center due to the pulsar's motion, and the reverse shock begins interacting with the nebula at ∼5000 years. At the age of ∼7500 years, the reverse shock reaches the pulsar, and as it continues toward the northeast, it sweeps the bulk of the PWN material away from the pulsar. A trail of PWN material begins to form behind the pulsar, along the direction of the passage of the reverse shock (9000 year panel of Figure 3). This supports the finding of Kolb et al. (2017) that the PWN trail reflects the velocity flow of the SN ejecta. Due to the lower ambient density to the east, the SN reverse shock has not yet reheated all of the SN ejecta in the northeastern part of the SNR and has not yet reached the PWN from this direction. The unshocked ejecta material can be seen in the light blue color in Figure 2 (left). The density enhancement at the northeastern edge of the PWN represents the innermost SN ejecta that have been swept up by the PWN.

Zoom In Zoom Out Reset image size

The right panel of Figure 2 shows the map of the age of the particles ejected by the pulsar at the final SNR age of 11,000 years. While we do not take into account the effects of any enhanced synchrotron cooling that may be caused by the PWN's interaction with the reverse shock, the map is consistent with the observed radio and X-ray morphologies of the SNR, assuming that the particles cool with age. As can be seen in Figure 1, the X-ray emitting region of the PWN in MSH 15-56 coincides with the youngest particles concentrated in the trail behind the pulsar, shown in blue. The particles with an age greater than 4000 years shown in orange are concentrated in the northeastern region of the PWN that coincides with the PWN's radio emission, as would be expected from synchrotron aging. The simulation also shows filamentary structures of younger particles that extend toward the northeast. These structures may explain the hard X-ray arc seen in the Chandra image in Figure 1 (right), but since the arc traces the edge of the radio PWN and the location of the reverse shock, it is more likely that it arises from a jet that has been swept back by the shock, similar to the structure inferred for the Geminga pulsar (Posselt et al. 2017).

In this paper, we measured the proper motion of the pulsar in SNR MSH 15-56 using Chandra X-ray images taken at two different epochs, resulting in a baseline of ∼15 years. We measure a pulsar velocity of ${720}_{-200}^{+250}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and a direction of motion of 10° ± 14° west of south. We used the measured pulsar velocity to constrain an HD model of the evolution of the SNR that invokes a moving pulsar and a density gradient in the surrounding ISM that cause the reverse shock to interact with the PWN asymmetrically. The model reproduces the morphology of the SNR at an age of 11,000 years, assuming an explosion energy of 10⁵¹ erg, an SN ejecta mass of 10 M_⊙, an average ambient density of 0.4 cm⁻³, and a distance of 4.1 kpc. The measured pulsar velocity places a constraint on the orientation of the ambient density gradient and indicates that the density increases from east to west. An increase in the average ambient density and the SN ejecta mass would produce a similar SNR morphology at an older SNR age. Constraints on the ambient density from upcoming XMM-Newton observations of the entire SNR shell will allow us to use the HD model to constrain the range of SN ejecta masses and ages that could produce the observed morphology of the SNR.

Support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number GO6-17057X issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060.