The Charge State of Polycyclic Aromatic Hydrocarbons across a Reflection Nebula, an H ii Region, and a Planetary Nebula

NGC 7023. NGC 7023 is the RN in Cepheus irradiated by the B2Ve (Rogers et al. 1995), massive spectroscopic binary (Millan-Gabet et al. 2001), Herbig Be star HD 200775. It is located some 430 pc from Earth (van den Ancker et al. 1997). The left panel in Figure 1 presents the morphology of NGC 7023's northwest  photodissociation region (PDR) for which Spitzer-IRS spectral map data are available. This PDR was previously studied in BBA13, BBA14, BBA15, and BBA16 employing the database-fitting technique used here. Other recent work studying the PAH emission from this source include Berné et al. (2007), Rosenberg et al. (2011), Pilleri et al. (2012), Stock et al. (2016), and Croiset et al. (2016).

Zoom In Zoom Out Reset image size

M17. M17 is the massive young star-forming region located in Sagittarius, some 1800 pc from Earth (Kharchenko et al. 2005). It is one of the most massive nearby star formation regions in the Galaxy (Povich et al. 2009). The most luminous of the seven ionizing O stars in M17 is the  binary system CEN1 (Chini et al. 1980). The middle panel of Figure 1 presents the morphology of M17  around the southwest photodissociation front—located in its southern bar—as seen by the SIRIUS (Nagashima et al. 1999) instrument at the IRSF 1.4 m telescope that covers the area for which Spitzer-IRS spectral map data are available.  The data cover an enormous 1' × 1' (∼1 pc²) area and trace emission from just beyond the hard H ii region–PDR boundary, located about 10'' (∼0.09 pc) from the northeastern corner of the field, to well into the giant molecular cloud (38; ∼2 pc). The spectral map data overlap with ISO-SWS observations studied by Verstraete et al. (1996), albeit avoiding the H ii region directly where PAHs are unlikely to survive the harsh environment. The well-known PAH emission in M17 has most recently been studied by Povich et al. (2007),  Doney et al. (2016), Stock et al. (2016), and Yamagishi et al. (2016).

NGC 40. NGC 40 is the PN located in Cepheus some 1250 pc from Earth (Stanghellini et al. 2008). NGC 40 is well studied and considered unusual when compared to the majority of Galactic PNe, notably because it shows spectral variability on a scale of hours (Balick et al. 1996; Grosdidier et al. 2001). Modeling the spectrum of the late-type WC8 Wolf–Rayet exciting central star suggests an effective temperature of up to 90,000 K (Bianchi & Grewing 1987; Leuenhagen et al. 1996; Marcolino et al. 2007). The right panel in Figure 1 presents the morphology of the region as seen by the WIYN/ESO 3.5 m telescope that covers the area for which Spitzer-IRS spectral map data are available.

The Spitzer-IRS low-resolution observations for the three targets were downloaded from the Spitzer Heritage Archive (SHA). The spectra cover 5.2–14.5 μm at a resolution of $R\equiv \lambda$/Δλ ∼ 64–128. The spectral maps were assembled using CUbe Builder for IRS Spectra Maps (CUBISM; Smith et al. 2007a) in its default configuration with the most recent calibration files. CUBISM's automatic bad-pixel-generation parameters were set to ${\sigma }_{{\mathtt{trim}}}=7$ and Minbad-fraction=0.5 and 0.75 for the global and recorded bad pixels, respectively. Spurious data at the ends of the SL slit were ignored by applying a wavsamp of 0.05–0.80. A rather large 0.2 cut at the high end was necessary to avoid striping artifacts that would appear for some of the lower signal-to-noise data, notably those with redundancies along the slit. No separate sky subtraction was performed due to the lack of such data. Since a local continuum is established and subtracted before isolating the PAH emission bands (see Section 3), and the actual sky is not expected to have much substructure; this does not pose a problem. The extracted SL orders SL1 and SL2 were written to file, after which remaining bad, NaN-valued pixels were interpolated in the spatial domain. With the help of the IDL Astronomy Library (Landsman 1993), the constructed spectral maps were resampled, combining 2 × 2 pixels, reaching an ultimate pixel size of 36 × 36. Note that no spatial corrections were applied for the wavelength-dependent point-spread function (PSF). The SL1 and SL2 data were spatially aligned and combined into a single data cube, where the SL2 order was multiplicatively spliced to the SL1 order by conserving integrated flux in the spectral region of overlap. Following this, the spectra were σ-clipped, where the flux contained in each resolution element was evaluated with respect to the mean of its two preceding and two following spectral elements. If the flux in the resolution element under consideration was removed three or more standard deviations from the surrounding mean, it was replaced with the linear interpolant across the adjacent elements. Lastly, a mask was constructed recording pixels with no spectrum in the sparse maps and pixels with a spectrum of low-quality and/or "abnormal" structure, e.g., spectra showing a strong hot dust continuum component associated with embedded sources.

Errors were propagated by taking the squared uncertainty maps from CUBISM and performing the same reduction steps as on the spectral data while assuming additive uncertainties (${\sigma }_{x}^{2}={\sum }_{i}{\sigma }_{i}^{2}$). Uncertainties from splicing and σ-clipping are ignored but are not expected to exceed 10%.

For each target, the spectra have been organized into bins using hierarchical clustering with a divisive scheme. A similarity matrix is constructed from the pairwise summed Euclidean distance norm of each two spectra, which were normalized to the total integrated intensity at each pixel.

As with k-means clustering, the established cluster bins divide the region into comparable morphological neighborhoods (zones) that can be used to help visualize spectroscopic and PAH population changes across them (see e.g., BBA14). In this study, two iterations, resulting in ${2}^{2}=4$ cluster bins, with a minimum of 10 spectra for a cluster bin, proved sufficient to map the salient spectral changes across each target.

Figure 2 presents the average 5.2–14.5 μm Spitzer-IRS spectrum of each bin for each target, and Figure 3 overlays the location of the spectra associated with each bin, hereafter referred to as zones, onto the overview image of each target.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The advantage of hierarchical over k-means clustering is that zone structures tend to better reflect discontinuities in the underlying morphology. This is because k-means clustering is based on average spectra, while hierarchical clustering, as used here, is based on spectral differences.

Keep in mind that the analysis and description that follow treat the data on a pixel-by-pixel basis and do not involve clustering techniques. However, the zones established through clustering are used to track comparable spectra—which are connected to morphologically similar regions—throughout the analysis and interpretation of the results.

To isolate the PAH emission spectrum, the stellar contribution, continuum emission, possible effects from extinction, and any molecular hydrogen and atomic emission lines must first be removed. Given the sizable amount of data, an automated approach is used. The narrow emission lines, continuum, stellar contribution, extinction, and broader PAH features are simultaneously modeled and fitted. The continuum is represented by a combination of several modified blackbodies at fixed temperatures and a power law, the stellar component by a blackbody at fixed temperature, the emission lines by Gaussian profiles, and the PAH features by Drude profiles. Initial line PAH centroids—together with their widths—and parameters describing the modified blackbodies have been taken from Smith et al. (2007b), augmented with other relevant atomic line positions. For the line profiles, the centroids and widths are allowed to vary by 5% and 10%, respectively. The initial line width is taken as 33% of the wavelength span in the observations taken up by the five resolution elements around the centroid position. For the PAH profiles, the centroids and widths are fixed. The power of all the emission profiles, continuum, and stellar component is forced strictly positive. A fully mixed model is chosen for the extinction, which leads to an overall diluting factor of $(1-\exp (-{\tau }_{\nu }))/{\tau }_{\nu }$ (e.g., Disney et al. 1989; Smith et al. 2007b). Here the optical depth ${\tau }_{\nu }$ is in terms of the visual extinction A_V (i.e., ${\tau }_{\nu }=1.8\times {10}^{21}{C}_{\mathrm{ext}}{A}_{{\rm{V}}}$). For the extinction cross section (${C}_{\mathrm{ext}}$), three curves are considered: the ${R}_{{\rm{V}}}=5.5$ and ${R}_{{\rm{V}}}=3.1$ curves from Weingartner & Draine (2001) and the Galactic center extinction curve from Chiar & Tielens (2006). For the latter, the visual extinction is related to the optical depth through ${A}_{{\rm{V}}}=19{\tau }_{9.7}$ (e.g., Roche & Aitken 1984; Mathis 1998). The analysis that follows uses the Galactic center extinction curve from Chiar & Tielens (2006). Section 5.5 discusses how the choice of extinction curve affects the results and their interpretation. Figure 4 illustrates the approach and depicts the breakdown of the 5.2–14.2 μm spectrum of a position in M17 into its different emission components. The figure shows that the spectrum is matched with high accuracy.

Zoom In Zoom Out Reset image size

The molecular hydrogen and atomic emission lines with a signal-to-noise ratio greater than 2 are subsequently removed from the spectra by linearly interpolating across the feature's instrumental width. If present, the ${{\rm{H}}}_{2}$ S(2) and/or [Ne ii] emission lines are treated differently, as they are heavily blended with the 12.7 μm PAH band at the resolution of the IRS's SL module. Separating these lines from the PAH feature follows the approach described by Shannon et al. (2015), in which an emission-line-free 12.7 μm PAH band is used as a template to fit the feature in the other spectra. Here an emission-line-free 12.7 μm PAH band from NGC 7023 is used for all targets. The fitting region is selected such that it excludes resolution elements affected by the emission lines. In the case of a strong [Ne ii] emission line, the entire 12.7 μm PAH band is replaced with the entire scaled template profile. In cases where there is only H₂ emission, the resolution elements of the H₂ (S2) line are only replaced with those from the scaled template profile. Shannon et al. (2015) reported that this method is highly reliable, with an error of no more than 2%–3%. Lastly, the complete PAH emission spectrum is isolated by correcting for extinction and subtracting the continuum.

The database-fitting approach follows that outlined in BBA16, with some adjustments and enhancements. Fitting an astronomical spectrum is a two-step process. First, density functional theory (DFT)–computed absorption spectra from PAHdb's libraries are converted into PAH emission spectra. Second, these modeled emission spectra are used to actually fit the observed astronomical PAH emission spectra. A brief overview of the approach taken in modeling the PAH emission process and the utilized fitting technique follows.

3.2.1. Converting PAH Absorption Spectra into PAH Emission Spectra

3.2.2. Fitting Astronomical PAH Spectra

Fitting technique. Spectroscopic fitting is performed using a nonnegative least-squares (NNLS) algorithm (Lawson & Hanson 1974) that has been adapted to minimize the ${\chi }^{2}$ (non-negative least-Chi-squares, NNLC; Désesquelles et al. 2009). Note that the plateaus underlying the distinct PAH features are included for the fit, and, due to the limited spectral resolution of the SL module, it is not necessary to resample the emission onto a uniformly sampled frequency grid, as was needed in BBA13.

The database analyses are performed using version 2.00 of the library of computed spectra and the AmesPAHdbIDLSuite² (dated 2015 August 27; Bauschlicher et al. 2010; Boersma et al. 2014a) from PAHdb. It is the AmesPAHdbIDLSuite suite of IDL³ -object classes that interacts with the spectral library, applies the PAH emission model to convert the absorption data into emission spectra, implements the NNLC fitting algorithm to fit each astronomical emission spectrum with the computed emission spectra, and does the bookkeeping for breaking the emission down into contributing PAH subclasses, e.g., charge and size. The database-fitting approach is illustrated in Figure 5, where the same spectrum of M17 shown in Figure 4 is broken down by PAH charge. The figure shows that the overall fit to the observed spectrum is good.

Zoom In Zoom Out Reset image size

Figure 6 presents three-color composite images for each of the targets representing charge and size, constructed from the PAHdb breakdown described above. One thing to keep in mind is that the anion contribution in the charge breakdown image is likely overestimated due to the overlap of the red wing arising from the anharmonic redshift with the possible anion contribution (red channel), which consistently falls to the red of the main PAH band (see also Section 5.6). For PAH size, a nonuniform histogram with three bins is constructed, binning the absolute contribution from PAHs with effective radii between 0 and 6.5 Å (${N}_{{\rm{C}}}\sim 20\mbox{--}52$), small; 6.5–10 Å (${N}_{{\rm{C}}}\sim 53\mbox{--}123$), large; and 10–12 Å (${N}_{{\rm{C}}}\sim 124\mbox{--}178$), very large. These correspond to the red, green, and blue channels, respectively. The effective radius ($\bar{a}$), in Å, is calculated as

$$\begin{eqnarray}&&\bar{a}\equiv \sqrt{\displaystyle \frac{A}{\pi }},\end{eqnarray} \tag{ 1 }$$

where A is the total area of the PAH in Å², calculated from counting the number of three to eight–membered rings (n_i) and multiplying that number by the area of an individual ring, i.e.,

$$\begin{eqnarray}A & = & {n}_{3}0.848+{n}_{4}1.96+{n}_{5}3.37\\ & & +{n}_{6}5.09+{n}_{7}7.12+{n}_{8}9.46\ [{\mathring{\rm A} }^{2}].\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Table 2 provides a few examples of PAHs and their effective radii.

Table 2.  Eight PAHs and Their Effective Radii from Equations (1) and (2)

Formula	Size Bin	Structure	Radius (Å)
C20H12	s		2.75
C24H12	s		3.37
C32H14	s		4.03
C66H20	s		6.25
C96H24	L		7.74
C112H26	L		8.44
C150H30	L		9.94
C190H34	VL		11.3

Note. s: small; L: Large; VL: very large.

Download table as:  ASCIITypeset image

The literature provides several "traditional" approaches to separate the PAH bands from the underlying continuum and their plateaus (e.g., Sellgren et al. 2007; Smith et al. 2007b). The approach adopted here is based on that described in van Kerckhoven et al. (2000), Hony et al. (2001), Peeters et al. (2002), and van Diedenhoven et al. (2004). Given the sizable amount of data, this procedure has been automated using a stiff (${\sigma }_{{\mathtt{tension}}}=5$) spline continuum, as demonstrated in Figure 7. Here anchor points are chosen at the minima within a five-element search window surrounding fixed positions between the observed features to isolate the 6.2, "7.7," 11.2, 12.0, 12.7, 13.5, and 14.2 μm PAH bands. Note that the 6.0 and 11.0 μm satellite features are included when isolating the 6.2 and 11.2 μm PAH bands, respectively. At the resolution of the SL data, these satellite features are heavily blended with the main features. The "7.7" μm complex was further decomposed by fitting a second-order polynomial between points near 8.0 and 9.0 μm, allowing extraction of the 8.6 μm PAH band. The 7.6 and 7.8 μm components were separated by integrating the emission falling shortward and longward of 7.72 μm for the 7.6 and 7.8 μm PAH bands, respectively. It is noted that, compared to BBA16, an additional anchor point has been introduced at 10.7 μm to avoid the continuum cutting off part of the 11.2 μm PAH band in spectra with a strong 10–15 μm plateau, notably for NGC 40.

Zoom In Zoom Out Reset image size

The analysis is performed on the extinction-corrected, emission line–removed, and continuum-subtracted spectra from Section 3.2. Uncertainties have been statistically propagated. While determination of the spline continuum will have added to the uncertainty, this should be small and systematic for the high-quality data. Thus, continuum subtraction should not significantly affect differences and ratios. Indeed, several studies have now established that PAH band trends are mostly independent of the way the bands are traditionally isolated (e.g., Galliano et al. 2008; Peeters et al. 2012).

PAH band strength ratios between bands in the 6–15 μm region have long been used as qualitative proxies to track variations in the properties of the emitting PAH population, e.g., charge and structure. Comparing the behavior of these ratios with the ratio determined from the database-fitting approach allows one to quantitatively calibrate these proxies.

Most commonly, the 6.2 μm PAH band is used as a proxy for PAH cations, while the 11.2 μm PAH band is used as a proxy for neutral PAHs, i.e.,

$$\begin{eqnarray}&&\displaystyle \frac{{n}_{{\mathrm{PAH}}^{+}}}{{n}_{{\mathrm{PAH}}^{0}}}\propto \displaystyle \frac{{I}_{6.2}}{{I}_{11.2}},\end{eqnarray} \tag{ 3 }$$

where ${n}_{{\mathrm{PAH}}^{+}}$ and ${n}_{{\mathrm{PAH}}^{0}}$ are the PAH cation and neutral densities determined using PAHdb, respectively, and ${I}_{6.2}$ and ${I}_{11.2}$ are the 6.2 and 11.2 μm PAH band strengths, respectively. Figure 5 shows that these proxies are not quantitative by themselves, as PAHs in all three charge states contribute to each band. Nonetheless, they have proven to be generally reliable tracers for the charge state of the PAH population  (e.g., Hony et al. 2001; Peeters et al. 2002; Bregman & Temi 2005; Compiègne et al. 2007; Galliano et al. 2008; BBA16).

Figure 8 presents the traditionally determined 6.2/11.2 μm PAH band strength ratio map for each of the targets, and the top panels in Figure 9 show the corresponding ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio map, where ${n}_{{\mathrm{PAH}}^{+}}$ and ${n}_{{\mathrm{PAH}}^{0}}$ are directly taken from the PAHdb charge decomposition. The bottom panels in Figure 9 show the PAHdb-determined large fractional size (${n}_{{\mathrm{PAH}}_{n{\rm{C}}}\gt 50}$/(${n}_{{\mathrm{PAH}}_{n{\rm{C}}}\gt 50}$ + ${n}_{{\mathrm{PAH}}_{n{\rm{C}}}\leqslant 50}$)) contribution maps. The black line in Figure 10 traces the  PAH charge, and the orange line traces the  PAH size changes along the line shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Knowledge of the charge state of the PAH population can be directly related to the physical conditions in the local environment through the PAH ionization parameter (Bregman & Temi 2005; Galliano et al. 2008; Fleming et al. 2010; BBA15; BBA16; Stock et al. 2016). The PAH ionization parameter is defined as ${G}_{0}{T}_{\mathrm{gas}}^{1/2}/{n}_{{\rm{e}}}$, where G₀ is the strength of the radiation field in terms of the Habing field (Habing 1968; ${G}_{0}=1\equiv 1.6\times {10}^{-3}$ erg cm⁻² s⁻¹), ${T}_{\mathrm{gas}}$ is the gas temperature (in K), and n_e is the electron density (in cm⁻³).

When assuming only two accessible ionization states and parameters applicable for circumcoronene⁴ (C₅₄H₁₂),  which is considered representative for an average interstellar PAH (${n}_{{\rm{C}}}=50\mbox{--}100$; e.g., Croiset et al. 2016), the PAH ionization parameter can be expressed as (see Boersma et al. 2014a)

$$\begin{eqnarray}&&2.66({n}_{{\mathrm{PAH}}^{+}}{/n}_{{\mathrm{PAH}}^{0}})\approx {G}_{0}{T}_{\ \mathrm{gas}}^{1/2}{/n}_{{\rm{e}}}\ [{10}^{4}\ {\mathrm{cm}}^{3}].\end{eqnarray} \tag{ 4 }$$

The reasons for only considering two PAH ionization states (PAH${}^{0}$, PAH${}^{+}$) is due to the unreliability associated with determining the anion contribution. In Section 5.6, among other things, the quantification of the PAH charge is further discussed.

Figure 2 clearly shows pure rotational ${{\rm{H}}}_{2}$ lines present in spectra of NGC 7023 and M17. Following an approach comparable to that outlined by Fleming et al. (2010), the pure rotational ${{\rm{H}}}_{2}$ lines are utilized to determine the gas temperature (${T}_{\mathrm{gas}}$), H₂ column density (${N}_{{{\rm{H}}}_{2}}$), and H₂ ortho-to-para ratio (${R}_{\mathrm{OP}}$). The ${{\rm{H}}}_{2}$ molecular parameters were obtained from the HITRAN Database (Rothman et al. 2013). Here ${T}_{\mathrm{gas}}$, ${N}_{{{\rm{H}}}_{2}}$, and ${R}_{\mathrm{OP}}$ are determined by minimizing the Euclidean distance norm between observations and model obtained from fitting a straight line with parameters ${T}_{\mathrm{gas}}$ and ${N}_{{{\rm{H}}}_{2}}$ while iterating on ${R}_{\mathrm{OP}}$. At least three ${{\rm{H}}}_{2}$ lines with a signal-to-noise ratio greater than three are required for a pixel to be evaluated. In addition, ${R}_{\mathrm{OP}}$ is forced strictly positive and needs to converge and not exceed three. Figure 11 demonstrates the approach by showing the population diagram constructed from the ${{\rm{H}}}_{2}$ lines present in a 5.2–14.2 μm SL Spitzer-IRS spectrum from M17 near the location of the spectrum shown in Figures 4, 5, and 7. Note that the ${{\rm{H}}}_{2}$ lines are corrected for extinction, and uncertainties are propagated. Figure 12 presents maps of ${T}_{\mathrm{gas}}$  and ${\mathrm{log}}^{10}{N}_{{{\rm{H}}}_{2}}$ for each target. Due to the lack of any apparent morphological substructure, the ortho-to-para-ratio maps have been omitted from the figure and are not further discussed. The reader is directed to Fleming et al. (2010) for a more detailed analysis and discussion of the ${{\rm{H}}}_{2}$ ortho-to-para-ratio across NGC 7023.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Results on the state of the PAH population are presented for PAH charge and size, looking at both the morphologies present in the different maps and possible PAH proxy calibrations.

The morphological structure of the PAH charge distribution across each target is shown in Figure 6, Figure 8, the top row of panels in Figure 9, and in Figure 10. Figure 6 captures the charge variations across each target using a three-color composite image, where the red, green, and blue channels depict the contribution from PAH anions, neutrals, and cations, respectively. The figure shows a morphological structure in good agreement with that present in Figure 1 for each target. That is, for NGC 7023, there is a clear delineation between the dense and diffuse medium to the north and south of the PDR, respectively, with emission from neutral PAHs peaking in the denser northwest region. Similar behavior holds for M17, where the "bright" and "dark" regions fall to the northeast and southwest, respectively, with emission from neutral PAHs peaking in the northeast. For NGC 40, the emission intensity peaks in the ring surrounding the central white dwarf. While Figure 6 does not apply a normalization, Figure 8 does, thus avoiding any morphological features simply appearing due to a more-is-more type of behavior. Figure 8 uses the 6.2/11.2 μm PAH band strength ratio as a proxy for the PAH charge. For NGC 7023, the ratio shows a gradual decrease when moving away from the irradiating star. However, the "ring" visible in Figures 1 and 6 is absent here. The striking feature in the 6.2/11.2 μm PAH band strength ratio map of M17 is the "hot spot" located slightly northeast of the center of the map. Apart from the obvious circular symmetry around the central star, NGC 40 does not show much substructure in its 6.2/11.2 μm PAH band strength ratio map. Overall, the measured ratios are low—perhaps the ratio peaks off the bright "ring" visible in Figures 1 and 6. Figure 9 uses the PAHdb-determined ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio to track the PAH charge morphology across each target. For NGC 7023, there is a general decrease in the ratio when moving away from the irradiating source, and the maps shown in Figures 8 and 9 are very comparable. For M17, the "hot spot" present in Figure 8 is absent in Figure 9. However, a smaller, less bright "hot spot" appears slightly offset to the east of the one present in Figure 8. For NGC 40, the morphology in Figure 9 shows a uniform brightening at the northeastern edge of the region that is choppy in Figure 8. However, the interior regions show more morphological agreement between the two figures. Finally, Figure 10 traces the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio along the lines shown in Figure 3 and provides a one-dimensional view of the PAH charge variations across each target as a function of distance from the irradiating source. Although the (somewhat noisy) trace across M17 does not show significant variation, the traces for NGC 7023 and NGC 40 do. For NGC 7023, the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio jumps from ∼0.5–0.6 in the vicinity of the exciting star to ∼0.9 some 25'' from HD 200775. The ratio remains at about 0.8 across the diffuse region until it crosses the PDR, where it starts a steady decline as it moves further into the dense region. For NGC 40, the behavior of the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio across the nebula clearly reflects a spherical morphology. It is symmetric with respect to the exciting star, peaking at the edge of the ring to ∼0.8–0.9, dropping between the exciting star and ring to ∼0.1–0.3, and rising to 0.4 as it crosses the Wolf–Rayet star at its center.

Moving away from morphology, all of the data from Figures 8 and 9 are combined to quantitatively calibrate the 6.2/11.2 μm PAH band strength ratio with the PAHdb-determined ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio and the PAH ionization parameter ${G}_{0}{T}^{1/2}/{n}_{{\rm{e}}}$. This calibration is done in Figure 13. The figure shows that the dynamic range in abscissa values is dominated by NGC 7023, which shows a reasonably good correlation with the 6.2/11.2 μm PAH band strength ratio. NGC 40 lines up well with NGC 7023, albeit only occupying low ordinate and abscissa values. M17 also follows the overall trend with a slightly larger range in abscissa values compared to NGC 40 and a coagulation of points along a vertical line around  ${G}_{0}{T}^{1/2}/{n}_{{\rm{e}}}\ [\times {10}^{4}\ {\mathrm{cm}}^{3}]$ of about three. The poor squared linear correlation coefficient (${R}^{2}=0.56$) is mostly due to M17, as its number of data points exceeds that on NGC 7023 and NGC 40 combined by almost a factor of five. Omitting the M17 data yields a linear correlation coefficient of ${R}^{2}=0.94$.

Zoom In Zoom Out Reset image size

In BBA15 and BBA16, the 8.6/11.2 μm PAH band strength ratio was shown to be very useful as a proxy for the PAH charge. Plotting the 8.6 versus 11.2 PAH band strength for NGC 7023 resulted in a bifurcated plot establishing two well-defined limiting lines: one associated with the dense medium, the other with the diffuse medium (BBA14). Figure 14 presents these relationships for the targets studied here. The left panel of the figure shows an overall good trend of the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio with the 8.6/11.2 μm PAH band strength ratio, where the range in the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio is dominated by NGC 7023. Compared to the 6.2/11.2 μm PAH band strength ratio shown in Figure 13, both M17 and NGC 40 are more in line with the established trend. The somewhat poor squared linear relation coefficient (${R}^{2}=0.62$) is again due to M17 dominating the coefficient because of its large number of data points when compared to the other two targets.

Zoom In Zoom Out Reset image size

The right panel in Figure 14 presents the 8.6/11.2 μm PAH band diagnostic template, where the data on M17 and NGC 40 have been scaled into the domain of NGC 7023 by multiplying both the ordinate and abscissa with a constant factor, i.e., preserving the slope and therefore their associated 8.6/11.2 μm PAH band strength ratio and ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio. Overplotted are lines with a constant ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio. The figure shows that all of the M17 data (green shaded area) fall inside the space occupied by that of NGC 7023 (purple shaded area), while little of the NGC 40 data (red shaded area) overlaps with that  of M17 or NGC 7023. The bulk of the NGC 40 data extends toward and abuts the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}=0$ limit.

4.2.1. PAH Size

From the pure rotational ${{\rm{H}}}_{2}$ lines present in the spectra, the ortho-to-para ratio, molecular hydrogen density, and gas temperature for NGC 7023 and M17 were derived from fits to population diagrams such as that presented in Figure 11. One thing to note in that figure is that the reported uncertainty for ${N}_{{{\rm{H}}}_{2}}$ is high, which is not immediately expected when looking at the figure. The large uncertainty comes from propagating the uncertainty on ${\mathrm{log}}^{10}{N}_{{{\rm{H}}}_{2}}$ (=19.96; ${\sigma }_{{\mathrm{log}}^{10}{N}_{{{\rm{H}}}_{2}}}=0.2905$), which is what is established by fitting the dashed straight line through the three data points shown in Figure 11. With x as ${\mathrm{log}}^{10}{N}_{{{\rm{H}}}_{2}}$, the propagated uncertainty is given by

$$\begin{eqnarray}&&{\sigma }_{{N}_{{{\rm{H}}}_{2}}}=\displaystyle \frac{d({10}^{x})}{{dx}}{\sigma }_{x}={10}^{x}{\mathrm{log}}^{e}(10){\sigma }_{x},\end{eqnarray} \tag{ 5 }$$

which equates to the reported ${10}^{19.96}\cdot 2.303\cdot 0.2905\,=6.1\times {10}^{19}$ cm⁻².

Figure 12  shows that the molecular hydrogen density ranges from about 19.75 to 22.6 and 19.5 to 22.2 for NGC 7023 and M17, respectively. The temperature of the gas ranges from about 450 to 1100 and 400 to 1100 K for NGC 7023 and M17, respectively.  Both the molecular hydrogen density and temperature maps for NGC 7023  show a clear separation across the PDR boundary, with higher densities and cooler temperatures to the north of the PDR. The M17 map for the molecular hydrogen density shows no real substructure. The map for the gas temperature does not either; however, the color scaling seems dominated by a few hot spots.

Croiset et al. (2016) traced the PAH size evolution in NGC 7023 along a line comparable to that shown in Figure 3. The top frame of their Figure 6 shows an overall steady increase in the 11.2/3.3 μm PAH band strength ratio from 2 to 3 with distance from the exciting star almost all the way to the PDR boundary. There it sharply drops and quickly recovers when moving across the brightest spot in the PDR. Arguably, the noisy and far less well spatially sampled orange trace in Figure 10 shows, in broad strokes, similar behavior for the ${n}_{{\mathrm{PAH}}_{{nC}\gt 50}}/{n}_{{\mathrm{PAH}}_{{nC}\leqslant 50}}$ ratio, with a small increase in the ratio from 1.2 to 1.3 with distance from HD 220775 up to the PDR interface (∼51''), where it sharply drops. However, the ratio fails to recover when moving deeper into the molecular cloud. However, it should be noted that the trace in Figure 10 extends some 20'' further into the molecular cloud and actually does not capture the brightest spot in the PDR, and, given the poor spatial sampling, it could well be that the sharp drop and its subsequent recovery are captured in a single pixel.

The PAH charge distribution is discussed in terms of the morphology of the targeted regions and the quantitative calibration of the PAH band strength ratio proxies with the PAHdb-determined charge breakdown.

5.2.1. The Morphology

5.2.2. The PAH Charge Proxy Calibrations

With the aim of developing diagnostic tools that are straightforward to use and that can be applied to all astronomical PAH spectra and to assess the presence of principal components (see, e.g., Berné et al. 2007; Pilleri et al. 2012), as discussed in BBA13 and BBA16, PAHdb template spectra for PAH charge and size have been derived. These template spectra are determined by normalizing the integrated flux for each of the PAH subclass spectra obtained from the database fit at each pixel to unity and taking the average for each target. The resulting five PAH subclass templates—i.e., three for PAH charge (anion, neutral, and cation) and two for size, large (${n}_{\mathrm{carbon}}\gt 50$) and small (${n}_{\mathrm{carbon}}\leqslant 50$)—are presented in Figure 15.⁵ Note that the templates are extended to cover the entire 2.5–20 μm wavelength range, while the fitted data cover only 5.2–14.2 μm. The figure shows that each of the PAH template spectra compares well for each different target. This is corroborated by Table 3, which presents the Pearson correlation coefficient between each of the template spectra derived for each target. The table shows that the template spectra derived from NGC 7023 match best with those derived from M17. In turn, those derived from NGC 40 match best with those derived from M17. However, the Pearson correlation coefficients between the template spectra derived in BBA16 for NGC 7023 and the RN NGC 2023—with coefficients of 0.99, 0.98, 1.00, 0.99, and 0.96 between their anion, neutral, cation, large, and small template spectra, respectively—surpass any of those reported here for NGC 7023.

Zoom In Zoom Out Reset image size

The validity of the template spectra is evaluated by fitting each target with the template spectra for the charge derived for the northwest PDR in NGC 7023. Here the observed spectrum at each position is fitted as the weighted sum of the three individual NGC 7023 charge template spectra:

$$\begin{eqnarray}&&{I}_{\mathrm{obs}}={{aI}}_{{\mathrm{PAH}}^{-}}+{{bI}}_{{\mathrm{PAH}}^{0}}+{{cI}}_{{\mathrm{PAH}}^{+}}.\end{eqnarray} \tag{ 6 }$$

Figure 16 presents the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio, i.e., c/b ratio, maps recovered for each target.

Zoom In Zoom Out Reset image size

Following Shannon et al. (2016), the Structural SIMilarity (SSIM) index is used to quantify morphological similarity between the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio maps determined from the template fit decomposition (i.e., Figure 16) and those from a direct database decomposition (i.e., Figure 9). The SSIM index is calculated on various windows, and for two images x and y, it is defined as (Wang et al. 2004)

$$\begin{eqnarray}&&\mathrm{SSIM}(x,y)=\displaystyle \frac{(2{\mu }_{x}{\mu }_{y}+{c}_{1})(2{\sigma }_{{xy}}+{c}_{2})}{({\mu }_{x}^{2}+{\mu }_{y}^{2}+{c}_{1})({\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{c}_{2})},\end{eqnarray} \tag{ 7 }$$

where ${\mu }_{x/y}$ and ${\sigma }_{x/y}$ are the average and variance of images x and y, respectively, and ${\sigma }_{{xy}}$ is their covariance. Here c₁ and c₂ are constants defined as ${c}_{1}={({k}_{1}L)}^{2}$ and ${c}_{2}={({k}_{2}L)}^{2}$, respectively. Here L is the dynamic range of the pixel values, i.e., 255 for an 8 bit image, and k₁ and k₂ are, respectively, by default 0.01 and 0.03. The convolution window is a 2D Gaussian with a width of 1.5 pixels, which is used for convolution. Two identical windows/images will have an SSIM index of 1. The original Matlab implementation⁶ from Wang et al. (2004) was ported to IDL⁷ and validated. It is noted that Shannon et al. (2016) utilized a uniform convolution window. Figure 17 presents the SSIM index maps and gives the averaged SSIM index. The figure shows that the results from the template fit and the direct database decomposition compare well for NGC 7023 and M17; NGC 40 fares far less well.

Zoom In Zoom Out Reset image size

For M17, the "hot spot" seen in the 6.2/11.2 μm PAH band strength ratio map, which is absent in the ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio map from the direct database decomposition, is reproduced by the template fitting, hence explaining the low SSIM index values near those positions in Figure 17.

While the template spectra for PAH charge derived from NGC 7023 produce overall acceptable fits, there are astronomical (local) environments in which they do not perform as well systematically. This could point to either a substantial difference in the underlying PAH population or that the ongoing changes to the PAH population are not well captured with variations in PAH charge alone.

The analyses of the ${{\rm{H}}}_{2}$ population diagrams  provide insight into the  ${{\rm{H}}}_{2}$ density and gas temperature.  A direct comparison of the results with those found by Fleming et al. (2010) on the northwest PDR in NGC 7023  show that the morphology of the  molecular hydrogen density and temperature maps is in good agreement. For M17,  many pixels did not meet the conditions set in Section 3.5 and have been left out of the maps. Possibly, the chosen extinction curve has a significant influence on these results (see Section 5.5).

Since the temperature of the gas (${T}_{\mathrm{gas}})$ shows up in the definition of the PAH ionization parameter (${G}_{0}{T}^{1/2}/{n}_{{\rm{e}}}$), it can be used in a simple model to separately (from that outlined in Section 4.2) try to calibrate the 6.2/11.2 μm PAH band strength ratio. The strength of the radiation field (G₀) is modeled by taking a blackbody at the effective temperature of the irradiating star, diluting it with the stellar radius over the projected distance squared, and integrating it from 91.2 to 111 nm. Determining the electron density (${n}_{{\rm{e}}}$) is far less straightforward, and therefore ${n}_{{\rm{e}}}\equiv 1$ cm⁻³ is taken. This value is more or less consistent with that at the PDR in NGC 7023 (${n}_{{\rm{e}}}\sim 0.56$ cm⁻³) when taking the total gas density as 4 × 10³ cm⁻³ (Chokshi et al. 1988), assuming all free electrons come from C ii, and using the diffuse ISM gas-phase carbon abundance of $1.4\times {10}^{-4}$ (Cardelli et al. 1996). Figure 18 presents the results. The figure shows broad-stroke similarities with Figure 13; that is, the PAH ionization parameter values are of similar magnitude, and the dynamic range is determined by NGC 7023. However, Figure 13 shows a much tighter relation, and, specifically for M17, each zone is much more coherent. The most likely reason for this is that the electron density is probably not the same and constant across each target and that using the projected distance is not sufficient to capture the actual morphological details of the region.

Zoom In Zoom Out Reset image size

The results and their interpretation are based on extinction-corrected data. The amount and spectral shape of this correction is highly dependent on the adopted extinction model (e.g., mixed versus screen), the shape of the used extinction curve, and its interplay with other components of the spectra (i.e., continuum, emission lines, etc.). Figure 19 compares the wavelength-dependent dilution factor derived from the three extinction curves considered here when assuming a mixed model (cf. Section 3.1) and a visual extinction A_V of 19, translating to ${\tau }_{9.7}=1$ in the case of the extinction curve toward the Galactic center from Chiar & Tielens (2006). The figure shows that there is a difference in the amount of actual extinction for the same A_V. While this is significant when comparing the extinction curve from Chiar & Tielens (2006) with those from Weingartner & Draine (2001), the biggest impact on any results presented here and their subsequent interpretation would come from differences in spectral shape, since they heavily rely on PAH band strength ratios.

Zoom In Zoom Out Reset image size

While the spectral shape of the core silicate absorption feature centered around 9.7 μm is practically identical between ∼8 and 12 μm, the features from Weingartner & Draine (2001) continue to broaden where the profile from Chiar & Tielens (2006) levels out. From Weingartner & Draine (2001), the ${R}_{{\rm{V}}}=5.5$ extinction curve shows relatively more extinction at shorter wavelengths compared to the ${R}_{{\rm{V}}}=3.1$ curve. From this, it is easily inferred that the choice of extinction curve can have a significant impact on the determined 6.2/11.2 μm PAH band strength ratio. Also, the strong interplay of the chosen extinction curve with the other components of the spectra proved to be, in some cases, highly unpredictable. Notably for M17, there appears to be a degeneracy between the arrived at total visual extinction (A_V) and continuum when using the ${R}_{{\rm{V}}}=5.5$ extinction curve from Weingartner & Draine (2001), causing an instability in the fitting that results in the "banded" pattern apparent in the top panel of Figure 20, where A_V jumps from zero to high values in the span of a single pixel. The extinction curve from Chiar & Tielens (2006), shown in the bottom panel of Figure 20, does not show this behavior.

Zoom In Zoom Out Reset image size

Figure 21 explores the impact of the choice of extinction curve on the measured 6.2/11.2 μm PAH band strength ratio, together with the associated PAHdb-derived ${n}_{{\mathrm{PAH}}^{+}}/{n}_{{\mathrm{PAH}}^{0}}$ ratio and the large PAH fraction. The figure shows that this impact can be as much as 40% when choosing either the ${R}_{{\rm{V}}}=5.5$ extinction curve from Weingartner & Draine (2001) or that from Chiar & Tielens (2006). Despite this, the overall results for M17 and their interpretation when using the extinction curve from Chiar & Tielens (2006) do not change. However, care needs to be taken when choosing an extinction curve, and, in the case of M17 here, that from Chiar & Tielens (2006) could be considered "better," as it avoids the "banded" behavior in the visual extinction map; cf. Figure 20.

Zoom In Zoom Out Reset image size

Besides the choice of extinction curve (Section 5.5), the results presented here also rely strongly on the contents of PAHdb's spectral libraries and, as described in BBA16, on the spectral resolution of the astronomical data, the details of the employed PAH emission model, and the fitting algorithm used.

As far as the contents of PAHdb's spectral libraries are concerned, it is very difficult to fully and quantitatively assess the influence of their (in)completeness on the accuracy of the PAH subclass breakdown. However, in BBA15, it was shown that, at this stage, the breakdown and calibration of the PAH charge using a two-level system is robust. In addition, Andrews et al. (2015) showed that the spectral breakdown into any PAH subclass can be done reliably. Nonetheless, while the overall quality of the PAHdb fits is considered good, they have difficulties matching some of the emission. Notable are the 6.2 μm PAH feature and the PAH emission originating between 10 and 15 μm. In addition, the contribution from PAH anions to the spectrum is unreasonable (see also BBA16). The following three paragraphs discuss each of these in turn.

Hudgins et al. (2005) suggested that the 6.2 μm band may be attributed to PANHs, PAHs containing nitrogen, and subsequent analysis of the emission from the northwest PDR in NGC 7023 required PANH cations to fit the bulk of the emission in the 6.2 and 11.0 μm bands (BBA13). However, the majority of large PANHs in the PAHdb spectroscopic libraries is limited to a collection of isomers in which the loci of the nitrogen atom are permuted rather than a collection of structurally different PANH molecules. The computation of the spectra of a large variety of (large) PANHs and adding them to PAHdb is required to further test this suggestion.

The emission between 10 and 15 μm is attributed to the CH_oop vibrational motion of peripheral hydrogen atoms. Each of the distinct bands in the 10–15 μm is associated with a particular hydrogen adjacency class—the number of neighboring hydrogen atoms protruding from a single aromatic ring, namely, the 11.2, 12.0, 12.7, and 13.5 μm PAH bands with the solo, duo, trio, and quartet adjacency classes, respectively (e.g., Hony et al. 2001; Bauschlicher et al. 2009). Currently, PAHdb's spectral libraries are limited in the number of large PAHs with varying adjacency classes. Most of the larger PAHs (${n}_{\mathrm{carbon}}\gt 50$) have condensed structures and are part of a series where one characteristic is permuted, e.g., the location of a seven-membered ring defect (Bauschlicher 2015). The computation of a large variety of spectra of large PAHs covering all hydrogen adjacency classes and adding them to PAHdb is required to overcome this limitation.

Although the presence of PAH anions in space has been suggested (e.g., Bregman & Temi 2005; Bauschlicher et al. 2009), their contribution in the PAHdb breakdown spectra here is highly uncertain. Bands with red wings and shoulders are intrinsic to emission from highly vibrationally excited molecules because of the anharmonicity of the vibrational potential, and the current emission models available do not treat anharmonic effects. This is the reason only two ionization states (neutral and positive), with singly and doubly charged PAH cations treated as a single class, are considered, with contributions from singly charged PAH anions ignored. Extending PAHdb to include computed PAH spectra that treat the effects of anharmonicity and/or enhance the employed PAH emission model to better treat anharmonicity should help overcome this issue. For recent efforts to understand PAH anharmonic effects on a fundamental spectroscopic level, see, e.g., Maltseva et al. (2015), Mackie et al. (2016), and Candian & Mackie (2017) and references therein. Inclusion of anharmonic effects in the PAHdb emission models is planned once they are better understood.