© ESO, 2015

1. Introduction

Understanding photodissociation regions (PDRs; also known as photon-dominated regions; for example Tielens & Hollenbach 1985b; Sternberg & Dalgarno 1989) is an essential part of comprehending the physics and chemistry of the interstellar medium (ISM). The classical definition of a PDR is the region in which the energy balance is dominated by far-UV photons with energies less than what is required to ionize hydrogen (13.6 eV). In practice this definition means that classical PDRs form the “skin” of the ionized H ii regions and represent the buffer zone between the fully ionized H ii region material and molecular clouds. Typically this means that they occur at A_V< 5 with respect to the cloud surface. In a broader sense, a large portion of all neutral ISM material exists under these conditions, underlining the importance of understanding the physics and chemistry of these regions.

The prototypical PDR, the Orion (Bar) region, has been extensively studied (e.g., Tielens & Hollenbach 1985a; Stacey et al. 1993; Hogerheijde et al. 1995; Jansen et al. 1995; Simon et al. 1997; Rubin et al. 2011; Boersma et al. 2012; Ossenkopf et al. 2013; Andree-Labsch et al. 2014; Nagy et al. 2015); and reviews of the whole complex by Genzel & Stutzki 1989; O’Dell 2001; Tielens et al. 1993; Habart et al. 2010). This has been largely due to two factors: a) its close proximity (450 pc, Genzel & Stutzki 1989) and b) the fact that the Orion Bar is viewed almost exactly edge on. This combination allows for the relatively easy separation of the emission from each region across the transition from ionized to PDR material and subsequent comparison with detailed PDR models (e.g., Tielens et al. 1993).

These models are highly sophisticated codes that include a detailed chemical network, as well as descriptions of relevant heating and cooling processes. As a function of the incident radiation field strength (usually denoted G₀¹) and density (or thermal pressure), the PDR models calculate the thermal equilibrium temperature and dominant atomic and molecular species in chemical balance. The models also predict the dominant emission line strengths that can be compared with observations. From this data the PDR models synthetic spectra are constructed and compared to those observed. It is now generally agreed that most PDRs are clumpy; that is, they take the form of dense clumps (n> 10⁵ cm^-3) embedded within an interclump medium 10 to 100 times less dense (Burton et al. 1990b; Meixner & Tielens 1993; Röllig et al. 2006). Such two-phase models can reproduce much of the observed atomic and molecular emission from PDRs.

Measurements of PDR diagnostics have been sparse in previous decades because the diagnostic features of PDR physical conditions predominantly emit in the far-infrared (FIR) spectral regime at wavelengths between ~30 and 1000 μm. This wavelength region is notoriously difficult to observe, because much of this range requires above-atmosphere observations, using space-based observatories such as ISO, and Spitzer, balloons like BICE (Nakagawa et al. 1998) or aircraft, such as the KAO (Kuiper Airborne Observatory) and SOFIA (Stratospheric Observatory for Infrared Astronomy) to cover the shortest wavelengths and ground-based sub-mm observatories to cover the longer wavelength emissions. In recent years, this has become much easier, because the Herschel Space Observatory (Pilbratt et al. 2010) was able to collect full spectra in the PDR-significant range.

In this paper we expand the sample of PDRs observed by Herschel using the PACS and SPIRE instruments (Poglitsch et al. 2010; Griffin et al. 2010), to include two Galactic PDRs generated by massive star-formation (IRAS 23133+6050 and S 106) in order to understand the physical conditions and compare their properties to those of other Galactic PDRs, including the prototypical Orion Bar. We obtained spectra of the important atomic PDR cooling lines ([O i] 63, 145 μm, [C ii] 158 μm, and [C i] 609 μm), as well as the CO and ¹³CO ladders from mid- to high-J. These observations allow us to probe various depths into the respective PDRs, with the atomic lines providing constrains on the warm, dense surface layers at low A_V (e.g., T = 100−1000 K, n = 10⁴−10⁵ cm^-3, A_V< 4) and the CO ladders probing the warm, low A_V regions and deeper, cooler regions at higher A_V (e.g., T = 10−100 K, n = 10⁴−10⁵ cm^-3, A_V> 4).

The nature of our targets and prior knowledge of them is outlined in Sect. 2. A description of the observations, data reduction procedures and the measurement of the spectra is given in Sect. 3. The results of this work, starting with the PDR diagnostics including the molecular tracers such as ¹²CO are discussed in Sect. 4 in terms of the physical conditions of these PDRs. Subsequently we compare these observations to the outputs of a sophisticated PDR model in Sect. 5. These results are discussed in terms of the geometry of the PDRs and the derived conditions with those of other Galactic PDRs in Sect. 6. Finally, we present our conclusions and a short summary in Sect. 7.

2. IRAS 23133+6050 and S 106

Both IRAS 23133+6050 and S 106 are well-studied regions of star-formation. In this section we review the pertinent literature in order to guide our later discussion. The properties drawn from the literature are summarized in Table 1. GLIMPSE survey images (Benjamin et al. 2003; Churchwell et al. 2009), obtained using the Spitzer/IRAC instrument (Fazio et al. 2004), of each source are shown in Fig. 1.

Fig. 1 Pointing for S 106 (left) and IRAS 23133+6050 (right). Overlaid on each GLIMPSE (Benjamin et al. 2003; Churchwell et al. 2009) IRAC [8.0] image is the corrected SPIRE beam aperture (43′′ and 37′′ respectively). We adopt the convention of coloring points relevant to S 106 and IRAS 23133+6070 red and blue respectively, beginning with the apertures in this figure.

2.1. IRAS 23133+6050

IRAS 23133+6050², appears to be a classical example of an isolated star-forming region within a molecular cloud (e.g., Kurtz et al. 1994, 1999 show it to be an isolated ultra-compact H ii region). From an astrometric perspective, there has been some debate regarding the distance of IRAS 23133+6050. Early estimates by Humphreys (1978) gave the distance to the Cassiopeia OB2 star-forming region, of which IRAS 23133+6050 is part, at 2.63 kpc. Peeters et al. (2002) found a distance of 5.5 kpc using kinematical methods and ISO-SWS spectroscopy, which agrees with an earlier calculation by Fich & Blitz (1984) who found d = 5.2 ± 0.9 kpc. However other authors, for example Lebrón et al. (2001), adopt a smaller value of 3.1 kpc, closer to initial estimates, based on spectrophotmetric observations of associated stars as derived by Brand & Blitz (1993). We adopt the higher value throughout this paper.

Kurtz et al. (1999) showed that the radio emission from IRAS 23133+6050 is consistent with the H ii region being powered by a late O star of spectral type O9.5. This star powers a cometary ultra compact H ii region with a diameter of 0.2 pc (Martín-Hernández et al. 2003) based on their calculated distance of 5.5 kpc. Using ISO observations of the MIR [S iii] lines Martín-Hernández et al. (2003) found the electron density of the H ii region to be about 2 × 10³ cm^-3. Comparatively few studies have discussed the PDR/molecular cloud envelope of the IRAS 23133+6050 region. Assuming pressure equilibrium, the expected density contrast between a typical H ii region and its PDR is expected to be around 30; which leads to an expected PDR density of the order of 6 × 10⁴ cm^-3. IRAS 23133+6050 was included in a sample of UC-H ii regions observed in the sub-mm by SCUBA, providing fluxes at 450 and 850 μm (Thompson et al. 2006).

2.2. S 106

S 106 (Sh 2-106; Sharpless 1959) is an extremely complex source, consisting of multiple components and IR sources. S 106 has been found to be at a distance of 1.2–1.8 kpc (Schneider et al. 2007), we adopt 1.5 kpc following Noel et al. (2005). S 106 is believed to be powered mainly by a single O-type star (S106-IR; Richer et al. 1993), which is surrounded by a very small optically thick disk (Hoare & Muxlow 1996). The double lobed structure of the H ii region is created by the shadows of this disk.

The PDR and molecular cloud components of S 106 have been the focus of prior studies (see Hodapp & Schneider 2008, for a thorough review). The overall structure of the S 106 system is very similar to bipolar nebulae formed by a combination of stellar outflows and a circumstellar disk. That the central star is undergoing significant mass loss has long been known (e.g., Drew et al. 1993; Gibb & Hoare 2007; Parkin et al. 2009) which is consistent with the idea that the large ionized lobes were sculpted in this way. However the origin of the central dark lane has been more controversial.

Early models of the S 106 system invoked a large optically thick torus of dense material as the explanation for the double lobed appearance of the system and central dark lane. Higher resolution radio/sub-mm observations (e.g., Richer et al. 1993) dispelled this notion though, and replaced it with the idea that the central regions of S 106 are dominated by a variety of clumps of denser material, which have in turn been suggested to be debris from the dissolution of a large dense torus (Little et al. 1995). In subsequent years it has been suggested that the dark lane is the projected shadow of an edge on accretion disk³ around the central young star (e.g., Noel et al. 2005), however recent SOFIA results seem to be inconsistent with this model (Simon et al. 2012).

Subsequently several authors have investigated the PDR environment of S 106 each using different datasets (e.g., van den Ancker et al. 2000; Schneider et al. 2003; Noel et al. 2005). The first of these, by van den Ancker et al. (2000), used large aperture ISO-SWS and LWS (de Graauw et al. 1996; Clegg et al. 1996) to find the global properties of the FIR emission from S 106. They conclude, on the basis of very strong [Si ii] emission, that the PDR associated with S 106 must have densities of around 10⁵ cm^-3 and UV radiation field strength of greater than 10⁵. van den Ancker et al. (2000) also conclude that this PDR must be 5−15′′ away from the central star. A subsequent study by Schneider et al. (2003) showed that the PDR of S 106 is much more extended than previously thought, extending up to several arcminutes from the central star. In addition, Schneider et al. (2003) used PDR diagnostics to show that a two phase PDR model was necessary to understand the line emission from S 106. They concluded that a very dense molecular component with an H₂ density greater than 3 × 10⁵ cm^-3 was necessary, along with a less dense interclump phase. This finding agreed with the previous study of CO emission by the same authors which found evidence of a very dense component (Schneider et al. 2002). Parallel efforts to those using higher spatial resolution mid-IR instruments have yielded further insight into the central structures. Noel et al. (2005) studied S 106 using near-IR spectro-imaging of a combination of hydrogen recombination lines and molecular hydrogen emission. From the perspective of this work, the Noel et al. (2005) study came to two key results: 1) the observations of H₂ lines were found to be consistent with shocks driven by outflows of the central star and; 2) the internal structure of the S 106 region based on their line measurements (presented in their Fig. 12). Shocks would be an expected part of any high-mass star-forming region, however previous studies found little evidence for their presence (e.g., van den Ancker et al. 2000).

3. Observations, data reduction, and analysis

3.1. Observations

A log of all Herschel observations used is presented in Table 2. The observations were performed as part of Herschel program OT1_epeeters_01. The objects were chosen to provide a sample to detect the long wavelength PAH features. This paper concerns the PDR properties of those objects which were not included in other PDR surveys.

The SPIRE observations each used the pointed spectroscopy instrument mode with no mapping to fill in the spatial gaps between the pixels. The PACS observations used the spectroscopy “Range Scan” mode, which provides a complete spectrum from 50–200 μm without any spectral gaps. These observations were combined to produce a complete spectrum for each object covering from 50−650 μm for the area on the sky covered by the central SPIRE spaxel beam after correcting for constant beam size. The pointings for each object are shown in Fig. 1 along with the corrected SPIRE beam size for each source (see Sect. 3.2). The spectral resolution varied as a function of wavelength for each instrument, with SPIRE providing resolutions between ~400 and ~1300 and PACS giving resolutions between 1000 and 5000. In practice this means that none of the detected lines were spectrally resolved.

3.2. Data reduction

Spire

The SPIRE data were reduced using the Herschel interactive processing environment (HIPE; version 11.0; Ott 2010). Each of our observations was taken in the “Bright Source” mode, which enables SPIRE to observe sources which would usually saturate the detectors. The reduction process used the “SPIRE Spectrometer Single Pointing User Reprocessing” HIPE script. This allows the user to calibrate bright source mode observations and use the newest version of the SPIRE calibration data which at the time of reduction was SPIRE_CAL_11_0.

Accurate calibration of our sources was achieved by using the HIPE semi-extended correction tool which implements the methodology described by Wu et al. (2013). The semi-extended correction tool calculates the degree to which the observed source differs from a point source using the break between the orders of the SPIRE detectors where the beam size changes. In order to use this tool it was necessary to reduce the observations using the “point source” calibration rather than the “extended” calibration. The semi-extended correction tool results in a spectrum which has been corrected to a constant beam size, 43′′ for S 106 and 37′′ for IRAS 23133+6050. This eliminates the flux jump between the two orders of the SPIRE instrument and allows easier cross calibration with the PACS data. The semi-extended correction tool only corrects the fluxes for the central pixel as the beams of the long and short wavelength detectors are centered with respect to each other. For each of the surrounding, off-axis, SPIRE pixels for both objects also contain semi-extended structure (e.g., Fig. 1) and as such they cannot be reliably calibrated at this time due to this software limitation.

PACS

The PACS data were reduced using HIPE version 11.0 and the PACS calibration set PACS_CAL_56_0. The observations of both S 106 and IRAS 23133+6050 were performed using the unchopped spectroscopy mode with corresponding background observations (see Table 2) because both sources are extended beyond the maximum PACS chop throw. Exposure times and ObsIDs are given in Table 2.

The PACS data are affected by several systematics, mainly order leakage and spectral ghosts⁴. The regions of order leakage were trimmed as necessary. This principally affected the longer wavelength spectra in band R1 (105–210 μm), where it proved impossible to have a well calibrated spectrum beyond around 190 μm. Regions affected by these problems were trimmed from the final data, resulting in gaps in the spectral coverage.

It was found that the continuum was offset slightly in each PACS spectral order, likely a reflection of the varying beam size, and as such the orders were scaled so that they agreed with the longest wavelength order. It was found that the offsets between the orders were additive, in that the line fluxes for lines in overlapping regions matched before scaling. For S 106, the blue orders agreed with the final red order (B2A, B2B, R1B. Hence the only scaling necessary was for R1A, which had to be increased by the addition of 850 Jy to agree with the continua of the other orders. For the IRAS 23133+6050 measurements, the first PACS spectral order was found to have a higher background than expected resulting in a strong offset of 950 Jy with respect to the continuum of the other bands. While these values sound large, they represent factors of only around 20%.

The final PACS spectra were created by coadding the pixels to match size of the beam produced by the SPIRE semi extended source correction. This was accomplished by calculating the relative contributions of each pixel to this beam and summing all of the pixels in the correct proportions. The final fluxes were achieved by scaling this coadded PACS spectrum so that the continuum slopes and intensities matched. This scaling factor was found to be multiplicative in order for the fluxes of the ¹²CO lines to be consistent across the PACS/SPIRE gap. For S 106 this factor was relatively small (1.5), while for IRAS 23133+6050 it was much greater (4).

3.3. Analysis

Fig. 2 Full PACS + SPIRE spectra for IRAS 23133+6050 (top, blue) and S 106 (bottom, red, multiplied by 0.01 for clarity). In addition we include the IRAS photometric points for each source (triangles), and for IRAS 23133+6050 also the WISE photometry (Wright et al. 2010) as squares and the SCUBA observations of Thompson et al. (2006) as diamonds. A modified blackbody fit for each spectrum is shown as a dashed line.

3.3.1. Line fluxes

In each of the coadded final spectra, the spectral lines were analyzed in the following ways. The line fluxes were extracted via fitting Gaussians and local continua to the data. It was found that better fits resulted from fitting a linear local continuum to regions of the spectrum on either side of the line, and then fixing the linear parameters fitted for these regions when fitting the full spectrum of a Gaussian line profile plus linear continuum. The signal-to-noise estimates for each line measurement were calculated by using the root mean square (rms) noise in the continuum spectral windows after the subtraction of the linear fit. This rms value can then be combined with the integrated line flux, width and the instrumental resolution to yield the final signal-to-noise ratio for the measured line. There is only one significant blended pair of lines present in the data, the 371.65 μm CO 7−6 and 370.42 μm [C i] 2−1 lines. Initially we performed a fit using two Gaussians to the complex. We have checked this fit by calculating the flux of the 7−6 CO line using the neighboring CO lines and subtracting this flux from the total flux of the blended complex. The residual of this process yielded the same [C i] flux as the initial fitting procedure. The measured line fluxes are given in Table 3.

For S 106 we can compare the measured PDR cooling line fluxes here to those found by other observations in the past. The most direct comparison can be drawn to those presented by van den Ancker et al. (2000) and Schneider et al. (2003) which are based on observations by ISO-LWS, KAO-FIFI and KOSMA-SMART (Clegg et al. 1996; Poglitsch et al. 1991; Graf et al. 2002). A direct comparison between the measured surface brightnesses is given in Table 4. Some care should be taken in comparing the presented quantities as the beam sizes differ by large factors from the Herschel measurements presented here. For the dominant cooling line, [O i] 63 μm, it is easy to see why the ISO-LWS surface brightness measurements would be much smaller as they were measured using a larger beam by a factor of around five, and hence the overall surface brightness is lower as the flux of S 106 is peaked strongly in the central regions near the exciting source IRS4. van den Ancker et al. (2000) measured surface brightnesses for a separate set of ISO-LWS observations which were also centered on S 106 and found that the [O i] 63 μm line had a surface brightness between the values found by Schneider et al. (2003) and our measurements. Evidently different pointings around the S 106 region can result in very different results, even considering the fact that the comparatively large LWS beam encompasses the majority of strong PDR cooling line emission from the S 106 complex (e.g., see the Sofia/GREAT map of [C ii] 158 μm presented by Simon et al. 2012).

Overall, our observations of the [O i] line agree within ~30% with the LWS measurements if we assume that the S 106 flux is concentrated at the center of the beam and compensate for the different beam sizes. The same process yields similar results for the other LWS cooling line measurements – for example for the [C ii] 158 μm line our results agree with the LWS observations almost exactly for the van den Ancker et al. (2000) pointing and to within 30% for the Schneider et al. (2003) pointing, assuming that the beam size quoted by Schneider et al. (2003) is correct (see footnote to Table 4). Overall then, we find agreement within the quoted absolute flux uncertainties of 30% for both ISO-LWS and PACS (LWS: Swinyard et al. 1996; PACS: PACS Spectrometer Calibration Document⁵). At first glance, the KAO surface brightness measurement of [O i] (from Schneider et al. 2003) that we include in Table 4 seems incompatible with our measurement, as it was for a beam smaller than the SPIRE beam by a factor of three and found a lower surface brightness. However, Schneider et al. (2003) discuss this problem as the discrepancy between KAO/Herschel is almost the same as between the KAO and LWS measurements. Schneider et al. (2003) conclude that the KAO and LWS measurements can be reconciled by taking into account the spatially varying emission intensities – specifically that if one assumed a Gaussian profile for the [O i] emission the LWS beam size, the LWS flux could be recovered from the KAO measurement.

3.3.2. Total FIR flux

The total flux of the FIR dust continuum was determined by fitting a modified blackbody to each spectrum. In order to better constrain the peak of this spectrum, we searched for additional observations of our sources in the 20−60 μm regime such that we might acquire extra data points via aperture photometry using the SPIRE beam. However, realistic measurements of the fluxes for these sources were not available as they are very bright and saturated in both MIPS and WISE data. Figure 2 includes IRAS photometric points, however these are integrated totals for the entire source and not scaled to the Herschel beam. In addition we investigated using MSX photometry, however the MSX fluxes agree well with the IRAS points. We therefore investigated using spectral data from ISO-SWS and LWS to help constrain the shorter wavelength regime. The ISO spectrometers used very different beam sizes though, with the LWS aperture being of ~120′′ diameter while the SWS aperture was 20″ × 33″ at its longest wavelengths. This corresponds to very different areas of each object dominating the beam. For SWS the bright central parts of the nebulae fill the beam, while for LWS the colder outer areas fill the beam. For S 106, the LWS spectra agree well with the final PACS spectrum, albeit with an offset of around 50%. This is likely because there is extended bright structure in the LWS beam and as such it gives a similar spectral slope. For IRAS 23133+6050, the picture is less clear, as the PACS spectrum has a much steeper gradient than the LWS spectrum. IRAS 23133+6050 is much more compact than S 106, with all of the hot material being confined to the center. The shallower gradient of the LWS continuum then may reflect the extra cold material in the beam.

For each object the peak of the spectral energy distribution falls just before the PACS wavelength range, therefore in order to measure the total infrared flux we fit a modified blackbody to the spectrum and integrated all of the flux of the fitted blackbody between 10 and 1000 μm. The aperture mismatches with the existing ISO-SWS and -LWS data and our observations proved impossible to reconcile and as such we have no strong constraints on the blue end of the spectrum. As a test of this uncertainty, we performed an additional fit in which the FIR flux was fit with two blackbodies. The single blackbody fit discussed previously was held constant, while an extra component was added to fit the IRAS points shortward of the Herschel data. The FIR flux found in this way was around 25% greater than that measured using only one blackbody. This figure is likely to be an overestimate of the true FIR flux, as the flux of the 12 μm IRAS filter is dominated by PAH emission. For repeatability and simplicity, we adopted the single blackbody measurement with an uncertainty of 25% on the derived total FIR fluxes. The total FIR fluxes for each object derived in this manner are listed in Table 3.

4. PDR diagnostics

It is clear from inspection of Table 3 that the [O i] 63 μm line dominates cooling for both PDRs, dominating both of the other important cooling lines ([O i] 146 μm and [C ii] 157 μm) by factors of between five and eight. This is typical of dense PDRs. In addition CO rotational lines extend to very high-J levels also indicating warm dense PDRs. On the other hand, the [C ii] and [C i] emission can also arise in low density gas illuminated by weak radiation fields. There is likely a mix of densities, temperatures, and filling factors that combine to produce the observed emission. We next use PDR line diagnostics to investigate the various phases.

4.1. FIR continuum and stellar spectral types

The FIR continuum can be used to directly probe the UV radiation field impinging on the PDR system. The total energy absorbed by the dust in the system must be the same as that re-emitted. Previous authors (e.g., Meixner et al. 1992) have used this relationship to calculate a value for G₀. This calculation is summarized in Appendix B of Meixner et al. (1992) where it is shown that the derived G₀ values weakly depend on the emitting structure. For a spherically symmetric source (e.g., IRAS 23133+6050), it can be assumed that the absorbing and radiating surface areas are equal, and hence the relationship can be summarized as: G₀ = ΣFIR/1 Habing field. For more complex sources the relationship between the radiating and absorbing surface areas is different, for example, Meixner et al. (1992) show that for a disk system that the derived G₀ will be different by a factor of two. S 106 would seem to be a mixture of these two paradigms, with a small disk around the central star being further encased within a molecular cloud/PDR, therefore there will exist a range of values of G₀ spanning these geometries. For each object the resulting values of G₀ assuming a spherically symmetric geometry are very similar at around 1.3 × 10⁵. For the disk region of S 106 though, G₀ could be up to 2.6 × 10⁵.

Fig. 3 PDR diagnostic diagrams for observed fine-structure lines – based on the PDR models of Wolfire et al. (2010) and Hollenbach et al. (2012). Top: [C ii] 158 μm/[O i] 63 μm versus ([O i] 63 μm + [C ii] 158 μm)/ΣFIR. Middle: [C ii] 158/[O i] 145 versus ([O i] 145 μm + [C ii] 158 μm)/ΣFIR. Bottom: [C i] 370/ [C i] 609 versus ([O i] 145 μm + [C ii] 158 μm)/ΣFIR with the density grid extended down to log n = 2. In each figure the blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively.

Fig. 4 PDR diagnostic diagrams for ratios of 12CO(J=15-14), 12CO(J = 7−6) and 12CO(J = 4−3) versus ([O i] 63 μm + [C ii] 158 μm)/ ΣFIR. In each figure the blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively. The data points indicate a much denser environment with a stronger UV field than that seen for the fine-structure line diagnostics.

An additional measure of G₀ can be calculated using the distance from the ionizing source to the PDR and the stellar spectral type of the exciting star. Tielens (2010) gives the relationship between G₀ and the stellar luminosity as ${\mathit{G}}_{\mathrm{0}}\mathrm{=}\mathrm{625}\mathrm{\times }\frac{{\mathit{L}}_{\mathrm{\ast }}\mathit{\chi }}{\mathrm{4}\mathit{\pi }{\mathit{d}}^{\mathrm{2}}}$(1)where L_∗ is the stellar luminosity (in erg s^-1), χ is the fraction of photon energy emerging from the H ii region between 6 and 13.6 eV, and d is the distance between the star and the ionization front in cm. Both IRAS 23133+6050 and S 106 are powered by late O8-9 type stars, which have far UV luminosities around 3.25 × 10⁴ times that of the sun⁶. For IRAS 23133+6050, the ionized region has a radius of 4.5′′, which at a distance of 5.5 kpc translates to 0.1 pc. This implies a G₀ value of around 5 × 10⁴. Similarly for S 106, the radius of the ionization front is around 0.1 pc (based on a distance of 1.5 kpc and a radius of 12.5′′; Noel et al. 2005), which implies a G₀ of about 8 × 10⁴ when using the same assumptions. Noel et al. (2005) found that the H ii region of S 106 was approximately cylindrical, and as such the G₀ we calculate using their radius is only appropriate for the parts of the ionization front closest to the exciting star and therefore is actually an upper limit. In good agreement with our analysis, the Habing field for S 106 determined from the Herschel PACS 70 and 160 micron fluxes is 3 × 10⁴G₀ (value at a distance of 0.1 pc from S106 IR in a 12′′ beam, Schneider et al., in prep).

4.2. Diagnostic diagrams

The measured line fluxes given in Table 3 can be combined in various ways to produce ratios which are sensitive to the physical conditions in the PDR. The variation of pairs of these ratios with respect to the fundamental PDR parameters G₀ and n⁷are shown in Figs. 3 and 4. To construct these diagrams we used the results of the PDR models of Wolfire et al. (2010) and Hollenbach et al. (2012). There are a wide variety of different PDR models available in the literature, and they have been extensively compared and validated against each other, for example see Röllig et al. (2007), in which the model we have chosen to employ here is referred to as the HTBKW model. These models have all been shown to reasonably agree for most PDR diagnostics and all reproduce the same trends when run on idential data (Röllig et al. 2007).

In each case we consider a range of G₀ values between 10³ and 10⁶ Habings and a range of densities between 10³ and 10⁷ cm^-3. In Table 5 we give an overview of the physical conditions indicated by the various diagnostics discussed in addition to those implied by the total FIR flux and the stellar spectral types of both sources.

4.2.1. [C II] and [O I]

The primary diagnostics usually employed for PDRs involve ratios between FIR fine-structure cooling lines. Usually the ratios involving [C ii] 158 μm and the two [O i] lines are employed reflecting that these lines are the best understood of those included in PDR models. Typically such ratios do not correspond to unique PDR parameters, as their dependence of the line ratios on G₀ and n is complex. In order to produce a single value, other diagnostics must be employed such that intersections at specific G₀ and n values are found. The simplest combination of diagnostics which yields such physical conditions occurs when considering [C ii] 158 μm/[O i] 63 μm and ([O i] 63 μm + [C ii] 158 μm)/ΣFIR diagnostic diagram (Fig. 3, top).

This diagram implies that both PDRs have similar physical conditions, with log G₀ around 4 and log  n around 3.75, well within typical PDR conditions. These values are compatible with those which were discussed in Sect. 2. For IRAS 23133+6050 the density of ~10⁴ agrees nicely with that expected given the density of the H ii region, while for S 106 this density likely represents the lower density material which makes up the majority of the PDR volume (e.g., Schneider et al. 2003).

The ratio of the [O i] lines measured for both sources is of the order of [O i] 145 μm/[O i] 63 μm ≃ 0.15, which agrees with the measurements of Schneider et al. (2003) who found 0.17 using ISO-LWS observations. This clearly indicates that one or both of the [O i] lines is optically thick (e.g., Tielens & Hollenbach 1985b; Fig. 2): if both lines were optically thin their ratio would be less than 0.1. More recently it has been shown that in many cases the [O i] line ratio indicates optically thick conditions, though other explanations for such ratios exist. It was concluded by Liseau et al. (2006) that self absorption of the [O i] 63 μm line along the line of sight with modest optical depths of around 1–2 could be responsible for high values of the [O i] 145 μm/[O i] 63 μm ratio⁸. Those results refer mainly to an ISO-LWS sample of dark clouds, for which the density required for both [O i] lines to be optically thick was considered very unlikely. For IRAS 23133+6050 and S 106 though, such densities may be plausible – however [O i] 63 μm self absorption would likely also play a role in shaping the flux ratio, rendering it difficult to disentangle the two.

4.2.2. [C I]

In addition, we also investigate the [C i] fine-structure lines by considering the ratio of the fluxes of the 370 μm line over the 609 μm line (bottom panel of Fig. 3). This ratio indicates very different physical properties to those mentioned previously, with the ratio of the [C i] lines indicating low density, weak radiation field PDR conditions. It should be noted though that it has been shown that the [C i] intensities given by models are generally uncertain (maximal variation of around a factor of ten; Röllig et al. 2007) and therefore the results for this line ratio are highly model dependent.

4.2.3. ¹²CO

Fig. 5 [C i] 609 μm/12CO (J = 4−3) versus ([O i] 63 μm + [C ii] 158 μm)/ ΣFIR PDR diagnostic diagram. The blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively.

In addition to the fine-structure cooling line diagnostics we can also employ the properties of the molecular lines in the spectra in the same way. Kaufman et al. (1999) provide the various diagrams for the different ¹²CO transitions, here we investigate ratios of three: 4−3, 7−6 and 15−14, but using the updated PDR models of Wolfire et al. (2010) and Hollenbach et al. (2012).

The ratio of the two mid-J lines: ¹²CO(J = 7−6)/¹²CO(J = 4−3) is shown in Fig. 4 (top), and indicates a much stronger UV field (log G₀ = 4.9) and denser environment (log  n = 4.8) for both sources than the conditions indicated by the fine-structure lines. For the highest-J¹²CO line, J = 15−14, even denser and more strongly illuminated conditions are found, with densities an order of magnitude higher and radiation fields around double that found from the lower excitation lines. Diagnostic diagrams for the ¹²CO (J = 15−14) line, as compared to the lower excitation CO lines, are shown in the lower two panels of Fig. 4.

In addition to the pure ¹²CO ratios, we also show the ratio of the [C i] fine-structure line with the ¹²CO(J = 4−3) line in Fig. 5. This diagnostic clearly shows a similar pattern to that of the other ¹²CO diagnostic ratios, i.e., high log G₀ (~5) and log  n (~5).

4.3. CO ladder fluxes

Fig. 6 12CO (squares) and 13CO (triangles) ladder fluxes for IRAS 23133+6050 and S 106 in terms of the excitation energy of that level.

In the previous section we exploited the simple relationships between different CO transitions resulting from PDR models to find the physical conditions. It is also possible to derive such conditions directly from the fluxes of the whole set of molecular lines. In Fig. 6 we show the ¹²CO and ¹³CO line fluxes against the energy of the upper level E_J. The detections of the ¹²CO lines for both IRAS 23133+6050 and S 106 are measured across the entire SPIRE/PACS range, that is, from J_u = 4 to 23, while for ¹³CO the range is smaller, from J_u = 5 to 12, owing to their weaker fluxes. The ¹²CO fluxes shown in Fig. 6a are remarkably similar for both objects, peaking at around E_j = 550 K, or J_u = 13 or 14. The only difference between IRAS 23133+6050 and S 106 comes from the higher-J lines, where for S 106 the four highest-J lines show additional excitation, whilst those for IRAS 23133+6050 continue to decrease. The possible implications of this will be discussed in the following section. The ¹³CO line fluxes seem to peak at around E_j = ~ 300 K, or J_u = 10, however the data are noisy and it is not clear that this apparent peak is real.

The peak of the ¹²CO fluxes is itself an indicator of the physical conditions of both sources. Using the density of the high density component predicted by the simple diagnostics (n ~ 10^5.6 cm^-3) and the radiative transfer code RADEX (van der Tak et al. 2007) we calculated the fluxes of ¹²CO lines for a variety of temperatures. The required temperature to fit the peak is then around 400 K. In Fig. 7 we show this model, combined with a model representing the low density component at the same temperature, which have been scaled to match the overall flux of the whole CO ladder. However, the combination of density and temperature which fits the peak is not unique, and higher or lower density models can fit the peak almost equally well with suitably reduced or increased temperatures. Our choice of n = 10^5.6 cm^-3 is driven by the simple diagnostics discussed in previous sections. In contrast, the highest and lowest “plausible” solutions considered have densities and temperatures of 10⁷ or 5 × 10⁴ cm^-3 and 200 or 1500 K respectively. Thus RADEX models lead us to conclude that the most plausible explanation for the mid-J (J_U ≃ 6−20) ¹²CO ladder is a (very) high density PDR with a temperature of around 200−400 K. We attempted to fit the high-J S 106 lines using RADEX, however it was not possible to find a good fit to those points without severely compromising the quality of the fit to the other points.

Fig. 7 Combination of RADEX models scaled to fit the S 106 12CO ladder. The low density component with log n = 4, T = 300 K is represented by the green line, while the dominant contribution arises from the dense component with log n = 5.6 and a temperature of 400 K which is represented by the purple line. The aggregate model (black line) matches the observed 12CO ladder from JU = 4 to 20.

4.3.1. CO rotation diagrams

To further illustrate the differences between our sources and those found in the literature we create rotation diagrams from the ¹²CO and ¹³CO fluxes (Fig. 8). Following Justtanont et al. (2000) we construct rotation diagrams showing the relationship between the line fluxes F_ji, the temperature T_x and the column density N₀ as ${\mathit{E}}_{\mathit{j}}\mathit{/}\mathit{k}{\mathit{T}}_{\mathit{x}}\mathrm{=}\ln \left(\frac{{\mathit{N}}_{\mathrm{0}}}{\mathrm{\Phi }\mathrm{\left(}{\mathit{T}}_{\mathit{x}}\mathrm{\right)}{\mathit{d}}^{\mathrm{2}}}\right)\mathrm{-}\ln \left(\frac{\mathrm{4}\mathit{\pi }{\mathit{F}}_{\mathit{ji}}}{{\mathit{A}}_{\mathit{ji}}\mathit{h}{\mathit{\nu }}_{\mathit{ji}}{\mathit{g}}_{\mathit{j}}}\right)\mathit{,}$(2)where g_j refers to the statistical weight of the upper level j, ν_ji is the frequency of the transition, Φ(T_x) is the partition function, A_ji is the oscillator strength for the particular transition and d is the distance to the source. Justtanont et al. (2000) refer to the second term as Y_j and we follow this convention. The temperature and column densities can then be found in the limit of optically thin emission in local thermodynamic equilibrium (LTE) by fitting a straight line to the E_j, Y_j coordinates. There are many potential problems with this method, mainly involving the aforementioned assumptions (see, for example, Goldsmith & Langer 1999). In light of that, it is not surprising that both IRAS 23133+6050 and S 106 do not follow exactly straight lines. Nevertheless, since it is a common procedure, we will analyze the observations using this technique and compare to other methods and models.

Fig. 8 12CO (diamonds) and 13CO (triangles) rotation diagrams for IRAS 23133+6050 and S 106. Two ultra-compact H ii regions from the Karska et al. (2014) sample which display the warm CO component are also shown as orange, purple lines. Simple fits to the lower-J CO lines are shown to demonstrate the deviation of the higher-J CO levels from this trend for S 106 and the Karska sources. A fit to the highest-J CO lines for S 106 is also included. The Karska sources have been offset in y for clarity. The first 12CO and 13CO points for both of our sources are 4−3 and 5−4 respectively.

The curvature exhibited at low E_j is expected and usually interpreted as requiring a second, warm, temperature component (cf. Fig. 7). These levels trace the warmish material on the surface of the (putative) molecular cloud, the colder molecular cloud material would only manifest itself in unobserved lower-J transitions (J< 3). However recent work has shown that this curvature can be explained in the optically thin limit by deviations from LTE at densities below 10⁵ cm^-3 and constant temperatures (Neufeld 2012). This effect arises because, in the low density limit, the population is dominated by the low lying transitions, with collisional excitation and spontaneous decay balanced for the higher level transitions.

From these best fits to the ¹²CO lines of S 106 and IRAS 23133+6050 we find excitation temperatures of around 190 K and around 10⁴⁹ molecules along the line of sight (summarized in Table 6). We have also included a linear fit to the high excitation component of the S 106 ¹²CO spectrum – this fit yields a much higher excitation temperature of around 330 K with a lower density. The derived excitation temperatures for ¹²CO are systematically lower than those found by the RADEX models discussed in the previous section. This discrepancy is likely due to the assumptions about optical depth which underpin the use of rotation diagrams – it seems unlikely that many of the low-J CO lines are optically thin, for example⁹. The RADEX models make no such assumptions and are therefore a somewhat more trustworthy indicator of the temperature of the various components, on the other hand the RADEX models assume either an arbitrary gas temperature distribution or constant temperature, and constant molecular abundances. However, we could not find a plausible solution for the very high excitation S 106 points with RADEX models, whereas the rotation diagrams tell us that the temperature of that component is higher than that of the medium density component by at least 100 K. However, even this argument relies on the assumption that our sources are in LTE, which is not altogether clear from inspection of Fig. 8. As such the parameters of the rotation diagram fits should be regarded with some caution, we will discuss how well these results compare to the results of the more detailed models in Sect. 6.3.

The linear component seen in the four highest-J¹²CO lines of S 106 betrays the possible existence of an additional physical component. There appear to be two possible explanations for this phenomenon: an even higher density PDR; or shocks in the S 106 environment. Such additional excitation in the upper CO levels has been observed commonly in a related class of objects: protostars. We have therefore included in the rotation diagram several prominent examples from a recent survey of high mass star-forming regions using PACS by Karska et al. (2014). The two objects (G34.26+0.15 and NGC 7538) were chosen for comparison on the grounds that they are ultra-compact H ii regions and therefore closest in structure to S 106 and IRAS 23133+6050. In fact, NGC 7538 is actually a bipolar H ii region similar to S 106. It is clear that these sources display an extra component, similar to that seen for S 106. In order to visualize the offset caused at high-J, we have included best fit lines to the mid-J linear portion of the data for each object in Fig. 8. In high mass protostars such as those studied by Karska et al. (2014) the additional excitation is thought to arise from a combination of shocks and PDR emission (for low mass protostars such excitation is usually thought to be dominated by shocks, for example van Kempen et al. 2010). These possibilities will be discussed in Sect. 6.3. In agreement with van den Ancker et al. (2000), we conclude that the majority of the CO emission observed towards S 106 is due to a dense PDR, particularly for the low to mid-J lines, however it is clear that the high-J CO lines may be indicative of a contribution by a shock.

Fig. 9 PDR model fits to 12CO ladders. a) Two component PDR model for S 106; b) three component PDR model for S 106; c) two component PDR fit for IRAS 23133+6050. In each panel the dashed line represents the low pressure component and the dot-dashed lines the high pressure component. In the centre panel the final, very high density PDR component for S 106 is represented by a triple-dot-dashed line.

4.4. Summary of diagnostics

The simple diagnostics reveal two very similar sources, with subtle differences possibly driven by the different physical structures of the two sources. Both show high G₀, as revealed by the total IR fluxes, but with lower radiation fields indicated by the various fine-structure and molecular line ratios. Thus a two component model is required, where the PDR is comprised of two distinct regions: i) a hot (T ~ 400 K), dense region characterized by high G₀ (>10⁵) and high density (>10⁵ cm^-3) traced by the various ¹²CO lines; and ii) a less dense region region dominating the emission of the atomic cooling lines characterized by moderate UV fields (G₀~ 10⁴) and densities (10⁴ cm^-3) with a temperature of around 300 K. The presence of such a distinction is not a surprise however, as models of PDRs have included density contrasts of this order since the early nineties (e.g., Burton et al. 1990b) and it has been observed in a number of objects, from Orion to M17 (e.g., Stutzki & Guesten 1990; Tauber et al. 1994). In addition the very high excitation CO lines for S 106 show an even hotter component, which may arise from a shock or a much hotter PDR than those considered thus far. A detailed comparison between our sources and typical galactic PDRs will be presented in Sect. 6.1.

5. Numerical PDR models

The densities and UV fields found in the previous section can be regarded as line of sight averages as they do not take into account the geometry of the system or the absolute fluxes of the lines. In order to extract extra information from the observations, we used the simple analysis as the input for numerical PDR codes. This process should help illuminate the structure along the line of sight producing the emission observed.

We use the PDR model described by Wolfire et al. (2010), but modified to use either constant density or constant thermal pressure. The model builds on that of Kaufman et al. (2006), which in turn combined the initial PDR models of Tielens & Hollenbach (1985b) and Wolfire et al. (1990). In addition the model includes both the molecular hydrogen processes incorporated into the Meudon PDR code (e.g., Le Petit et al. 2006) and molecular freeze-out and oxygen chemistry as described by Hollenbach et al. (2012). The model assumes a plane parallel slab geometry illuminated from one direction with a specified 1D UV field and calculates the thermal balance and abundances of the dominant atomic and molecular species as well as their line emission. These models have been used to interpret a variety of observations of molecular and atomic species in a wide range of physical conditions including OH⁺, H₂O⁺, and H₃O⁺ in diffuse gas (Hollenbach et al. 2012), H₂O and HF in diffuse gas (Sonnentrucker et al. 2015), [C i] and CO in low UV field molecular clouds (Burton et al. 2014; Lee et al. 2014), H₂ in dense PDRs (Sheffer et al. 2011) and [C ii] and [O i] in dense PDRs and extragalactic nuclei (Kaufman et al. 2006; Croxall et al. 2012).

In order to get good fits to the observations of both PDRs constant thermal pressure models were used, with one exception which will be discussed shortly. In each case the final model shown in Fig. 9 is a combination of high pressure and low pressure models with the parameters guided by the properties of each regime from the simple diagnostics. Because the simple diagnostics report only average properties within the beam, it is necessary to scale the high and low density models for each source to match the overall fluxes. Hence we scaled the models using the CO line fluxes and found good agreement with the exception of the highest-J lines for S 106 which display additional excitation as discussed earlier. In order to attempt to fit these lines, a third, extremely high density, constant density PDR model component was introduced¹⁰. The degree to which each component needs to be scaled can be interpreted as the filling factor for the high density model, and a rough estimate of the number of PDRs on the line of sight for the low density model.

Since our models are for constant thermal pressure, both the temperature and density vary with depth within the cloud. We characterize the physical conditions where the line emission arises at the depth where one half of the intensity is produced. For IRAS 23133+6050, the peak CO line emission (CO 13−12) comes from the high pressure component at T = 217 K, n = 1.3 × 10⁶ cm^-3. The [O i] 63 μm line comes from a region closer to the cloud surface where temperatures are higher and densities are lower. We find T = 2030 K, and n = 7.6 × 10⁴ cm^-3. The low pressure component provides about half of the observed CO 5−4 emission at T = 28 K, and n = 9.0 × 10³ cm^-3. It also provides half of the [C ii] emission at T = 530 K, n = 260 cm^-3. At depths where these lines arise, including the high-J CO lines, grain photoelectric heating dominates and we find no need for additional and exotic heating processes. For S 106 we consider only the (a) model which neglects the extra, high-J component. For CO 13−12 we find T = 214 K, n = 1.4 × 10⁶ cm^-3, and [O i] 63 μm we find T = 2070 K, n = 7.9 × 10⁴ cm^-3. For the CO 6−5 emission and [C ii] we find T = 40 K, n = 2.7 × 10⁴ cm^-3 and T = 370 K, n = 1.2 × 10³ cm^-3 respectively. Again, we find that photoelectric heating dominates in regions where these are formed.

The high density component dominates the overall CO ladder fluxes for upper levels greater than J_u ≃ 9, and therefore the filling factor of the high density component is well constrained by our fits at ~0.62 and ~0.57 for S 106 and IRAS 23133+6050 respectively. For the low density components factors of 4.6 and 5.1 were found for S 106 and IRAS 23133+6050 respectively, strongly suggesting that additional low density PDRs are present along the line of sight. Finally, the filling factor of the extra very high density PDR component included to fit the highest-J S 106 lines was 0.006. The parameters of the various components which make up the composite models are given in Table 7.

In addition to the CO line fluxes, the models also report the fluxes of each of the fine-structure lines in the SPIRE/PACS range. In Table 8 we show the ratio of the measured line fluxes and those predicted by the numerical models. For most lines the model does a good job of replicating the observed line flux within a factor of 2−3. As the [O i] 63 μm and [C ii] 157 μm lines are commonly employed as tracers of shocks, we investigated whether the underprediction of the fluxes of these lines could be due to a contribution of a possible shock. We have estimated the possible contribution to these lines from shocks in S 106 by measuring the fraction of the total ¹²CO flux present in the high excitation component which possibly represents a shock compared to the dominant, cooler, components which are well fit by PDR models. This ratio is around 5%, showing that even in S 106 where there is evidence for a shock in the CO lines, it is over an order of magnitude smaller in flux than the PDR emission in the CO lines. Assuming that the effect on the ionic lines would be of similar magnitude, shocks cannot explain the discrepancy between our models and observations. In fact, the differences in line flux estimates between different PDR models are of the same order of magnitude as the discrepancies we observe between our numerical model fluxes and the observed line fluxes: therefore we conclude that the extra emission in these lines is unlikely to be due to shocks. In particular the IRAS 23133+6050 [C ii] 157 μm line provided the most discrepant line flux ratio of around four, which is somewhat strange given both sources were fit with similar parameters and only the [C ii] in IRAS 23133+6050 shows a large deviation from the model. We suspect that additional low density [C ii] line emission lies along the line of sight to the source, or alternatively this could reflect complexity in the geometry of the source which could perhaps be resolved by higher spectral resolution measurements which would reveal the velocity structure.

6. Discussion

6.1. Comparison with Galactic PDRs

While the Orion region remains the prototypical example of a PDR, many other examples have been observed and interpreted using PDR models. In Table 9 we show the properties of the regions studied here alongside a sample of directly comparable PDRs associated with H ii regions, as well as the properties of some PDRs generated by cooler stars¹¹. If we consider initially the properties of the interclump phase, our sources lie somewhere between these two classes of sources, with lower overall radiation fields than M17, W49A or Orion – as one might expect given their relatively cooler exciting stars and the fact that M17 and W49A are illuminated by clusters of O stars. In fact, the central clusters of both M17 and W49A both possess several examples of stars with the spectral type listed in Table 9. The G₀ values we find for our sources are actually much closer to those found for the reflection nebulae, particularly that of NGC 2023. In terms of density, the interclump medium for our PDRs is consistent with the other examples shown in Table 9, with the exception of the Orion bar which presents a density around an order of magnitude higher. For the dense clump component the density derived for S 106 is very similar to that found before for S 106 by Schneider et al. (2003), and also close to the value quoted by Meixner et al. (1992) for M17. These values are typical of those quoted for the clumps in a dense PDR, however the value of ~10⁷ cm^-3 for the dense component of the Orion bar PDR by Tielens et al. (1993) is actually very similar to the very high density PDR component we invoked in our numerical models to attempt to fit the highest-J CO lines in S 106. This would suggest that the high-J component might also be present for the Orion bar region. Hints of this phenomenon may be visible in the SPIRE observations of the Orion bar presented by Habart et al. (2010), where the presented RADEX fits to the ¹²CO lines begin to deviate in the mid-J regime (J ~ 12), however the region was not observed using the PACS spectroscopy range scan mode employed for our objects, so observations with SOFIA or other FIR observatories would be required to confirm this hypothesis.

Several groups have shown that there is a relationship between the incident radiation field and the PDR density for a variety of objects. Studies have examined this relationship for a sample of reflection nebula PDRs (Young Owl et al. 2002), while a more recent study found a slightly different relationship for spatially resolved regions of the Carina nebula (Oberst et al. 2011). If we consider the upper bound values for the inter-clump density in Table 9, our two sources fall almost perfectly on both of these trends, as is shown in Fig. 10. The dashed line represents the expected conditions for systems in a state of equilibrium between the pressure and temperature of the H ii region and PDR respectively under the assumption that the H ii region is a Stromgren sphere and has an electron temperature of 8000 K (Young Owl et al. 2002; referred to as “Young-Owl et al. theoretical” in Fig. 10).

Fig. 10 Correlations between G0 and density n0 found by Young Owl et al. (2002) and Oberst et al. (2011). Reflection nebula points from Young Owl et al. (2002), Galactic PDRs as stated in Table 9. IRAS 23133+6050 and S 106 fit into the trend between the larger H ii regions with hotter exciting stars and the RNe PDRs with B type exciting stars. Dashed and solid black lines: theoretical prediction for the relationship between radiation field and density for a pressure bound H ii region and the best fit to the Young Owl et al. (2002) data. Green dot-dashed line: best fit to the Carina nebula data of Oberst et al. (2011). For a selection of sources we have also included the points representing the dense PDR component (denoted “dc”) as opposed to the default values representing the interclump medium (denoted “ic”).

Conversely, Oberst et al. (2011) suggest that their best fit line represents systems which are density bounded as opposed to the ionization bounded RNe, and hence they found a different relationship. We have also included the H ii region PDR sources from Table 9 while the RNe sources in Table 9 form part of the Young Owl et al. (2002) sample and are represented by the black diamonds in Fig. 10. Our sources are consistent with the H ii region PDRs with late O type exciting stars and the RNe formed by early B stars. Amongst the H ii region PDRs with hotter exciting stars there seems to be a greater spread in properties away from the trends, with W49A and W3 IRS 5 in particular having either higher G₀ or lower n_H by at least an order of magnitude than predicted by either relationship. We have also included the properties of the dense phase for several sources, denoted as “dc” in Fig. 10. For our sources and NGC 2023, the properties of the dense phases fall very near either of the relationships discussed, the dense component of M17 however falls well above even the matter bounded trend given by Oberst et al. (2011) by several orders of magnitude, this may be due to the dense clumps being independently gravitationally bound or transient objects – that would therefore have much higher densities than would be expected for pure pressure equilibrium between the clumps and their surroundings. Overall it seems that the simpler sources with a dominant exciting star (e.g., IRAS 23133+6050, S 106, NGC 2023, S 125, S 140) agree well with the trends, while the much more complex sources (e.g., W49A, M17, W3) have much more variance. However one of the most complex sources, 30 Dor, is consistent with both predictions, so the possible trend between simple/complex sources could easily be a manifestation of small number statistics.

It is expected that the flux ratio of the 158 μm [C ii] and the 63 μm [O i] lines should decrease with increasing G₀ and n as the [O i] line becomes the dominant source of cooling. We see this progression between our two sources, with IRAS 23133+6050 having a cooler ionizing star and lower G₀ and a correspondingly higher ${}^{\mathrm{[}\mathrm{C}{ii}\mathrm{]}}\mathit{/}_{\mathrm{[}\mathrm{O}{i}\mathrm{]}}$ ratio than S 106 by around a factor of two. The two values we find for this ratio fit nicely into the progression seen in Table 9, with lower values for the H ii regions, (with the puzzling exception of W49A although this may be due to the large physical size of W49A in comparison to the other objects) and higher values for the RNe PDRs. As with the comparisons between the radiation field and densities of the regions, the closest Galactic analog of IRAS 21333+6050 we measured is actually the reflection nebula NGC 2023. This likely arises because the O8.5-9.5 star illuminating IRAS 21333+6050 has a surface temperature relatively close to that of the B1.5 central star of NGC 2023. On the other hand, S 106 is closer to the H ii region PDRs listed, and actually shows the same line ratio as the Orion bar region.

In context of other Galactic PDRs then, on the basis of the spectral types of the exciting stars alone the PDRs studied here might be expected to have properties between those of the massive star-forming regions and clusters and the RNe PDRs, and this is indeed what we find. With both regions having properties similar to the PDRs generated by the hottest B stars, especially IRAS 23133+6050, which is very similar to the PDR of the RNe NGC 2023. The general consistency between our two objects can also be explained using similar logic. Both sources have similar H ii region electron densities and sizes, and similar types of exciting stars, around O8-9, although from the literature presented in Sect. 2 it seems likely that the star illuminating S 106 is slightly hotter by about one subtype. Despite this, both sources have very different morphologies, with S 106 having a very complex dual lobed H ii region and a small stellar disk region surrounded by a large extended PDR, while IRAS 23133+6050 is still an ultra-compact H ii region yet to break out of its natal molecular cloud. This dichotomy poses an interesting problem and might constrain models for the evolution and dynamics of H ii regions and PDRs.

6.2. S 106

S 106 in particular has a rich vein of literature and prior studies, of which a few discuss its PDR specifically (for example van den Ancker et al. (2000) and Schneider et al. (2003) listed in Table 9). We discussed these studies in terms of comparing the measured surface brightnesses for the various PDR cooling lines in Sect. 3.3.1, here we discuss the overall results of those studies and how they compare to our own. In general there is broad agreement between the literature findings and our results for S 106, in that there is a clumpy PDR consisting of dense clumps irradiated by a strong UV field producing the high-J CO lines, along with a less dense PDR phase permeating the region. We have extended these findings by presenting the full CO ladder available to the PACS and SPIRE Herschel instruments, which have revealed the presence of a further, even higher density, component. The key difference between the prior studies and our results here is that we are probing only the innermost region of S 106 (a circular aperture of radius 21.5′′), whereas van den Ancker et al. (2000) used ISO-LWS spectra of the central region with a beam some three to five times larger than our Herschel observations. The van den Ancker et al. (2000) concluded that the S 106 PDR was of very high density and experienced an extremely intense UV radiation field (G₀ > 10⁵) despite measuring a similar set of the common PDR cooling lines as we presented earlier from the Herschel observations and found evidence of a more moderate environment. The key difference arises because van den Ancker et al. (2000) also considered ISO-SWS observations including the [Si ii] 34.8 μm line, the strength of which could only be explained by invoking stronger radiation fields and higher densities.

A subsequent study by Schneider et al. (2003) spatially resolved the larger scale structure of S 106 PDR cooling line emission with pixel scales of 20−40′′ using the FIFI instrument on the KAO, and also including ISO-LWS spectra. Schneider et al. (2003) found density variations within the PDR in a similar manner to those described in this work, and very similar parameters for the two phases were recovered. The only slight difference arises in the density of the dense component, in which we found a value slightly above the range quoted by Schneider et al. (2003). We also find a value for the average UV radiation field experienced by the interclump medium of G₀ = 1.6 × 10⁴ while Schneider et al. (2003) give a lower limit of G₀ > 10³. The fact that we recover very similar values for the density and incident radiation field for the interclump medium in the region immediately adjacent to the star as was found by previous studies as averages for the entire region can then be easily explained by inspection of the maps presented by Schneider et al. (2003) which show that the emission from S 106 is very strongly peaked near the exciting star in the center of the region. By using the full Herschel CO ladder we have extended this picture to include a further excitation component, which may represent a further, yet higher density, PDR component or a shock – a conundrum we will discuss in the following section.

6.3. The diagnostic value of the CO ladder

We have used the observed CO fluxes in order to constrain the PDR properties of our sources. The shape of the CO ladder (as well as the curves in the rotation diagrams) provides the number of components with distinct physical conditions present in each source. In combination with the simple diagnostics which use the standard PDR lines, the ladder was then used to find the temperature of each component using RADEX models. Subsequently we used the CO fluxes to guide the application of numerical PDR models – leading to a reasonable fit to all of the other observables. However, there are several complexities to this picture, firstly: to what extent do the results of the rotation diagram analysis make sense in terms of our more complex models? It is not clear that the level populations are in LTE while the lower-J lines may be optically thick and the temperature distribution is not constant, so it might be expected that the results we find for the rotation diagrams compare poorly to the more detailed approaches. Secondly: we have detected the presence of an extra excitation component in the high-J S 106 CO lines which could be ascribed to either very-high density/thermal pressure PDR, or a strong shock.

The rotation diagram fits to both sources revealed very similar rotation temperatures of around 190 K (188 and 190 K for IRAS 23133+6050 and S 106 respectively) for the mid-J¹²CO lines. In Sect. 4.3.1, we compared these temperatures to those found using RADEX fits to the respective ¹²CO ladders. Using a density of 10^5.6 cm^-3 as indicated by the diagnostic diagrams gave temperatures of around 100 K warmer. We also found reasonable fits to the CO ladders using RADEX with higher densities and a temperature of 200 K. In light of the results of our numerical models, these results are actually reasonable estimates. For example, for IRAS 23133+6050 we found that the peak CO line emission emanated from a region of T = 217 K at a density roughly three times higher than that predicted by the simple diagnostic diagrams. This concurrency suggests that the skepticism with which we referred to the rotation diagram figures was perhaps unwarranted; in fact, for these data one could use the rotational temperature found using the rotation diagrams as a guide for using RADEX models to find the PDR density in the vicinity of the CO molecules.

In a more general sense the CO ladders we find are typical, and are routinely found in extragalactic contexts. This is because any observation of a star-forming galaxy will inevitably include star-forming complexes like S 106 and IRAS 23133+6050. However, such observations may well include regions dominated by other physical processes, including analogues of PDRs dominated instead by X-rays or cosmic rays (XDRs or CDRs), PDRs with enhanced mechanical heating (mPDRs) or an environment dominated by shocks (e.g., van der Werf et al. 2010; Rangwala et al. 2011; Rosenberg et al. 2014 and references therein). Each of these processes has some effect on both the relative strengths of the CO ladder and the fine-structure lines. The challenge in interpreting such observations then is to tease apart the different effects and attempt to create a coherent model of the environment. In contrast our sources are relatively simple and do not seem to be dominated by the wide variety of other physical processes that can be important for star-forming galaxies or AGN. Therefore, the extreme conditions we infer in our sample represent a baseline for the more complex sources, in that such conditions would be expected in the PDRs associated with massive star-formation in these extragalactic sources.

There appear to be two plausible explanations for the extra heating in the very high-J (J> 20) CO lines seen in S 106; shocks or a very high density PDR. Previous work on S 106 with respect to shocks is somewhat conflicting. van den Ancker et al. (2000) found that the IR emission from S 106 was dominated by PDR emission on the basis of the [S i] 25.2 μm/[Si ii] 34.8 μm line ratio, which is usually thought to be representative of shocks. However, the van den Ancker et al. (2000) study used ISO-SWS observations with a comparatively large beam - possibly implying that the environment is dominated by PDR material, rather than exclusively PDR material. A later spatially resolved study by Noel et al. (2005) found evidence of thermalization in near-IR H₂ line ratios, which could be evidence of shocks in the S 106 region, which would be included in the beam of our observations. Such ratios could also be produced by a very high density PDR (e.g., Burton et al. 1990a), although Noel et al. (2005) conclude that this cannot be the case for one of the thermalized regions, as it lacks detections of Brγ, PAHs and continuum emission. As we noted earlier, shocks are common in star-formation regions of most masses (e.g., Santangelo et al. 2013; Karska et al. 2014). The bipolar morphology of S 106 lends further credence to the idea that shocks are present in S 106 and are contributing to the extra excitation observed in the CO ladder.

In contrast, the evidence supporting the alternative view – that the extra excitation arises from the high density PDR itself – is somewhat limited and circumstantial. If taken at face value, the center panel of Fig. 9 and the negligible differences in terms of the fine-structure line fluxes between the two S 106 models in might appear to support the idea that the emission arises from a very high density PDR. However, we found that a good fit to the high-J lines could be only reached using an additional PDR model with a density of 10⁸ cm^-3 and G₀ of 1.6 × 10⁵. Such densities could perhaps be plausible in two circumstances. Firstly, the small circumstellar disk around IRS4 in the core of S 106; secondly, small dense clumps interspersed within the main PDR. However, both of these scenarios could be problematic. The first possibility would require a much higher values of G₀ than used in our model by virtue of being much closer to the primary star, however the second circumstance would require small clumps of PDR material with densities several orders of magnitude beyond any previously observed. Clearly further studies will be required to decisively settle this issue. Nevertheless, the general shape of the CO ladder fluxes is very similar to Fig. 5 of Burton et al. (1990b), in which at high densities and radiation fields, the high-J CO lines can be excited in a thin layer on the surface of PDRs at temperatures of several hundred K. The transition between the two regimes is not at precisely the same J_U, however the shape of the curve is almost identical.

7. Conclusions and summary

We have obtained full Herschel PACS and SPIRE spectra of the star-forming regions S 106 and IRAS 23133+6050 and analyzed them in terms of the emission from their PDRs. Simple PDR diagnostic diagrams were employed to find the physical conditions and hints of the physical structure of the two PDRs. These results were tested using a detailed numerical model of PDR physics. The study reveals that the PDRs around S 106 and IRAS 23133+6050 are remarkably similar in some respects, despite being dramatically different types of objects – S 106 is a large, complex star-formation region, with many components, while IRAS 23133+6050 is a relatively simple cometary ultra-compact H ii region. Despite these differences, both objects were found to have a similar range of physical conditions and structures, in common with many other PDR observations and models. The structure of each region, implied by our observations, is that of a classical PDR with intermediate properties (G₀ = 10⁴, n = 10⁴ cm^-3) and with embedded clumps of much denser material (G₀ = 10⁵, n = 10⁶ cm^-3). We employed numerical models which were fitted to the ¹²CO ladder fluxes, initially to determine the temperature of the PDR environments from the peak of the CO ladder, as well as to investigate whether the conditions indicated by the simple diagnostics were consistent with those indicated by a full numerical treatment. During this process it was noticed that the S 106 CO ladder possesses an extra excitation component at the high-J lines. This component is only reconcilable with PDR models upon invoking extreme densities, so we tentatively ascribe shocks as a probable cause. Ultimately though, both PDRs studied possess very similar properties to those of the Orion complex, which is surprising as Orion and the two objects studied here are in very different evolutionary stages. The youngest object, IRAS 23133+6050, is still embedded in its natal cloud, while in Orion the cluster of exciting stars have emerged from this phase and are steadily eroding their environment. The consistency of properties on various PDRs between objects and in time underlines their importance as a key phase of the ISM.

Acknowledgments

D.J.S. thanks the referee, N. Schneider, for a very thorough and thoughtful referee report which has greatly improved the discussion sections of the paper. D.J.S. also thanks E. Polehampton at Rutherford Appleton Laboratories for his assistance in reducing the SPIRE data and understanding the Herschel SPIRE beam profiles. D.J.S. and E.P. acknowledge support from an NSERC Discovery Grant and an NSERC Discovery Accelerator Grant. M.G.W. was supported in part by NSF grant AST-1411827. Studies of interstellar chemistry at Leiden Observatory are supported through advanced- ERC grant 246976 from the European Research Council, through a grant by the Dutch Science Agency, NWO, as part of the Dutch Astrochemistry Network, and through the Spinoza premie from the Dutch Science Agency, NWO. C.B. is grateful for an appointment at NASA’s Ames Research Center through San José State University Research Foundation (NNX14AG80A). J.C. acknowledges support from an NSERC Discovery Grant. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

All Tables

All Figures

	Fig. 1 Pointing for S 106 (left) and IRAS 23133+6050 (right). Overlaid on each GLIMPSE (Benjamin et al. 2003; Churchwell et al. 2009) IRAC [8.0] image is the corrected SPIRE beam aperture (43′′ and 37′′ respectively). We adopt the convention of coloring points relevant to S 106 and IRAS 23133+6070 red and blue respectively, beginning with the apertures in this figure.
In the text

Fig. 2 Full PACS + SPIRE spectra for IRAS 23133+6050 (top, blue) and S 106 (bottom, red, multiplied by 0.01 for clarity). In addition we include the IRAS photometric points for each source (triangles), and for IRAS 23133+6050 also the WISE photometry (Wright et al. 2010) as squares and the SCUBA observations of Thompson et al. (2006) as diamonds. A modified blackbody fit for each spectrum is shown as a dashed line.

In the text

Fig. 3 PDR diagnostic diagrams for observed fine-structure lines – based on the PDR models of Wolfire et al. (2010) and Hollenbach et al. (2012). Top: [C ii] 158 μm/[O i] 63 μm versus ([O i] 63 μm + [C ii] 158 μm)/ΣFIR. Middle: [C ii] 158/[O i] 145 versus ([O i] 145 μm + [C ii] 158 μm)/ΣFIR. Bottom: [C i] 370/ [C i] 609 versus ([O i] 145 μm + [C ii] 158 μm)/ΣFIR with the density grid extended down to log n = 2. In each figure the blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively.

In the text

	Fig. 4 PDR diagnostic diagrams for ratios of 12CO(J=15-14), 12CO(J = 7−6) and 12CO(J = 4−3) versus ([O i] 63 μm + [C ii] 158 μm)/ ΣFIR. In each figure the blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively. The data points indicate a much denser environment with a stronger UV field than that seen for the fine-structure line diagnostics.
In the text

	Fig. 5 [C i] 609 μm/12CO (J = 4−3) versus ([O i] 63 μm + [C ii] 158 μm)/ ΣFIR PDR diagnostic diagram. The blue and red points represent the values found for IRAS 23133+6050 and S 106 respectively.
In the text

	Fig. 6 12CO (squares) and 13CO (triangles) ladder fluxes for IRAS 23133+6050 and S 106 in terms of the excitation energy of that level.
In the text

	Fig. 7 Combination of RADEX models scaled to fit the S 106 12CO ladder. The low density component with log n = 4, T = 300 K is represented by the green line, while the dominant contribution arises from the dense component with log n = 5.6 and a temperature of 400 K which is represented by the purple line. The aggregate model (black line) matches the observed 12CO ladder from JU = 4 to 20.
In the text

Fig. 8 12CO (diamonds) and 13CO (triangles) rotation diagrams for IRAS 23133+6050 and S 106. Two ultra-compact H ii regions from the Karska et al. (2014) sample which display the warm CO component are also shown as orange, purple lines. Simple fits to the lower-J CO lines are shown to demonstrate the deviation of the higher-J CO levels from this trend for S 106 and the Karska sources. A fit to the highest-J CO lines for S 106 is also included. The Karska sources have been offset in y for clarity. The first 12CO and 13CO points for both of our sources are 4−3 and 5−4 respectively.

In the text

Fig. 9 PDR model fits to 12CO ladders. a) Two component PDR model for S 106; b) three component PDR model for S 106; c) two component PDR fit for IRAS 23133+6050. In each panel the dashed line represents the low pressure component and the dot-dashed lines the high pressure component. In the centre panel the final, very high density PDR component for S 106 is represented by a triple-dot-dashed line.

In the text

Fig. 10 Correlations between G0 and density n0 found by Young Owl et al. (2002) and Oberst et al. (2011). Reflection nebula points from Young Owl et al. (2002), Galactic PDRs as stated in Table 9. IRAS 23133+6050 and S 106 fit into the trend between the larger H ii regions with hotter exciting stars and the RNe PDRs with B type exciting stars. Dashed and solid black lines: theoretical prediction for the relationship between radiation field and density for a pressure bound H ii region and the best fit to the Young Owl et al. (2002) data. Green dot-dashed line: best fit to the Carina nebula data of Oberst et al. (2011). For a selection of sources we have also included the points representing the dense PDR component (denoted “dc”) as opposed to the default values representing the interclump medium (denoted “ic”).

In the text