Discovery of Radial Spectral Hardening in the Hot Bubble of Planetary Nebula BD+30°3639 with Median Energy Imaging

We apply an imaging analysis technique to study the spatial distribution of the X-ray spectral hardness across the hot bubble of planetary nebula BD+30°3639. Hot bubble emission is typically photon starved, thus limiting the methods for spatial–spectral analysis; however, this technique uses the statistics of photon energies across the nebula to identify spatial variations in the X-ray spectral hardness. Using the median energy value of the X-ray photons, we identified a rise in median energy values toward the projected edge of the nebula, which we refer to as radial spectral hardening. We explored the origin of this radial spectral hardening with an X-ray spectral analysis of distinct regions of high and low median energy values. Given the results of the X-ray spectral fits and the fact that the hot bubble is embedded within a young, dense, planetary nebula, we conclude that the radial spectral hardening is due to an increased column density at the projected nebular edge.


INTRODUCTION
A planetary nebula (PN) forms during the late evolutionary stages of low-to intermediate-mass star (1-8 M ) when a nascent fast stellar wind collides with the older, denser, and slower asymptotic giant branch (AGB) stellar wind (Kwok et al. 1978).According to this interacting stellar wind theory, material from the AGB wind is swept into the nebular shell and ionized by the hot central star.This wind-wind interaction generates a shocked region that reaches X-ray emitting temperatures (> 10 6 K) and is called the "hot bubble".Hot bubbles are expected to scale with the square of the fast wind velocity, however, Chandra and XMM-Newton observations indicate that plasma temperatures are much cooler than expected (e.g., Kastner et al. 2008;Ruiz et al. 2013).Understanding the origin of these cooler plasma temperatures is a unresolved area of research being explored with a variety of approaches (Soker & Kastner 2003;Akashi et al. 2007;Steffen et al. 2008;Heller et al. 2018;Toalá & Arthur 2016, 2018).
The Chandra X-ray Observatory unambiguously resolved hot bubble emission for the first time with observations of the PN BD+30 • 3639.The plasma temperature was found to be 3 × 10 6 K by studying the low energy resolution CCD X-ray spectrum (Kastner et al. 2000).High-resolution grating X-ray spectroscopy revealed that the hot bubble emission was carbon-rich and iron-deficient and provided evidence for a nonisothermal plasma with plasma temperatures ranging from 1.7 to 2.9 × 10 6 K (Yu et al. 2009).Additionally, hot bubbles are embedded within the optically-and infrared-bright nebular shells and variations of this overlying nebular material are believed to influence the detected X-ray emission (Kastner et al. 2002;Maness et al. 2003).The spatial and spectral decomposition of the X-ray emission from the hot bubble in NGC 7027 (Montez & Kastner 2018) qualitatively demonstrated the degree of influence that the surrounding nebular material can have on the observed hot bubble plasma.The spatial resolution of Chandra is sufficient for the mapping planetary nebulae hot bubbles, but most hot bubbles studied by Chandra are photon-starved (Kastner et al. 2012;Freeman et al. 2014).In this study, we explore a novel methodology designed to glean spatial information Montez across the hot bubble of BD+30 • 3639 using the photon statistics of the X-ray emission.
2. OBSERVATIONS 2.1.Chandra Observations BD+30 • 3639 has been observed for a total of nearly 400 ks with Chandra ACIS-S imaging and zeroth-order ACIS-S gratings observations.Although it is tempting to compile all of the observations, the time-dependent sensitivity due to build-up of a contaminant on the ACIS optical blocking filter (Marshall et al. 2004) would limit the utility of the methodology employed in this work.This is especially problematic because BD+30 • 3639 primarily emits soft photons ( 1 keV), where the most change in detector sensitivity has occurred.Instead, we focus only on the observations acquired in Cycle 9 (Ob-sIDs 9932 and 10821 with exposures of 38.1 and 39.1 ks and acquired on 2009-01-27 and2009-01-22, respectively;PI: Yu).Instrumental changes in sensitivity are negligible since these two observations occurred within five days of each other.
We analyzed the observations with CIAO (version 4.11 Fruscione et al. 2006) and the calibration database CALDB 4.8.3.The chandra repro CIAO script was used to reprocess the data with sub-pixel event repositioning (Li et al. 2003(Li et al. , 2004) ) and the merge obs CIAO script was used to merge the two observations.The 77 ks merged observation yields 9485 counts in the 0.15 to 1.5 keV energy range from a region consistent with the nebula (an ellipse with a semi-major and semi-minor radii of 3 and 2.5 , respectively).The average count density is ≈ 400 counts per square arcsecond or ≈ 100 counts per ACIS-S native pixel (0. 492 × 0. 492).
We have also extracted source and background spectra from ObsID 10821 with the CIAO script specextract following the analysis thread for extended sources1 .We modeled the emission using Astrophysical Plasma Emission Code (APEC) models with variable abundances (Smith et al. 2001) and the Tuebingen-Boulder ISM absorption model (Wilms et al. 2000).In §5 we provide further details on the spectral analysis.

Hubble Observations
We obtained the Hα λ656.3 nm narrow-band images of BD+30 • 3639 from the Hubble Space Telescope.The Hα narrow-band filter (F656N) image from Hubble proposal 11122 (PI: Balick) was acquired from the Hubble Legacy Archive2 .These observations were obtained on 2008-03-02, within a year of the Cycle 9 Chandra observations.The F656N images were flux scaled according to Dudziak & Walsh (1997), however, we did not correct for contamination from [NII] (F658N), nor the continuum from the F547M filter (not available).The optical image is dominated by Hα emission.We performed a cosmetic cleaning of hot pixels, saturation regions, and cosmic rays with negligible impact on our analysis.The central star has been masked for display purposes.

MEDIAN ENERGY IMAGE
We introduce an imaging product called the median energy image, which consists of an image wherein each pixel is generated from the statistical properties of photons that reside within the pixel boundaries.We generate a number of evenly distributed spatial bins based on the desired pixel size, w, and overall image size.Each photon is collected into its appropriate bin based on its physical coordinates and we determine the median value of the photon energies, Ẽ, the standard deviation of the photon energies, σ E , and the total number of photons, N , for the collection of photons in each pixel.
To estimate the variance of the median value of the photon energies, σ 2 Ẽ , we adopt the standard error on the median, where N is the number of photons in a given pixel.We note that this formula is intended for normal distributions and typical hot bubble plasma emission is not normally distributed.As an alternative, we considered a boot-strapping error estimation derived from thousands of draws of N distinct photons from the observed plasma distribution.We found that the standard error on the median provided realistic uncertainty estimates and boot-strapping typically resulted in overly optimistic error estimates (smaller errors).
The median values of the photon energies, hereafter referred to as the "median energy" values, Ẽ, median energy error values, σ Ẽ , and total count values, N , are used to populate the pixels in 2D arrays that form the median energy image, median energy error image, and the count image, respectively.In Figure 1, we depict the count and the median energy images.
We calculated the radius from the central star of BD+30 • 3639 to each pixel.These radii are color-coded in the left panel of Figure 2 in the count-median energy plane.In the right panel of Figure 2, we plot the radial distribution of the pixels in the median energy images color-coded by the number of X-ray counts in each pixel.Error bars for the median energy values are included in both panels based on σ Ẽ .The radial extent of the X-ray emission from BD+30 • 3639 is 3. 0 (see Figure 1).In Figure 3, we depict the Hubble Hα emission, the X-ray count, and the median energy images.In the top row, contours from the smoothed Hα image are overlaid on all the images, while in the bottom row the spectral extraction regions discussed in §5 are overlaid on the images.The central star in the Hα image is masked and filled in with the local average flux.The X-ray count and median energy images are filtered to only show those pixels with at least 5 counts.

RESULTS
The median energy imaging of BD+30 • 3639 reveals a distinct morphology compared to the morphology seen in the X-ray count image.The median energy image exhibits a partial ring-like structure of higher median energies ( Ẽ ∼ 0.8 − 0.9 keV) and a central region with lower median energies ( Ẽ ∼ 0.6 − 0.7 keV).Based on Figure 2, there is no apparent trend between the number of counts and the median energy values.Specifically, higher median energy values exist for both low and high count pixels.The partial ring-like structure of higher median energy is also apparent in the radial profile shown in Figure 2, where the rise in median energy begins at ≈ 1. 0 and rises till the edge of the nebula at ≈ 2. 5 to 3 .We refer to this apparent spectral hardening of the hot bubble emission as "radial spectral hardening" and now discuss its possible origins.

DISCUSSION
The median value is known as a robust statistic, meaning that the median value provides a good representation of the distribution even in the presence of strong outliers and non-optimal distributions (e.g., non-normal and/or low count distributions).The utility of median energy values have been shown in various studies.For example, in the Chandra Orion Ultradeep Project (Feigelson et al. 2005), the median energy values are used to estimate N H for point sources with a few to tens of counts.In Hong et al. (2004), the median and quartile energy val-  ues are used to classify spectral properties of X-ray point sources with limited statistics.In the lower-temperature regime, the Chandra Planetary Nebula Survey (Chan-PlaNS) used the median energy values and proxies for the column densities to characterize both point sources and extended sources of X-ray emission from PNe and their central stars (Kastner et al. 2012;Freeman et al. 2014).Our analysis presents the first attempt to use the median energy values as an imaging tool.
For well-behaved and well-characterized instrumental sensitivity and calibration, changes in the median energy values for a source modeled by an absorbed thermal plasma could be caused by variations in the elemental abundances, plasma temperature, and/or absorbing column.These three potential origins are not neces-

Note-Model descriptions:
Model I: oxygen (O) and neon (Ne) abundances allowed to vary in two regions independently.Column densities and temperatures are tied together (N H,R2 = N H,R1 and T X,R2 = T X,R1 , respectively) and allowed to vary.
Model II: Abundances fixed, column densities are tied together (N H,R2 = N H,R1 and allowed to vary, and temperatures are allowed to vary in two regions independently. Model III: Abundances fixed, temperatures are tied together (T X,R2 = T X,R1 ) and allowed to vary, and column densities are allowed to vary in two regions independently.Errors are derived from 90% confidence ranges.The abundances of carbon (C), nitrogen (N), magnesium (Mg), and iron (Fe) are set to those reported in Yu et al. (2009), 19.5, 0.17, 0.6, and 0.13, respectively, as are O, and Ne for Models II and III.See Figure 3 for depiction of regions and Figure 4 for plots of the best-fit spectra.
sarily mutually-exclusive and have a degree of interdependence.As the plasma temperature changes so does the ionization stage and thus the strength of specific emission lines in the X-ray spectrum.Emission lines at lower energy values are more sensitive to absorption than those at higher energy values, so photons from Cvi at ≈0.37 keV will be absorbed more readily than Neix at ≈0.92 keV.Additionally, absorption becomes less efficient as the plasma temperature rises since a hotter plasma emits fewer soft photons.We independently consider each of these three potential origins in the case of BD+30 • 3639 by modeling X-ray spectra from the high median energy and low median energy regions detailed in Figure 3.

Varying Elemental Abundances
The elemental abundances of hot bubbles in PNe is an open question.The central star often has a distinct chemical composition compared to the nebular chemical composition.When modeling the hot bubble emission, one can assume the chemical composition of the central star or the nebula.However, when there is sufficient signal to constrain hot bubble abundances, the chemical composition is often found to be discrepant with both the central star and nebular abundances (e.g., Toalá et al. 2019).Mixing between central star and nebular abundances are included in the theoretical efforts to understand hot bubble X-ray emission (Steffen et al. 2008;Heller et al. 2018;Toalá & Arthur 2016, 2018).In the X-ray emitting temperature ranges of hot bubbles, the abundances of carbon, nitrogen, oxygen, neon, and iron have the strongest influence on the spectral shape of the X-ray emission emission.At CCD spectral resolution, it is often difficult to obtain meaningful constraints on many of these elements, especially for carbon and iron (see Montez & Kastner 2018, for a detailed discussion).However, the high energy resolution grating spectrum of BD+30 • 3639, which resolves the major ionized lines responsible for the X-ray emission, provides strong constraints on the important elemental abundances (Yu et al. 2009).
The apparent radial spectral hardening in BD+30 • 3639 could be explained if the neon abundance increased radially outward (and/or if oxygen decreased radially inward).To test this scenario, we assumed a single absorbing column and plasma temperature (fit as free parameters) and allowed the abundances for oxygen and neon independently vary in the two regions (Model I in Table 1).As expected, the neon abundance was higher in the high-median energy region and lower in the low-median energy region, while oxygen trended towards a value that is a factor of two lower than that reported by Yu et al. (2009).In Figure 4 we plot the spectra, best-fit model, and residuals for the low and high median energy regions.The fit to the spectrum of the low-median energy region is poor for energies 0.5 keV.We also attempted spectral fits to the carbon, nitrogen, and iron abundances but the abundances and plasma properties became unconstrained, suggesting that the spectrum is relatively insensitive to these parameters.Instead, we fixed the abundances of these elements to the single-temperature model in Yu et al. (2009) 3 .

Varying the Plasma Temperatures
Assuming a constant column density and a homogeneous chemical composition in the hot bubble, the radial spectral hardening could suggest a radially-increasing temperature.To explore this scenario, we fit a single column density to the low-and high-median energy regions with abundance values fixed to that from Yu et al. (2009) and allowed two independent temperatures, one for each of the two regions (Model II in Table 1).As expected, the high-median energy region results in a higher plasma temperature (T X = 2.4 ± 0.2 MK) compared to the low-median energy region (1.8 ± 0.1 MK).For comparison, the spectral fits of the high-resolution grating X-ray spectroscopy of BD+30 • 3639 led Yu et al. (2009) to argue for a two-component model with temperatures of 1.7 and 2.9 × 10 6 K.The spectral fit for this scenario is well-behaved and its goodness of fit is acceptable and better than that of Model I.It is worthwhile to note that heat conduction and nebular mixing studies predict radially-decreasing temperatures, in contrast to the prediction of Model II of radially-increasing temperatures.

Varying the Absorbing Columns
For the cooler plasma temperatures measured from hot bubbles, intervening material can be efficient at absorbing the soft (< 1 keV) X-ray photons.Assuming the hot bubble is an isothermal plasma with homogeneous chemical composition, we performed spectral fitting on the two median energy regions assuming a single plasma temperature and two independent column densities (see Model III in Table 1).Similar to Model II, the spectral fit for this scenario is well-behaved and its goodness of fit is acceptable and better than that of Model I.The best fit parameters suggest that the column density in the high-median energy region is a factor of two higher than that in the low-median energy region.
Since we view the hot bubble through the nebula, line of sights near the projected nebular edge will pass through more nebular material than those in the center, so higher column densities of nebular material at projected edge of nebula should be expected.The radius of the nebula projected onto the sky is 0.02 pc (Kastner et al. 2012) and the thickness of the nebular shell is likely an order of magnitude smaller, so the ionized thin shell of a spherically-symmetric nebula is unlikely to provide sufficient absorbing column to explain the radial spectral hardening.However, there is compelling evidence that the three-dimensional nebular morphology of BD+30 • 3639 is an elongated ellipsoid with the long axis tilted in our direction (Masson 1989;Bryce & Mellema 1999;Lee & Kwok 2005;Freeman & Kastner 2016).Such a morphology implies that we view the hot bubble through more nebular material than in a spherically-symmetric case.Additionally, dust and molecular shells (Bernard et al. 1994) and a large scattering halo comprised of ionized and neutral hydrogen gas (Harrington et al. 1997) -also see Figure 3 -exist around BD+30 • 3639, indicating the presence of additional material that could absorb the soft X-ray emission.Harrington et al. (1997) presented selective and total extinction maps of BD+30 • 3639 using the ratio of optical and radio imaging observations.These maps provide clear evidence for increased absorption at the projected edge of the nebula.The extinction behavior is consistent with the radial spectral hardening behavior seen in the median energy image, specifically, the extinction decreases as the radial distance from the central star decreases.Harrington et al. (1997) argued that a majority of the absorption is caused by the nebula itself and not the interstellar medium.The extinction maps of Harrington et al. (1997) and Lee & Kwok (2005) show that the extinction is higher on the western side of the nebula than the eastern side, which is consistent with the median energy image, which is harder on the western side versus the eastern side.This asymmetry suggests that our spectral fit in the high median energy regions is weighted towards the eastern side of the hot bubble emission where more counts are detected.As a result, the range of column densities across the hot bubble are probably larger than that indicated by the best-fit values (1.53 × 10 21 cm −2 ).
The extinction studies (Harrington et al. 1997;Lee & Kwok 2005) have argued that dust is present in the ionized nebular gas in addition to the dust and molecular shells.Comparing the morphology of the median energy image to that of the distribution of polycyclic aromatic hydrocarbons (PAHs) in the 3.3 µm image presented in Bernard et al. (1994) reveals a strong correspondence between the PAH emission and the radial spectral hardening morphology.Specifically, both images show partial ring-like morphologies with a break in similar northern locations.The 3.3 µm emitting region and the apparent radial spectral hardening of the X-ray plasma may be casually related.There are three models fit to the spectra.These three models are described in §5 and Table 1.The three smaller panels in each row show the residuals produced by each model fit.

CONCLUSIONS
Median energy imaging of the hot bubble in PN BD+30 • 3639 provides a new perspective on the X-ray photon distribution of this well-studied planetary nebula.We find that the median energy image exhibits a unique morphology that is distinct from that shown in the count image.Specifically, a ring-like structure of higher median energy values surrounds a central region of lower median energy values that we call radial spectral hardening.The radial profile confirms the increase to higher median energy values towards the nebular rim.
We explored a variety of origins to explain the observed radial spectral hardening using X-ray spectra that isolated the high-median energy and low-median energy regions.Three scenarios focused on the abundances, temperatures, and column densities in each region.We studied the behavior of the X-ray spectra under these scenarios and found that all three scenarios produced acceptable spectral fits, suggesting a degeneracy in the origin of the radial spectral hardening.Combinations of these scenarios, which were not considered here, could also potentially explain the radial spectral hardening.
Given that the hot bubble is embedded within a young, dense planetary nebula with an elongated nebular shell aligned with our line of sight, we suggest that the radial spectral hardening of the hot bubble X-ray emission in BD+30 • 3639 is due to an increased column density at the projected nebular edge.Similar to the case of NGC 7027 (Montez & Kastner 2018), in BD+30 • 3639 the impact of overlying nebular material on the hot bubble characteristics is important.In Xray spectral analysis with existing CCD-resolution spectroscopy, assuming a constant column density across the extended emission introduces a bias that can obscure or confuse the hot bubble physical parameters.The median energy imaging technique provides a promising avenue for exploring such variations in relatively low count sources.

Montez
The median energy imaging for planetary nebulae was originally developed by the author while attending the Astro Hack Week unconference in 2017 at the University of Washington and the author expresses gratitude to the organizers and participants of that workshop.The author is grateful to Annie Blackwell, Roel Olvera, Daniel Castro, and Regina Jorgenson who have worked with the author on applying the median energy imaging technique to supernova remnants in projects for the NSF-funded Maria Mitchell Observatory Research Experiences for Undergraduates program.This research has made use of data obtained from the Chandra Data Archive and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa.This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2018, 2013).RMJ acknowledges support from NASA contract NAS8-03060.In addition to the median energy image, we also considered images derived from the mean and standard deviation of the photon energy distribution in each pixel.From the mean image, I µ E , the standard deviation image, I σ E , and median energy image, I Ẽ , we calculated an image of the non-parametric skew, I NP , using The four images are shown in Figure 5.The mean energy image shares similar morphological properties as the median energy image, however, the mean energy image has less contrast between the lower and higher mean energy region compared to the median energy image.Note that the colorbar scales in the mean and median energy images in Figure 5 are equal to aid in the comparison of their contrast.This indicates the robustness of the median value as it is less susceptible to outliers.The values of the standard deviation peak at ∼ 0.23 keV with the majority failing below 0.27 keV.The standard deviation values appear uncorrelated with the various morphologies of the counts, mean, and median energy images.
A strong imaging contrast arises in the non-parametric skew image with the central region showing a positive skew (mean is greater than the median) while the outer region showing a zero to slightly negative skew (mean is equal or less than median).We plot the pixel values in Figure 6.There is no clear correlation between non-parametric skew and the number of counts.There is a clear indication that the central region (small radii) has a predominately positive skew, whereas above 1. 4 the non-parametric skew has a broader range of values.While the contrast of the non-parametric skew is high, perhaps more than that seen in the median energy image, it is difficult to interpret the skew since the energy distribution, i.e., the X-ray spectrum, is multi-modal.Nevertheless, the non-parametric skew supports the notion that the distribution of photons changes across the nebula.

Figure 1 .
Figure 1.Count (top row) and median energy (bottom row) images of BD+30 • 3639 derived from the Cycle 9 observations obtained by Chandra.All images have bins comparable to the native ACIS-S pixel size (0. 492 × 0. 492).The images in the left column depict the full X-ray observation, while the images in the right column are filtered to only show pixels with 5 or more counts.The median energy image is derived from the energy distributions in each count image pixel.

Figure 2 .
Figure2.Pixel values from the count and median energy images.In the left panel, the median energy pixel values are plotted as a function of the number of X-ray counts and the symbols are color-coded by the radial distance from the central star.In the right panel, the radial profile of the median energy pixels is plotted and the symbols are color-coded by the number of X-ray counts.In each panel, the error bars are standard errors on median energy values as described in the text and we only plot pixels with at least 5 counts.

Figure 3 .
Figure 3. From left to right we depict the Hubble Hα image, the Chandra X-ray count image, and the median energy image.In the top row we overlay contours from the smoothed Hubble image and in the bottom row we overlay spectral extraction regions from the high-median energy region (outer wedge-shaped regions) and a low-median energy region (inner circular region).

Figure 4 .
Figure4.X-ray spectra and spectral fits from the high (left) and low (right) median energy regions.There are three models fit to the spectra.These three models are described in §5 and Table1.The three smaller panels in each row show the residuals produced by each model fit.

Figure 5 .
Figure 5. Statistical imaging of BD+30 • 3639.The mean energy image, standard deviation image, median energy image, and non-parametric skew image are shown in the four panels and distinguished by the colorbar label.

Figure 6 .
Figure 6.Pixel values from the statistical images.In the left panel, each non-zero pixel from the count and non-parametric skew images are plotted against each other and color-coded by the radial distance from the pixel centroid to the central star.In the right panel, the radial profile of the non-parametric skew pixels is plotted and color-coded by the median energy values.

Table 1 .
Spectral fits for median energy-defined regions.