PROBING X-RAY ABSORPTION AND OPTICAL EXTINCTION IN THE INTERSTELLAR MEDIUM USING CHANDRA OBSERVATIONS OF SUPERNOVA REMNANTS

The linear relationship between optical extinction (${A}_{{\rm{V}}}$) and hydrogen column density (N_H) has long been observed and utilized to estimate the X-ray or optical brightness for new sources or to fit the broadband spectrum of X-ray sources. It is also used to obtain distance estimates for X-ray sources using their measured column densities (see, e.g., Durant & van Kerkwijk 2006; Güver et al. 2010; Ratti et al. 2010; Nielsen et al. 2012; Soria et al. 2012).

The photoelectric absorption by interstellar material causes rapid attenuation of soft X-rays in the spectra of Galactic sources. Measuring the extent of this attenuation yields information about the total column density along the line of sight to the source. Although this is commonly expressed in terms of the equivalent hydrogen column density N_H in the soft X-ray band (0.1–10 keV), it is predominantly caused by abundant heavier elements such as O, Ne, Si, Mg, and Fe. Optical extinction is caused by grains of the same elements. Because dust tends to follow the metal distribution in the interstellar medium (ISM), when averaging over many different lines of sight and over distances that are larger than the clumping scale of the ISM, it is reasonable to assume an approximate linear relationship between between A_V and N_H.

There have been a number of different studies that have sought to accurately determine the relation between optical extinction and the hydrogen column density empirically. The methods vary across these studies, and the results show discrepancies greater than the statistical errors for each. Reina & Tarenghi (1973) used X-ray binaries and extended sources and found a relation of N_H = 1.85 × 10²¹ × A_V cm⁻² (hereafter N_H is in units of cm⁻² and A_V is in magnitudes); Gorenstein (1975) used supernova remnants (SNRs) to find N_H = (2.22 ± 0.14) × 10²¹ × A_V; while Predehl & Schmitt (1995) used a combination of ROSAT point sources and SNRs and measured N_H = (1.79 ± 0.03) × 10²¹ × A_V. Recently, Güver & Özel (2009) collected a sample of 22 SNRs, for which the hydrogen column density and optical extinction were previously measured, and found N_H = (2.21 ± 0.09) × 10²¹ × A_V. This had the advantage of using the high-quality data from modern X-ray telescopes such as Chandraand XMM-Newton and focusing on sources with little to no intrinsic absorption,⁴ but had the disadvantage of being unable to account for systematic errors that may vary for each published value. For example, when using values from the literature, there is no way to account for the variety of choices that are made during the data processing pipeline. The selected regions or the particular emission models may be interesting for the objectives of a particular study, but less ideal for the determination of the hydrogen column density and the comparison to the optical extinction measurements to determine the slope of the relation.

In this study, we take advantage of the wealth of SNR data available in the Chandra archive to investigate and quantify the systematic errors present in the determination of the hydrogen column density from the analysis of the X-ray data. We analyze all of the observations from the Chandra archive using standardized procedures so that we can quantify the systematic errors on each measurement. Using only observations from the Chandra archive also ensures a completely uniform data set: all observations were performed using the ACIS detector, and each had spectra generated and analyzed using the same treatment for background subtraction and model fitting routines using the spectral analysis software xspec (version 12.8.1; Arnaud 1996, with NEIVERS 1.1). This consistent treatment of a uniform data set gives us the opportunity to quantify the existing systematic errors in a way that had not been possible before.

The uniformity of our data set also allows us to explore the uncertainties associated with fitting spectra with a variety of models. SNRs can be home to a wide range of plasma properties due to the wide range of ages, host environments, and composition of the sampled gas. In general, however, SNRs exhibit continuum emission driven by thermal bremsstrahlung accompanied by emission lines produced by ejecta material or swept up ISM. In the 0.5–5.0 keV range, where we predominantly perform the spectral fits, these features manifest most visibly in magnesium, silicon, and sulfur lines (though features due to iron, oxygen, and neon are also possible within this range). Non-thermal features can also be present due to particle acceleration in SNR shocks. Due to the wide range in spectral features and plasma conditions, there also exist a range of models that can be used to fit X-ray spectra, which will be detailed in Section 2.

Making use of the high spatial resolution and large number of counts in the Chandra archival data, we also investigate the systematic uncertainties associated with differing lines of sight toward larger remnants. We assume that all absorption is due to interstellar gas and dust along a line of sight, with no intrinsic absorption from the optically thin remnant. This assumption allows there to be differing values for N_H for different regions in a given SNR due to real differences in gas along a line of sight, which is supported by deep observations of large remnants such as W49B (Sasaki et al. 2013) and CTB 109 (also known as SNR G109.1−1.0; Keohane et al. 2007), for which many regions can be fit. When the archival observations have sufficient duration to have multiple high-count regions, we are able to fit multiple regions to understand how the value of N_H varies over the remnant. By fitting multiple regions with a range of models, we can determine the magnitude of the total systematic error on N_H measurements in the direction of the SNR. We note that any dust that is intrinsic to the SNR contributes negligible extinction. The integrated column density of swept up dust is small, and that of dust condensed from the stellar ejecta is smaller still. Moreover, SNR shocks are very efficient dust destroyers (Temim et al. 2015), making any contribution to the overall extinction exceedingly small.

Finally, and possibly most importantly, many of the previous analyses of Chandra and other SNR data used the solar abundances of Anders & Grevesse (1989) derived from observations of the Sun and meteorites. Wilms et al. (2000) showed that the abundances in the ISM can differ from solar values, which can affect the fit value of N_H. It is also important to consistently use an ISM absorption model that features improved calculations of the absorption cross sections when seeking to improve systematic errors in the hydrogen column density determination. We use the wilms abundance table and tbabs absorption model within xspec in order to accurately model the interstellar absorption to measure the hydrogen column density. The tbabs model in xspec adopts the cross sections from Verner & Yakovlev (1995). This standard that we adopt for the entire sample leads to larger (∼30%) column densities than those reported in the earlier studies.

In Section 2, we discuss the variety of models used to fit our data, as well as the results of simulated data sets. Sections 3 and 4 detail the processes we used to process the data, as well as our approach to identifying systematic errors and determining total uncertainties. Finally, in Section 5 we use our new data set to derive the slope of the N_H − A_V relation.

The models we used to fit each thermal region (raymond, mekal, nei, and pshock) were chosen to span the range of complexity one might expect from a region of plasma in an SNR. The raymond model (Raymond & Smith 1977) is the simplest, modeling a hot, diffuse gas in collisional equilibrium. Similarly, the mekal model also calculates spectra of a plasma in collisional equilibrium, but with improved handling of the Fe-L complex (Kaastra et al. 1996). The nei and pshock models (Borkowski et al. 2001) are more complicated, providing spectra for a non-equilibrium ionization (NEI) collisional plasma (nei), and constant temperature NEI plasma heated by a plane-parallel shock (pshock). For non-thermal regions, the models powerlaw and/or srcut were used. Both models produce a spectrum we would expect from synchrotron emission from a power-law distribution of shocked electrons interacting with an SNR's magnetic field, with an additional exponential cutoff of the distribution for the srcut model (Reynolds 1998; Reynolds & Keohane 1999).

As discussed previously, differing lines of sight can sample different paths through the ISM, so there is a chance for real variations in N_H measurements across an SNR. There can also be small systematic variations introduced by fitting the regions with the models discussed above. Each region may be composed of multiple plasma components, so fitting to a single plasma model may introduce a systematic error or a bias. To show this, we performed fits to simulated spectra produced with a range of N_H. We produced the simulated spectra using the sedov model in xspec (version 12.8.1 using NEIVERS 1.1), which models the total emission from an SNR undergoing adiabatic expansion (Borkowski et al. 2001). The complete set of simulated spectra cover a grid of inputs for shock temperature (T_s = 0.3–0.9 keV), electron temperature (T_es = 0.15–0.45 keV), ionization age (${\tau }_{0}=1-10\times {10}^{11}$ s cm⁻³), and N_H (0.5–1.0 × 10²² cm⁻²). For each parameter set, we simulated 200 spectra and then binned these to have at least 50 counts per bin. We fit each spectrum with a pshock model, using the tbabs absorption model with the wilms abundance table. We present the average of the best-fit values for each parameter set in Table 1, along with the standard deviation of all fit values.

The results in Table 1 show that while a single plasma model approximately recovers the assumed N_H value from a complex plasma, it can introduce significant statistical errors and suggests a possible bias. In particular, the N_H values we measure from the simulated data are close to but often lower than the assumed values. This, in fact, is not surprising: fitting a single-temperature model to a simulation that is meant to represent the entire remnant (i.e., a Sedov spectrum which integrates over a range of temperatures) cannot adequately describe the spectrum or yield a very accurate value of the hydrogen column density.

To minimize these effects in the analysis of the actual data, we take a three-pronged approach. First, we extract spectra from as small regions of the remnant as possible to avoid creating complicated, multi-temperature regions with multiple plasma components. Second, we perform fits with thermal, non-thermal, and mixed thermal/non-thermal models to capture the spectral characteristics of the regions correctly. Third, we use Bayesian statistical tools (discussed in Section 4.2) to combine measurements from different spatial regions as well as from different spectral model results to assess any systematic uncertainties in the N_H measurement for a given remnant. It is important to note, however, that the simulated spectra results on their own are not the sole motivation for the Bayesian analysis; if we believed that the simulated results were immediately comparable to the Chandra data results, then we could use the simulated results to establish the magnitude of a bias term. In Section 5, we will use the total uncertainties we determine for the N_H toward each remnant to more accurately constrain the N_H/A_V relationship.

3.1. Pre-processing

3.2. Region Selection and Data Processing

We used the new level = 2 event file to select regions for spectral extraction. We selected regions using ds9 and chose areas that had ≥10,000 counts, as well as sampling a variety of lines of sight across the remnant. In some remnants, it was preferable or necessary to select regions with non-thermal emission. We processed these using procedures identical to the thermal regions, but fit them with non-thermal models. In general, regions were selected using the morphology of the remnant as a guide, such that the regions would contain one variety of plasma (the importance of which is highlighted in Section 2). We extracted a spectrum from each region and then grouped it to a minimum of 25 counts per bin.

As a typical example of the regions selected for a remnant, we show SNR G109.1−1.0 in Figure 1, with a number of the selected regions marked with circles and the regions selected for local background shown with a dashed circle and rectangle. In the full analysis of G109.1−1.0, we used a total of 11 regions. However, for clarity and simplicity, we chose to highlight only six in this figure. We show the resultant binned spectra from two of these regions in Figure 2, again as two representative examples of the types of spectra we encountered in the analyses. Note that small contributions from model components (such as the power-law component in Region 3 of Figure 2) may appear insignificant, but are important for producing acceptable fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We fit a variety of models (discussed in Section 2) to each grouped spectrum. We fit each model with default ISM abundances, but it was often necessary to allow some of the plasma abundances to vary in order to properly fit regions with ejecta-enriched plasma. The number and variety of the free abundances were allowed to change from model to model for a given region, and only the plasma models' abundances were allowed to vary. For non-thermal regions, this process was the same, but simplified by the fact that there were no free abundance parameters possible for those models.

Despite our best efforts to select regions that would be composed of only one variety of plasma (and thus could be fit using one model), sometimes this was infeasible (for example, if the regions had to be large in order to contain an adequate number of counts). In such cases, we needed to add other components to the models. Most commonly, this required adding a non-thermal model to the thermal model to compensate for regions near non-thermal filaments or a central pulsar wind nebula. However, in some cases we added Gaussian features to fit spectral lines that were either not properly fit or non-existent in the thermal model, or to account for ejecta mixing into a non-thermal region, causing emission lines to appear on top of a power-law spectrum. For the most part, the data were good enough that we could find regions that allowed for a single model to produce an acceptable fit, but when included, these additional components did not make statistically significant changes in the best-fit value of N_H.

We determined if a spectral fit was good both by visually checking the final fit, as well as by using the χ²/ν fit statistic. We show in Figure 3 an example of an acceptable fit compared to a poor fit. We did not place a hard upper limit on the fit statistic because some large χ²/ν values were dominated by a single feature that was not captured by the model but did not affect the inferred N_H. Nevertheless, nearly all of the fits used in our final analysis had ${\chi }^{2}/\nu \leqslant 2$.

Zoom In Zoom Out Reset image size

4.1. Results of The Spectral Analysis

For each SNR region, we included all models that produced acceptable fits. As a typical example, we show in Table 2 the complete fit results for six of the regions of G109.1−1.0. It is evident from this table that vnei and vpshock models with thawed abundances were commonly the best models for remnants with thermal regions, with a powerlaw component added in when necessary. The raymond and mekal models ranged from fairly successful to very poor, which presumably reflects significant NEI contributions in the data. For SNRs with completely non-thermal regions, the powerlaw model was the default model with consistent success. (Two of the remnants, G120.1+1.4 and G04.5+6.8, were better fit by an srcut model.) All errors presented on the fit N_H values in Table 2 are 90% confidence ranges. The complete fit results for all remnants used in our analysis have been released as a Zenodo data set (http://zenodo.org/record/17183#.V22YjpMrKqA).

Table 2.  Example Fit Values for SNR G109.1−0.0

Region	Model	Thawed Elements	Fit NH	χ2/ν
			(1022 cm2)
1 (NE)	vnei + powerlaw	Mg, Si	${1.07}_{-0.07}^{+0.07}$	146.94/97
	vpshock + powerlaw	Mg, Si, S	${1.05}_{-0.07}^{+0.08}$	143.81/96
3 (N)	nei + powerlaw	None	${0.95}_{-0.09}^{+0.09}$	132.71/111
	vnei + powerlaw	Mg, Si	${0.93}_{-0.16}^{+0.10}$	130.55/109
	pshock + powerlaw	None	${0.95}_{-0.04}^{+0.09}$	131.83/111
	vpshock + powerlaw	Mg, Si	${0.89}_{-0.06}^{+0.09}$	126.88/109
4 (E)	vnei + powerlaw	Mg, Si	${0.90}_{-0.07}^{+0.07}$	160.49/124
	vpshock + powerlaw	Mg, Si	${0.93}_{-0.07}^{+0.08}$	153.14/124
6 (Center)	vnei	Mg, Si, S	${0.86}_{-0.07}^{+0.08}$	138.73/110
	vpshock	Mg, Si, S	${0.71}_{-0.05}^{+0.08}$	130.71/110
9 (SE Inner)	vnei	Mg, Si, S	${1.24}_{-0.07}^{+0.06}$	139.21/117
	vpshock	Mg, Si, S	${1.23}_{-0.05}^{+0.05}$	137.94/117
11 (SE Outer)	nei	None	${0.87}_{-0.06}^{+0.06}$	142.52/111
	vnei	Mg, Si	${0.88}_{-0.07}^{+0.08}$	142.15/109
	pshock	None	${0.88}_{-0.06}^{+0.08}$	130.27/111
	vpshock	Mg, Si	${0.89}_{-0.08}^{+0.08}$	129.88/109

Download table as:  ASCIITypeset image

4.2. Bayesian Analysis to Determine N_H and Its Uncertainty

Each of the fits within a region have formal (and often asymmetric) uncertainties, as we showed in Table 2 for SNR G109.1−1.0. From these fits, the dispersion in the measurements arising from spatial sampling or from different continuum models can be calculated. For a number of remnants, this contribution is comparable to the formal uncertainties. In others, however, the difference between the various measurements of N_H is significantly larger than the formal uncertainties, pointing to systematic uncertainties originating from spatial sampling, from the choice of model, or both. For the purposes of our analysis, we do not distinguish between the variance of the ISM along different lines of sight and the systematic error introduced by the incomplete plasma models and we combine these different types of systematic error into a single term (see, e.g., Sinervo 2003). Our aim in this section is to determine the most likely N_H value for each remnant, as well as a measure of the combined formal and systematic uncertainties in this quantity. We accomplish this by finding the parameters of the underlying N_H distribution for each remnant that is consistent with our sample of measurements. We use this distribution to find the most likely value of the hydrogen column density and its uncertainty.

We start with a parametric form of the N_H distribution and use the set of measurements from the individually fitted regions in order to estimate its parameters. We take the assumed underlying distribution to be a Gaussian

$$\begin{eqnarray}&&P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}})=C\mathrm{exp}\left[-\displaystyle \frac{{({N}_{{\rm{H}}}-{N}_{\mathrm{Hc}})}^{2}}{2{\sigma }^{2}}\right],\end{eqnarray} \tag{ 1 }$$

with a mean N_Hc and a standard deviation σ that can be different for each remnant. In this and the following expressions, C is a proper normalization constant such that

$$\begin{eqnarray}&&{\int }_{0}^{\infty }P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}}){{dN}}_{{\rm{H}}}=1.\end{eqnarray} \tag{ 2 }$$

We also need to model the individual measurement uncertainties ${P}_{i}({N}_{{\rm{H}}})$, where i represents a particular SNR region/model combination that yields a single measurement of N_H for that region. In general, the χ² surface for each measurement is not simple or symmetric around the minimum, leading to asymmetric formal errors in these individual measurements. However, around the minimum, it is accurate to represent the likelihood using two half Gaussians with different dispersions that smoothly connect at the most likely column density ${N}_{{\rm{H}}0,i}$ for each measurement; i.e.,

$$\begin{eqnarray}&&{P}_{i}({\rm{data}}| {N}_{{\rm{H}}})\\ \quad =\;\left\{\begin{array}{ll}{C}_{i}\mathrm{exp}\left[-\displaystyle \frac{{({N}_{{\rm{H}}}-{N}_{{\rm{H}}0,i})}^{2}}{2{\sigma }_{-,{\rm{NH,}}i}^{2}}\right], & {N}_{{\rm{H}}}\lt {N}_{{\rm{H}}0,i}\\ {C}_{i}\mathrm{exp}\left[-\displaystyle \frac{{({N}_{{\rm{H}}}-{N}_{{\rm{H}}0,i})}^{2}}{2{\sigma }_{+,{\rm{NH,}}i}^{2}}\right], & {N}_{{\rm{H}}}\gt {N}_{{\rm{H}}0,i}\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where "data" stands for the most likely column density ${N}_{{\rm{H}}0,i}$ and the two associated uncertainties ${\sigma }_{-,{N}_{{\rm{H}}},i}$ and ${\sigma }_{+,{N}_{{\rm{H}}},i}$ for the ith SNR region/model combination.

We want to calculate the quantity $P(\sigma ,{N}_{\mathrm{Hc}}| {\rm{data}})$, which measures the posterior likelihood of the parameters of the N_H distribution, given the observations. Using Bayes' theorem, we can write this as

$$\begin{eqnarray}&&P(\sigma ,{N}_{\mathrm{Hc}}| {\rm{data}})={C}_{2}P({\rm{data}}| \sigma ,{N}_{\mathrm{Hc}})P(\sigma )P({N}_{\mathrm{Hc}}),\end{eqnarray} \tag{ 4 }$$

where C₂ is a normalization constant and P(σ) and $P({N}_{\mathrm{Hc}})$ are the priors over the values of the Gaussian dispersion, σ, and the peak of the N_H distribution, N_Hc. Here, "data" stands for the ensemble of the ${N}_{{\rm{H}}0,i},{\sigma }_{-,{N}_{{\rm{H}}},i}$ and ${\sigma }_{+,{N}_{{\rm{H}}},i}$ values for a particular SNR. We take a flat prior over the Gaussian dispersion σ between σ_min that is equal to 0.1 times the smallest formal uncertainty obtained from a spectral fit for each remnant and σ_max that is equal to 10 times the largest difference between two N_H measurements for each remnant:

$$\begin{eqnarray}P(\sigma )=\left\{\begin{array}{ll}0, & \sigma \leqslant {\sigma }_{\mathrm{min}}\\ \displaystyle \frac{1}{{\sigma }_{\mathrm{max}}-{\sigma }_{\mathrm{min}}}, & {\sigma }_{\mathrm{min}}\lt \sigma \lt {\sigma }_{\mathrm{max}}\\ 0, & \sigma \gt {\sigma }_{\mathrm{max}}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

Similarly, we take a flat prior over the centroid of the N_H distribution N_Hc that spans the range from 0.1 times the smallest N_H measurement to 10 times the largest N_H measurement per source. These limits ensure that the particular minimum and maximum values of the prior distribution do not affect the results.

In Equation (4), the quantity $P({\rm{data}}| \sigma ,{N}_{\mathrm{Hc}})$ measures the likelihood that we will make a particular set of measurements for the column density given the values of the parameters of the column density distribution. We need to estimate this quantity, given the measurement likelihoods given in Equation (3). We will assume that each measurement is independent, so that

$$\begin{eqnarray}&&P({\rm{data}}| \sigma ,{N}_{\mathrm{Hc}})\\ &&\quad =\displaystyle \prod _{i}\displaystyle \int {{dN}}_{{\rm{H}}}{P}_{i}({\rm{data}}| {N}_{H})P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}}).\end{eqnarray} \tag{ 6 }$$

Combining this last equation with Equation (4), we obtain the posterior likelihood

$$\begin{eqnarray}&&P(\sigma ,{N}_{\mathrm{Hc}}| {\rm{data}})={CP}(\sigma )P({N}_{\mathrm{Hc}})\\ &&\quad \times \displaystyle \prod _{i}\displaystyle \int {{dN}}_{{\rm{H}}}{P}_{i}({\rm{data}}| {N}_{{\rm{H}}})P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}}),\end{eqnarray} \tag{ 7 }$$

where C is the overall normalization constant. We can use Equation (7) to determine the parameters σ and N_Hc given the individual region/model fits for each remnant. The dispersion (σ) is then a measure of the systematic uncertainty associated with the remnant. If a remnant has a significant systematic uncertainty, that contribution is considered when we define the most likely value of the N_H for the SNR and its uncertainty.

In Figure 4, we show as examples the Bayesian credible regions over the parameters of the N_H distribution that correspond to the measurements for G109.1−1.0 and G4.5+6.8. We show with a dot the most likely values of the centroid and dispersion; i.e., the peak of the posterior likelihood given by Equation (7). In the case of G4.5+6.8, the best-fit value of the dispersion is very close to zero and this parameter is consistent with being zero within the p = 0.05 credible region, indicating a negligible level of systematic uncertainty. On the other hand, for remnants such as G109.1−1.0, the best-fit value of the dispersion is not consistent with being zero, indicating the presence of systematic uncertainties arising from spatial sampling across the remnant or the differences among spectral models.

Zoom In Zoom Out Reset image size

The previous Bayesian analysis to quantify the systematic error term yields two categories of remnants. In the first category the systematic dispersion (σ) is consistent with zero, so we assign a zero systematic uncertainty, i.e., set $P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}})=\delta ({N}_{{\rm{H}}}-{N}_{\mathrm{Hc}})$ in Equation (6). This allows us to find a properly weighted average of the N_H measurements using their formal uncertainties and report the best-fit value and its formal uncertainty in Table 3. In the second category of remnants, we compute the posterior likelihood over their N_H by weighing each Gaussian distribution with a given centroid N_Hc and dispersion σ with the likelihood calculated in Equation (7) that those pair of parameters represent the observed N_H values:

$$\begin{eqnarray}&&P({N}_{{\rm{H}}})\\ &&\quad =\displaystyle \int \displaystyle \int P({N}_{{\rm{H}}};\sigma ,{N}_{\mathrm{Hc}})P(\sigma ,{N}_{\mathrm{Hc}}| {\rm{data}}){{dN}}_{\mathrm{Hc}}d\sigma .\end{eqnarray} \tag{ 8 }$$

We use this distribution to infer the best-fit values of N_H and its uncertainty for each source and report this in Table 3.

Table 3.  Measured Column Densities for All SNRs

SNR	Number of	NH	Error
	Regions/Models	(1022 cm2)
G0.0+0.0	3	9.86	0.34
G04.5+6.8	4	0.54	0.011
G06.4−0.1	22	0.63	0.19
G53.6−2.2	2	0.67	0.03
G54.1+0.3	1	2.55	0.04
G69.0+2.7	1	0.45	0.02
G109.1−1.0	30	0.90	0.14
G111.7−2.1	8	1.17	0.16
G116.9+0.2	2	0.92	0.07
G119.5+10.2	1	0.38	0.11
G120.1+1.4	14	0.80	0.10
G130.7+3.1	6	0.54	0.01
G184.6−5.8	1	0.30	0.02
G260.4−3.4	3	0.35	0.15
G263.2−3.3	1	0.03	0.01
G327.6+14.6	2	0.16	0.01
G332.4−0.4	9	0.87	0.35

Download table as:  ASCIITypeset image

We presented in the previous section the hydrogen column densities toward a large sample of SNRs that we measured using the Chandra data archive. We now combine these N_H measurements with the optical extinction measurements presented in Güver & Özel (2009) for our sample of SNRs to determine a relation between these two quantities that accounts for systematic errors. The optical extinction toward the remnants in our sample was determined through several different methods, which we list in Table 4. Most of these methods involve measuring the intensity ratio of two emission lines, for which the intrinsic ratio is known. Comparing the observed intensity ratio to the intrinsic ratio yields a reddening, which is then converted into a measurement of the optical extinction. These methods are challenging, in general, because a high signal-to-noise ratio spectrum is necessary to obtain the intensity ratios.

Several emission line pairs are frequently used for this purpose. The H_α (6563 Å) to H_β (4861 Å) line ratio, referred to as the Balmer decrement method, is one of the most well known and reliable ones among these pairs. Other emission line ratios from the SII multiplet (Miller 1968) and FE[II] IR transitions (Oliva et al. 1989) are also used in a similar way. For all of these pairs of lines, the ratio depends only very weakly on the temperature and density of the emitting plasma, leading to minimal uncertainties in the calculation of the ratio (Osterbrock 1989; Lequeux 2005).

Another common method is to use stars near an SNR to estimate the extinction toward that remnant. This method relies on having an existing estimate of the distance to the remnant, and then identifying stars at similar distances. If suitable stars can be observed, the extinction measurements of those stars can be applied to the remnant (Ruiz-Lapuente 2004; Koo et al. 2008).

In Figure 5, we plot the hydrogen column density measurements and their uncertainties presented in Table 3 against the measurements of the optical extinction summarized in Table 4. Following Güver & Özel (2009), we assign a 15% error to the optical extinction measurements for which uncertainties have not been reported (denoted by "-" in Table 4). In addition, upon inspecting the uncertainties of individual measurements in Table 3, we note that a best-fit line between these two quantities will be heavily influenced by a small number of remnants where N_H measurements have very small (<5%) formal uncertainties. Even though the sample of regions and the range of models we considered in the spectral fits did not allow us to measure a systematic uncertainty for these remnants, considering the many possible sources of systematic uncertainty that are usually present in the determination of the hydrogen column densities, we assign a 5% error to the N_H measurement of SNRG130.7+3.1, SNRG69.0+2.7, SNRG004.5+06.8, SNRG53.6−2.2, SNRG54.1+0.3, and SNRG00.0+0.0 to avoid biasing the results.

Zoom In Zoom Out Reset image size

Using these data, we obtain a best-fit linear relation between the hydrogen column density and the optical extinction that is described by

$$\begin{eqnarray}&&{N}_{{\rm{H}}}=2.87\pm 0.12\times {10}^{21}\ {A}_{{\rm{V}}}\;{{\rm{cm}}}^{-2},\end{eqnarray} \tag{ 9 }$$

where A_V is in magnitudes and the error represents the 1σ uncertainty. We show the best-fit line in Figure 5.

We can test the effect of imposing a floor to the N_H uncertainties by also fitting to the raw values. This increases the best-fit value of the N_H/A_Vto (2.92 ± 0.11) × 10²¹ cm⁻², which is within 1σ of the best-fit value presented above, but dominated by remnants with small errors. Finally, we note that the remnants with either very low or high N_H will have the greatest leverage on the final fit. To test the magnitude of this effect, we remove the points with the highest and lowest values of the hydrogen column density. This results in a best-fit value of (2.76 ± 0.13) × 10²¹ cm⁻², again within 1σ of the best-fit value using all of the data and a minimum error imposed on the N_H measurements.

These results indicate that a consistently analyzed sample of SNRs using ISM abundances for spectral fitting produces an N_H/A_V relation that is significantly higher than previous estimates. This has the effect of increasing the N_H derived from existing E(B−V) measurements. Conversely, using the new relation to estimate the brightness of an optical counterpart of an X-ray source will result in less extinction compared to estimates calculated with previous values of the N_H/A_V relation. As anticipated by Güver & Özel (2009) and further discussed by Watson (2011), the primary reason for this difference is the change in the abundances that are used for the ISM. The default metallicity library of the frequently used spectral analysis package xspec uses the solar abundances given by Anders & Grevesse (1989) and this is indeed the main library that is utilized in the studies of the spectral properties of X-ray sources. However, these values are known to be on average 45% higher than the updated values by Asplund et al. (2009) using the solar spectrum and by Wilms et al. (2000) for the ISM, and therefore result in a smaller N_H value for a given amount of total absorption due to the matter in the ISM. We showed here in our systematic analysis that a change in the abundance table results in ≈30% change in the coefficient of the linear relation between the hydrogen column density and optical extinction.

As a final caveat, we note that when the hydrogen column density is measured using the ISM abundances (e.g., from Wilms et al. 2000), the relation we presented here should be used to predict or compare with the extinction in the optical band. On the other hand, if an N_H is found using solar abundances (e.g., from Anders & Grevesse 1989), then the relation reported in Güver & Özel (2009) should be employed for self consistency.

In this paper, we presented a comprehensive study of interstellar X-ray extinction using the Chandra SNR archive. We used standardized procedures and made use of the high energy and spatial resolution of the data set to assess the uncertainties in the measurement of the hydrogen column density from X-ray spectra. In contrast with earlier work, we modeled the interstellar extinction using the latest ISM model as well as interstellar abundances of Wilms et al. (2000). We also modeled SNR spectra with a variety of thermal, non-thermal, and mixed models to evaluate the effects of the continuum models on the measured hydrogen column density. In addition to assessing the uncertainties in this measurement arising from the range of models, we explored uncertainties due to the spatial variations within individual remnants and the different physical conditions of the remnants such as their compositions, temperatures, and non-equilibrium regions.

We used a Bayesian statistical analysis tools to determine the systematic uncertainties in the hydrogen column density for each remnant. We used the hydrogen column density measurements and their uncertainties in combination with the measurements of the optical extinction toward these same remnants to determine the relation between these quantities. We found a best-fit linear relation described by ${N}_{{\rm{H}}}=(2.87\pm 0.12)\times {10}^{21}\;{A}_{{\rm{V}}}\;{\mathrm{cm}}^{-2}$.

A couple of different avenues could lead to further progress in the determination of the empirical relation between the optical extinction and the hydrogen column density, as also discussed in Güver & Özel (2009) and Watson (2011). First, a larger number of high-quality optical or near-infrared spectra of SNRs need be obtained to increase the sample of precise optical extinction measurements. This is because while the extensive archives of X-ray satellites allow us to determine the hydrogen column density and its uncertainty accurately, optical/NIR SNR data, especially for strongly reddened regions in the Galaxy, are scarce. The second strategy targets different populations of sources that are bright in both optical and X-ray bands, such as blazars. The measurement of the hydrogen column density and the optical extinction in such a population would provide an independent measurement of the relation between these two quantities and help cross check the results.

This work was supported in part by Chandra award No. AR0-11007X. T.G. has been supported by Istanbul University: project Nos. 49429 and 57321. F.Ö. gratefully acknowledges partial support from NSF grant AST 1108753. P.S. acknowledges partial support from NASA contract NAS8-03060. This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO, ChIPS, and Sherpa. The authors thank K. Arnaud for helpful discussions regarding the tbabs model in xspec.