© ESO, 2013

1. Introduction

Photoelectric heating is a major heating process in photodissociation regions (PDRs), and its efficiency (ϵ_pe) is one of the key parameters to understanding the energy balance there. Theoretical investigations suggest that small grains, in particular polycyclic aromatic hydrocarbons (PAHs), play a dominant role in photoelectric heating (Bakes & Tielens 1994). ϵ_pe is defined as the fraction of energy absorbed by dust that is converted into kinetic energy of the ejected electrons and therefore into gas heating. Observationally it can be estimated by measuring the ratio of the energy emitted in the gas cooling lines ([O i] 63 μm and [C ii] 158 μm are the strongest ones from Av ≲ 5 in PDRs) against the total UV energy absorbed by dust grains, which is probed by the dust infrared (IR) emission. Previous observations show a wide variation of ϵ_pe in different sources, ranging from 10^-4 in W49N, which is illuminated by an intense UV field (Vastel et al. 2001) to 1−2% in the relatively low-UV excited PDR of the Horsehead Nebula (Goicoechea et al. 2009), although mechanical heating also plays a role in W49 (Nagy et al. 2012). Mizutani et al. (2004) showed a variation of ϵ_pe from 0.06 to 1.2% across the 40′ × 20′ area of the Carina Nebula, and an anti-correlation between ϵ_pe and the intensity of the local UV radiation field. For extragalactic sources, Mookerjea et al. (2011) presented ϵ_pe of 0.3−1.2% toward one H ii region in M 33 at a resolution of 50 pc. These variations have been attributed to differences in the mean charge state of the grains, i.e., a positive grain charge results in a decreased efficiency, and the correlation with the intensity of the UV radiation field supports this interpretation. However, considering two properties of different dust populations together, large grains that are the main carriers of the dust IR emission and small grains that are expected to be the main contributors of the photoelectric effect, limits the interpretation. Habart et al. (2001) investigated the photoelectric heating efficiency on small grains and indicated that the observed spatial distribution of [O i] 63 μm and [C ii] 158 μm in L1721 can be better explained by a model with a varying abundance of small grains across the cloud. Joblin et al. (2010) concluded that the evolution of PAHs and very small grains at the border of PDRs should be considered to model the gas energetics. Before Herschel, estimating the gas-cooling energy in spatially resolved PDRs was only rarely possible because of the poor spatial resolution in the far-infrared (FIR) wavelength range. To investigate the charge state of PAHs, mid-infrared (MIR) spectroscopic observations are needed, which can be provided by the Infrared Spectrograph (IRS; Houck et al. 2004) onboard the Spitzer Space Telescope (Werner et al. 2004). The most evident spectral difference between neutral PAHs (PAH⁰) and ionized PAHs (PAH⁺) is the strength of the 6−9 μm complex relative to the 11.3 μm feature, which is supported by laboratory experiments (DeFrees et al. 1993; Pauzat et al. 1994; Langhoff 1996) and theoretical calculations (Szczepanski & Vala 1993b,a; Hudgins & Allamandola 1995a,b). Strong variations of the observed intensity of the individual bands have been seen in different environments, reflecting an evolution of the charge state of PAHs in different physical conditions (Galliano et al. 2008; Sakon et al. 2004; Peeters et al. 2002; Joblin et al. 1996). Recently, Joblin et al. (2008) and Pilleri et al. (2012b) proposed an alternative method of deriving the fraction of PAH⁺ using a spectral decomposition approach based on a few template spectra (Berné et al. 2007; Rapacioli et al. 2005) that include PAH⁰, PAH⁺, and evaporating very small grains (eVSGs).

The combination of the observations by Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) onboard Herschel (Pilbratt et al. 2010) and the IRS onboard Spitzer enables us to investigate the relation between ϵ_pe and the fraction of ionized PAHs in spatially resolved PDRs. In this paper, we report the results for six PDRs studied in the WADI (Warm And Dense Interstellar medium; Ossenkopf et al. 2011) guaranteed time key program of Herschel. WADI is aimed to investigate the physics and chemistry of PDRs and shocked regions with a wide range of physical properties. In this study, we investigate six PDRs located in five different regions that show clear detections of [C ii] and [O i] with PACS and were observed with the IRS.

2. Observations and data reduction

2.1. Targets

Fig. 1 PACS footprints of the blue camera ([O i] 63 μm) and the red camera ([O i] 145 μm and [C ii] 158 μm) overlaid on IRAC 8 μm, except for f) Mon R2. (There IRAC 8 μm is saturated, and the total intensity of AIB bands derived from fitting of IRS spectra (see text) is shown.) Green circles or boxes show the areas from which our PACS and IRS spectra are extracted (subregions in Table 3). For Mon R2 f), the positions of four IR sources (IRS1−4) and PDR 1 from Berné et al. (2009) are marked.

Table 1 shows the global properties of the investigated regions. They cover three orders of magnitude in far-ultraviolet (FUV; hν = 6−13.6 eV) flux, G₀ in units of the Habing field (1.6    ×    10^-6 W m^-2, Habing 1968). Figure 1 shows the PACS field-of-view and the extracted area for each PDR. Some of these sources have been studied and modeled in detail based on the observations with the Herschel-Heterodyne Instrument for the Far-Infrared (HIFI; de Graauw et al. 2010) in WADI (see Joblin et al.2010 for NGC 7023 and Fuente et al.2010 and Pilleri et al.2012a, for Mon R2).

2.1.1. Horsehead

The Horsehead Nebula emerges from the western edge of L1630 as a dark cloud at visible wavelengths, and its outer edge is delineated by a bright and narrow filament in the MIR emission (Abergel et al. 2003, Fig. 1a). The exciting source σ Ori illuminates the filament from the west, and a steep change in column density marks the western edge of the filament, while on the eastern side the MIR emission decreases because of the extinction of the incident radiation by dense material (Abergel et al. 2003). Using observations of the H₂ 1-0 S(1) line emission, Habart et al. (2005) proposed a model with a density gradient from 10⁴ to 10⁵ cm^-3 with a scale length of about 10′′ of this filament. Goicoechea et al. (2009) mapped the [O i] 63 μm emission at a low spectral resolution and concluded that a nonlocal and non-LTE treatment can be important to model the [O i] emission. They derived a value of ϵ_pe of 1–2%.

2.1.2. Ced 201

Ced 201 is a reflection nebula, illuminated by the B9.5 star BD + 69° 1231. Based on radial velocity measurements, Witt et al. (1987) suggested that this reflection nebula is the result of an accidental encounter of a small dense molecular cloud with an unrelated star. Kemper et al. (1999) modeled cooling lines of this PDR and concluded that heating by PAHs and VSGs is required to explain the observed cooling-line emissions. Young Owl et al. (2002) estimated ϵ_pe to be <0.3% based on low spatial resolution KAO observations of [C ii] 158 μm and the upper limit of [O i] 63 μm.

2.1.3. NGC 7023

NGC 7023 is a prototype PDR, which has been widely studied at many wavelengths. The UV radiation from the exciting B star, HD 200775, creates three main PDRs in this nebula; the brightest NGC 7023 North-West (hereafter NGC 7023 NW) is located about 40′′ northwest of the star, another PDR lies about 70′′ south (NGC 7023 S), and NGC 7023 East (hereafter NGC 7023 E) is located about 170′′ east of the star (Berné et al. 2007). Pilleri et al. (2012b) showed a clear difference in the spatial distribution of PAH⁰ and PAH⁺ in these three PDRs. Joblin et al. (2010) presented the first results from HIFI observations along a cut through NGC 7023 NW and S, suggesting that both the [C ii] emission and the aromatic infrared-band (AIB) emissions in the MIR arise from the regions located in the transition zone between atomic and molecular gas, providing new insights into the importance of the PAH charge evolution in the energetic studies of PDRs.

2.1.4. Carina

The Carina Nebula is a massive star-forming region complex with 65 O-type stars at a distance of 2.3 kpc (Smith 2006), which provides the nearest example of a very massive star-forming region and has been observed at many wavelengths (Smith & Brooks 2007). Trumpler 14 and 16 are the most prominent clusters. The PDR properties have been investigated using FIR and submilimeter emission lines (Kramer et al. 2008; Mizutani et al. 2004; Brooks et al. 2003). The region excited by Trumpler 14 has a FUV flux of G₀ = 7    ×    10³ and a high density of 2    ×    10⁵ cm^-3 (Kramer et al. 2008). Preibisch et al. (2011) present a large deep 870 μm continuum map, showing that the total mass is about 2 × 10⁵   M_⊙. Gaczkowski et al. (2013) recently published PACS and SPIRE continuum maps of the Carina Nebula complex to study embedded young stellar objects. In this study, we analyzed a PACS position in Carina North (hereafter Carina N), excited by Trumpler 14.

2.1.5. Mon R2

Mon R2 is a close-by ultracompact H ii (UCH ii) region at a distance of 850 pc, which comprises several PDRs that can be spatially resolved at both mm and IR wavelengths. The most intense UV source, which coincides with the infrared source IRS1, is located at the center of the cometary shaped UCH ii. The brightest PDR of Mon R2 is illuminated by an extremely intense UV field (G₀ ~ 5 × 10⁵; Rizzo et al. 2003), and its PAH emission peaks at about 20′′ northwest of IRS1 (Fig. 1). Several other PDRs lie around IRS1 and span different physical conditions such as temperature, density, column density, and UV field (e.g. Rizzo et al. 2003; Berné et al. 2009; Pilleri et al. 2012a). All PDRs have very bright MIR spectra, consisting of the emission from AIBs, H₂ rotational lines, and the continuum. The spatial extent of these PDRs (~0.03 pc) yields an angular size of ~8′′, as shown by PAH and H₂ emissions in the MIR (Berné et al. 2009). The shape of the MIR spectra varies at different PDRs, reflecting the photo-processing of the AIB carriers.

2.2. Far-infrared spectroscopy with Herschel

The FIR spectroscopic observations were performed with PACS as part of the WADI program. The observational parameters are listed in Table 2. All measurements are single pointed observations with the PACS spectrometer consisting of a 5 × 5 array of spatial pixels (spaxels) with a size of 9.4′′ each (see Fig. 1). We pipelined the data from level 0 to 2 using the Herschel data processing system (HIPE; Ott 2010) version 8. We converted the unit of the pipelined spectra, Jy/spaxel, to W m^-2 μm^-1 sr^-1 using the field of view of a spaxel (9.4′′ × 9.4′′). For the Chop/Nod observations, we checked the contamination at the OFF positions for [O i] 63 μm, 145 μm, and [C ii] 158 μm using the script provided in HIPE (ChopNodSplitOnOff.py), which provides separate spectra for the ON and OFF positions. Only the [C ii] 158 μm emission in the Horsehead shows significant OFF contamination. Since the northern side of the OFF positions shows stronger [C ii] 158 μm, and the emission from the southern side can be attributed to a diffuse component that is not associated with the region, we used only the OFF measurement of the southern side to obtain the flux from the source. For the unchopped grating-scan observations, some emission is detected at the OFF positions in [C ii] 158 μm (Carina N and Mon R2) and [O i] 63 μm (Carina N). The line intensities at the OFF positions are 4−7%, 2−5%, and 1−4% of that at the ON positions. We subtracted the OFF emission in Mon R2 but not in Carina N, because the detection of [O i] 63 μm in the Carina N OFF position indicates a real contamination from local dense clouds, while the [C ii] 158 μm emission traces more diffuse region and can be attributed to the large-scale diffuse Galactic [C ii] emission. To determine the continuum levels, we subtracted the OFF measurements in all observations to remove the telescope background.

The uncertainty of the absolute flux calibration is 11−12% for the line and continuum emission¹. A larger uncertainty is expected in weak sources for unchopped grating scans. The reproductivity in total absolute flux (telescope + source) of the unchopped grating scan is 4% (peak-to-peak). The relative uncertainty against the source flux can be expressed as 0.04 × F(telescope + source)/F(source), where F is the flux. For Mon R2, this means an uncertainty of <13% and <5% for the continuum at 105−180 μm and at the peak of [C ii] 158 μm. For Carina N, it is <10% and <5% at the peak of the [O i] 63 μm and [C ii] 158 μm emission lines and 20−50% for the underlying continuum level.

Fig. 2 Examples of [O i] 63 μm, 145 μm, and [C ii] 158 μm spectra in each PDR. The numbers labeled in figures after the name of the object show subregions (see Table 3).

We defined the region to be used in our study depending on the morphology of each source so that the widest range of physical conditions is included. Since our PACS observations are not fully sampled, we cannot exploit the entire spatial information. Instead, we selected a few typical areas in each source and extracted the [O i] and [C ii] line intensities and the continuum flux from PACS observations and combined them with the IRS results. These areas are listed in Table 3, shown as green boxes or circles in Fig. 1, and are explained in detail in the following. We extracted spectra based on the geometrical area and did not apply a beam-size correction. The FWHM of the PACS spectrometer beam at 63 μm and 158 μm is ~9′′ and ~11.5′′, respectively, and the uncertainty from the difference of the beam size is estimated to be ~15% (see Appendix A).

For Horsehead, Carina N, and Mon R2, all observations were made with almost the same position angle, i.e., the PACS spaxels observed almost the same area in the sky for different wavelengths (Fig. 1). In this case, we extracted the areas to be examined on the basis of the PACS spaxels. In the Horsehead, the spaxels are aligned to the ridge. We defined three areas; the first one lies at the western side of the ridge in the ionized gas (an average over four spaxels except for the north-end spaxel), the second area is along the ridge (over five spaxels), and the third area is in molecular gas, at the eastern side of the ridge (five spaxels). In Carina N, we selected three regions; a single spaxel toward a clump seen in the IRAC 8 μm map, a ridge-like east-west structure in the 8 μm map including this clump (five spaxels), and the southern 5 × 3 spaxels where the PAH emissions in the IRS spectra are prominent. In Mon R2, we selected three spaxels around PDR1 (following the definition in Berné et al. 2009), where the PAH⁺ fraction is suggested to be low, and toward the southwest and northeast inner edges of the PAH⁺ distributions shown in Berné et al. (2009). They correspond to the three green boxes in Fig. 1f from southwest to northeast.

For Ced 201, NGC 7023 E, and NGC 7023 NW, the [O i] 63 μm line was observed in different seasons of the year than the [O i] 145 μm and [C ii] 158 μm observations, and the difference of the position angle is significant. Therefore, we defined circles as areas to be studied, and took a weighted mean of PACS spectra based on the geometrical overlap between the defined circles and the PACS spaxels. For Ced 201 and NGC 7023 E, two circles were defined with a common center for both PACS footprints and diameters of 9.4′′ and 28.2′′. For NGC 7023 NW, we selected three typical regions that trace different ionization fractions of PAHs based on the analysis in Pilleri et al. (2012b). The areas were defined by circles with 12′′ diameter, at the interface with strong PAH emission, toward the molecular cloud, and in the cavity toward the exciting star (Fig. 1).

After extracting the spectra in each region as described above, we obtained the [O i] 63 μm and [C ii] 158 μm intensities, as well as [O i] 145 μm when available, by a Gaussian fit after linear baseline subtraction (Fig. 2, Table 3). In general, the ratio of [O i] 63 μm/[C ii] 158 μm traces the density of the PDRs (Röllig et al. 2006). Among our PDRs, NGC 7023 E and Ced 201 show a low ratio (<1), indicating a low density, and Mon R2 has the highest ratio (≳4). This trend is roughly consistent with the density estimate in previous studies (Table 4, Appendix B). The [O i] 145 μm/63 μm is >0.09 for all targets and 0.3 at maximum, which exceeds the PDR model prediction for the relevant physical parameter ranges (Röllig et al. 2006; Kaufman et al. 1999). Possible reasons are high optical depths of the [O i] 63 μm emission, a suprathermal population of [O i] by collision with H₂, and foreground absorption in [O i] 63 μm (Liseau et al. 2006). Although the optical depth of the [O i] 63 μm line is taken into account by the PDR models for simple geometries, the overlapped several PDR clumps along the line-of-sight or the edge-on geometry can cause more significant self-absorption of [O i] 63 μm (Habart et al. 2003; Okada et al. 2003), and the non-local calculation suggests a higher [O i] 145 μm/63 μm ratio (Elitzur & Asensio Ramos 2006). The foreground absorption cannot be quantified in our PDRs because of the lack of velocity-resolved [O i] 63 μm observations. However, the HIFI observations of the velocity-resolved [C ii] 158 μm imply only minor effects toward our PDRs.

2.3. Mid-infrared spectroscopy with Spitzer

We analyzed the MIR spectra observed with the Short-Low (SL1 and SL2) module of IRS onboard Spitzer except for NGC 7023 E, for which we used the ISOCAM highly-processed data product (Boulanger et al. 2005). We used the same data cubes as in Pilleri et al. (2012b) for the Horsehead, Ced 201, and NGC 7023, and the data reduction is described in that paper as well. For Mon R2, the IRS data and their reduction are shown in Berné et al. (2009). For Carina N, we collected the IRS observations from the Spitzer data archive and analyzed them using CUBISM (Smith et al. 2007). To estimate the continuum, Short-High (SH) and Long-Low (LL) spectra were also reduced and extracted when available. Their intensities are scaled to match the SL spectra, then the whole spectra are scaled to match the photometric data; IRAC 8 μm and MIPS 24 μm for Horsehead, Ced 201, and NGC7023 E, and IRAC 8 μm for NGC7023 NW and Carina N. For Mon R2, the absolute flux scaling was not applied because IRAC 8 μm is saturated, and the available spectral range does not cover the MIPS 24 μm band. The correction should be less significant in Mon R2 because it is the brightest source in our samples. Applying this correction does not affect any trend in the following results and changes none of our conclusions. We extracted a spectrum from the same area as for the PACS spectra described in the previous subsection.

3. Analysis

3.1. Total infrared intensity

Fig. 3 Continuum fit at the same positions as in Fig. 2. The black lines are the observed spectra, while the red and blue lines are the models used to estimate the TIR. The blue asterisks are the data points that are used for fitting the FIR thermal dust emission (see text). The blue curve in the FIR is the fit with the thermal dust model, that in the MIR indicates either the direct integration, the linear fit, or the assumed constant value. The red curve below 14 μm shows the fit described in Sect. 3.2. The black asterisk in NGC 7023 NW is the MIPS 24 μm flux.

Observationally, ϵ_pe is measured as ([O i] + [C ii]) / TIR, where TIR is the total infrared flux (3−1100 μm; Dale & Helou 2002), representing the fraction of the input energy that is converted into the cooling lines. The total FIR flux, 42.5−122.5 μm (Helou et al. 1988) or the integrated flux of the thermal emission by large dust grains, is often used as a tracer of the input energy. Here we used the TIR instead, which is the emission from all dust grains, including PAHs and VSGs, to represent the total input energy. As dominant cooling lines, we summed the intensities of [O i] 63 μm, 145 μm, and [C ii] 158 μm, except for Carina N, where the [O i] 145 μm was not observed, and we used ([O i] 63 μm + [C ii] 158 μm). A second PACS observation 1.3′ away in Carina N, which is not used in this study because it does not overlap with the IRS observations, shows an [O i] 145 μm/63 μm intensity ratio of 0.05−0.15. If our position has the same ratio, neglecting the [O i] 145 μm emission underestimates ϵ_pe by 8% at most.

To estimate TIR as the energy input tracer we need to integrate the grain emission of 3−1100 μm. The contribution from 3−5.5 μm is negligible because of the low intensity at 5.5 μm in all regions. To obtain the energy in the range of 5.5−14 μm, we fit the IRS spectra using the procedure described in Pilleri et al. (2012b) and integrated the resulting fit (see Sect. 3.2).

For λ > 14 μm, the available data and the quality are not uniform, and we fine-tuned approaches from region to region to estimate the integrated flux (see Fig. 3), although the basic approach was the same; we fit the FIR spectra from PACS observations with the thermal dust model, extrapolated the fit in the MIR to longer wavelengths, and connected the two fits at the intersection point.

In Horsehead, Ced 201, and NGC 7023 E, MIR spectra are available up to 35 μm. We estimated the integrated flux for 14−28 μm by direct integration excluding the strong emission lines for Ced 201 and NGC 7023 E, and by a linear fit to the spectra in the Horsehead because of the lower signal-to-noise ratio (S/N). The 28−32 μm range is fitted by a linear function, and the results are extrapolated to longer wavelengths. In NGC 7023 NW and Carina N, IRS/SH data are available for 14−19 μm, which are directly integrated, and we assumed a flat spectrum from 19 μm up to the wavelength where the big-grain emission, peaking in the FIR, exceeds this level. For Mon R2, we assumed that the flux at λ > 14 μm is constant (Fig. 3).

For the PACS Chop/Nod LineSpec observations (Table 2), the continuum levels are available around the emission lines. They are indicated as blue asterisks in Fig. 3. In NGC 7023 NW and Mon R2, the full spectral data from the PACS RangeSpec observations at >70 μm, after excluding strong emission lines, were used to fit the FIR continuum. In Carina N, the uncertainty of the absolute continuum flux is too large (see Sect. 2.2). Therefore we obtained the photometry data at blue (70 μm) and red (160 μm) bands from the Herschel data archive. We fit the FIR continuum using the dust emissivity of Ossenkopf & Henning (1994), for a gas density of 10⁶ cm^-3 and grains with thin ice mantles, considering the temperature and the column density as free parameters. Since we discuss only the integrated IR flux, the choice of the dust emissivity is not critical. Where this fit was exceeded by the extrapolation from the MIR spectrum, we connected these two at the intersection point. The final fits are shown as blue lines in Fig. 3.

Since the PACS observations are not fully sampled and the IRS spectra were extracted just using the geometrical area, the difference of the PSF could affect the resulting fit of the spectral energy distributions (SEDs). However, the uncertainty decreases for larger regions (average over several spaxels) and it does not change the general trend of ϵ_pe (Sect. 4; Fig. 5).

In NGC 7023 NW, the MIPS 24 μm flux is available but was not used to scale the IRS spectra because the spectral range of the IRS does not cover the MIPS 24 μm band. Figure 3 shows that the MIPS 24 μm flux matches the assumed constant continuum flux.

We translated the obtained TIR into the impinging FUV flux. By assuming that the entire FUV energy is absorbed by dust grains and reradiated as infrared flux, we estimate G₀(TIR) as 4π × TIR/(1.6 × 10^-6) (Table 4). The systematic uncertainty of TIR is ~15%, which comes from difference of the beam sizes at different wavelengths (Sect. 2.2). In this formula we assumed a spherical geometry. When we resolve an edge-on PDR with a shorter thickness compared to the length of the line-of-sight, G₀(TIR) is overestimated (Meixner et al. 1992). Another uncertainty of G₀(TIR) is the contribution of photons with an energy outside of 6 eV < hν < 13.6 eV to the thermal heating of the dust. For later B-type stars, the contribution of <6 eV becomes significant, which causes the overestimate of G₀(TIR). The SED of a B9 star by Castelli & Kurucz (2004) indicates that only 1/4−1/5 of the total energy is emitted in the range of 6 eV < hν < 13.6 eV from the star. We did not apply the correction and treated G₀(TIR) as an upper limit. Detailed comparisons of G₀(TIR) with previous studies and other diagnostics are described in Appendix B, and the final G₀ estimates, which are used in the following analysis, are also listed in Table 4.

3.2. Ionization of PAHs

Fig. 4 Fit of the MIR spectra using the PAHTAT procedure (Pilleri et al. 2012b) at the same positions as in Fig. 2, assuming the mixed extinction model and RV = 3.1. The black line shows the observed spectrum, the red line is the fitted spectrum, the blue line represents the continuum, green, orange, light-green and light-blue show the PAH0, PAH+, PAHx, and eVSG components.

To derive the contribution of PAH⁺ in each region, we applied the PAHTAT (PAH Toulouse Astronomical Templates) procedure described in Pilleri et al. (2012b). This procedure fits MIR spectra using a minimal set of template spectra. The PAH-related templates comprise PAH⁰, PAH⁺, larger ionized PAHs (named PAH^x, see Joblin et al. 2008), and evaporating very small grains (eVSGs). In the MIR the PAH⁰ and PAH⁺ templates are characterized by a strong difference in the relative strength of their 7.7 μm and 11.3 μm band features; PAH⁺ has a stronger 7.7 μm band. The PAH^x population consists of large ionized PAHs, which was introduced to provide a better fit to planetary nebula spectra: its MIR spectrum is similar to PAH⁺, but the 7.7 μm band is shifted to longer wavelengths (Joblin et al. 2008; Pilleri et al. 2012b). eVSGs are an intermediate population between PAHs and classical VSGs. They present both a broad band and continuum emission in the MIR range. The PAHTAT procedure also allows simultaneous fitting of the gaseous emission lines, underlying continuum, and the extinction by dust grains along the line-of-sight. More details of the fitting tool can be found in Pilleri et al. (2012b). We assumed a linear continuum except for Mon R2, where a continuum of two slopes connected at 10 μm was adopted since big grains contribute significantly to the MIR continuum because of the strong UV field. To explore the parameter space of the fitting, we used all combinations of the following different assumptions: the dust extinction curves with R_V = 3.1 and 5.5, geometries in which the absorbing materials are placed in the foreground or are fully mixed with the emitting materials, and since the presence of PAH^x in PDRs with an intermediate UV radiation field strength is doubtful, we also tested the results including or excluding the PAH^x component from the fit. The difference between these different assumptions is included into the uncertainty of the results, although they are typically smaller than the errors estimated by the discrepancy between the observed spectra and the model (see below). The fitted spectra are shown in Fig. 4, and the derived fraction of each PAH and eVSG contribution to the total band flux is shown in Table 5.

In the following we define the fraction of positively ionized PAHs by the integrated intensity ratio as , where the integration is made over the wavelength from 5.5 μm to 14 μm. The uncertainty of this fraction can be attributed to two different effects: (1) fitting errors and (2) systematic errors in the assumptions within PAHTAT, i.e., due to intrinsic uncertainties in the templates. To estimate the first contribution, we integrated the absolute discrepancy of the observed spectra from the model. We assumed that this residual consists of the uncertainty of the band emissions of small grains (PAH⁰, PAH⁺, PAH^x, and eVSGs) and the continuum emissions; ^∫(residual) = Δ(band) + Δ(continuum), and the ratio of the uncertainty Δ(band)/Δ(continuum) is the same as the ratio of the integrated flux . Thus, the uncertainty of the integrated band emission is defined as . Then we attribute it to the uncertainty of the PAH⁺ flux, i.e., the final uncertainty of f(+) is expressed as Δ(band)/. This is conservative, because it assumes that the entire uncertainty in the band emission comes from that of PAH⁺.

The second uncertainty in the obtained f(+) comes from the templates for PAH⁰ and PAH⁺. The templates are constructed based on observational data and a mathematical blind signal separation (BSS) method (Pilleri et al. 2012b). The spectral properties of the BSS-extracted spectra, templates, and their assignment to PAH⁺ and PAH⁰ have been discussed in Rapacioli et al. (2005), Berné et al. (2007), Joblin et al. (2008), and Berné et al. (2009). Recently, Rosenberg et al. (2011) compared the BSS extraction to the theoretical spectra with the density functional theory and showed good agreement, although this is for a spectral range of 10−19.5 μm. A precise determination of the uncertainty is difficult, however, given the nature of the extraction method and our limited knowledge of the PAH populations that exist in space. Using a Monte Carlo approach, Rosenberg et al. (2011) showed in their BSS extraction a 1σ uncertainty of typically 10% and up to 30% in some parts of the spectrum (see their Fig. 3). However, this uncertainty in the templates systematically propagates to the uncertainty of f(+) and does not change the trend and our conclusion discussed below. Therefore, we present our errorbars in the following figures based on the fitting error described above.

4. Results and discussion

4.1. Relation between ϵ_pe and f(+)

Fig. 5 ϵpe = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR versus the fraction of ionized PAHs (f(+)).

Figure 5 shows the relation between the ϵ_pe = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR and the fraction of positively ionized PAHs f(+). Both ϵ_pe and f(+) strongly vary among the regions; ϵ_pe varies from 0.1% to 0,9%, and f(+) from 0 (+11)% to 87 (±10)%. All positions with a high f(+) show a low ϵ_pe, and all positions with a high ϵ_pe show a low f(+). This trend supports the theoretical expectation in which ϵ_pe decreases when grains are positively ionized, because the energy required to eject electrons from positively charged grains is higher than that from neutral grains. Here we take the ionization of PAHs as a gross indicator of the positive charging of grains in general, which is generally quantified by the charging parameter (γ) as shown in Sect. 4.2. Figure 5 directly compares of the observationally derived ϵ_pe and f(+), independent of the use of PDR models to quantify the physical conditions, over a wide range of the physical properties of PDRs.

Using ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR as a tracer of ϵ_pe contains several assumptions. We neglect the heating by the collisional de-excitation of vibrationally excited H₂. The H₂ de-excitation heating becomes a dominant heating process in a dense PDR with a low UV field (Röllig et al. 2006). Among our PDRs, only G₀ and the upper value of the density of the Horsehead are close to the regime where the contribution from the photoelectric heating and H₂ heating is comparable. Those in other PDRs indicate that the photoelectric heating is dominant. We also neglect the cooling by H₂ emission. This is partly justified because the H₂ emissions as a consequence of the excitation by the UV-pumping and the H₂ formation do not contribute to the estimate of the photoelectric heating efficiency. Habart et al. (2011) showed intense H₂ lines in the Horsehead (0−0 S(0) to S(3) and 1−0 S(1)), which are, as a sum, comparable with the ([O i] + [C ii]) intensity. In NGC 7023 E and Ced 201, the sum of the H₂ pure rotational emissions S(0) to S(3) (Habart et al.2011 for NGC 7023 E, and from the IRS spectra in the Spitzer archive for Ced 201) is ~80% and ~60% of the ([O i] + [C ii]) intensity, respectively. Neglecting these H₂ emissions may underestimate ϵ_pe by a factor of ≲2. In other regions, the H₂ pure rotational emissions from the IRS spectra in the Spitzer archive show the intensity of ≲20% of ([O i] + [C ii]).

As mentioned in Sect. 2.2, the [O i] 63 μm emission is indicated to be optically thick in most targets. Therefore, the optical-depth-corrected [O i] 63 μm intensity is higher than the observed values. However, we observe the emission that leaves the cloud, and this contributes to the net cooling of the cloud. This justifies the direct use of the observed ([O i] + [C ii]) to express the gas cooling, except if the foreground absorption is significant, which is unlikely for our PDRs, as mentioned above.

As discussed for G₀(TIR) (Sect. 3.1), photons with an energy of <6 eV contribute to TIR, but not to the photoelectric heating. Therefore, estimating ϵ_pe with the observed TIR provides a lower limit. For the exciting source(s) of spectral types earlier than early B, the effect is not significant in terms of the absolute factor and the variation between different spectral types. For Ced 201, which has an exciting source of B9.5, we may underestimate ϵ_pe by a factor of a few. This uncertainty, however, does not change the overall trend in Fig. 5 and our conclusions.

Our definition of ϵ_pe assumes a spherical geometry, where both line and continuum emissions radiate isotropically. If one assumes an edge-on geometry, where the continuum emission can escape through the molecular side while the line emissions, especially [O i] 63 μm, becomes completely optically thick and can be emitted only from the front side, ϵ_pe should be divided by 2 (Tielens 2005).

4.2. Comparison of ϵ_pe with theory

Fig. 6 ϵpe = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR versus the charging parameter (γ = G0T1/2/ne). Solid line is the theoretical calculation from Bakes & Tielens (1994), dashed lines are that from Weingartner & Draine (2001b) for RV = 3.1 (105bC = 0.0 and 6.0 for the lower and upper line), dotted lines are for RV = 5.5 with Case A (105bC = 0.0 and 3.0 for the lower and upper line; see text).

Fig. 7 ϵPAH = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/(PAH + [O i] 63 μm + [O i] 145 μm + [C ii] 158 μm) versus the charging parameter (γ). The black curve depicts our calculations using 25 different size distributions from Weingartner & Draine (2001a), which looks degenerated (see text).

Fig. 8 Fraction of PAH+ (f(+)) versus the charging parameter (γ).

Theoretically, ϵ_pe can be expressed as a function of the charging parameter, defined as γ = G₀T^1/2/n_e, where T is the gas temperature and n_e is the electron density. γ is proportional to the ratio of the ionization and recombination rate (Bakes & Tielens 1994). We calculated γ from G₀, n_H, and T listed in Table 4. They were carefully estimated in each subregion (Appendix B) to be independent of the existing PDR models. We estimated n_e by assuming that most electrons are provided by carbon ionization and all carbon atoms are singly ionized, i.e., n_e = x(C)n_H, where x(C) = 1.6 × 10^-4 (Sofia et al. 2004) is the elemental carbon abundance. In Fig. 6 we show the relation between ϵ_pe and γ. There is a trend that ϵ_pe decreases when γ increases. Bakes & Tielens (1994) theoretically estimated ϵ_pe based on the physics of the photoelectric effect on dust grains. Weingartner & Draine (2001b) modeled ϵ_pe with some improved physical parameters and various size distributions of dust grains in Weingartner & Draine (2001a). Dashed and solid lines in Fig. 6 show some of their estimates, representing the strongest variation of ϵ_pe depending on the dust size distribution. Dense regions characterized by R_v = 5.5 (dotted lines) have a lower heating efficiency. The parameter b_C is the total C abundance in the log-normal populations of the size distribution, and the model with a larger b_C shows a stronger enhancement at the size of <10 Å (Weingartner & Draine 2001a). Case A is the case without the constraint on the total grain volume (Weingartner & Draine 2001a). Although the observed trend that ϵ_pe decreases when γ increases is reproduced by models, the observation indicates a weaker dependence of ϵ_pe on γ, and all models overestimate ϵ_pe at low γ. In the following, we focus on the photoelectric heating only on PAHs.

When we consider ϵ_pe on PAHs, ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR is not an appropriate definition, because the TIR expresses the total energy emitted by all dust grains. The total energy that PAHs absorb is converted into either MIR AIB emission or cooling line emission through the photoelectric effect on PAHs. A fraction of the energy (generally taken to be 0.5) remains behind as electronic excitation energy after the photoelectric effect (Tielens 2005), which also results in MIR AIB emission. Therefore ϵ_PAH is defined as ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/(PAH band emission + [O i] 63 μm + [O i] 145 μm + [C ii] 158 μm). The trend of ϵ_PAH against γ (Fig. 7) is very similar to that reported in Fig. 6.

Following Tielens (2005), we computed a simple theoretical estimate of ϵ_PAH. For a given PAH containing the number of carbon atoms of N_c, (1)where f(Z = 0) is the neutral fraction and IP is the ionization potential. With a typical photon energy of 10 eV and an ionization potential of 7 eV, this equation becomes (2)The neutral fraction is given by (3)where γ₀ is the ratio of the ionization rate over the recombination rate, which for small PAHs is given by (4)We substitute Eq. (4) into Eq. (3) and then into Eq. (2), and integrate over the size distribution from Weingartner & Draine (2001a) up to N_c = 100. In contrast to ϵ_pe, in which a different size distribution gives a large difference (Fig. 6), ϵ_PAH is insensitive to the adopted size distribution. In Fig. 7, the results with 25 different size distributions in Weingartner & Draine (2001a) are plotted, but they look degenerated. Although the match between the theoretical estimates and the observed values is not excellent, it is much better than in Fig. 6. This result suggests that the photoelectric heating is dominated by PAHs.

In Fig. 8, the relation between the fraction of PAH⁺ (f(+)) and γ is shown. PAHs are almost neutral at γ < 10³ and almost fully ionized at γ > 10⁶, and there is a transition in-between. Since γ is defined by the environment properties and it is proportional to the ratio of the grain ionization and recombination rate, Fig. 8 confirms that the fraction of PAH⁺ is a good indicator for the positive charging of grains in general. This trend is also consistent with the correlation between the intensity ratio of PAH(6.2 μm)/PAH(11.3 μm) and γ presented in Galliano et al. (2008).

5. Summary

We analyzed Herschel/PACS and Spitzer/IRS spectroscopic observations in six PDRs and showed that the photoelectric heating efficiency (ϵ_pe) is lower in regions with a large fraction of positively ionized PAHs (f(+)). Based on examining the photoelectric heating efficiency on PAHs, we found a dominant contribution of PAHs to the photoelectric heating.

Our PDR sample covers a wide range of physical conditions (100 ≲ G₀ ≲ 10⁵) and provides a good test case for investigating the relation between ϵ_pe and the charge state of PAHs. We estimated ϵ_pe as ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR. ϵ_pe varies in our PDRs between 0.1 and 0.9%. f(+) was obtained from a fit of the MIR spectra with a set of template spectra representing PAH-related species and varied from 0% (+11%) to 87% (±10%). All positions with a high f(+) show a low ϵ_pe, and all positions with a high ϵ_pe show a low f(+). This trend supports a scenario in which a positive grain charge results in a decreased heating efficiency. The theoretical estimate of ϵ_pe shows a stronger dependence on the charging parameter (γ) than the observed ϵ_pe reported in this study, and overestimates ϵ_pe at low γ. The photoelectric heating efficiency on PAHs, ϵ_PAH = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm) / (PAH band emission + [O i] 63 μm + [O i] 145 μm + [C ii] 158 μm), shows a much better match between the observations and the theoretical estimates, indicating a dominant contribution of PAHs on the photoelectric heating. PDR models that fully account for the relative contribution of different PAH and eVSGs populations are needed. Velocity-resolved observations of [O i] 63 μm in the future with for instance SOFIA/upGREAT (Heyminck et al. 2012) will enable us to investigate the foreground absorption and the optical depth effect in detail (cf. Boreiko & Betz 1996).

Acknowledgments

PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). HIFI has been designed and built by a consortium of institutes and university departments from across Europe, Canada and the United States under the leadership of SRON Netherlands Institute for Space Research, Groningen, The Netherlands and with major contributions from Germany, France and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronómico Nacional (IGN), Centro de Astrobiología (CSIC-INTA); Sweden: Chalmers University of Technology − MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University − Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC. We thank the Herschel helpdesk and the Spitzer helpdesk for their support in analyzing data. We thank the referee for useful suggestions that greatly improved the paper. Part of this work was supported by the German Deutsche Forschungsgemeinschaft, DFG, project number SFB956 C1, by the Spanish program CONSOLIDER INGENIO 2010, under grant CSD2009-00038 Molecular Astrophysics: The Herschel and ALMA Era (ASTROMOL), by CNES, by a Ramón y Cajal research contract, and by the Spanish MICINN through grants AYA2009-07304 and CSD2009-00038.

References

Appendix A: Relative flux uncertainty between 63 μm and 158 μm caused by the difference of the beam size

Fig. A.1 Black curves show the simulated r, defined as Eq. (A.1), as a function of the position of the point source for the beam size of a) 9′′ and b) 11.5′′. The relative contribution of the extended source compared to the height of the point source at 63 μm is 0.0, 0.02, 0.1, 0.3, and 0.7 from upper to lower curves. Shaded areas represent the observed values with errors of ±3σ (blue lines in Fig. A.2). c) The flux ratio of 9′′ beam and 11.5′′ beam at spaxel 2 as a function of the position of the point source with the same model. Shaded areas show ranges when r in a) or b) matches the observed value.

Fig. A.2 Histograms of the observed r for [O i] 63 μm and [C ii] 158 μm. Red lines show the fit with a Gaussian profile, and blue lines show the center and ±3σ of the fitted Gaussian.

We estimate the uncertainty of the relative flux in one spaxel at different wavelengths that arises from the different PSF and the unknown spatial structure of the sources by considering the combination of two extremes, pure point sources and flat extended emission, and by comparing the contrast between neighboring spaxels for these combinations with the contrast actually measured in the observations.

We consider three spaxels 1−3 along a line whose centers are located at (0.5, 0), (1.5, 0), and (2.5, 0) in units of a spaxel size (9.4′′). We put a point source represented by a Gaussian profile with the FWHM of 9′′ and 11.5′′ for 63 μm and 158 μm, respectively, centered at various positions between (0.5, 0) and (2.5, 0), and calculate the mean flux ratio of the adjacent spaxels as (A.1)where F is the flux falling in one spaxel. r is lowest when the center of the point source is located at (1.5,0), which is the center of spaxel 2, and the minimum r is 7.1 and 4.0 for 63 μm and 158 μm, respectively (upper lines in Figs. A.1a,b). On the other hand, we derive the observed r of the [O i] 63 μm and [C ii] 158 μm line intensities from the original 5 × 5 PACS spaxels for all regions. The histogram of the distribution of r is well represented by a Gaussian with a small number of outliers at high r values (Fig. A.2), and the fitted Gaussian has a center of 1.1 for both lines and σ of 0.2 and 0.1 for 63 μm and 158 μm. Even the upper 3σ (1.7 and 1.4) is smaller compared to the above model with a point source, indicating that the contribution of extended emission is significant in our regions. Therefore, we modify the model by adding a constant extended component to all spaxels, and estimate r, which is shown as black lines in Figs. A.1a and b with varying constant values. Then we estimate F(spaxel 2, 63 μm)/F(spaxel 2, 158 μm), which expresses the uncertainty of the relative flux at 63 μm and 158 μm, for each model (black lines in Fig. A.1c). The observed r with 3σ errors are shown as shaded areas in Figs. A.1a and b, and the corresponding ranges are shaded in Fig. A.1c, which is numerically calculated by surveying the constant emission value and determining the position range of the point source, where r is in the observed range for each constant emission value. It reaches to within 15% of unity. Therefore, we take 15% as the uncertainty of the relative flux between 63 μm and 158 μm caused by the beam size difference.

Appendix B: Estimate of the charging parameter

Here we describe the estimate of G₀, the gas density n_H, and the gas temperature T to be used for deriving the charging parameter γ in individual subregions (Table 4). The rotational temperature derived from the low-J H₂ pure rotational emissions up to S(4) is adopted as T in all regions. G₀ and n_H affect γ linearly and we estimate them independent of an existing PDR model. These properties are representative of the warm PDR surface where most of the H₂ rotational emission and PAH emission comes from.

Appendix B.1: Horsehead

Habart et al. (2005) estimated G₀(star) ~ 100 at the PDR interface based on the effective temperature of the exciting star and assuming geometrical dilution at the projected distance. Under this assumption, our three subregions are close enough to each other to make no difference for the G₀(star) between them. G₀ derived from TIR, G₀(TIR), is listed in Table 4. We adopt both G₀(star) and G₀(TIR) for the possible range in individual subregions. For the gas temperature T, we use the range of the estimate from the intensity ratio of H₂ S(2)/S(0), S(3)/S(1), and S(4)/S(2) in Habart et al. (2011).

Abergel et al. (2003) showed that the lower limit of the density behind the filament is ~2    ×    10⁴ cm^-3 from the infrared brightness profile. Habart et al. (2005) modeled the spatial distribution of H₂, PAH, CO and 1 mm continuum emissions using a PDR model and suggested the density gradient from n_H = 10⁴ cm^-3 in the H₂ emitting region to n_H = 2 × 10⁵ cm^-3 in the inner cold molecular layers. Since the region we are interested in is the H₂-emitting gas, we still include the n_H value in the H₂-emitting region even for subregion 3, i.e. we adopt n_H = (0.1−2) × 10⁵ cm^-3 for subregions 2 and 3. n_H ~ 10⁴ cm^-3 at the interface is also confirmed by the pressure equilibrium with the ionized gas. In the ionized gas, n_e is estimated to be 100−350 cm^-3 with T_e = 7500 K (Compiègne et al. 2007). Assuming the pressure equilibrium, 2n_eT_e = n_HT, gives n_H = (0.3−2) × 10⁴ cm^-3. Therefore, we adopt n_H = (0.3−2) × 10⁴ cm^-3 for subregion 1.

Appendix B.2: Ced 201

Kemper et al. (1999) estimated G₀ = 200 by comparing the observations with a PDR model in units of the average interstellar radiation between 2 eV and 13.6 eV. On the other hand, Young Owl et al. (2002) derived G₀ = 300 by the infrared continuum emission, with a geometry correction factor of 2. We adopt G₀ = 200 as a lower limit and G₀(FIR) as an upper limit for each subregion. The fraction of eVSGs determined by PAHTAT is correlated to G₀ (Pilleri et al. 2012b) and the adopted G₀ for each subregion is also consistent with it.

For n_H, Young Owl et al. (2002) derived n_H = 4 × 10² cm^-3 from the lower limit of [C ii] 158 μm/[O i] 63 μm using a PDR model. Kemper et al. (1999) estimated n_H of (5 ± 1) × 10³ cm^-3 by simple excitation models of CO and ¹³CO emissions. They also modeled the emission of CO, C, C⁺, CS, and HCO⁺ using a PDR model, which suggests n_H of 1.2 × 10⁴ cm^-3. We use all these ranges for the two subregions. T is estimated as ~330 K from the ratio of H₂ S(1) and S(3) emissions (Kemper et al. 1999).

Appendix B.3: NGC 7023 E

Pilleri et al. (2012b) modeled the spatial distribution of the MIR AIB emission, assuming a spherical shell geometry and the energy balance taking into account, and derived the spatial variation of n_H and G₀ along a cut in NGC 7023 E. At the PDR front, G₀ is calculated to be 250 based on the spectral type of the star, geometrical dilution and assuming an attenuation of A_v = 1.5 around the star (Pilleri et al. 2012b). The density at the PDR front is derived as 1.4 × 10³ cm^-3. Around the MIR AIB emission peak, which corresponds to the regions in this study, G₀ is estimated to be 120−170 and n_H = (0.5−1) × 10⁴ cm^-3. On the other hand, G₀(TIR) is 970 and 900 for two subregions. We adopt 120 as a lower limit of G₀, and the corresponding G₀(TIR) as an upper limit for each subregion, and n_H = (0.1−1) × 10⁴ cm^-3 for both subregions. For the gas temperature T, we use 258−370 K from the intensity ratio of H₂ S(2)/S(0) and S(3)/S(1) (Habart et al. 2011).

Appendix B.4: NGC 7023 NW

The same modeling as for NGC 7023 E was made by Pilleri et al. (2012b). Our three subregions are not exactly on their cut because we chose the subregions to maximize the variation of the PAH⁺ fraction. Subregion 2 is located at the PDR interface. Pilleri et al. (2012b) estimate G₀(star) =   2600 at the PDR front, whereas G₀(TIR) = 7700 at this subregion. We adopt G₀ = 2600−7700 for subregion 2. Subregion 1 is located closer to the star than the PDR interface. The scaling of G₀(star) by the projected distance gives G₀ = 10⁴, higher than G₀(TIR), possibly because either the real distance to the star is larger than the projected distance, and/or the assumption that all UV radiation is absorbed and converted into the IR emission underestimates G₀ in such a cavity because many UV photons pass the region. Nevertheless, we conservatively cover both values, i.e., G₀ = (0.5−1) × 10⁴, for subregion 1. The adopted G₀ in subregions 1 and 2 is also consistent with the correlation between the fraction of eVSGs and G₀ from Pilleri et al. (2012b). For subregion 3, we use this correlation to derive G₀ (see Eq. (5) in Pilleri et al. 2012b), which results in G₀ = 500−3100. Together with G₀(TIR) = 5400, we adopt G₀ = 500−5400 for subregion 3.

In Pilleri et al. (2012b), the modeled density profile sharply increases from 1.1 × 10³ cm^-3 to 2    ×    10⁴ cm^-3 within ~8′′. Fuente et al. (1996) showed high-density filaments of a few 10⁵ cm^-3 based on HCO⁺ observations. Therefore, we adopt 10³−2 × 10⁵ cm^-3 for subregion 2, and (0.2−2) × 10⁵ cm^-3 for subregion 3. Berné & Tielens (2012) examined n_H in the cavity in detail from several diagnostics and suggested 150 ± 100 cm^-3. We take this as n_H for our subregion 1. For the gas temperature T, we use 430−450 K from the intensity ratio of H₂ S(3)/S(1) and S(4)/S(2) (Fuente et al. 2000).

Appendix B.5: Carina N

We estimate G₀(star) by computing the contribution from all OB stars of Trumpler 14 listed in Smith (2006), which gives (7−8) × 10³ for the three subregions. On the other hand, G₀(TIR) is (1.5−1.6) × 10⁴. The point-source catalog of the Wide Field Infrared Survey Explorer (WISE; Wright et al. 2010) contains several sources with a flux stronger at 4.6 μm than at 3.4 μm close to our region, which indicates a contribution from embedded protostars. We consider both G₀(star) and G₀(TIR) as a possible range in individual subregions.

Kramer et al. (2008) modeled CO and [C i] emissions in the Carina N region and suggested n_H of 2 × 10⁵ cm^-3. This is consistent with a pressure equilibrium with the H ii region. Using the emission lines ratio from ions with similar ionization potential observed with IRS, we can estimate n_e assuming the elemental solar abundance and isothermal thin emission. [Ne ii] 12.8 μm/[S iii] 18.7 μm, [P iii] 17.9 μm/[Ne ii] 12.8 μm, and [Ar iii] 8.99 μm/[S iii] 18.7 μm give n_e of 10⁴, 2 × 10⁴, and 9 × 10³−10⁴ cm^-3, respectively, when T_e = 10⁴ K. Since the gas temperature T in PDRs is derived as T = 290−529 K from H₂ emission lines, n_H can be estimated as 3 × 10⁵−10⁶ cm^-3. We use 2 × 10⁵−10⁶ cm^-3 for all subregions.

Appendix B.6: Mon R2

We calculate G₀(star) from the contribution of IRS1-4 (see Fig. 1) with the luminosity listed in Henning et al. (1992), as (5.2−5.9) × 10⁴, 10⁵, and (1.1−1.2) × 10⁵ in the three subregions. While G₀ = 5 × 10⁵ is used to characterize the radiation field at the ionization front (Rizzo et al. 2003), our subregions do not cover the nearest regions of infrared sources to avoid unresolved complex spatial distributions. G₀(star) matches G₀(TIR) well. We consider both G₀(star) and G₀(TIR) as a possible range in individual subregions.

Rizzo et al. (2005) and Pilleri et al. (2013) performed an LVG analysis with C₃H₂ and C₂H at several positions in Mon R2 and derived n_H> 10⁶ cm^-3. Ginard et al. (2012) showed similar results from different molecules; the molecular hydrogen density of a few 10⁵ cm^-3 to 10⁶ cm^-3. These molecules traces the cold molecular gas with a temperature of T ~ 45 K (Giannakopoulou et al. 1997). On the other hand, Berné et al. (2009) derived n_H of (0.4−4) × 10⁵ from PDR model calculations using H₂ rotational emission lines, and the temperature derived from the H₂ rotational emission lines is T = 314−574 K. These n_H and T are consistent with a pressure equilibrium with the cold molecular gas. In the H ii region, the emission measure of 1.9 × 10⁷ pc cm^-6 and the geometrical mean diameter of 23.7′′ (Takahashi et al. 2000; Wood & Churchwell 1989) give an estimate of n_e = 1.4 × 10⁴ cm^-3. With T_e = 7600 K (Downes et al. 1975) and T = 314−574 K, the pressure equilibrium gives n_H = (4−7) × 10⁵ cm^-3, which is also compatible with the estimate from H₂ rotational emissions, although the pressure equilibrium is an inadequate assumption in an UCH ii like Mon R2. We adopt the estimate range from H₂ emissions, i.e., (0.4−4) × 10⁵ for all subregions.

All Tables

All Figures

Fig. 1 PACS footprints of the blue camera ([O i] 63 μm) and the red camera ([O i] 145 μm and [C ii] 158 μm) overlaid on IRAC 8 μm, except for f) Mon R2. (There IRAC 8 μm is saturated, and the total intensity of AIB bands derived from fitting of IRS spectra (see text) is shown.) Green circles or boxes show the areas from which our PACS and IRS spectra are extracted (subregions in Table 3). For Mon R2 f), the positions of four IR sources (IRS1−4) and PDR 1 from Berné et al. (2009) are marked.

In the text

	Fig. 2 Examples of [O i] 63 μm, 145 μm, and [C ii] 158 μm spectra in each PDR. The numbers labeled in figures after the name of the object show subregions (see Table 3).
In the text

Fig. 3 Continuum fit at the same positions as in Fig. 2. The black lines are the observed spectra, while the red and blue lines are the models used to estimate the TIR. The blue asterisks are the data points that are used for fitting the FIR thermal dust emission (see text). The blue curve in the FIR is the fit with the thermal dust model, that in the MIR indicates either the direct integration, the linear fit, or the assumed constant value. The red curve below 14 μm shows the fit described in Sect. 3.2. The black asterisk in NGC 7023 NW is the MIPS 24 μm flux.

In the text

	Fig. 4 Fit of the MIR spectra using the PAHTAT procedure (Pilleri et al. 2012b) at the same positions as in Fig. 2, assuming the mixed extinction model and RV = 3.1. The black line shows the observed spectrum, the red line is the fitted spectrum, the blue line represents the continuum, green, orange, light-green and light-blue show the PAH0, PAH+, PAHx, and eVSG components.
In the text

	Fig. 5 ϵpe = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR versus the fraction of ionized PAHs (f(+)).
In the text

Fig. 6 ϵpe = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/TIR versus the charging parameter (γ = G0T1/2/ne). Solid line is the theoretical calculation from Bakes & Tielens (1994), dashed lines are that from Weingartner & Draine (2001b) for RV = 3.1 (105bC = 0.0 and 6.0 for the lower and upper line), dotted lines are for RV = 5.5 with Case A (105bC = 0.0 and 3.0 for the lower and upper line; see text).

In the text

	Fig. 7 ϵPAH = ([O i] 63 μm + [O i] 145 μm + [C ii] 158 μm)/(PAH + [O i] 63 μm + [O i] 145 μm + [C ii] 158 μm) versus the charging parameter (γ). The black curve depicts our calculations using 25 different size distributions from Weingartner & Draine (2001a), which looks degenerated (see text).
In the text

	Fig. 8 Fraction of PAH+ (f(+)) versus the charging parameter (γ).
In the text

Fig. A.1 Black curves show the simulated r, defined as Eq. (A.1), as a function of the position of the point source for the beam size of a) 9′′ and b) 11.5′′. The relative contribution of the extended source compared to the height of the point source at 63 μm is 0.0, 0.02, 0.1, 0.3, and 0.7 from upper to lower curves. Shaded areas represent the observed values with errors of ±3σ (blue lines in Fig. A.2). c) The flux ratio of 9′′ beam and 11.5′′ beam at spaxel 2 as a function of the position of the point source with the same model. Shaded areas show ranges when r in a) or b) matches the observed value.

In the text

	Fig. A.2 Histograms of the observed r for [O i] 63 μm and [C ii] 158 μm. Red lines show the fit with a Gaussian profile, and blue lines show the center and ±3σ of the fitted Gaussian.
In the text