HIGH-ENERGY X-RAY IMAGING OF THE PULSAR WIND NEBULA MSH 15−52: CONSTRAINTS ON PARTICLE ACCELERATION AND TRANSPORT

Since the pulsar B1509−58 is very bright, we needed to minimize its contamination in the PWN imaging and spectral analysis. For the Chandra data, the pulsar contamination was removed by image filtering. However, we were not able to do the image filtering for NuSTAR because its point-spread function (PSF) is broad. Therefore, we selected the off-pulse interval in the NuSTAR data for the image and spectral analyses below.

We extracted the pulsar events in the NuSTAR observations in a 30'' radius circle in the 3–79 keV band, barycenter-corrected the events using R.A. = 15^h13^m5552, Decl. = −59°08'088 (J2000; Caraveo et al. 1994) to produce event lists, and divided each observation into three sub-intervals, yielding 12 sub-intervals for the four observations. We performed H test (de Jager et al. 1989) on the event lists to measure the period for the subintervals and fit the period to a linear function to find the spin period and the spin-down rate. The pulsations were measured with very high significance in each subinterval, and the measured period and the spin-down rate were 0.1517290191(14) s and 1.5281(4) × 10⁻¹² s s⁻¹ for 56450 MJD, respectively. We folded the light curves using the measured period, and show the resulting pulse profile in Figure 1. We used phases 0.7–1.0 for the PWN to minimize the pulsar contamination in all the subsequent NuSTAR data analyses in this paper. The other phase interval was used for the pulsar analysis which will be presented elsewhere. We note that there is contamination from the DC emission of the pulsar even in the off-pulse interval. For example, ∼1.6% of the DC counts are expected in a circle of R = 30'' at a distance of 60'' from the pulsar, but much less at larger distances.

Zoom In Zoom Out Reset image size

In order to produce energy-resolved PWN images, we first produced a merged Chandra image of the five observations in Table 1 using the merge_obs tool of CIAO 4.5 in the 0.5–2 keV, 2–4 keV, and 4–7 keV bands with bin size 4 pixels (Figure 2). Note that the central 4'' ×4'' corresponding to the pulsar emission was removed in these images.

Zoom In Zoom Out Reset image size

For NuSTAR observations, we extracted events in the energy bands 3–7 keV, 7–12 keV, 12–25 keV, and 25–40 keV for the off-pulse phase. After the phase selection, the pulsar component is expected to be reduced significantly. The NuSTAR absolute aspect reconstruction accuracy on long timescales is ∼8'' (90% confidence), which can blur the resulting merged image obtained with the three observations and two modules. Therefore, we aligned the images by registering the pulsar to the known position before phase filtering. Since only one point source (the pulsar) was significantly detected in each observation, we were not able to fully correct the position (e.g., for translation, rotation and scale). We note that the rotational misalignments are measured and corrected with high accuracy (Harp et al. 2010), and a small residual change in scale is not a concern for the spatial scales of our analyses. Therefore, we assumed that the position offsets were caused by pure translations.

In order to produce deblurred images of the NuSTAR observations for comparison with the low-energy high-resolution Chandra images, we corrected for the exposure and deconvolved the NuSTAR images with the PSF using the arestore tool of CIAO 4.5. We then merged the deconvolved images (see Figure 2). The number of iterations in the deconvolution process was determined by comparing the deconvolved NuSTAR image to the Chandra image in a similar band. We chose the energy bands so that the average photon energy weighted by the response and the spectrum in a NuSTAR band is similar to that in a Chandra band, and used the 3–7 keV and 4–7 keV bands for NuSTAR and Chandra, respectively.

Using the 2D images, we produced projected profiles along the jet axis (the southeast to northwest direction) in order to compare the deconvolved 3–7 keV NuSTAR profile with the 4–7 keV Chandra profile. Here, we filled the pulsar region in the Chandra data which were removed above with the average counts of the surrounding pixels. We rotated the images 60° clockwise with the origin being the pulsar position, so that the jet structure lies in the horizontal direction (x axis). We projected the images in Figure 2 onto the axis along the jet, subtracted background, smoothed the profile over a 25'' scale, and normalized the scale with respect to the brightest point at the center. The backgrounds were assumed to be flat over the detector chips. The background normalization factor was first determined by taking a box in a source-free region, and then further adjusted by matching the y-projected profiles of the source and the background at large distance from the center for each energy band. We found that the results presented below are not sensitive to the background subtraction since background accounts for only small fraction of the intensity. We find that NuSTAR-measured profile in the 3–7 keV band is similar to that measured with Chandra in the 4–7 keV band (see dashed and dot–dashed curves in Figure 3), and that the results of the deconvolution are not very sensitive to the number of iterations (e.g., 15–50), and we used 20 iterations.

Zoom In Zoom Out Reset image size

While the deconvolved NuSTAR image (Figure 2(e)) shows similar overall morphology to the high-resolution Chandra image in the similar band (Figure 2(d)), there are differences. Most notably, the small arc-like structure and the elongation in the central region (R ≲ 30'') are not resolved in the NuSTAR images. This is because the structures are smaller than the FWHM (∼18'') of the NuSTAR PSF. Also note that RCW 89, ∼6–7' north, is not clearly visible in the NuSTAR data. This is mainly because the NuSTAR observations did not have much exposure in that region; most of RCW 89 fell outside the FoV during the observations.

To measure the size of the narrow structures around the pulsar, we projected events between −50'' and 50'' in the y axis direction onto the x axis in several energy bands. The profile is very asymmetric and is not smooth on large scales (R ≳ 100''). Furthermore, the deconvolution produces artificial structures in the outer regions due to the paucity of counts. Therefore, defining a size (e.g., full width 1% maximum) is impractical on a large scale. However, the source images are smooth on smaller scales (R ≲ 50''), allowing us to measure the FWHM and HWHM in the northern and the southern directions of the projected profiles in several bands without smoothing the images. We measured the sizes by calculating the relative brightness with respect to the peak, and show them in Figures 4(a)–(c). Although there are different structures in the northern and the southern directions, the HWHM's are similar to each other and to half the FWHM.

Zoom In Zoom Out Reset image size

We calculated the spectrum- and response-weighted average energy for each energy band, fit the widths to a power-law function R(E) = R₀E^m as suggested by Reynolds (2009), and measured the decay index m for various y-integration widths (e.g., ∼40–130''). The measured decay index was stable over this range, as shown in Figure 4(d). Note, however, that our measurement is based on deconvolved images, and our uncertainties are therefore approximate.

The image analysis shows a spectral change with radius in the PWN and we therefore tried to see differences in spectra at different radii. We extracted spectra in various regions as described, and backgrounds from source-free regions, and fit them with an absorbed power-law model.

Since spatial blurring due to the PSF size is much more significant for NuSTAR than Chandra, we did not attempt to fit the spectrum jointly except for one case of using a large aperture (Section 3.3.1), where PSF "blurring" is not a large effect. However, we jointly fit the spectra taken with single telescope at different epochs.

We used the χ² and the lstat statistics in XSPEC 12.8.1 to fit the spectra (Arnaud 1996). Results from the two methods were consistent, and we primarily report the results obtained with the χ² statistics. Since the NuSTAR data are not sensitive to the hydrogen column density (N_H) and the results are not affected by small change of N_H, we froze it at a previously reported value (0.95 × 10²² cm⁻²; Gaensler et al. 2002). We used a cross-normalization factor to account for a slight difference between NuSTAR module A and B, and between observations. For the Chandra data fitting, we let N_H vary and introduced a cross-normalization factor between observations.

3.3.1. Spectrum of the Entire Nebula

We first measured the total spectrum of the PWN using a source extraction aperture of R = 5' centered at the pulsar position for a phase interval 0.7–1.0 in the NuSTAR data. Photons with energies up to ∼20–30 keV were detected above background for each observation. We extracted a spectrum for this aperture in each of the three NuSTAR observations, N1, N3 and N4 in Table 1, and jointly fit the spectra to a power law. The measured power-law index is 2.06 and the absorption-corrected 3–10 keV flux is 5.9× 10⁻¹¹ erg cm⁻² s⁻¹ (see Table 2). For the Chandra data, we extracted source spectra using the same 5' aperture, ignoring the central 5'' in order to minimize the pulsar contamination, and jointly fit the five Chandra spectra of each region to a power-law model in the 0.5–7 keV band because background dominates above 7 keV. We note that the Chandra data fits were not acceptable with χ²/dof of 2632/2214 (p = 1 × 10⁻⁹), having large residuals in the low-energy band below 2 keV. This is perhaps because the large regions are a mixture of subregions with different spectra (e.g., see Section 3.3.4). We therefore fit the data above 2 keV only with frozen N_H. When removing photons below 2 keV, the remaining data were fit to a single power-law model with a slightly smaller photon index, having χ²/dof = 1757/1704 (p = 0.18). The cross-normalization factors for the five Chandra observations are all within 1%. Note that letting N_H vary also yields an acceptable fit (χ²/dof = 1756/1703) with N_H = 0.91(4) and Γ = 1.90(1).

Table 2. Best-fit Parameters for the Total PWN Emission Spectrum

Dataa	Modelb	Radius	Energy	NHc	Γs	FPLd	Ebreak	Γh	χ2/dof
		'	(keV)	(1022 cm−2)			(keV)
C	PL	5	2–7	0.95	1.912(5)	5.98(2)	...	...	1757/1704
C	PL	5	3–7	0.95	1.90(1)	6.00(3)	...	...	1408/1359
N	PL	5	3–7	0.95	2.03(2)	5.93(7)	...	...	543/587
N	PL	5	3–20	0.95	2.06(1)	5.91(5)	...	...	1895/1887
N	PL	5	15–30	0.95	2.1(1)	6.6(1)	...	...	576/579
N + C	PL	5	2–20	0.95	1.950(4)	5.77(5)	...	...	3878/3592
N + C	BPL	5	2–20	0.95	1.918(5)	5.88(5)	6.3(3)	2.12(2)	3691/3590

Notes. 1σ uncertainties are given in parentheses at the same decimal place as the last digit. ^aN: NuSTAR, C: Chandra. ^bPL: powerlaw model, and BPL: bknpower model in XSPEC. ^cFrozen. ^dAbsorption corrected 3–10 keV flux in units of 10⁻¹¹ erg s⁻¹cm⁻².

Download table as:  ASCIITypeset image

The Chandra-measured spectrum in the 2–7 keV band is significantly harder than that measured with NuSTAR in the 3–20 keV band. We note that our results are consistent with the previous measurements made with BeppoSAX and INTEGRAL (Mineo et al. 2001; Forot et al. 2006). Mineo et al. (2001) reported a photon index of 1.90(2) in the 1.6–10 keV band for a 4' aperture, which is consistent with our Chandra measurement in the 2–7 keV band. In the hard band, the reported photon indices were 2.1(2) and 2.12(5) for BeppoSAX (20–200 keV) and INTEGRAL (15–100 keV), respectively. Note that the large apertures used for BeppoSAX and INTEGRAL include the RCW 89 region, but the effect of RCW 89 is negligible because the emission is very soft (Yatsu et al. 2005) and the telescopes operate only above 15 keV. The photon index we measure with NuSTAR in the 15–30 keV band is 2.1(1) for the 5' aperture, which agrees with the previous measurements. The results of our measurements are summarized in Table 2.

Since the large apertures include many subregions with different spectral properties as we show below (see Sections 3.3.1, 3.3.3, and 3.3.4), a single power law may not properly represent the combined spectrum. In particular, we find that the best-fit photon index for the Chandra data becomes smaller as we ignore lower energy spectral channels, that is, the spectrum appears to harden (is concave up) as we move to higher energies. However, this is the opposite to what we see with NuSTAR (Table 2). While this may imply a spectral break in the X-ray band, some other effects such as contamination from the pulsar and/or RCW 89, or cross-calibration systematics between the two instruments may have some impact. Therefore, we investigate some possibilities below in order to see if the discrepancy in the spectral index measurements of NuSTAR and Chandra is caused by a spectral break.

First, we note that the pulsar contamination was not completely removed in the Chandra data. Excising 5'' leaves 2%–5% pulsar emission in the R = 5' aperture, which may bias the PWN spectrum. In order to see the effects quantitatively, we fit the pulsed spectrum of the pulsar in the NuSTAR data (total spectrum minus the DC level in Figure 1) with a power-law model and find that the photon index is 1.36(1) and the 3–10 keV flux is 2.24(3) × 10⁻¹¹ erg s⁻¹ cm⁻², which broadly agree with the previous measurements (Cusumano et al. 2001; Ge et al. 2012). We added this pulsar component to the Chandra fit assuming 5% of the pulsar emission is in the 5' aperture after the 5'' excision. Thus, the spectral model was a double power-law model, one for the pulsar emission and the other for the PWN emission. We froze the pulsar component, fit the PWN spectrum, and find that the spectral index of the PWN does not change. Since the spectral index of the pulsed spectrum may be different in the soft band, we changed the spectral index of the pulsar component to 1.19 as reported by Cusumano et al. (2001) in the 1.6–10 keV band, and find that the photon index of the entire PWN softens only by ΔΓ = 0.01. We further increased the pulsar flux by 10% and find no change in the spectral index of the PWN. We verified the results by increasing the excision region to 10''. Note that the DC component of the pulsar is not included in this study. However, the unmodeled DC component is much smaller than the pulsed component as seen in Figure 1, so the effect would be negligible in the Chandra data fit.

Second, we consider the effect of the pulsar DC component in the NuSTAR data. Although the DC component is negligible in the Chandra data due to the image filtering, the DC component presents in the NuSTAR data because time filtering does not remove the DC emission. Although it is not possible to measure the DC spectrum accurately, we estimated 3–10 keV DC flux using a 30'' aperture as follows. With NuSTAR, we measure the total (pulsed+DC+nebula) and the pulsed flux to be 3.29 × 10⁻¹¹ erg s⁻¹ cm⁻² and 2.24 × 10⁻¹¹ erg s⁻¹ cm⁻², respectively. The nebula flux is measured to be 7.5 × 10⁻¹² erg s⁻¹ cm⁻² with Chandra (Section 3.3.3). By subtracting the pulsed and the nebula flux from the total flux, the 3–10 keV DC flux is estimated to be 3.1 × 10⁻¹² erg s⁻¹ cm⁻². We assumed that the photon index is 1.7, similar to the NuSTAR-measured value for the 30'' aperture (Section 3.3.3). We included the DC emission in the NuSTAR fit of the 5' aperture spectrum, and followed the procedure described above for the pulsar contamination estimation in the Chandra data. This procedure effectively removes the DC component from the PWN spectrum. However, note that removing such a hard spectrum only softens the spectrum of the entire PWN, making the discrepancy larger. We therefore arbitrarily changed the photon index of the DC component to 2.5 to mitigate the possibility of having very soft DC emission and find that the photon index of the entire PWN hardens only by ΔΓ = 0.02.

Third, we varied the background level by ±30%, and found that spectral indices change only by 0.02 and 0.01 for Chandra and NuSTAR, respectively. We also used different background regions and found that the spectral indices do not change significantly.

Finally, we estimate the contamination of the RCW 89 emission in the NuSTAR data. Using the NuSTAR PSF, we estimated the contamination of a structure at ∼6' into the 5' circle to be ∼14%. We extracted the RCW 89 spectrum from the NuSTAR observation N1 which sampled the RCW 89 best among the NuSTAR observations. We added the spectrum as additional background in the PWN spectral fit. Since the other observations N3 and N4 sampled only a small fraction of the RCW 89 region, we used the RCW 89 spectrum extracted from N1 for these observations as well. We varied the normalization of the RCW 89 background from 0.14 to 0.28 in order to account for the fact that the actual RCW 89 region may be larger than what we sampled with NuSTAR and found that the photon index changes by ≲0.02.

The large discrepancy in the spectral index measurements between NuSTAR and Chandra cannot be explained by a combination of the above effects. We therefore consider alternatives below. We note that there may be cross-calibration systematics between NuSTAR and Chandra. For example, Kirsch et al. (2005) showed that systematic uncertainties between X-ray observatories caused by cross calibration are significant using observations of the Crab nebula, for which the authors found that Chandra and BeppoSAX/LECS measured smaller spectral indices than BeppoSAX/MECS and INTEGRAL did. Since NuSTAR is calibrated so that the spectrum of the Crab nebula is a simple power law with a photon index of 2.1, larger than the Chandra-measured value of 1.95 (Kirsch et al. 2005), the spectral break we see for MSH 15−52 could be explained by imperfect cross-calibration.

However, the cross-calibration effects may be different from source to source depending on the spectral shape, and the measurements made for the Crab nebula may not be directly translated into the case of MSH 15−52. We therefore consider an alternative; the discrepancy of spectral indices between the Chandra and the NuSTAR measurements is caused by a spectral break. We jointly fit the 2–7 keV Chandra and the 3–20 keV NuSTAR data with a single power law and a broken power law. We find that a single power law does not describe the data well having χ²/dof = 3878/3592 (p = 5 × 10⁻⁴), but a broken power law with a break energy of 6.3 keV does (χ²/dof = 3691/3590, p = 0.12). The results are presented in Table 2 and Figure 5.

Zoom In Zoom Out Reset image size

We compared the NuSTAR and the Chandra spectra in the same energy range below the possible break at 6 keV in order to see if the break is a single sharp break. If so, we expect the NuSTAR and the Chandra spectra to be same regardless of the difference in the effective area shape. We fit the NuSTAR and the Chandra spectra in the common energy band, defined as 3–E_max, where we vary E_max from 4.5 keV to 7 keV. In the lowest energy band (3–4.5 keV), the NuSTAR and Chandra results agreed with photon indices of 1.99(7) and 1.98(2). However, as the upper energy range increased to 5 keV and above, the NuSTAR results were significantly softer than that of Chandra. For example, the photon indices are 2.03(2) and 1.90(1) in the 3–7 keV band for NuSTAR and Chandra, respectively, inconsistent with each other with 90% confidence. This suggests that the spectral break is likely to be caused by the cross-calibration effect. However, if the spectral break is real, the observational discrepancy in the common energy band may imply that the broadband spectrum is not sharply broken at 6.3 keV but slowly curves over a energy range (e.g., 4–7 keV) probably because different regions in the 300'' aperture have different break energies. In this case, NuSTAR collects relatively more photons above the break than Chandra does since it has rising effective area in that band, yielding a softer spectrum.

3.3.2. Azimuthal Variation of the Spectrum

In order to see if the PWN spectrum varies azimuthally, we first extracted NuSTAR events in 30'' radius circles for six, twelve, and eighteen azimuth angles for three radial distances, 60'', 120'', and 180'' from the pulsar, respectively. The regions do not overlap. For each region, backgrounds were extracted from an aperture of R = 45'' in a source-free region on the same detector chip. We jointly fit the NuSTAR spectra for the three observations N1, N2, and N4. The energy ranges for the fit were 3–20 keV, where the source events were detected above the background. We performed these analysis for the same regions with the Chandra data in the 0.5–7 keV band.

We show the results in Figure 6, where the azimuth angle ϕ is defined from east in a clockwise direction. A sinusoidal variation of the spectral hardness is clearly visible for each radial group, and the spectrum is hardest in the jet region (ϕ ∼ 300°). The flux values also peak in the jet regions but do not seem to vary sinusoidally. Note there is a small discrepancy between the NuSTAR and Chandra measurements in Figure 6; the NuSTAR-measured spectral indices are larger than those measured with Chandra in general. While this may suggest the spectral break we see in the spectrum of entire nebula (Section 3.3.1), it could be due to PSF mixing in the NuSTAR data; when there is a sharp spatial contrast such as the jet or image edge, it is convolved with the PSF in the NuSTAR data.

Zoom In Zoom Out Reset image size

It appears that the photon index covaries with N_H in Figure 6. We checked if there is a correlation between the two quantities using Pearson's product moments, and found no clear correlation in any radial group.

3.3.3. Spectral Variation in the Northern Nebula and in the Jet

Since we observe significant spectral variation in the azimuthal direction, we analyzed the northern nebula and the jet separately. For the northern nebula we used annular regions, ignoring the southern part. Hence, each region covers the upper ≳180° in azimuth angle. The innermost region was a circle with radius 30'' centered at the pulsar, and annuli with width 30'' or 60'', and boxes were used for the outer regions (see Figure 7(a) white). We used the off-pulse phase only for NuSTAR, and ignored the central pulsar using a circle with radius 5'' for Chandra.

Zoom In Zoom Out Reset image size

We separately fit the NuSTAR (N1, N2, and N4) and the Chandra (C1–C5) spectra of each region with an absorbed power-law model. Since there is significant thermal contamination from RCW 89 in the Chandra data at large distances (see the two upper white box regions in Figure 7(a)), we ignore the low-energy data below ∼3 keV for the two rectangular regions for the Chandra fits. We also tried to model the RCW 89 regions using the vnei plus a power law in XSPEC, and fit the Chandra data in the 0.5–7 keV band. The result for the photon index was sensitive to the remnant model but broadly agree with what we found by fitting data above 3 keV only (see also Yatsu et al. 2005). Our results for the Chandra data are consistent with, but more accurate than those obtained by DeLaney et al. (2006) who used ∼60 ks of observations taken from 2000 August to 2003 October. We present our measurements in Figures 7(b) and (d).

While our results show that N_H increases with radius, we note that it is possible to force the N_H to be constant and allow only the photon index to vary. For example, an N_H value of 0.957 ± 0.005 × 10²² cm⁻² constant over the field with photon indices of 1.58–2.16 fits the data out to R = 200'' in Figure 7 (χ²/dof = 12289/12363), which implies no radial variation of N_H and smaller variation of the photon index with radius. However, the N_H profile given in Figure 7(d) provides better fit with χ²/dof = 12223/12353 corresponding to F-test probability of 2 × 10⁻¹⁰. We also verified that the results do not change if we fit the spectra only in the 0.5–6 keV (below the spectral break, see Section 3.3.1) in order to reduce the effect of the complex continuum. We note that better constraining Γ and N_H by jointly fitting the NuSTAR data is not possible because of the PSF mixing in the NuSTAR data and the spectral break (Section 3.3.1).

There is a hint of a possible break in the linear slope of the radial profile of the photon index Γ(R) (Figure 7(b)). We measured the location of the break in the northern nebula using a broken line fit. We first fit the Chandra measured photon index profile, and found that the break occurs at R_break = 71 ± 3''. We note that using a constant N_H over the field changes Γ only slightly and gives a consistent result (R_break = 68 ± 2''). The NuSTAR profile gives a larger R_break = 150 ± 10'' because of the large photon index at smaller radii which might be biased by mixing from outer regions. We also find that a single power-law model Γ(R) = Γ₀R^η with η = 0.149 ± 0.003 broadly agrees with the data (see Figure 7(b)).

Since spectral softening is expected in the jet direction as well, we measured the spectral variation along the southern jet. To do this, we extracted source spectra using non-overlapping circular apertures with radii 10'', 10'', 15'', 20'', 25'', 35'', and 25'' along the jet (see yellow circles in Figure 7(a)), which we refer to as regions J1–J7. Note that the center of J1 is R ∼ 15'' from the pulsar, and all the NuSTAR observations (N1–N4) were used for this analysis. We fit the spectra in each region with an absorbed power law, and measured the photon index and flux. The results are presented in Figures 7(c) and (e).

The spectral indices of the J2 region measured with Chandra and NuSTAR are very different, which may suggest that there is a strong spectral break. However, we note that measuring the spectral parameters with NuSTAR was difficult for regions with sharp spectral changes because the NuSTAR PSF changes from a circular shape to an elliptical shape with off-axis angle (An et al. 2014), and thus regions with different off-axis angles have different degrees of azimuthal mixing. The four NuSTAR observations had different pointings and thus different off-axis and azimuthal angles. In particular, in the regions J1–J2 where we use small apertures and the source spectrum strongly varies, spatial mixing has significant impact on the NuSTAR results. Therefore, the discrepancies in the spectral index between the NuSTAR and Chandra measurements, and even between the NuSTAR observations are expected. The mixing was not a concern in the analysis of the northern nebula in which spectral variation is not severe.

We also measured the location of the break in the radial profile of the photon index in the jet direction using the Chandra measurement in Figure 7(c). Here we ignored the first data point for the reason described below. A fit to a broken line gave a break location of R_break = 110 ± 30''. A single power-law model also fits the data with a power index η = 0.12 ± 0.01 (Figure 7(c)).

We note that the first Chandra data point, corresponding to region J1, shows a very soft spectrum compared to the next one in J2, unexpected in synchrotron cooling models (e.g., Reynolds 2003; Tang & Chevalier 2012). The steady-state solutions may not be applicable to the inner region within ∼1' of the pulsar for this source, as DeLaney et al. (2006) found strong variability in the brightness and morphology in that region.

We tried to see if the spectral hardness in J1 varied over time. We first jointly fit the Chandra spectra of the region taken from the five observations C1–C5 with a common N_H, but separate photon index and cross normalization for each observation, and found that the photon indices are all within the 1σ uncertainty of the value in Figure 7. We carried out the same analysis for the J2 region, and found that the spectrum of one observation (C1, Obs. ID 754) was slightly softer than the others (Γ = 1.62 ± 0.12) but not significantly. It is probably because the region in this observation fell on the detector chip gap. Therefore, we conclude that the spectral index did not change significantly over time in this region.

Since spectral hardness covaries with N_H, the spectral difference between the J1 and J2 regions may be less significant if we consider the covariance. In order to investigate the effect of covariance, we ignored Obs. ID 754 because the J2 region in this observation was on the chip gap. We then jointly fit the spectra of each region, varied both N_H and Γ using the steppar tool in XSPEC and found that the 99% contours do not overlap, which suggests that the difference is significant with the covariance as well, and the spectrum of the J1 region is significantly softer than that of the J2 region. If we take the best-fit values, the N_H variations of ∼3 × 10²¹ cm⁻² imply extremely high densities of n ∼ 2000 cm⁻³ for an assumed line-of-sight distance of 0.5 pc (similar to the transverse distance for the assumed distance to the source of 5.2 kpc) in the regions with high N_H.

We further spatially resolved the J1–J2 regions using overlapping circular regions with radius 5''. We fit the spectrum of each region with a power-law model, and found spectral softening in the innermost regions (R ≲ 20''). We show the photon indices in Figure 8.

Zoom In Zoom Out Reset image size

Note that N_H increases with radius in the northern nebula. In the jet direction, we used finer spatial resolution, and see a more complicated change; there is a dip at R = 30–70''. At large distances, we find that N_H is large. This structure is visible in the 2D N_H map as well (see Figure 9). Note also that the power-law index (η) of the photon index profile is larger in the northern nebula than in the jet, that is, the spectral steepening is more rapid, which was also implied by the imaging analysis above (e.g., Figure 3).

Zoom In Zoom Out Reset image size

3.3.4. 2D Maps of the Spectral Parameters

The NuSTAR images in the hard band (≳7 keV; Figures 2(e)–(h)) show that the source shrinks with energy in the 2D projection on the sky, which is attributed to the synchrotron burn-off effect. While the effect has been shown for this source in a previous study of azimuthally integrated spectra (Schöck et al. 2010), this is the first time it is shown with a 2D imaging analysis over a broad X-ray band. In particular, we found that the burn-off effect is stronger in the northern direction than in the jet direction (Figure 3). Although it is very difficult to build a full 2D model that can reproduce the measured images, matching overall morphological changes in energy may give us important clues to understand particle outflow in the PWN.

The change of the size with energy can be used for inferring properties of the particle flow in the PWN using theoretical models (Kennel & Coroniti 1984; Reynolds 2009). For example, the results of Reynolds (2009) can be rewritten in a form appropriate for sources whose spectral break Δ between radio and X-ray power laws, and size index m can both be measured:

$$\begin{equation*} {\Delta \over (-m)} = \left(1 \over \epsilon \right) (1 + 2\epsilon + (3 + 2\alpha _r) m_\rho /3 + (1 + \alpha _r) m_b), \end{equation*}$$

where m_x is the index of an assumed power-law function for a quantity $x\propto R^{m_x}$ (x = ρ or b for mass density and magnetic-field strength, respectively), α_r is the energy index of the radio spectrum, and is a confinement parameter (e.g., = 1 for conical jets, Reynolds 2009). For the values we derive for MSH 15−52 of α_r = 0.2, Δ = 0.86, and m = −0.2 (Section 3.2), the formula gives

$$\begin{equation*} 1 + 1.1\, m_\rho + 1.2\, m_b = 2.3 \epsilon. \end{equation*}$$

This condition requires either that the mass in the flow is not constant (for instance, due to mass evaporated from thermal gas filaments joining the outflow; Lyutikov 2003), or that magnetic flux is not conserved (for instance, due to turbulent amplification or reconnection): either m_b or m_ρ, or both, must be positive. Various combinations of gradients can reproduce our results. For instance, a conical flow = 1 requires 1.1 m_ρ + 1.2 m_b = 1.3, approximately met if density is constant and magnetic field rises as radius—or vice versa. A confined flow with = 0.44 (roughly parabolic) would have m_ρ≅ − m_b, which could be satisfied with both constant density and constant magnetic field, or with one dropping as fast as the other increases.

We find that the integrated spectrum of the entire nebula measured with NuSTAR is a simple power law with photon index 2.06 (Table 2). This is similar to values previously reported based on BeppoSAX and INTEGRAL data (Γ = 2.08 ± 0.01, 2.12  ±  0.05; Mineo et al. 2001; Forot et al. 2006). The Chandra-measured parameters for a single power-law model above 2 keV imply a harder spectrum than that measured with NuSTAR (see Table 2), which is also seen in the spectra extracted for other inner regions (e.g., Figures 6 and 7). We note that our Chandra spectral fit results are consistent with that measured previously with BeppoSAX (Γ = 1.90 ± 0.02 in the 1.6–10 keV band for a 4' aperture, Mineo et al. 2001).

Since the large aperture include many subregions with different spectral parameters, we expect to see a harder spectrum at high energies if the spectra of the subregions are simple power laws in the 0.5–20 keV. However, the 3–20 keV NuSTAR spectra are much softer than the Chandra spectra, which is the opposite to what is expected from a sum of simple power-law spectra. We find that contamination of the pulsar and the backgrounds can explain only ∼0.01 of the photon index discrepancy of the NuSTAR and Chandra measurements, while the measured difference is 0.15. The correction for the pulsar, the RCW 89 contamination and the background variation in the data seem not to explain the discrepancy between the NuSTAR and the Chandra results.

We find that the discrepancy between the NuSTAR and the Chandra measurements is likely to be caused, at least in part, by imperfect cross calibration of the instruments. However, if the break is real, it implies a break in the energy distribution of the shock accelerated electrons at ∼200 TeV (having the peak synchrotron power at 6 keV) for a magnetic field strength of 10 μG, assuming synchrotron emission (e.g., Equation (5) in TC12).

We note that a spectral break in the X-ray band has recently been reported for G21.5−0.9 (Nynka et al. 2014) and the Crab Nebula (Madsen et al. 2014) and may be common to other young PWNe as well. If true, this has important implications on the particle acceleration mechanism in PWNe. Sensitive broadband X-ray observations of other PWNe will be helpful.

Using the broadband X-ray data obtained with NuSTAR and Chandra, we find an azimuthal variation of the spectral index. For the PWN 3C 58, Slane et al. (2004) suggested a possible azimuthal variation of the spectral hardness based on a scenario where current flows out from the pulsar's pole and returns in the equator (Blandford 2002). However, Slane et al. (2004) did not find an obvious azimuthal variation in the PWN 3C 58 within R ∼ 2 pc for a distance 3.2 kpc we find for MSH 15−52 at R = 1.5–4.5 pc. Furthermore, we find the azimuthal spectral variation in MSH 15−52 is likely sinusoidal and different from that of the flux. This azimuthal spectral variation of the emission may hint at a large-scale current flow, however, it could also be due to azimuthal diffusion of jet particles in MSH 15−52.

A radial change of the spectral index of MSH 15−52 was reported by Schöck et al. (2010). While they integrated the spectrum over the full azimuthal angle, we measured the profiles for the northern and the jet directions separately because the two regions are different. We found that significant softening with radius is seen in both directions, more significantly in the northern region. Interestingly, the radial profile of the photon index (rate of spectral steepening) flattens with radius as is also seen in 3C 58 (Slane et al. 2004).

An outflow model considering both diffusion and advection was developed by Tang & Chevalier (2012, TC12 hereafter), where they calculated the change of the spectral index with distance from the central pulsar with an assumed electron injection spectrum and diffusion coefficient, and were able to reproduce the radial variation of the spectral index for three compact PWNe, the Crab nebula, G21.5−0.9, and 3C 58. The model has been applied to PWNe where the particle escape times (the times for particles to diffuse a distance R in the Bohm limit) are longer than their ages (Equation (2) in TC12):

$$\begin{equation*} t_{\rm esc}\approx 16{,}000\left(\frac{R_{\rm PWN}}{2\ \rm pc}\right)^2\left(\frac{E_e}{100\ \rm TeV}\right)^{-1} \left(\frac{B}{100\ \mu \rm G}\right)\ \rm yr, \end{equation*}$$

where R_PWN is the radius of the PWN, E_e is the energy of synchrotron emitting particles, and B is the magnetic-field strength in the PWN. Using the size R ∼ 10 pc, E_e = 100–600 TeV (Forot et al. 2006; Nakamori et al. 2008), and an estimation of the magnetic-field strength of 8–17 μG (Gaensler et al. 2002; Aharonian et al. 2005), we find that t_esc is 5000–7000 yr, greater than the spin-down estimated age of τ_c = 1700 yr. Since the photon spectrum we are using in this work corresponds to a smaller E_e, t_esc can be larger than the above estimation. However, we note that there have been suggestions that the true age of the PWN is ≳ 6000 yr based on the association with RCW 89, larger than the spin-down estimation (e.g., Seward et al. 1983). If so, there may be particles escaping the PWN and the particle spectrum in the outer parts of the PWN becomes steeper, which may be the case for old (≳ 10⁵ yr) PWNe. For young PWNe, TC12 uses a reflecting boundary condition at the outer edge of the PWN. We note that the TC12 model is for spherically symmetric PWNe and may not be optimal for MSH 15−52. However, the azimuthal variation in the northern nebula is not large (see Figure 6), and thus the model may provide a reasonable description of the source in that region.

In this model, the angular size of the "flat" region where the radial profile of the spectral index is flat can be used to estimate the diffusion coefficient (Equation (14) of TC12). This is calculated using the following equations:

$$\begin{equation*} \theta _{\rm flat}\approx \theta \left(\frac{6^{1/2}}{2}\left[\frac{\nu _R}{\nu }\right]^{1/4} -1\right) \end{equation*}$$

and

$$\begin{equation*} \nu _{R}=1\times 10^{17}\left(\frac{D}{10^{27}\ \rm cm^2\ s^{-1}}\right)^2 \left(\frac{1\ \rm pc}{R_{\rm PWN}}\right)^4\left(\frac{100\, \mu \rm G}{B}\right)^3\ \rm Hz, \end{equation*}$$

where θ_flat is the angular size of the region that has a flat photon index profile, θ is the angular size of the PWN, ν is the photon frequency, and D is the diffusion coefficient.

We find that photon index profile steepens more slowly beyond R = 71'' and R = 110'' in the northern and the jet regions, respectively. Note that the radial profile of the spectral index in the northern nebula shows a flat region between R = 70''–200'' although the profile seems not to show any flat region in the southern nebula. We use the value for the northern nebula for the size of the flat region of TC12. Assuming the size of the source is R_PWN ∼ 300'' and using the above formulae with ν = 2.4 × 10¹⁷ Hz (1 keV), we estimated the diffusion coefficient to be 4–13× 10²⁷ cm² s⁻¹ for B = 8–17 μG, which is slightly larger than that estimated for 3C 58 by numerical modeling (2.9 × 10²⁷ cm² s⁻¹; TC12).

Using the diffusion coefficient we estimated above, we calculate the critical particle energy E_R for which the diffusion distance is equal to the size of the PWN (see Table 3 of TC12) and where the electron distribution has a break (Gratton 1972), using formulae given by TC12: R = (4D/QE_R)^1/2 and Q = 1.58 × 10⁻³B² erg s⁻¹. For MSH 15−52, we find E_R to be 130–190 TeV for B = 8–17 μG. It is interesting to note that this is similar to the maximum electron energies inferred from broadband spectral-energy-distribution modeling (130 or 250 TeV; Nakamori et al. 2008), and that inferred from the possible spectral break at 6 keV we measured in Section 3.3.1.

We note that the spectrum in the jet direction is significantly softer in the innermost J1 region compared to that in the farther J2 region. Such behavior is not expected in simple advection and/or diffusion models, since the synchrotron emitting particle spectrum only softens with distance. This simple picture may not be appropriate in the regions where the particle flow may be more complicated due to magnetic hoop stress as suggested for this source by Yatsu et al. (2009). The authors found a ring-like structure with R ∼ 10'' using Chandra data and interpreted the structure as the termination shock for this PWN. Based on the morphology, the authors further suggested that the shock accelerated particles are diverted and squeezed toward the poloidal direction right below the ring due to magnetic hoop stress (e.g., Lyubarsky 2002). We also find that the jet structure becomes narrower to R ∼ 35'' and then broader. Furthermore, the spectral hardness turn-over, non-monotonic variation of the spectral index (see Figure 7), happens near the location where the jet is narrowest, which might be occurring because of the compression of the magnetic fields and particles.

We presented 2D maps of the spectral parameters. The maps visualize the properties of the source very well, and can be compared with 3D PWN models.

We showed that the 2D map of N_H has a shell-like structure. The density is low near the central pulsar, increasing out to R ∼ 3' (see also Figures 6 and 7). We note that the fit value of N_H could in principle be degenerate with other spectral parameters. However, we do not find clear evidence of correlation between N_H and photon index or flux from our analysis (see Section 3.3.4), and using a constant N_H degrades the fit significantly.

A higher column density is observed in the south and the east directions (see Figure 9). If the material responsible for N_H was produced by the supernova, one would expect the pulsar to have a kick in the opposite direction of the material, toward the northeast direction, which is consistent with the direction of the kick velocity for PSR B1509−58 estimated by Rots (2004) based on 2800 days of timing. However, Livingstone & Kaspi (2011) found no evidence of proper motion using 28 yr of timing data. Nevertheless, we do not see any enhanced emission in the shell-like structure, which makes the supernova ejecta scenario less plausible.

Alternatively, the structure may be an interstellar bubble produced by the stellar wind of the supernova progenitor (e.g., Castor et al. 1975). In the wind model, the size of the bubble is given by a simple formula:

$$\begin{equation*} R_s(t)=28\left(\frac{\dot{M}_6 V_{2000}^2}{n_0}\right)^{1/5} t_6^{3/5}\ \rm pc, \end{equation*}$$

where $\dot{M}_6$ is the mass loss rate of the progenitor in units of 10⁻⁶M_☉ yr⁻¹, V₂₀₀₀ is the speed of the wind in units of 2000 km s⁻¹, n₀ is the number density (cm⁻³) of the interstellar medium, and t₆ is the time in units of 10⁶ yr. The radius of the ring structure we observe is ∼5 pc, much smaller than the calculated value for a typical O6 star. However, the value can vary significantly for different input parameters, and due to spatial non-uniformity of the interstellar medium or radiative loss (e.g., Weaver et al. 1977).

The asymmetric shell structure may be explained by outbursts of massive stars. The massive star progenitors of core-collapse supernovae undergo a variety of instabilities that drive episodic mass loss: e.g., pulsations driven by bumps in the continuum opacity (Lamers & Nugis 2002; Fryer et al. 2006; Paxton et al. 2013) and explosive shell burning (Quataert & Shiode 2012; Arnett et al. 2014). These outbursts occur up until the collapse of the star and are believed to have large asymmetries. An outburst a few thousand years prior to collapse could explain the features we observe at 5 pc in the remnant.

If, as suggested above, the shell-like structure at R ∼ 3' was formed by the stellar wind, the structure would have to avoid being swept up or destroyed by the SN ejecta. If the supernova ejecta did not fill a full spherical shell, a part of the wind-produced shell can be left over. In this case, density of the shell is expected to be higher in the direction where the supernova ejecta were less dense. We see such a trend when comparing our N_H map with the radio image of the SNR (Figures 2 and 3 of Gaensler et al. 1999); there are more ejecta in the northern region than in the southern region.

We have estimated the mass of hydrogen contained in the observed shell-like structure. Using the measured radial profile of N_H shown in Figure 7(d), the excess mass compared to the central region is ∼460 M_☉, large compared to the ∼12 M_☉ one would estimate for a sphere with R = 5 pc for typical interstellar density of 1 cm⁻³. Furthermore, the large amount of material in the structure should produce HI emission, which we do not see in the 20 cm map (Figure 4 of Gaensler et al. 1999). This may be because the radio continuum emission was not subtracted in the radio map and/or because the X-ray measurement is sensitive only to foreground material while the radio observations are sensitive to both foreground and background structures.

We note that we cannot unambiguously rule out the possibility of a constant N_H over the field; the observed N_H being an artifact of a more complex underlying continuum. Thus, it is very difficult to clearly interpret the structure using X-ray observations only. Nevertheless, if the shell-like structure in the N_H map is intrinsic to the source, it may support the idea of the existence of an underdense region around the supernova progenitor, which was suggested to explain the discrepancy between the pulsar's characteristic age of ∼1700 yr (Kaspi et al. 1994) and the SNR age of >10,000 yr based on the RCW 89 association (Seward et al. 1983). This requires further confirmation by observations in other bands.