SPATIALLY RESOLVED SPECTROSCOPY OF A PULSAR WIND NEBULA IN MSH 15−56

The southwest vicinity of MSH 15−56 was observed with XMM-Newton on 2004 August 11 (Observation ID: 0204270101). The aim point of the telescope was $\alpha = 15^{\mathrm{h}} 52^{\mathrm{m}} 14\buildrel{\mathrm{s}}\over{.}99$, δ = −56°18'120 (J2000.0), about 10' southwest of the center of the SNR. EPIC-MOS1/2 and PN operated in the prime full window mode with the medium filter. The data mode was set to "imaging."

Although the duration of the whole observation was 57.6 ks, half of the total exposure time was seriously damaged by background flares. Therefore, we created a Good Time Interval (GTI) file that can be used to create a filtered event file that only accepts good times where the background was low. To distinguish usable exposures from noisy time slots, we defined selection criteria based on the event rate. For this purpose, we made histograms of the event rate for 100 s intervals in the 10–12 keV energy band and fitted the histogram with a function consisting of a Gaussian and a power-law function. While the background flares show high variability, the good times are dominated by the steady cosmic X-ray background (CXB) event and correspond to the time intervals in the Gaussian regime on the histogram.

The best-fit parameters of the Gaussian function were μ = 0.075 ± 0.03 and σ = 0.03 ± 0.01 counts s⁻¹, respectively, for both MOS1 and MOS2, which corresponds to the X-ray events from the target object and CXB. We also obtained the best-fit Gaussian function for PN: μ = 0.483 ± 0.009 and σ = 0.076 ± 0.009 counts s⁻¹.

Consequently, we defined the upper threshold of the event rate as 2.706σ, which is equivalent to 0.15 and 0.69 counts s⁻¹ for MOS and PN, respectively. The data sets were then filtered with the GTI file based on the criteria defined above using Science Analysis System (SAS) version 10.0. Finally, we obtained a usable data set with a 28.5 ks exposure.

The field of MSH 15−56 was also observed with the Chandra X-Ray Observatory (Observation ID: 1965) on 2001 August 18 using the back-illuminated CCD S3 of the ACIS (Burke et al. 1997). The detector was operated in the VFAINT, timed-exposure mode. We therefore applied a standard data reduction process to generate a level 2 event file with VFAINT-mode background cleaning using CIAO version 4.1 with CALDB version 4.1.3. In the process, we disabled pixel randomization to resolve fine structures just around the pulsar. The resulting usable exposure time was 56.1 ks.

2.2.1. Large-scale Structures (XMM)

To perform a morphological analysis and spatially resolved spectroscopy, we combined X-ray images of EPIC MOS and PN utilizing "Image_script," which is an analyzing script of SAS. The resulting X-ray images of the southwest rim of MSH 15−56 in four energy bands (0.2–1.2, 1.2–2.0, 2.0–10.0, and 0.2–10.0 keV) are shown in Figure 1. To clarify the physical relationship between the SNR and PWN, Figure 2 compares the X-ray and radio structures in MSH 15−56.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the 0.2–1.2 keV energy band, only a diffuse structure was observed. The soft component seems to fill the inside of the SNR and shows a sharp edge at its southwest rim, while the harder energy bands revealed a compact nebula with a coma and a diffuse tail. Those morphological features distinguish the object from other background point sources. The hard X-ray nebula extends about 2' in the direction of the SNR center (see Figure 2).

2.2.2. Fine Structures (Chandra)

Chandra reveals more fine structures. Figure 3 shows a close-up view of the PWN in the 0.4–8.0 keV energy band with a binning size of 4 × 4 pixels (∼2'' × 2''). Although the PWN apex appears to be a point source, we tried to search for the extended component correlating with the wind termination shock, as seen in the Mouse nebula (Gaensler et al. 2004). To constrain the size of the emitting region, we analyzed the Chandra image using Sherpa in CIAO (Freeman et al. 2001). Since the size of the brightest part of the PWN appears to be comparable with the point-spread function (PSF) of the HRMA, we evaluated the intrinsic size of the nebula. For this purpose, we made a PSF of the HRMA as a convolution kernel with the ray trace simulator ChaRT (Carter et al. 2003) for the target position and the observed spectrum.

Zoom In Zoom Out Reset image size

The close-up X-ray view of the apex compared with the HRMA's PSF is shown in Figure 4, and Table 1 shows the results of model fitting with a two-dimensional Gaussian function indicating that the angular size of the coma exceeds the PSF to some extent. The angular size of the coma is 17 ± 04 in FWHM. The model fitting also constrains the precise coordinate of the point source,

$$\begin{eqnarray*} &&\alpha = 15^{\rm h} 52^{\rm m} 12 \buildrel{\mathrm{s}}\over{.}672, \quad \quad \delta = -56^\circ 18^\prime 58{\buildrel{\prime\prime}\over{.}} 26, \end{eqnarray*}$$

though that still contains a systematic error of about ±06 (with a 90% confidence level).⁴

Zoom In Zoom Out Reset image size

Table 1. Results of Two-dimensional Model Fitting of the Chandra Image

Parameter	Best Fit	Lower Bound	Upper Bound
Gaussian component
X-position (J2000)	$15 ^{\rm h} 52 ^{\rm m} 12 \buildrel{\mathrm{s}}\over{.}672$	$-0 \buildrel{\mathrm{s}}\over{.}013$	$+0 \buildrel{\mathrm{s}}\over{.}012$
Y-position (J2000)	−56°18'5826	−018	+019
FWHM (arcsec)	171	−040	+042
Amplitude (counts)	12.5	−4.5	+8.9
Constant component
Amplitude (counts)	0.068	−0.002	+0.002

Note. Errors are within a 90% confidence level.

Download table as:  ASCIITypeset image

The tail was resolved into a wedge shape that is brightening at the edges. The west side of the wedge is clear and shows a mildly curved line extending ∼1' from the point source, while the opposite side is relatively unclear due to the superimposed arc-like structure on the south edge of the SNR (see below). The apparent half opening angle of the wedge is θ_app ∼ 25°.

Near the PWN apex, an arc-like faint filament extends eastward, labeled "ARC" in Figure 3. This feature is also observed with XMM-Newton but only in the hard X-ray band. Figure 5 describes the surface brightness distribution along the line normal to the arc. While the soft X-ray component appears to fill the SNR as seen in Figure 2, the hard X-ray component was only detected at the edge of the SNR. The hard X-ray filament has a width of about 20'', which may be due to the three-dimensional geometry projection of the SNR. Nevertheless, the sharply rising outer edge implies the existence of a very thin (∼10'') structure that may be associated with a non-thermal X-ray filament as seen in young SNRs (Bamba et al. 2003).

Zoom In Zoom Out Reset image size

2.3.1. SNR Component

Because the PWN appears in the SNR on the celestial sphere, we first evaluated the SNR component before investigating the PWN in detail. For this purpose, we selected regions "A" (a = 32, b = 12) and "B" (r = 3') denoted in Figure 2. Region A is the brightest part of the southwest rim in the X-ray energy band and contains the PWN. Conversely, we chose region B, about 4' apart from the PWN, to avoid contamination from non-thermal structures such as the PWN and arc. From these regions, source photons were accumulated using "evselect" and screened with the GTI filter mentioned above. For background subtraction, we selected another region labeled "BG" outside the SNR in Figure 2. After the background subtraction, we obtained 4613 events (MOS1)/4753 events (MOS2) and 4273 events (MOS1)/3943 events (MOS2) from regions A and B, respectively. Figure 6 shows the obtained X-ray spectra, both of which show several line features below 2 keV, while region A also shows a hard X-ray continuum above 2 keV that may have originated in the PWN and arc.

Zoom In Zoom Out Reset image size

In the spectral analysis of region B, we tried three thin-thermal plasma models: a non-equilibrium ionization (NEI) model and two equilibrium ionization models of MEKAL and APEC. To account for the line features, the metal abundances of Ne, Mg, Si, and Fe were set to be variable. The resulting best-fit parameters are summarized in Table 2. The three plasma models converged under similar conditions. The obtained high n_et infers that the diffuse plasma filling the SNR is in an ionization equilibrium state.

Table 2. Results of Spectral Fitting of Regions A and B

Model	NH	kT	Metal Abundancea				net	EMb	Photon Index	Flux0.5-10 keVc	$\chi _{\nu }^2(\nu)$
	(1022 cm−2)	(keV)	Ne	Mg	Si	Fe	(cm−3 s)			(10−12erg s−1cm−2)
Region A
NEI+PL	0.35 ± 0.03	0.61 ± 0.03	$0.58_{ -0.19}^{+ 0.25}$	$0.59_{ -0.18}^{+ 0.12}$	$1.02_{ -0.24}^{+ 0.27}$	$0.38_{ -0.06}^{+ 0.05}$	>4.4 × 1011	$1.00_{ -0.47}^{+ 0.07}$	$1.87_{ -0.12}^{+ 0.13}$	1.46 ± 0.08	1.115 (139)
MEKAL+PL	0.35 ± 0.05	0.61 ± 0.03	$0.63_{ -0.32}^{+ 0.58}$	$0.94_{ -0.34}^{+ 0.64}$	$1.09_{ -0.36}^{+ 0.67}$	$0.34_{ -0.09}^{+ 0.18}$	...	$1.06_{ -0.45}^{+ 0.53}$	$1.79_{ -0.22}^{+ 0.19}$	$1.40_{ -0.12}^{+ 0.13}$	1.167 (140)
APEC+PL	$0.36_{ -0.04}^{+ 0.08}$	0.62 ± 0.03	$0.71_{ -0.40}^{+ 0.45}$	$0.79_{ -0.37}^{+ 0.43}$	$0.95_{ -0.41}^{+ 0.48}$	$0.36_{ -0.12}^{+ 0.13}$	...	$1.19_{ -0.42}^{+ 0.96}$	$1.74_{ -0.34}^{+ 0.17}$	$1.36_{ -0.16}^{+ 0.11}$	1.128 (140)
Region B
NEI	$0.39_{ -0.06}^{+ 0.04}$	$0.62_{ -0.01}^{+ 0.03}$	$0.92_{ -0.20}^{+ 0.19}$	$0.46_{ -0.11}^{+ 0.12}$	$0.71_{ -0.18}^{+ 0.16}$	$0.32_{ -0.08}^{+ 0.07}$	>3.5 × 1011	$3.04_{ -0.84}^{+ 0.58}$	...	...	0.903 (61)
MEKAL	0.34 ± 0.05	0.65 ± 0.03	$0.87_{ -0.30}^{+ 0.37}$	$0.92_{ -0.27}^{+ 0.36}$	$0.90_{ -0.25}^{+ 0.32}$	$0.33_{ -0.06}^{+ 0.08}$	...	$2.44_{ -0.47}^{+ 0.54}$	...	...	1.034 (61)
APEC	0.36 ± 0.05	0.65 ± 0.03	$0.99_{ -0.30}^{+ 0.41}$	$0.81_{ -0.25}^{+ 0.36}$	$0.86_{ -0.24}^{+ 0.34}$	$0.38_{ -0.08}^{+ 0.11}$	...	$2.52_{ -0.54}^{+ 0.59}$	...	...	0.891 (60)

Notes. Errors are within a 90% confidence level. ^aThe abundance of each metal is normalized by the solar abundance. ^bEmission measure $10^{-17}({1}/{4\pi D_{\rm A}^2}) \int n_{\rm e} n_{\rm H} dV$, where D_A is the angular diameter distance to the source (cm), and n_e and n_H (cm⁻³) are the electron and hydrogen densities, respectively. ^cUnabsorbed energy flux of the power-law component over the energy range of 0.5–10 keV.

Download table as:  ASCIITypeset image

For region A, we employed two-component models consisting of a thermal plasma model and a power-law function that explains the high-energy continuum. Moreover, we found that the parameters of the thermal component were very similar to those of region B: a column density of N_H = 3.5 × 10²¹ cm⁻², a plasma temperature of kT = 0.61 keV, a relatively low Fe abundance below 0.4 in the solar abundance unit, and a high ionization parameter of n_et > 4 × 10¹¹ cm⁻³ s. Conversely, the power-law function has a photon index of ∼1.8, which seems to be a reasonable value for a PWN (Kargaltsev & Pavlov 2008).

2.3.2. Point Source

Next, we tried to investigate the spectrum of the point source at the PWN apex with the Chandra data. To prevent photons from the SNR or PWN contaminating the spectrum of the point source, we selected a circular region centered on the point source with a radius of 3'', in which more than 90% of photons from a point source would be accumulated (Weisskopf et al. 2000b), while the background events were sampled from an annulus just around the source region (3'' < R < 20''). Subsequently, we extracted photons using "specextract" from the selected regions. After the background subtraction, we obtained 143 photons for a 0.5–10.0 keV energy band (the background contribution was 2.2%). The poor photon statistics meant that the column density could not be constrained by a simultaneous fitting. We therefore employed a fixed value of N_H = 3.5 × 10²¹ cm⁻² for the column density, based on the spectroscopy of region B shown above. In the spectral fitting of the point source, we examined the model functions of a power-law function and a blackbody, and obtained a photon index of Γ = 1.3 ± 0.3 and a temperature of kT = 0.8 ± 0.1 keV, respectively. In addition, we also tried to fit the spectrum with a mixture model of a blackbody and a power-law function, but could not attain reliable results for the limited photon statistics. The obtained X-ray spectrum is shown in Figure 7 with model functions, and the resultant parameters are summarized in Table 3. The results of the spectral analysis are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Results of Spectral Fitting of the Point Source

Model	NH	kTBB	Photon Index	Normalization	$\chi _{\nu }^2(\nu)$
	(1022 cm−2)	(keV)
Power law	0.35 (fixed)	...	$1.33^{+0.26}_{-0.25}$	$4.12^{+0.75}_{-0.70}$ a	1.036 (12)
Blackbody	0.35 (fixed)	$0.83^{+0.12}_{-0.10}$	...	$2.95^{+0.86}_{-0.58}\times 10^{-2}$ b	1.152 (12)

Notes. Errors are within a 90% confidence level. ^aEnergy flux over the energy 0.5–10.0 keV in units of 10⁻¹⁴ erg s⁻¹ cm⁻². ^bSource radius in kilometers assuming the distance to the source to be 4.1 kpc.

Download table as:  ASCIITypeset image

2.3.3. Spatially Resolved Spectroscopy of the PWN

To investigate the spectral evolution of the pulsar wind along the flow line, we performed a spatially resolved spectroscopy using the XMM data. For this purpose, we divided the PWN into seven rectangular regions with dimensions of 20'' × 100'', as shown in Figure 8. For each region, we then conducted a spectral fitting with a two-component model of an NEI model and a power-law function. The limited photon statistics prevent us from setting all the parameters free. Therefore, we employed fixed parameters: a column density of N_H = 3.5 × 10²¹ cm⁻², metal abundances (Ne = 0.58, Mg = 0.59, Si = 1.0, and Fe = 0.38), and an ionization parameter of n_et = 4.4 × 10¹¹ s cm⁻³, according to the spectroscopy of region A.

Zoom In Zoom Out Reset image size

Figure 9 describes the photon index of the power-law component, the energy fluxes of the power-law component, and the energy flux of the thermal components as functions of distance from the PWN apex. Compared with the thermal component showing constant surface brightness, the power-law component seems to have originated in the PWN, and its photon index obviously steepens with increasing distance from the point source.

Zoom In Zoom Out Reset image size

PWN structures (i.e., a bow shock and a tail) are sensitive to the ambient conditions. Accordingly, we first clarify the environmental condition around the PWN in MSH 15−56. Based on the spectroscopy of region A, the hot gas in the SNR seems almost to be in an ionization equilibrium condition with an ionization parameter of n_et ⩾ 4.4 × 10¹¹ s cm⁻³. Assuming that the emitting region is an oblate ellipsoid, the emitting volume amounts to V = (4/3)π d³(32 × 32 × 12) ∼ 2.6 × 10⁵⁷ cm³. Adopting the emission measure listed in Table 2, the number density of the plasma can be estimated, and we have n_p ∼ 0.28f^−1/2(d/4.1 kpc)^−1/2 cm⁻³, where f is a filling factor. Comparing the lower limit of n_et with this number density, we obtained a lower limit of the ionization timescale, τ_ion > 50f^1/2 kyr, which exceeded the dynamic SNR age τ ∼ 15 kyr to some extent (Kassim et al. 1993). If that is due to the overestimation of the emitting volume, then the filling factor f should be smaller than 0.1.

Evolutions of SNRs depend on the number density of the ambient medium. Blondin et al. (1998) calculated the transition time from the Sedov phase to the pressure-driven snowplow (PSD) phase:

$$\begin{equation} t_{\rm PSD} = 29 \times E_{51}^{-4/17} n_{\rm 0}^{-9/17} \; kyr, \end{equation} \tag{ 1 }$$

where E₅₁ is the explosion energy in units of 10⁵¹ erg. At that moment, the SNR radius achieves $R_{\rm SNR} = 6.3 \times 10^{19} E_{51}^{5/17} n_0^{-7/17}$ cm, which corresponds to $17^\prime E_{51}^{5/17} n_0^{-7/17} (d/4.1\,{\rm kpc})^{-1}$. Our subsequent analysis for the PWN will show that the pulsar remains inside the SNR, though its velocity exceeds 1000 km s⁻¹, which implies that the age should be less than 21 kyr. To end the Sedov phase within 21 kyr, we need n₀ > 1.8 cm⁻³, which contradicts the large SNR radius observed. Therefore, we conclude that the SNR MSH 15−56 is presumed to be still in the Sedov phase.

The temperature of the hot gas filling the SNR is 0.61 keV, and hence the sonic velocity in the gas is $c_{\rm s} = \sqrt{(2\gamma kT)/(\mu m_{\rm p})} \sim 440\mu ^{-1/2}$ km s⁻¹, where γ is the specific heat ratio (γ = 5/3 for mono-atomic gas), m_p is the proton mass, and μ is the mean molecular weight.

For the point source at the PWN apex, we performed spectroscopy using the Chandra data. However, the poor photon statistics prevented us from distinguishing the radiative processes, blackbody, synchrotron, or a mixture thereof. With a single blackbody model, we obtained a temperature of ∼0.8 keV, which appears unfeasibly high compared with the standard NS cooling model that predicts that a 10 kyr old isolated neutron star could not radiate in an X-ray energy band (Shapiro & Teukolsky 1986). It should also be noted that there are several counterexamples of the standard model, such as the Vela pulsar (τ ∼ 11 kyr) and Geminga (τ ∼ 500 kyr), that emit thermal X-rays. It is therefore difficult to reject the possibility of thermal emission from the pulsar in MSH 15−56 at this moment. Dickel et al. (2000) proposed that the point source should be a magnetar to explain the position gap between the hard X-ray PWN and the radio plerion. The known magnetar candidates show high-temperature blackbody components up to 0.7 keV (Woods & Thompson 2006), but the obtained source radius of ∼0.02 km, which corresponds to $L^{\rm BB}_{\rm PSR} = 1.8 \times 10^{33}$ erg s⁻¹, is uncomfortably small for a magnetar candidate such as the anomalous X-ray pulsars.

In contrast, the power-law model provided parameters. The obtained power-law index, Γ_PSR ∼ 1.33 ± 0.26, is consistent with the previous analysis reported in Plucinsky et al. (2002) and is a typical value for a pulsar accompanied by a PWN. In addition, the pair of photon indices of the point source and the entire PWN match the empirical correlation, Γ_PWN = 0.86 + 0.72Γ_PSR, as reported by Gotthelf (2003). It should be noted that the point source appears slightly extended with a size of 17 in FWHM. Therefore, the spectral data of the point source also contain X-rays from the extended component just around the pulsar, and it is difficult to eliminate this contamination even with Chandra.

The X-ray luminosity of the point source amounts to $L^{PL}_{\rm PSR} = 8.3\times 10^{31}$ erg s⁻¹ cm⁻² for the 0.5–10.0 keV energy band. Assuming that the non-thermal emission from region A is accounted for by the sum of the point source and the PWN, the flux ratio between the point source and the PWN is about 3%. Kargaltsev et al. (2007) claimed that the luminosities of PWNe and pulsars showed a power-law correlation. Particularly for middle-aged pulsars with non-thermal luminosities within the range 10³¹ erg s⁻¹ < L_PSR < 10³⁴ erg s⁻¹, the relation can be effectively described by a power-law function, log L_PWN = 0.89log L_PSR + 4.2 (or simply L_PWN = 5L_PSR). Compared with those samples in Kargaltsev et al. (2007), the point source in MSH 15−56 is about one order of magnitude fainter than the expected luminosity for its PWN. It should be noted that the luminosities of some PWNe listed in Kargaltsev et al. (2007) can be underestimated because their large-scale components are too faint to be seen. If that is the case, then the luminosity ratio L_PWN/L_PSR may be more reasonable.

Nonetheless, the observed L_PWN/L_PSR ratio resembles that of the Vela pulsar, which is also unusually underluminous in X-rays (Pavlov et al. 2001b). Considering previous pulsar studies, the non-thermal model provided the relatively reliable and appropriate parameters expected for a middle-aged Vela-like pulsar.

Despite expectations (Kaspi et al. 1996), no pulse emission from the point source has been detected in any observable energy bands. We therefore try to estimate the spin-down luminosity $\dot{E}$ using an empirical relation between the spin-down luminosity and the X-ray luminosity of the entire nebula (Seward & Wang 1988; Becker & Truemper 1997; Kawai et al. 1998; Kargaltsev et al. 2007). In this study, we employed one of the most recent model functions for pulsar–PWN systems based on Chandra and XMM data, $L_{\rm PWN}/(10^{32}\: {\rm erg} {\rm \: s}^{-1}) = 3.0 \times (\dot{E}/10^{36}\: {\rm erg}\: {\rm s}^{-1})^{1.1}$ for 0.5–10.0 keV (Kanai 2010), and found a spin-down luminosity of $\dot{E} = 7.9 \times 10^{36}$ erg s⁻¹. The obtained $\dot{E}$ is almost equivalent to that of the Vela pulsar, $\dot{E}_{\rm Vela}=6.7\times 10^{36}$ erg s⁻¹ (Pavlov et al. 2001a and references therein). Moreover, Gotthelf (2003) reported another relation between the spin-down luminosity and the spectral index of PWN, $\Gamma _{\rm PWN} = 2.36 - 0.021\dot{E}_{40}^{-1/2}$, where $\dot{E}_{40}$ is the spin-down energy in units of 10⁴⁰ erg s⁻¹. Substituting Γ_PWN = 1.7 into the formula yields $\dot{E} \sim 1 \times 10^{37}$ erg s⁻¹. However, the recent studies on PWNe accompanied by pulsars with known $\dot{E}$ revealed that the relation of $L_{\rm PWN}\hbox{--}\dot{E}$ shows very large scatter up to four orders of magnitude with respect to the best-fit relations (Cheng et al. 2004; Li et al. 2008; Kargaltsev & Pavlov 2010). Moreover, Arzoumanian et al. (2011) discovered a PWN around PSR J2022+3842 that has an unusually low conversion efficiency of spin-down power to X-ray luminosity and concluded that the $L_{\rm PWN}\hbox{--}\dot{E}$ relation proposed earlier may not apply for many young pulsars. Therefore, discovery of the pulse emission is required for more quantitative discussions.

Chandra and XMM-Newton revealed detailed structures of the PWN in MSH 15−56: a point source accompanied by a wedge-shaped broad tail extending ∼2' from the apex. The tail is extended and directed at the center of the SNR. Obviously, the observed PWN is analogous to the bow shock PWN running supersonically in the ISM, as reported in dozens of previous studies.

Currently, we know the precise position of the point source that can be the pulsar driving the PWN. The pulsar is ∼13' apart from the SNR centroid and the offset angle corresponds to 2.4 × 10¹⁹ cm for d = 4.1 kpc. If the pulsar was kicked by the SN explosion that generated MSH 15−56, then the proper motion velocity of the pulsar would be easily calculated by dividing the travel distance by the SNR age. Incidentally, we must take into account the inclination angle i between the direction of proper motion and the line of sight as shown in Figure 10, because the observed image is only a projection of three-dimensional structures on the celestial plane. Subsequently, we have a formula for the proper motion velocity⁵ of the pulsar,

$$\begin{eqnarray} v_{\rm PSR} &=& \frac{13^\prime d}{\tau _{\rm SNR}}\frac{1}{\sin {i}} \nonumber \\ &=& 1.5\times 10^3 \: \frac{1}{\sin {i}} \left(\frac{\tau _{\rm SNR}}{10\: {\rm kyr}} \right)^{-1} \left(\frac{d}{4.1\: {\rm kpc}} \right) \; km s^{-1}. \nonumber \\ \end{eqnarray} \tag{ 2 }$$

Since the above relation depends on the parameters of distance d, the SNR age τ_SNR, and the inclination angle i, which may contain large errors, the resulting velocity seems rather uncertain.

Zoom In Zoom Out Reset image size

We can also approach the proper motion of the pulsar based on the observed structures of the PWN. While simple analyses (Chen et al. 1996; Wilkin 1996) indicate that the bow shock structure should be collimated toward the tail, the observed tail seems to have a constant opening angle θ_app. Accordingly, we interpret this structure as a Mach cone; plasmas left behind from the PWN subsequently expand at the speed of sound, whereupon the surface of such plasma forms a conical structure as discussed in Chevalier & Reynolds (2011). If the surface of the PWN propagates in a medium with sonic speed c_s, then the opening angle of the PWN θ should follow

$$\begin{equation} \sin \theta = \frac{c_{\rm s}}{v_{\rm rel}}, \end{equation} \tag{ 3 }$$

where v_rel is the proper motion velocity with respect to the surrounding material. Here, we have two options for c_s: one is c_{s, ISM} ∼ 15 km s⁻¹ for the PWN already outside the SNR, and the other is c_{s, SNR} ∼ 440 km s⁻¹ for the PWN remaining in the SNR.

In the former case, we can set v_rel ≃ v_PSR. The inclination angle i should exceed the apparent opening angle θ_app ∼ 25°, because we can see the wedge-shaped bow shock (see Figure 3), and it also be smaller than ≃ 44° to have the PWN outside the SNR, as shown in Figure 10(b). Taking into account the effect of the inclination angle i, Equation (3) can be rewritten as

$$\begin{equation} v_{\rm PSR} = c_{\rm s, ISM} \sqrt{1 + 1/(\sin {i}\: \tan {\theta _{\rm app}})^2}. \end{equation} \tag{ 4 }$$

The constraint of the inclination angle, i > 25°, requires a proper motion of v_PSR < 80 km s⁻¹. If this is the case, then the travel time of the pulsar from the birthplace to the current location will be unfeasibly long, ∼260 kyr, which is one order of magnitude longer than the SNR age reported in past studies. Therefore, the PWN is likely still inside the SNR and the inclination angle should be within the range 44° < i < 90°, as shown in Figure 10(a). Even if we neglect the motion of gas inside the SNR and use Equation (4), i = 90° gives a lower limit, v_PSR > 1000 km s⁻¹. Taking into account the flow velocity of the downstream behind the SNR forward shock, the relative velocity is set at v_rel = v_PSR − v_HG, where v_HG is the velocity of the shock-heated gas behind the SNR forward shock. If the SNR is currently in the Sedov phase, meaning that the expansion velocity follows v_SNR = (2/5)R_SNR/τ_SNR, then we determine the velocity of the shock-heated gas behind the SNR forward shock to be v_HG = (3/4)v_SNR by assuming the strong shock condition (van der Swaluw et al. 2001). Consequently, the proper motion velocity can be written as

$$\begin{equation} v_{\rm PSR} = c_{\rm s,SNR} \sqrt{1 + 1/(\sin {i}\:\tan {\theta _{\rm app}})^2} + \frac{3}{10}\frac{R_{\rm SNR}}{\tau _{\rm SNR}}. \end{equation} \tag{ 5 }$$

Note that the above equation also depends on uncertain parameters of τ_SNR and i as in Equation (2). By combining Equations (2) and (5), we can separate the variables v_PSR and τ_SNR as follows:

$$\begin{eqnarray} v_{\rm PSR} &=& c_{\rm s,SNR} \left(1 - \frac{3}{10} \frac{\sin i}{\zeta }\right)^{-1} \sqrt{1 + 1/(\sin i \tan \theta _{\rm app})^2}, \nonumber \\ \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \tau _{\rm SNR} &=& \frac{R_{\rm SNR}}{c_{\rm s, SNR}} \left(\frac{\zeta }{\sin i} - \frac{3}{10}\right) \frac{1}{\sqrt{1 + 1/(\sin i \tan \theta _{\rm app})^2}},\nonumber \\ \end{eqnarray} \tag{ 7 }$$

where ζ ≡ 13'/18' is the distance from the SNR centroid to the pulsar in units of the SNR radius. Figure 11 shows the resulting pulsar's proper motion velocity and the age of the SNR as functions of the inclination angle i for various θ_app. The proper motion of the pulsar is within the range v_PSR = 1900 ± 100 km s⁻¹, which is about five times faster than the typical kick velocity of known pulsars (Frail et al. 1994). Meanwhile, the SNR age should be 9800  ±   ∼ 1200 yr for θ_app ∼ 25°. By substituting the current status of the SNR, R_SNR = 18' and τ_SNR = 9.8 kyr, into the Sedov solution, the density of the ISM surrounding the SNR can be constrained as

$$\begin{eqnarray} &&n_0 = \left(\frac{1.15}{R_{\rm SNR}}\right)^5 \frac{E_{51}}{m_{\rm p}}\: \tau _{\rm SNR}^2 \simeq 9.1 \times 10^{-2} \left(\frac{d}{4.1\:{\rm kpc}} \right)^{-3} \; cm^{-3}. \nonumber \\ \end{eqnarray} \tag{ 8 }$$

The number density of the surrounding material and that of the compressed ISM in the SNR (n_p = 0.28 cm⁻³) seem to follow the jump condition of the strong shock. Although the distance to the object remains unclear, the SNR age obtained seems reasonable when compared with past estimations. Moreover, although the obtained pulsar velocity is relatively fast, it remains in an acceptable range according to past pulsar studies.

Zoom In Zoom Out Reset image size

Assuming an isotropic pulsar wind in the pulsar rest frame, the ram pressure should balance the pulsar wind,

$$\begin{equation} \frac{ \dot{E} }{4 \pi r_{\rm TS}^{\rm F\: 2} c} = \rho v_{\rm rel}^2, \end{equation} \tag{ 9 }$$

where $r_{\rm TS}^{\rm F}$ is the radius of the forward termination shock from the pulsar, ρ is the mass density of the ambient matter, and v_rel is the pulsar's velocity with respect to the surrounding material. As shown in Figure 12, the standoff distance of the backward termination shock $r_{\rm TS}^{\rm B}$ is thought to be elongated because the bow shock becomes rather weak on this side, the ratio of the standoff distances between the forward and backward termination shocks follows

$$\begin{equation} \frac{r_{\rm TS}^{\rm B}}{r_{\rm TS}^{\rm F}} \approx \left(\frac{p_{\rm BS}}{p_{\rm SNR}} \right)^{1/2} =\left(\frac{\rho v_{\rm rel}^2}{2 n_{\rm p} kT} \right)^{1/2} \sim 3.7, \end{equation} \tag{ 10 }$$

where p_BS is the pressure at the bow shock and p_SNR is the internal pressure of the SNR (Bucciantini 2002; van der Swaluw et al. 2001). The current expanding velocity of the SNR can be derived with the Sedov solution, v_SNR = 2/5(R_SNR/τ_SNR) ∼ 860 km s⁻¹ for τ_SNR = 9.8 kyr. Therefore, the relative velocity between the pulsar and its surroundings is v_rel = v_PSR − (3/4)v_SNR ∼ 1300 km s⁻¹. Adopting the number density inside the SNR n_p = 0.28 cm⁻³, the mass density becomes ρ ∼ 4.7 × 10⁻²⁵ g cm⁻³. Therefore, the entire length of the termination shock becomes $4.7\: r_{\rm TS}^{\rm F} \sim 4{\buildrel{\prime\prime}\over{.}} 1$, which can be resolved by Chandra. In fact, the obtained image of the point source seems to exceed the PSF slightly, although the ellipticity is uncertain. The measured size of the coma is about 17 ± 04 in FWHM, which reveals another constraint on the spin-down luminosity of the pulsar, $\dot{E} \sim 1.5\times 10^{36}$ erg s⁻¹, which seems to be one order of magnitude smaller than the inferred value based on the empirical relation shown above.

Zoom In Zoom Out Reset image size

As shown in Figure 9, the spectral index Γ_PWN steepens from 1.7 to 2.5 with increasing distance from the point source. The spectral evolution can be explained by synchrotron cooling or adiabatic cooling, hence we can conclude that the point source at the PWN apex drives the entire PWN, although no pulsation has been detected at this moment.

Let us assume that the shocked pulsar wind balances the pressure of the gases within the SNR, whereupon the flow velocity of the pulsar wind will be constant. Here, we applied a theoretical model proposed by Tanaka & Takahara (2010), who assumed that the PWN is expanding at constant velocity, to the bow shock nebula in which relativistic gas is flowing at v_rel and expanding within the bow shock. For simplicity, we also assumed a constant opening angle of the bow shock that bounds the shocked pulsar wind. In this model, a volume element at r(t) = v_relt from the pulsar is expressed as dV = Ωr(t)²dr, where dr is the thickness of the volume element and Ω is the solid angle of the bow shock. So the energy loss of the relativistic electrons is

$$\begin{equation} \dot{\gamma }_{\rm ad} = \left(\frac{\partial \gamma }{\partial V}\right) \left(\frac{\partial V}{\partial t}\right) = -\frac{3}{2}\frac{v_{\rm rel}}{r(t)}\gamma. \end{equation} \tag{ 11 }$$

Meanwhile the relativistic electrons also lose their energy via synchrotron radiation. The radiation loss is

$$\begin{equation} \dot{\gamma }_{\rm syn} = -\frac{4}{3} \frac{\sigma _{\rm T}}{m_{\rm e}c} U_{\rm B} \gamma ^2, \end{equation} \tag{ 12 }$$

where γ is the Lorentz factor of relativistic electrons, σ_T is the Thomson cross-section, U_B is the magnetic energy density, and m_e is the electron mass.

Subsequently, we compare the cooling times of these processes focusing on the rectangular region "C" (Figure 8), just upstream of the point where the spectral index changes. Here, we assumed that the shocked pulsar wind in region C could be in an equipartition condition and the emitting volume is a cylinder of length 20'' and radius 20'' that corresponds to ∼5.8 × 10⁵⁴ cm³ for a distance of 4.1 kpc. The spectral fitting (Figure 9) yielded a photon index, Γ_PWN = 1.7, and an unabsorbed energy flux of the power-law component, 4.2 × 10⁻¹³ erg s⁻¹ cm⁻², for this region. We therefore infer that the spectral luminosity at 1 keV is L_1 keV = 8.9 × 10¹⁴ erg s⁻¹ Hz⁻¹. Based on the equipartition argument (Longair 1994), we can constrain the minimum energy of the electrons and magnetic field in this region,

$$\begin{equation} E_{\rm min} = 3.5 \times 10^{43}(1+\beta)^{4/7}f^{3/7} \; erg, \end{equation} \tag{ 13 }$$

where f is the filling factor and β describes the energy contribution from protons and heavy ions. The obtained energy density corresponds to a magnetic field of B ∼ 26 μG in the equipartition condition. The cooling timescale for a relativistic electron emitting synchrotron X-rays with energy is therefore

$$\begin{equation} \tau _{\rm syn} =\frac{\gamma }{\dot{\gamma }_{\rm syn}} = 280\left(\frac{B}{26\:\mu {\rm G}}\right)^{-3/2} \left(\frac{\epsilon }{1\:{\rm keV}}\right)^{-1/2} \; yr. \end{equation} \tag{ 14 }$$

The adiabatic expansion also cools the relativistic electrons, although it does not change the spectral index since the process does not depend on the electrons' energy. Based on Equation (11), we can obtain the timescale of adiabatic cooling,

$$\begin{equation} \tau _{\rm ad} = \frac{3}{2}\frac{\gamma }{ \dot{\gamma }_{\rm ad}} = \frac{3}{2}\frac{r(t)}{v_{\rm rel}}. \end{equation} \tag{ 15 }$$

The cooling timescale is therefore about 900(1/sin i) yr for r(t) = 40'' and v_rel = 1300 km s⁻¹. The obtained timescale is longer than the cooling timescale of synchrotron estimated above, but is still not negligible.

Meanwhile, in terms of dynamics, the relativistic electrons in region C travel ∼30'' while cooling, hence the photon index steepens by ΔΓ = 0.5. Adopting a transverse pulsar velocity with respect to the surrounding SN ejecta v_rel ∼ 1300 km s⁻¹, a we obtain travel time of τ_dyn = 450 yr, which is consistent with the synchrotron cooling time obtained above.

The magnetization parameter σ of the pulsar wind just upstream of the termination shock can be estimated from the ratio of the termination shock radius to the size of the PWN. Since the shock transforms the kinetic energy _kin of the pulsar wind into the internal energy _int, σ can be rewritten as follows:

$$\begin{equation} \sigma = \frac{B_{1}^{2}/4\pi }{\epsilon _{\rm kin}} \sim \frac{B_{1}^{2}/4\pi }{\epsilon _{\rm int}} \sim \frac{B_{1}^{2}/4\pi }{B_{Eq}^{2}/4\pi }, \end{equation} \tag{ 16 }$$

where B₁ is the toroidal magnetic field just upstream of the shock and B_Eq is the equipartition magnetic field after the shock (Shibata et al. 1998). With the assumption of a toroidal magnetic field around the pulsar, combining the equations of continuity ($r_{\rm TS}^{\rm B\: 2}\; c/3 = r_{\rm n}^{2} v_{\rm n}$) and the conservation law of magnetic flux (r_TS cB₁ = r_nv_nB_Eq) yields

$$\begin{equation} \frac{r_{\rm n}}{r_{\rm TS}^{\rm B}} \sim \frac{1}{3\sqrt{\sigma }}. \end{equation} \tag{ 17 }$$

Substituting the backward shock radius, $r_{\rm TS}^{\rm B} = 3{\buildrel{\prime\prime}\over{.}} 2$, and the radius of the nebula r_n = 20'' that corresponds to the distance from the pulsar to the brightest point at which the pulsar wind achieves equipartition, we found σ ∼ 3 × 10⁻³. The obtained magnetization is about one order of magnitude smaller than that for the Crab pulsar σ_crab = 0.03 (Mori 2002) and that for PSR B1509−58 σ_B1509 = 0.01 (Yatsu et al. 2009), which were calculated in the same manner, though it effectively matches that for the Crab pulsar based on traditional estimation using expansion velocity (Kennel & Coroniti 1984).

We found a non-thermal X-ray filament nearby the PWN, despite its old age and hence relatively slow shock velocity. Here, we confirm whether or not the emission from the filament can be explained by shock acceleration. The limited photon statistics mean that we could not investigate precise spectral parameters, but the filament seems to emit X-rays up to a few keV (Figure 2). We therefore assumed a power-law spectrum of accelerated electrons with an exponential cutoff at E_max, N(E)∝E^−pexp (E/E_max) based on Aharonian & Atoyan (1999). Moreover, these electrons emit synchrotron radiation with a peak energy of

$$\begin{eqnarray} E_{\rm peak} &=& \frac{3}{2}\gamma ^2 \frac{e B_{\rm shell}}{2\pi m_{\rm e}} h = 6.9 \times 10^{-5} \left(\frac{B_{\rm shell}}{1\:\mu {\rm G}}\right)\nonumber \\ &&\times \left(\frac{E_{\rm max}}{1 {\rm \:TeV}}\right)^2 \; keV, \end{eqnarray} \tag{ 18 }$$

where B_shell is the magnetic field in the emitting region. If we assume E_peak = 0.1 keV, then the pair of E_max and B_shell can be constrained as

$$\begin{equation} E_{\rm max} \: B_{\rm shell}^{1/2} = 39 \; TeV \mu G^{1/2}. \end{equation} \tag{ 19 }$$

To solve E_max and B_shell independently, we took into account the shell width.

First, we assumed that the non-thermal emission originates upstream of the shock where the electrons are trapped and accelerated to relativistic velocity. According to Bamba et al. (2003), the magnetic field and maximum energy of the electrons can be estimated using the width of the rising edge. In this study, we employed a width of the outer edge Δr = 10'' that corresponds to ∼6.1 × 10¹⁷ cm and a shock velocity of v_SNR = 860 km s⁻¹. It should be noted that the photon statistics of the X-ray image were too poor to resolve the accelerated particles spreading upstream of the shock, therefore the Δr in this argument is only an upper limit. If the shell width is determined by competition between the advection and diffusion of electrons upstream of the shock, then the diffusion coefficient D can be written as follows:

$$\begin{equation} D \equiv \frac{1}{3} \eta c r_{\rm g} = v_{\rm SNR} \Delta r, \end{equation} \tag{ 20 }$$

where η is the ratio of the mean free path to the gyro radius, the so-called gyrofactor, and r_g = (E_max/eB_shell) is the gyroradius. Combined with Equation (19), we found the maximum energy of the electron and the upstream magnetic field of the shock,

$$\begin{eqnarray} E_{\rm max} &=& 13 \: \eta ^{-1/3} \left(\frac{\Delta r}{10^{\prime \prime }}\right)^{1/3} \left(\frac{d}{4.1\:{\rm kpc}}\right)^{1/3}\nonumber \\ && \times \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s^{-1}}}\right)^{1/3} \; TeV, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} B_{\rm shell} &=& 8.4\: \eta ^{2/3} \left(\frac{\Delta r}{10^{\prime \prime }}\right)^{-2/3} \left(\frac{d}{4.1\:{\rm kpc}}\right)^{-2/3}\nonumber \\ && \times \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s^{-1}}}\right)^{-2/3} \; \mu \mbox{G}. \end{eqnarray} \tag{ 22 }$$

Based on the theory of diffusive shock acceleration, we can estimate the time to accelerate electrons up to the energy required to emit X-ray photons,

$$\begin{equation} \tau _{\rm acc} = \frac{3}{v_{\rm u} - v_{\rm d}} \left(\frac{D_{\rm u}}{v_{\rm u}} + \frac{D_{\rm d}}{v_{\rm d}}\right), \end{equation} \tag{ 23 }$$

where D_u and D_d are the diffusion coefficients defined in Equation (20) for upstream and downstream, respectively, and v_u is the downstream flow velocity in the shock rest frame. Since we do not know the magnetic field just above the shock where the magnetic field may be amplified by turbulence, we therefore roughly assumed D_u = D_d = E_maxcη/(3eB_shell), while employing the jump condition for the flow velocity, v_u = 4v_d(= v_SNR). Adopting those assumptions, we obtain

$$\begin{eqnarray} \tau _{\rm acc} &\simeq & 20\: \frac{D_{\rm u}}{v_{\rm u}^2} = 2.1 \times 10^3 \: \eta \: \left(\frac{E_{\rm max}}{1\: {\rm TeV}}\right) \left(\frac{B_{\rm shell}}{1\: {\rm \mu G}}\right)^{-1}\nonumber \\ &&\times \left(\frac{v_{\rm SNR}}{1000\:{\rm km\:s^{-1}}}\right)^{-2} \; \mbox{yr}. \end{eqnarray} \tag{ 24 }$$

For E_max = 13 TeV and B_shell = 8.4 μG, the acceleration time becomes

$$\begin{eqnarray*} &&\tau _{\rm acc} \sim 4.4 \eta \left(\frac{\Delta r}{10^{\prime \prime }}\right) \left(\frac{d}{4.1\:{\rm kpc}}\right) \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s}^{-1}}\right)^{-1} \; \mbox{kyr}. \end{eqnarray*}$$

Conversely, the cooling time due to the synchrotron radiation limits the maximum energy of the electron, and thus should exceed the acceleration time shown above. The cooling time for the electrons with given energy E_max and a magnetic field B_shell can be written as

$$\begin{eqnarray} &&\tau _{\rm cool} = \frac{3 m_{\rm e} c}{4 U_{\rm B} \sigma _{\rm T} \gamma } \simeq 1.3 \times 10^7 \: \left(\frac{E_{\rm max}}{1 \: {\rm TeV}}\right)^{-1} \left(\frac{B_{\rm shell}}{1 \: {\rm \mu G}}\right)^{-2} \; \mbox{yr}. \nonumber \\ \end{eqnarray} \tag{ 25 }$$

Substituting Equations (21) and (22) into the above formula, we obtain a cooling time of

$$\begin{eqnarray*} &&\tau _{\rm cool} \sim 14\: \eta ^{-1} \: \left(\frac{\Delta r}{10^{\prime \prime }}\right) \left(\frac{d}{4.1\:{\rm kpc}}\right) \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s}^{-1}}\right) \; \mbox{kyr}. \end{eqnarray*}$$

The above estimations are based on the width of the non-thermal shell shown in Figure 3 that is relatively unclear, therefore the obtained τ_acc and τ_cool may be lower and upper limits, respectively. Nevertheless, the obtained cooling time τ_cool seems longer than τ_acc, and the result thus supports the fact that the discovered thin X-ray filament can be explained by synchrotron emission from the shock front of the SNR.

The width of the non-thermal filament can also be interpreted as competition between the flow velocity and the cooling time of the electrons downstream of the shock, τ_cool ≃ Δr/v_d ∼ 900 yr. Substituting 4v_d = v_SNR and Equation (25), the possible pair of E_max and B_shell is constrained as

$$\begin{eqnarray} E_{\rm max} B_{\rm shell}^2 &=& 1.4 \times 10^4\: \left(\frac{\Delta r}{10^{\prime \prime }}\right)^{-1} \left(\frac{d}{4.1\:{\rm kpc}}\right)^{-1} \nonumber \\ && \times \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s^{-1}}}\right) \; TeV\:\mu G^2. \end{eqnarray} \tag{ 26 }$$

Subsequently, combining the above equation with Equation (19), we obtain

$$\begin{eqnarray} E_{\rm max} &=& 5.5\: \left(\frac{\Delta r}{10^{\prime \prime }}\right)^{1/3} \left(\frac{d}{4.1\:{\rm kpc}}\right)^{1/3} \nonumber \\ && \times \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s^{-1}}}\right)^{-1/3} \; \mbox{TeV}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} B_{\rm shell} &=& 50\: \left(\frac{\Delta r}{10^{\prime \prime }}\right)^{-2/3} \left(\frac{d}{4.1\:{\rm kpc}}\right)^{-2/3} \nonumber \\ &&\times \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s^{-1}}}\right)^{2/3} \; \mu \mbox{G}. \end{eqnarray} \tag{ 28 }$$

The obtained magnetic field B_shell is somewhat large for an SNR but still seems acceptable because the magnetic field at the brightest vicinity in MSH 15−56 at a radio energy band amounts to 45 μG (Dickel et al. 2000). If we take an upstream magnetic field of B_ISM = B_shell/4 ∼ 13 μG using a conservative estimate, then the acceleration time can be evaluated from Equation (24),

$$\begin{eqnarray*} &&\tau _{\rm acc} \sim 1.2 \: \eta \: \left(\frac{\Delta r}{10^{\prime \prime }}\right) \left(\frac{d}{4.1\:{\rm kpc}}\right) \left(\frac{v_{\rm SNR}}{860\:{\rm km\:s}^{-1}}\right)^{-3} \; \mbox{kyr}. \end{eqnarray*}$$

The obtained value is comparable to a cooling time of τ_cool ∼ 900 yr, although that is one order of magnitude shorter than the estimation based on Bamba et al. (2003). Consequently, these estimations based on two different scenarios both support the existence of synchrotron electrons around the shock front.