NuSTAR STUDY OF HARD X-RAY MORPHOLOGY AND SPECTROSCOPY OF PWN G21.5−0.9

We extracted NuSTAR data using a 165'' radius circular region centered at the pulsar PSR J1833−1034 position. This region represents the largest circle within the same detector chip (to ensure that NuSTAR detector background is the same and predictable), and it encloses over 93% of the PWN photons.

We first fit the NuSTAR spectra with a single absorbed power-law model in the 3–45 keV energy band. We obtained a best-fit power-law index of Γ = 2.04 ± 0.01 and an X-ray flux of F_2–8 keV = 5.54 ± 0.03 erg s⁻¹ cm⁻². These results are consistent within their uncertainties to the parameters reported by Tsujimoto et al. (2011) (Γ = 2.05 ± 0.04, F_2–8 keV = 5.7 ± 0.5 erg s⁻¹ cm⁻²), who analyzed G21.5−0.9 for the purposes of cross-calibrating X-ray instruments operational at that time. All uncertainties (90% confidence level) quoted in the text and tables include both systematic and statistical errors. A more detailed comparison of G21.5−0.9 spectral fitting between NuSTAR and other X-ray telescopes as well as systematic error estimates will be addressed in a separate calibration paper.

A single power-law fit is not satisfactory, with reduced χ² = 1.33, with clear residuals evident between 5 and 10 keV (see the left panel of Figure 1). We then fit a broken power-law model (bknpower) to the NuSTAR spectra, yielding reduced χ² = 1.09. Statistical tests using the F-distribution (ftest in XSPEC) show the broken power-law model is statistically favored over the single power-law model with high significance, so that a spectral break is unambiguously required. The best-fit power-law indices were Γ₁ = 1.996 ± 0.013 and Γ₂ = 2.093 ± 0.013, with the best-fit break energy of E_break = 9.7 ± 1.3 keV.

Zoom In Zoom Out Reset image size

In order to investigate any instrumental effects that mimic a break around 9 keV, we analyzed NuSTAR observations of the known power-law source 3C273. We followed the same procedures used for G21.5−0.9, and found no systematic residuals around 9 keV in the NuSTAR data. A spectral break in another PWN, the Crab Nebula, has recently been detected by NuSTAR at similar energies of 8–11 keV (K. C. Madsen, in preparation).

Dust scattering can also mimic the low-energy spectral behavior observed in the NuSTAR spectra. Recent analysis of two low-mass X-ray binaries with similar absorption columns of ∼3 × 10²² cm⁻², GX5−1 and GX13+1, present radial profiles of models of the dust scattering halos. These profiles indicate that 0.1% of the source photons at 2.5 keV are enclosed between 50'' < r < 600'' (Smith et al. 2002). Including the energy-dependant halo intensity profile relationship I∝E⁻² (Predehl & Schmitt 1995; Smith 2008), a radius of r < 160'' encloses over 99% of the source photons in both 3–5 keV and 5–8 keV. Similar analysis reveals that the r < 30'' extraction region, presented in Section 3.2, has a scattering loss of ∼0.2% and ∼0.1% of source photons in 3–5 keV and 5–8 keV, respectively. Loss of photons by dust scattering therefore has a negligible effect, and cannot explain the spectral break in the NuSTAR spectrum. Additionally, because our scattering estimates are based on point sources, we cannot conclusively rule out the possibility of a very small effect of dust scattering on the radial behavior of the low-energy spectrum presented in Section 3.2.

An X-ray spectral break at ∼12 keV was first suggested by Tanaka & Takahara (2011). They extrapolated the spectral fits from the hard and soft X-ray bands reported by Tsujimoto et al. (2011), and noted that a break at ∼12 keV was likely. NuSTAR not only confirms the break, but provides the first measurement of the ∼9 keV break in broadband 3–45 keV spectra of G21.5−0.9 in a single X-ray telescope.

Obtaining a well-calibrated X-ray spectrum from G21.5−0.9 is particularly valuable when creating PWN theoretical models. Specifically, Tanaka & Takahara (2011) presented synchrotron emission from radio to GeV energies by modeling the electron spectrum emanating from the pulsar of G21.5−0.9. They incorporated adiabatic and energy losses in a one-zone model. Figure 2 shows the NuSTAR best-fit broken power-law data overplotted on Figure 3 from Tanaka & Takahara (2011). More details are presented in Section 6.1.

Zoom In Zoom Out Reset image size

A number of young PWNe, including G21.5−0.9, exhibit spectral softening from the center of the PWN outward. This is the classic signature of synchrotron burn-off. Slane et al. (2000) and later Safi-Harb et al. (2001) extracted Chandra spectra at various annuli of G21.5−0.9 and showed spectral softening from Γ = 1.43 ± 0.02 in the central 5'' radius circle to Γ = 2.13 ± 0.06 in the outer annulus at a radius of ∼40''. NuSTAR has the ability to probe the synchrotron burn-off effects in the hard X-ray band above 10 keV.

We extracted NuSTAR spectra from the central 30'' region as well as from four nested annuli, each 30'' in width. After experimenting with various annulus widths, we found that the annulus width of 30'' is the best compromise between the NuSTAR angular resolution and the expected spatial variation from the Chandra results.

The broadband 3–45 keV fit results from both the absorbed powerlaw and bknpower models are shown in Table 1. The power-law index fit with a single power-law model increases from Γ = 1.97 ± 0.01 in the 30'' inner region to Γ = 2.25 ± 0.02 in the outer 90''–120'' radius annulus, confirming the spectral softening discovered by Chandra. However, the residuals for the inner 30'' and 30''–60'' regions also show a spectral break around 9 keV. We fit a broken power-law model, and the results are shown in Table 1 and plotted in Figure 3. We conclude that the spectral break at 9 keV is detected with high significance in the inner 30'' radius circular region and the 30''–60'' annulus. Interestingly, as is shown in Figure 3, the second (high energy) power-law component shows no spatial variation. The spectral softening is observed only below 9 keV.

Zoom In Zoom Out Reset image size

Table 1. Spectral Fitting of the NuSTAR G21.5−0.9 Data

Parameter	r ⩽ 165''		r ⩽ 30''		r = 30''–60''		r = 60''–90''
	Power Law	Broken P.L.	Power Law	Broken P.L.	Power Law	Broken P.L.	Power Law
Γ1	2.039 ± 0.011	$1.996^{+0.013}_{-0.012}$	$1.964^{+0.011}_{-0.012}$	1.852 ± 0.011	2.051 ± 0.012	$1.98^{+0.02}_{-0.03}$	2.09 ± 0.014
Γ2	⋅⋅⋅	$2.093^{+0.013}_{-0.012}$	⋅⋅⋅	$2.099^{+0.019}_{-0.017}$	⋅⋅⋅	$2.14^{+0.02}_{-0.03}$	⋅⋅⋅
Ebreak	⋅⋅⋅	$9.7^{+1.2}_{-1.4}$	⋅⋅⋅	$9.0^{+0.6}_{-0.4}$	⋅⋅⋅	9.74 ± 1.0	⋅⋅⋅
Fx(2–8 keV)	5.45 ± 0.03	5.27 ± 0.08	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
Fx(15–50 keV)	5.47 ± 0.03	5.11 ± 0.08	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
χ2/dof	1.33(5117)	1.09(5112)	1.32(3564)	0.99 (3557)	1.08(2257)	0.960(2255)	0.942(1924)

Notes. Spectral fitting of G21.5−0.9 with various extraction regions. The column header indicates the region from which the data was extracted. N_H was frozen to 2.99 × 10²² cm⁻² for all fits. Flux is listed in units of 10⁻¹¹ erg s⁻¹ cm⁻². The goodness of fit is evaluated by the reduced χ². The errors are 90% confidence level.

Download table as:  ASCIITypeset image

A measurement of the PWN radius in different energy bands is complementary to the radially dependent spectral analysis. The angular resolution of NuSTAR is comparable to the PWN size (∼80''). Therefore, it is essential to take into account the blurring effects by the PSF convolution.

We adopted a forward-folding method to measure the size of the PWN using Sherpa, CIAO's modeling and fitting package (Fruscione et al. 2006). The use of Sherpa allows us to use the goodness-of-fit test (C-statistics) to find the best-fit parameters. We convolved an assumed source model, a circular 2D Gaussian profile, with the effective PSF described in Section 3, and fit the NuSTAR flux images. We also ran the conf command in Sherpa to determine 90% confidence intervals for all fit parameters.

We confined our fitting range to within a radius of 160''. In general, a circular 2D Gaussian provides a good fit to the entire PWN, leaving only small residuals. The best-fit Gaussian FWHMs in five energy bands are plotted in Figure 4. The trend of decreasing PWN size with energy is evident, confirming the synchrotron burn-off effect at energies above 10 keV. In Section 6.2 we employ the FWHM radius to measure the synchrotron cooling length. We fit the data with a power-law L(E)∝E^m. The best fit of L(E) is shown in Figure 4 and yields an index of m = −0.21 ± 0.01. These results are used in Section 6.2 to constrain physical conditions in the nebula.

Zoom In Zoom Out Reset image size

The Eastern Limb and North Spur are only partially resolved in the raw NuSTAR mosaics due to the small size of G21.5−0.9 compared to the NuSTAR PSF, combined with the low surface brightness of these features. To detect these features with high confidence, we applied two different methods. First, we confirmed the detection of the Eastern Limb and Northern Spur in the flux images in each energy band (Section 4.2.1). Second, we sharpened the NuSTAR flux images using the Lucy–Richardson deconvolution algorithm (Richardson 1972; Lucy 1974).

4.2.1. 1D Profile Analysis

In order to detect the Eastern Limb and North Spur in the raw NuSTAR image, we analyzed line intensity profiles of 6 image pixels (a total of ∼9'') in width. We chose projection axis lines along different orientation angles around the center of the PWN. The left image in Figure 5 contains two lines, colored red and blue, that indicate the positions along which the flux profiles were taken. The red line was chosen to fall along the regions of brightest Eastern Limb emission, while the blue line overlays a region of G21.5−0.9 that does not contain any emission from either the Eastern Limb or the North Spur. The intensity profiles on the right of Figure 5 correspond to the lines drawn in the left of Figure 5. The red profile shows clear and significant excess emission when compared to the blue profile, confirming that the Eastern Limb is detected in the 3–6 keV band. This also holds for the 6–10 keV and 10–15 keV. While the two profiles in the 15–20 keV band are very similar, the lines are separate and distinct between 60''–40'', when accounting for statistical uncertainty. With these intensity profiles we confirm the detection of the North Spur and Eastern Limb at energies as high as 20 keV.

Zoom In Zoom Out Reset image size

We repeated this process with the North Spur, orienting the red line along the north/south axis. Lastly, we also confirmed the detection of the Eastern Limb and North spur in residual maps from fitting the PWN, as mentioned in Section 4.1. The residual maps had areas of faint excess, indicating that there exists emission aside from the PWN. However, the broad NuSTAR PSF prevents us from studying any further morphology of the faint emission. We therefore turn to image deconvolution in an attempt to remove the effects of the PSF.

4.2.2. Image Deconvolution: Method and Verification

We applied an iterative deconvolution technique to the NuSTAR images using arestore (CIAO v4.4) and the effective PSFs described in Section 4. The iterative image deconvolution can produce artificial features if the process is over-iterated and/or the background region is deconvolved. Great care was taken to ensure that the deconvolved image was both stable and reproducible. We deconvolved the 3–6 keV mosaiced NuSTAR image with several iterations, from 20 up to 200 in increments of 20. Each of these images exhibited the same features of the North Spur and Eastern Limb that are characterized by areas of brightened emission to the north and northeast of the PWN. Additionally, we confirmed that the features visible through deconvolution were not dependent on telescope rotation, detector module, or observation. The data sets were first grouped by NuSTAR module, summing all the FPMA images before deconvolving, and likewise for all the FPMB images. The deconvolved images have identical features. We also grouped the data by epoch, summing the 2012 and 2013 data separately. As before, the two deconvolved images are very similar.

As a final verification of our deconvolution process, we compared the output images to the 3–6 keV image from Chandra. Figure 6(a) shows the Chandra image at the top left corner. One can see the Eastern Limb and North Spur, which is identified as "knot" in the image. The North Spur is visible as the excess of emission 100'' north of the PWN. The Eastern Limb is an arced feature that begins at the southern edge and wraps clockwise around to the northern edge, located approximately ∼150'' from the center of the plerion. For reference, the outer green circle has a radius of 165''. Figure 6(b) shows the final deconvolved 3–6 keV NuSTAR image, created with 20 iterations. The NuSTAR features correspond well to those seen by Chandra, allowing us to confidently extend this deconvolution technique to higher energies.

Zoom In Zoom Out Reset image size

4.2.3. Image Deconvolution: Analysis

Once we established that the NuSTAR results reproduce the Chandra images stably for a range of Lucy–Richardson iterations, we applied the same deconvolution to the higher energy images.

Figures 6(b)–(f) show the NuSTAR deconvolved images at 3–6 keV, 6–10 keV, 10–15 keV, 15–20 keV, and 20–25 keV, respectively. One can visually identify the Eastern Limb and North Spur at energies as high as 15 keV, as seen in Figure 6(d). This is not surprising, since all of the previously reported spectral fitting of these two regions have required a non-thermal component. This is, however, the first direct measurement showing that these features have emission above 10 keV.

At energies above 15 keV the Limb and Spur become very faint. One can still see emission along the eastern edge of the Limb in Figure 6(e). We can, however, verify the existence of these features above 15 keV by producing intensity profiles of the deconvolved images at various azimuthal angles. Applying the same one-dimensional (1D) profile method as in Section 4.2.1, we obtained profiles along three different radial lines, each bisecting G21.5−0.9 at angles in increments of 45° as shown in the left image of Figure 7. The arrows overlaid on the lines indicate the direction along which the profiles were obtained.

Zoom In Zoom Out Reset image size

The intensity profiles themselves for the 3–6 keV image are shown on the right in Figure 7. Each profile corresponds to the line of matching color overlaid on the deconvolved image on the left. One can clearly see a sharp increase corresponding to the Eastern Limb in the blue and yellow profiles, while the red excess indicates the existence of the North Spur. Similarly, the Eastern Limb and North Spur are clearly visible in the profiles taken from the images up to energies of 15 keV (see Figure 8). The intensity profiles from the 15–20 keV image also show statistically significant photon excess from the Limb and Spur, confirming the existence of these features to energies as high as 20 keV.

Zoom In Zoom Out Reset image size

Chandra observations of G21.5−0.9 show spectral softening over the PWN, with the photon spectral index ranging from Γ ∼ 1.4 at the inner (radius < 5'') region to Γ ∼ 2.1 at the outer 35''–40'' radius annulus. NuSTAR observations confirm this by showing spectral softening below 9 keV, as shown in Table 1. The broad NuSTAR PSF causes some mixing of the spectra in different annuli, and as a result the variation of Γ with radius obtained by NuSTAR is less pronounced than that seen by Chandra. NuSTAR finds a spectral break at 9.7 keV in the integrated PWN emission. The spectral break is also observed in the inner regions with radius <30'' and in an annulus from 35''–40''. The break is statistically significant at radii >40''. Above 10 keV the photon spectral index remains constant.

There are several possible origins for the spectral break seen by NuSTAR. One possibility is that the NuSTAR PSF mixes the radially dependent power-law indices seen by Chandra (Slane et al. 2000; Safi-Harb et al. 2001), and softens the spectrum in such a way as to cause a sharp break. However, simulations that fold Chandra maps through the NuSTAR response indicate that this is not the case. The effects of the large interstellar medium (ISM) extinction (Tsujimoto et al. 2011) again would not cause such a sharp, defined energy break. Similarly, the loss effect of dust scattering at lower energies is negligible (below 2%) in the NuSTAR band because of the relatively low column absorption. Contributions from the pulsar are also not likely to be responsible. Most pulsars have a harder spectrum than the PWN they power. While this could in fact cause a sharp transition between photon indices, it would cause spectral hardening with energy, not the softening seen by NuSTAR. Finally, spectral breaks are often attributed to synchrotron cooling, as proposed by Tanaka & Takahara (2011). However, the break energy of 9 keV would require an unreasonably low magnetic field strength of ∼6 μG, compared to the ∼300 μG derived from equipartition arguments (Camilo et al. 2006). The prominent spectral steepening between radio and X-rays, if interpreted with simple cooling models, yields magnetic field strengths ranging from 25 μG (de Jager et al. 2008) to 64 μG (Tanaka & Takahara 2011).

It is likely that the spectral break results from physical effects, either due to a break in the injected electron energy spectrum or due to energy losses due to particle transport in the PWN. Pulsars emit pairs of relativistic electrons and positrons, which are accelerated at a termination shock near the pulsar itself. Downstream of the termination shock the accelerated electrons interact with the magnetic field, also produced by the pulsar, and subsequently emit synchrotron radiation from the radio through gamma energy bands. The injection spectrum can therefore shape and influence the spectrum of the synchrotron radiation, as noted by Tanaka & Takahara (2011). These authors propagate a broken power-law electron spectrum through a time-dependent model that includes energy losses due to synchrotron radiation, IC scattering, and adiabatic cooling. NuSTAR provides confirmation that the relatively simple injection spectrum used in Tanaka & Takahara (2011) is not adequate to fit the observed X-ray data, as seen in Figure 2. The model, noted as the solid black line, has a steep negative slope in the X-ray band and does not fit the X-ray spectra obtained from Chandra, XMM-Newton, INTEGRAL/IBIS, and now NuSTAR.

It is therefore reasonable to explore whether a more complex model can explain the 9 keV spectral break. Vorster et al. (2013) extend the models with the addition of diffusive losses as well as a broken injection spectrum with a discontinuity at the break energy. However, both the aforementioned SED models do not include the spatial dependence of parameters such as the magnetic field of the PWN, which provides additional complexity and can perhaps better fit the X-ray data. We explore these effects in the following sections.

Particle transport in PWNe has long been a matter of discussion and study. Unfortunately, few celestial objects are both bright enough and large enough to allow distinguishing radially dependant features to fit to the various existing models. The Crab Nebula, G21.5−0.9  and 3C 58 are such PWNe. They have been frequently observed and analyzed, and have provided valuable insights into the physics of magnetohydrodynamic (MHD) flows and relativistic shocks that govern the appearance of PWNe (e.g., Kennel & Coroniti 1984; Komissarov & Lyubarsky 2004; Chevalier 2005; Matheson & Safi-Harb 2010). Extending the energy range of analysis will allow for more detailed probing of radially dependant physics. This information is important both to understand PWNe and to isolate the properties of the relativistic shock at the PWN inner edge from the subsequent spectral evolution downstream. We shall ask: can the observable parameters of G21.5−0.9 from radio to X-ray be reproduced assuming the shock injects only a single structureless power-law particle spectrum?

One such parameter used as a spectral fingerprint of various models is the radial dependency of the power-law photon index. Numerous observations have confirmed softening of the G21.5−0.9 PWN spectrum with increasing projected radius (Safi-Harb et al. 2001; Warwick et al. 2001; Matheson & Safi-Harb 2005, 2010; Slane et al. 2000). This is shown at higher energies for the NuSTAR observations in Figure 4. The softening is associated with synchrotron cooling of the electrons, and the corresponding decrease in the maximum emitted X-ray energy in the bulk velocity flow downstream of the termination shock.

An alternative approach is to characterize the variation of source size with photon energy using a characterization of the "cooling length" L(E), such as the source FWHM. Since only the energy dependence of this length is important for modeling the spectral steepening, its precise definition is not important. This scale depends only on mapping the total number of photons as a function of radius and energy to determine the scale length L(E). This is more straightforward than spectral modeling because counting statistics are almost always limited, and length scale measurements do not require determining a spectrum at each radius. We therefore use the cooling length scale, with FWHM as its surrogate, as the fundamental diagnostic for extracting information about physical conditions.

We shall also make use of the observed steepening from radio to X-rays of G21.5−0.9  through the parameter Δ ≡ α_x − α_r. Most discussions of PWN physics (e.g., Chevalier 2005) simply take this as an intrinsic property, but we shall attempt to explain it through evolutionary effects.

There exist two prominent mechanisms that have been invoked to explain PWNe particle transport: advection and diffusion. Beginning with Wilson (1972) and Gratton (1972), diffusion has long been investigated as a cause for the characteristics of a PWN. However, the early models of particles propagating outward from a central source purely by diffusion, applied exclusively to the Crab Nebula, have been unable to account for the detailed X-ray properties of the Crab, such as its change of size with frequency (Ku et al. 1976). Pure advection models have also been proposed, such as Rees & Gunn (1974). The canonical theory involving pure advection as the method of particle transport was presented by Kennel & Coroniti (1984, hereafter KC84). When KC84 was used to predict the radial behavior of the spectrum, the resultant spectral photon index has little to no variation from the center of the nebula outward, and begins to vary only toward the PWN periphery (Tang & Chevalier 2012, hereafter TC12, Figure 2). This also does not match the observed behavior of a slowly steepening X-ray spectrum Matheson & Safi-Harb (2005).

TC12 provide a nuanced approach to particle transport by diffusion and present two updated models. They claim that the magnetic field is not predominately toroidal far from the termination shock, as is often assumed in PWN theory, but has a more complex geometry with cross-field transport best described by diffusion. Their first model consequently incorporates pure diffusion only, with both the magnetic field and diffusion coefficient constant with radius, and synchrotron emission as the only loss of energy. TC12 argue that such a model better explains the radial dependence of Γ as seen in the Crab, 3C58, and G21.5−0.9, with proper adjustment of the diffusion coefficient. Involving complexities such as an energy-dependent diffusion coefficient might be more physically reasonable, but the data are not sufficiently constraining to distinguish these cases from a simple diffusion model.

While the pure diffusion model of TC12 appears to describe observations of the Crab and 3C 58 relatively well, we argue that it is less appropriate for G21.5−0.9. It does provide a good description of Γ(r) for the Crab and 3C 58 (χ² ∼ 1); however, the fit of this model to the G21.5−0.9 >10 keV data is poor (χ² ∼ 3), as seen in Figure 6 in TC12. In addition, the ratio of advective to diffusive timescales that determines what process dominates is radially dependant. Using v∝1/r² downstream of the PWN termination shock (KC84), and the know relationships t_adv∝r/v and t_diff∝r, we obtain t_adv/t_diff∝r. G21.5−0.9 is very compact, with a size more than two times smaller than that of the other two PWNe, and thus advection is likely to dominate.

Finally, the advective model provides a good fit to the energy dependent cooling length scale (see below). TC12 presented the radially dependant spectral index Γ(r) rather than using the cooling length scale L(E) as characteristic of their model, which makes fitting the model to higher energy NuSTAR data too difficult.

The second model presented by TC12 involves a Monte Carlo simulation that includes both diffusion and advection transport close to and farther away from the pulsar, respectively. This allows a more complex treatment than previous analytical models (Massaro 1985). While likely a more appropriate approach for G21.5−0.9, TC12 only applied this simulation to the Crab and 3C 58.

KC84 provide a complete advective solution for a steady, spherically symmetric wind terminated by a MHD shock. Although KC84 is routinely applied to determine quantities such as the mean downstream magnetic field, its range of applicability is in fact more narrow. KC84 represents an idealized theory, suited to the case of constant injection of electrons in a spherical geometry, transported outward by pure advection in an ideal MHD flow with an ordered, toroidal magnetic field. The model does not attempt to reproduce the Crab spectrum from radio to X-rays, but only the optical to X-ray portion, with an injection spectral index of optically emitting electrons of α_o = 0.6. The predicted value of steepening of α_x − α_o of 0.51 is roughly appropriate (so their cooling break is at UV wavelengths). This value is fortuitously close to the value of 0.5 for a stationary, homogeneous source. Thus this model also cannot reproduce radio-to-X-ray SEDs of many other PWNe, including G21.5−0.9, which have a larger Δ. As mentioned above, the Γ(r) relationship predicted by the KC84 model is flatter in the PWN interior and softens toward the edges, as shown by Reynolds (2003) and Tang & Chevalier (2012), while generally a gradual variation in Γ(r) is observed. This has motivated several generalizations of KC84.

The values we obtain from the NuSTAR data for L(E) and α_x are clearly inconsistent with the model presented in KC84. The gradients implied by the assumptions of KC84, predict L(E)∝E^−1/9, independent of spectral index, and Δ = (4 + α)/9 for the physically important inner flow region. If we attempt to apply the KC84 formalism to describe the radio-to-X-ray spectrum of G21.5−0.9 with α_r = 0.0 and α_x = 0.9, we fail on both counts. First, we find m = −0.21 ± 0.01 (Section 4.1). Second, we find the observed spectral indices in the radio and X-ray bands produce Δ = α_x − α_r = 0.9 ± 0.1, instead of 5/9 from a KC model using the injected spectral index of α = 1.

Reynolds (2009, hereafter R09) noted that there are a number of physical effects not accounted for in KC84 that could produce a steeper spectral break than Δ = 0.5. For example, the magnetic field may have a significant radial or turbulent component, or it may not satisfy mass conservation due to cloud evaporation, or it may not have magnetic flux conservation (e.g., due to magnetic reconnection or turbulent amplification). R09 constructed a simple model that includes these effects. The model involves generating simple scaling relations for the downstream magnetic field B, fluid velocity field v, and fluid density ρ in terms of the dimensionless length scale L = r/r_o, where r_o is the inner injection radius. The non-spherical geometrical effects could also be parameterized in terms of jet width w = w_oL, where = 1 corresponds to a conical jet or a section of spherical outflow; a confined jet has < 1, while a flaring jet would have > 1. By assuming B, v and ρ all vary as power laws in L with indices m_b, m_v, and m_ρ, R09 obtains a series of general consistency relations that these indices must satisfy with each other and with observable parameters. These variables of m, α, and Δ, the energy scaling of the cooling length, particle injection index and spectral index break, respectively, provide constraints on the allowed values of , m_b, m_v, and m_ρ.

The results of R09 can be rewritten in terms of the functional form of m, and assuming α_r = 0 as observed for G21.5−0.9:

$$\begin{equation*} \Delta = (-m)(1 + 2\epsilon + m_\rho + m_b)/\epsilon \end{equation*}$$

independent of m_v. Since we observe Δ = 0.9 and m = −0.21, this gives 1 + m_ρ + m_b = 2.29. The possible solutions are thus very restrictive for our measured value of Δ and m. Additionally, if the flow conserves mass (disallowing, for instance, mass loading by evaporation of thermal material), then the density index m_ρ is linked to the velocity index m_v. In this case, some solutions are unphysical, such as those with m_v > 0, corresponding to an accelerating downstream fluid. With some judicious rejection of such solutions we can draw some interesting general conclusions, based on the possible values of the observables, which may guide further investigations.

For a conical or spherical flow, = 1, so either density or magnetic-field strength, or both, must rise with radius (since either m_ρ or m_b is positive, or both). Mass conservation links m_ρ to m_v, but if that assumption is abandoned, there is no relation and no constraint on the velocity profile—only the density profile matters for observable quantities. A steeply decelerating flow can produce m_ρ > 0 with mass conservation, but mass loading can also do this. Similarly, m_b > 0 can be a reasonable outcome of flux nonconservation through processes such as reconnection. Our observations require one or both of these effects: addition of mass to the flow through some kind of evaporation, and increase in magnetic-field strength beyond flux freezing.

A strongly confined jet ( < 1) can relax some constraints; for = 0.3, we only require m_ρ + m_b = −0.32. However, even here either mass or flux nonconservation is necessary. One can construct constant density solutions (m_ρ = 0) but such a geometry is disfavored due to the high symmetry of G21.5−0.9 as evidenced in both the broad axial symmetry observed in the radio (Furst et al. 1998) and in the soft X-ray band (Safi-Harb et al. 2001).

While Reynolds (2009) expands on the treatment of KC84, both fail to reproduce the steadily steepening spectrum with radius shown by G21.5−0.9 and other PWNe. This shortcoming is characteristic of models in which particles are transported outwardly purely by advection, so that all particles at a given radius have similar ages. To produce steady spectral steepening probably requires a mixture of particles of different ages at each radius. This could be caused by more complex fluid flow such as the back flows found in simulations by Komissarov & Lyubarsky (2004), or by particle diffusion.

We conclude that a model describing both the radio-to-X-ray spectrum of G21.5−0.9  and the size shrinkage with X-ray energy we observe can be accommodated in a pure advection model requiring the injection of only a straight power-law spectrum of electrons, N(E)∝E⁻¹. However, as with all pure advection models, the gradual rather than sudden steepening of the spectrum with radius is not reproduced. The viability of this explanation for the observed properties of G21.5−0.9 will depend on whether the addition of diffusion can reproduce the gradual steepening while preserving the successes of the advection model.

NuSTAR has detected, for the first time, the North Spur and Eastern Limb above 10 keV. Three main theories have been proposed to explain the nature of these features: they are extensions of the PWN itself, they are limb-brightened shock fronts propagating into surrounding ISM and accelerating cosmic rays, or they result from an interaction of ejecta with the envelope of the progenitor SNR, presumably a Type IIP SN (Bocchino et al. 2005). Since the North Spur and Eastern Limb have different spectral and spatial properties, we discuss them separately.

A multi-wavelength analysis is required to get a full understanding of the North Spur. This feature was observed in the radio and the soft X-ray, most recently by Bietenholz et al. (2011) and Matheson & Safi-Harb (2010), respectively. Bietenholz et al. reported a radio detection of the North Spur, with a 1.43 GHz flux density of 20.2 ± 1.8 mJy, and a FWHM size of 18'' × 8''. This is notable because it is the only feature, other than the PWN itself that is detected in the radio band. The Eastern Limb has no radio emission detected to date.

Matheson & Safi-Harb (2010) obtained ∼580 ks of Chandra data, and found that the North Spur has a spectrum comprised of non-thermal and thermal components. The thermal component, represented by the pshock model, is best fit to temperature of kT ∼ 0.2 keV and contributes only ∼6%–7% to the overall 0.5–8 keV flux. The non-thermal component, however, is equally well fit by either the srcut or powerlaw model.

The model srcut describes the synchrotron emission from a homogeneous source consisting of a power-law energy distribution of electrons with an exponential cutoff, radiating in a uniform magnetic field. The emitted spectrum is a power law that steepens slowly above the photon energy corresponding to the electron cutoff energy. This slow curvature can mimic a steeper power law in a limited energy band. However, in a broader energy band the srcut model can fall well below the extrapolation of a power law with the same slope at lower energies. Since X-ray emission is visible from the North Spur at energies as high as 20 keV, the correct spectral model must provide a photon flux from 15–20 keV that is statistically higher than the background.

Due to signal-to-noise limitations, we cannot use NuSTAR data to spectrally fit the Spur, however, we can use the deconvolved images to distinguish between the srcut and powerlaw models. We simulated spectra with these models using the fakeit command in XSPEC using the model parameters reported by Matheson & Safi-Harb (2010). A photon spectral index of Γ = 2.21 was applied to the powerlaw model, while the srcut model used a radio index of α = 0.8 and a rolloff frequency of ν_rolloff = 18 × 10¹⁷ Hz. The NuSTAR response files were based on a point source extraction with r = 30''. We were thus able to obtain fluxes from each respective model from 10–15 keV, 15–20 keV, and from 20–25 keV.

Both models have count rates higher than that of the background within the 10–15 keV energy band. In the 15–20 keV band, the powerlaw model has a count rate seven times higher than the background, while the srcut model is only three times higher. Finally, within the 15–20 keV band, the count rates are 4 and 1.5 higher than the background for the powerlaw and srcut models, respectively. This implies that the srcut model should be marginally visible up to 20 keV, and above 20 keV should have a count rate equal to that of the background. This matches well with what is seen in the NuSTAR images. The powerlaw model, however, should be detectable at energies as high as 25 keV.

If the North Spur were an extension of the PWN, it would have a spectral photon index similar to that of the PWN itself. While this is true, the analysis above indicates that the North Spur is not described by a powerlaw model extending to higher energies. If it were, the North Spur would easily be detectable at energies as high as 25 keV. However, the NuSTAR images with a combined exposure of 281 ks do not detect any statistical emission above 20 keV, indicating that srcut is the more plausible spectral model for the North Spur.

It is possible that the North Spur results from the interaction of the inner SN ejecta with the H-envelope of the progenitor (Bocchino et al. 2005; Matheson & Safi-Harb 2010). This is supported by the thermal component in the spectral fit of Matheson & Safi-Harb (2010), a pshock model with solar abundances, low temperature of kT ∼ 0.21 ± 0.4 keV, and low ionization timescales. This is also supported by the morphology of the North Spur itself. With the inclusion of projection effects, the North Spur is located between ∼75'' and the edge of the SN shell at ∼120''.

X-ray emission from the shell of G21.5−0.9 is clearly visible in the NuSTAR image. Our results confirm the existence of this shell, revealed by Chandra (Matheson & Safi-Harb 2010) in both their image and in an extracted shell spectrum (up to ∼6 keV), but also hinted at in earlier XMM-Newton data (Bocchino et al. 2005). The XMM-Newton observations reveal evidence of the shell in a ∼2–8 keV energy band image, after careful subtraction of a modeled dust scattering component below ∼5 keV. The NuSTAR detection extends up to much a higher energy of ∼20 keV. The morphology of the NuSTAR emission is striking in its similarity to the Chandra and especially the XMM-Newton image. Emission is detected from position angle (P.A.) ∼180° to ∼300° in the 6–10 keV image (P.A. = 0 at north, positive clockwise), with the extent shrinking as the energy increases until it is visible mainly in the east and north at the highest energies. This is consistent with the intensity distribution with position angle seen in the lower energy image.

A question unresolved by previous X-ray observations is whether the shell emission is thermal or non-thermal. The extended energy response of NuSTAR can be exploited to answer this question. Matheson & Safi-Harb (2010) found that the spectrum of the Eastern Limb could be characterized equally by four distinct models. A thermal fit to the pshock model provided a temperature of kT ∼ 7.5 keV while a non-thermal powerlaw model obtained a spectral photon index of Γ = 2.13. Two srcut models were also well fit to the Eastern Limb spectrum. Ideally srcut requires both a radio flux density and radio spectral index for the shell, but the Eastern Limb has not been detected in the radio. Therefore, two values of the radio index α that are reasonable for a SN shell were chosen (α = 0.3/0.5) while all other parameters were left free to vary. Care was taken to ensure the best-fit radio flux was below the upper limit reported by Bietenholz et al. (2011).

As with the North Spur, we extrapolated the four Chandra spectra for the Eastern Limb into the NuSTAR band. We created an effective area file for an extended source shaped like the Eastern Limb, then used the fakeit command in XSPEC to simulate spectral data. The thermal pshock model predicted a shell flux which would not produce the X-ray emission seen by NuSTAR at ∼15–20 keV. In contrast, the three non-thermal models produced X-ray fluxes consistent with imaging of the Eastern Limb by NuSTAR, although the srcut and powerlaw models could not be distinguished from each other. Nonetheless, the NuSTAR observations firmly establish the non-thermal nature of the shell X-ray emission.

The detection of a non-thermal shell up to quite high X-ray energy in G21.5−0.9 is interesting, and is in contrast to observations of other Crab-like SNRs. Recently shells have been detected in 3C58 (Gotthelf et al. 2007) and G54.1+0.3 (Bocchino et al. 2010), and a clear, detached shell of emission in Kes 75 (Helfand et al. 2003). However, the 3C58 shell is clearly thermal, with no sign of a non-thermal component. The shell of G54.1 + 0.3 can be fit with both thermal and non-thermal models; however, a thermally emitting shell seems much more likely. Assuming thermal emission, Bocchino et al. (2010) were able to use PWN–SNR evolutionary models to obtain an age consistent with the characteristic age derived from pulsar observations, obtain the proper SN–PWN radius ratio, and predict that the reverse shock has not encountered the PWN yet, consistent with other observations. The ages of 3C58 and G54.1+0.3 are ∼3000–5000 yr and ∼1800–3300 yr, respectively. G21.5−0.9 is thus unique in that it is much younger (∼290–1000 yr), and potentially has higher forward shock speed, both of which could lead to the detectable non-thermal shell. The Crab Nebula itself, of comparable age to G21.5−0.9 but much closer, still shows no non-thermal shell, presumably due to a very low interstellar medium density. The detection of this G21.5−0.9 non-thermal shell at quite high X-ray energy is thus likely due to its younger age compared to these other Crab-like SNRs.