Open Access
Issue
A&A
Volume 615, July 2018
Article Number A129
Number of page(s) 20
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201832611
Published online 27 July 2018

© ESO 2018

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Photodissociation regions (PDR) are key regions in the study of the interstellar medium. PDRs are the interfaces between molecular gas (where stars form), and the surrounding galactic medium (see review by Hollenbach & Tielens 1999). Dissociating UV photons produced by young stars are absorbed in PDRs by dust and gas allowing a transition from the atomic phase to the molecular phase. Understanding the structure, and physical and chemical processes in PDRs is necessary to constrain stellar feedback on molecular clouds, but also to be able to interpret observations of distant galaxies where the contribution of unresolved PDRs dominate the IR spectrum. Intense emission from fine-structure lines of C+, O, and C, as well as H2 rotational and rovibrational transitions and CO rotational transitions can be observed in PDRs. It is now admitted that emission in these atomic and molecular lines is mainly induced by the heating of the gas by UV photons from nearby massive stars involving the photo-electric effect on grains (Bakes & Tielens 1994; Weingartner & Draine 2001) as well as the collisional deexcitation of H2 excited by UV pumping (Sternberg & Dalgarno 1989; Burton et al. 1990; Röllig et al. 2006). Several stationary PDR models have been developed to analyse the emission in these lines (Tielens & Hollenbach 1985; Sternberg & Dalgarno 1989; Kaufman et al. 1999; Le Petit et al. 2006) and to constrain the chemical and physical processes that take place in PDRs. Over the years, these models simulating radiative transfer, chemistry, and thermal processes of the gas and dust have progressed in their description of the microphysics and chemical processes at play in these regions (Röllig et al. 2007). Today, some of these codes, e.g. the Meudon PDR code, can simulate very detailed micro-physical processes such as input and output line, and continuum radiative transfer, non local pumping by infrared photons, and detailed surface chemistry with stochastic processes (Goicoechea & Le Bourlot 2007; Gonzalez Garcia et al. 2008; Le Bourlot et al. 2012; Bron et al. 2014, 2016).

The analysis of line emission in PDRs is intimately connected to considerations of the morphology of these regions. Earlier observations of a few mid- and high-J CO lines from ground-based submillimeter and airborne far-IR observations have suggested that PDR interfaces reach high temperatures and densities both in the atomic and molecular gas, which requires the interface to be clumpy or filamentary (Stutzki et al. 1988). This led Burton et al. (1990) to propose a clumpy PDR model in which dense clumps (nH ~ 106–107 cm−3) are embedded in a less dense interclump medium (nH ~ 103–105 cm−3). In following studies on the bright PDRs M17 SW, Orion Bar and Carina, clumpy PDR models were used to analyse the observations quite successfully (Meixner et al. 1992; Hogerheijde et al. 1995; Kramer et al. 2008; Andree-Labsch et al. 2017).

The Herschel Space Observatory, with its three instruments HIFI, SPIRE, and PACS (Pilbratt et al. 2010), has opened the possibility to observe, more systematically and continuously in wavelengths, the warm molecular gas in galactic and extragalactic sources by covering all CO excitation lines from Jup = 4 to Jup = 50. This allows us to build full CO spectral line energy distributions (SLED) including high-J levels. Such CO SLEDs have been particularly studied in star-forming regions in order to provide information on the energetic processes acting in these objects. This includes low- and high-mass protostars (Yıldız et al. 2010; Visser et al. 2012; Manoj et al. 2013; Karska et al. 2013; Goicoechea et al. 2015a) and a couple of PDRs associated with young stars of high or intermediate mass (Köhler et al. 2014; Stock et al. 2015). Studies were also performed in the Galactic centre with observations towards Sgr A* (Goicoechea et al. 2013), as well as in the Sgr B2 cores, B2(M), and B2(N) (Etxaluze et al. 2013). In external galaxies, CO SLEDs have been obtained for Seyfert galaxies, starburst galaxies, (ultra)luminous infrared galaxies, and active galactic nuclei (Rangwala et al. 2011; Hailey-Dunsheath et al. 2012; Greve et al. 2014; Kamenetzky et al. 2014; Mashian et al. 2015; Rosenberg et al. 2015; Wu et al. 2015).

Modelling CO SLEDs in protostars or in active regions of galaxies is complicated by the fact that mechanical heating due to shocks (supernova explosions or stellar winds) is likely to be a source of energy powering this CO emission (see in particular Kazandjian et al. 2012, 2015). In PDRs, UV photons are expected to be the major actor and these CO lines offer the possibility of further constraining the gradients in gas density and temperature, as well as the underlying excitation processes. Using a PDR model, Stock et al. (2015) analysed the full CO SLEDs of two PDRs generated by massive star formation, S 106 and IRAS 23133+6050. Their best results were obtained with a two-component model comprising high-density clumps (nH ~ 106 cm−3) immersed in a strong far-UV radiation field (G0 ~ 105 in units of the Habing field; Habing 1968) and an interclump medium at lower density and irradiation field (nH ~ 104 cm−3, G0 ~ 104). However, the high value derived for the UV field on the clumps (a factor of 10 higher compared to the interclump medium) is striking.

In this work, we have tried to bring some new insights to this subject. We present a study of two prototypical PDRs, NGC 7023 NW and the Orion Bar. The Orion Bar is probably the most studied PDR (Tielens et al. 1993). NGC 7023 NW is well known to exhibit bright H2 filaments at the interface with the illuminating star (Lemaire et al. 1996). The two objects have been extensively studied in the past but often using a limited set of tracers. Here, we benefit from the Herschel HIFI, SPIRE, and PACS data. In order to include further tracers of the warm or hot molecular gas, we complete this dataset with ancillary data coming from ground-based instruments as well as by ISO and Spitzer space missions. The full data sets include mid- to high-J CO, H2, CH+, [CII], [OI], [CI], HD, OH and HCO+ lines. We have analysed these observations using the latest version of the Meudon PDR code (Le Petit et al. 2006). Our goal is to achieve better insights on the structure of the irradiated interface at the border of PDR and to determine whether UV photons alone can explain the high-J excited CO lines observed in PDRs.

The article is organised as follows. In Sect. 2, we compile information from the literature on the two PDRs of interest. The observations are described in Sect. 3, including new Herschel data as well as complementary data gathered from the literature or archives. The CO observations are described in Sect. 4, including both line profiles and intensity ladders. In Sect. 5, we describe the Meudon PDR model and the procedure used to fit the observational data, and present our modelling results. We used isobaric models to mimick the density gradient at the interface. In Sect. 6, we revisit the impact of our analysis on the interpretation of the emission in the mid- to high-J CO lines and conclude on the presence of sharp PDR interfaces at high thermal pressure that were likely shaped by the UV radiation field.

2. Prototypical PDRs: NGC 7023 and Orion Bar

2.1. NGC 7023

The source NGC 7023 is a reflection nebula in the Cepheus constellation illuminated by HD 200775 [RA(2000) = 21h01m36.9s; Dec(2000) = +68°09′47.8″], a spectroscopic binary system (B3Ve and B5, see Alecian et al. 2008). Its distance from the Sun was measured by HIPPARCOS at pc in the 1997 catalogue (van den Ancker et al. 1997). The new reduction by van Leeuwen (2007) gave 520 ± 150 pc, whereas Benisty et al. (2013) have proposed a distance of 320 ± 51 pc based on a study of the orbital parameters of the spectroscopic binary system HD 200775, which we adopt in the following. At this distance, 1″ corresponds to a physical length of 1.6 × 10−3 pc.

Chokshi et al. (1988) observed the [CII] 158 µm and [OI] 63 µm lines in NGC 7023 and derived at the emission peak located 50″ NW from the star a UV field of G0 = 2600 and a density of nH ~ 4 × 103 cm−3. Later observations have shown that the star formation process has shaped a cavity inside the molecular cloud with walls of dense PDRs at the north-west (NW), south, and east (Fuente et al., 1992, 1998; Rogers et al. 1995; Gerin et al. 1998). Very bright thin filaments were revealed by high angular resolution images in the extended red emission (ERE), vibrationally excited H2 emission lines (Sellgren et al. 1992; Lemaire et al. 1996; Witt et al. 2006) and in HCO+ millimeter lines (Fuente et al. 1996a). These structures also correspond to enhanced emission in polycyclic aromatic hydrocarbon (PAH) emission (An & Sellgren 2003; Berné et al. 2007), in [O I] (63 and 145 µm) and [CII] (158 µm) lines (Bernard-Salas et al. 2015), and in warm CO and dust emission (Rogers et al. 1995; Gerin et al. 1998; Köhler et al. 2014; Bernard-Salas et al. 2015). From these observations emerges a picture of the NGC 7023 NW morphology in which the PDR interface is made of filamentary structures at high-density, nH ~ 105 − 106 cm−3 (Martini et al. 1997; Lemaire et al. 1996; Fuente et al. 1996a), which are embedded in a more diffuse gas with nH ~ 103–104 cm−3 (Chokshi et al. 1988; Rogers et al. 1995). This filamentary structure is observed at small spatial scales of 0.004 pc or less. Whether it is composed of detached filaments or the result of compressed sheets is still unclear (Lemaire et al. 1996; Fuente et al. 1996a). Because of its geometry, brightness, and proximity, the NW PDR of NGC 7023 turns out to be one of the best sites to study the physical and chemical processes taking place in a PDR.

2.2. Orion Bar

The Orion Bar PDR lies ~2′ SE of the Trapezium stars cluster: θ1 Orionis C, A and E (Simón-Díaz et al. 2006; Allers et al. 2005), a cluster of O and B stars that creates a HII region that penetrates into the parent molecular cloud. The UV intensity impinging on the PDR has been estimated to be G0 = 1–4 × 104 in Habing units (Tielens & Hollenbach 1985; Marconi et al. 1998). The distance of the Orion nebula has been measured with great precision by Meixner et al. (2007) from trigonometric parallax, yielding a value of 414 ± 7 pc. An angular distance of 1 arcsec corresponds therefore to 2 × 10−3 pc. Because of its proximity and edge-on geometry, the Bar is one of the most studied PDRs. It is the prototype PDR associated with massive-star formation, which can be used as a template for more distant regions including extragalactic studies. Hogerheijde et al. (1995) reported spatial observations of the Bar in rotational transitions from a variety of molecules and concluded that the morphology of the molecular emission is mainly due to the geometry of the Bar, which changes from face-on to almost edge-on and then to face-on again. The bright PDR corresponds to the edge-on part. The overall observed spatial stratification of the Bar was found to require an average density of at least 5 × 104 cm−3 in order to account for the observed offsets of ionisation front and molecular lines (see also Wyrowski et al. 1997). Most models of the molecular emission in the Orion Bar used high-density clumps (nH ~ 106–107 cm−3) embedded in a lower density gas (nH = 5 × 104–105 cm−3). Clumps were introduced first to explain the emission of excited lines of CO and warm H2 emission (Parmar et al. 1991; Tauber et al. 1994; van der Werf et al. 1996) but were also found to be necessary to model the HCO+, HCN (Young Owl et al. 2000) and OH (Goicoechea et al. 2011) emission at the PDR interface.

The interface with the HII region is of special interest if one wants to study the feedback of star formation on the molecular cloud. Omodaka et al. (1994) have discussed that the Bar has been shaped by shock compression related to the expansion of the HII region. Giard et al. (1994) mapped the 3.3 µm PAH emission and revealed a clumpy structure of the gas down to scales of a few 10−3 pc behind the ionisation front and in front of the molecular interface traced by vibrational H2 (Tielens et al. 1993). This is also well illustrated in Fig. 1, which combines the IRAC 8 μm map dominated by PAH emission with the H2 map of van der Werf et al. (1996). Fuente et al. (1996b) observed the molecular interface and concluded that it consists of a corrugated dense layer (nH ~ 106 cm−3). Walmsley et al. (2000) proposed that it can be described by a filament using a cylindrical shape. Andree-Labsch et al. (2017) showed that this shape can be excluded due to the visible shadowing pattern seen in the maps. Recently, the Orion Bar has been observed in CO and HCO+ by ALMA (Goicoechea et al. 2016). The sharp edge (~2″) at which the CO and HCO+ emission start to occur coincides with the bright H2 vibrational emission. The key chemical transitions of the PDR including the conversion from H to H2 as well as that from C+ to C and then CO, which are referred in the following as the H/H2 and C+/C/CO transitions, should therefore be very close at this interface. This is incompatible with stationary PDR models at nH ~ 5 × 104 cm−3.

thumbnail Fig. 1.

Composite images of the Orion Bar (top panel) and NGC 7023 (bottom panel). Red indicates the 8 μm emission observed with Spitzer. Green shows the vibrationally excited emission of H2 from Lemaire et al. (1996) and van der Werf et al. (1996) for NGC 7023 and Orion Bar, respectively. Blue shows the 12CO emission, J = 6 − 5 for the Orion Bar (Lis & Schilke 2003) and J = 1−0 for NGC 7023 (Gerin et al. 1998). The circles represent the HPBW of Herschel at 550 GHz and 1900 GHz. The square indicates the position of the central spaxel of the PACS observations.

Open with DEXTER

3. Observations

The observations of NGC 7023 NW were programmed in the framework of the WADI (Ossenkopf et al. 2011) Guaranteed Time Key Program (GTKP) and were centred on the bright H2 filaments, a position referred to as the H2 peak (Joblin et al. 2010, α2000 = 21h01m32.4s, δ2000 = +68°10′25.0″). The Orion Bar was observed at the CO+ peak position (Stoerzer et al. 1995, α2000 = 5h35m20.61s, δ2000 = −5°25′14.0″; cf. Fig. 1) as part of the HEXOS GTKP (Bergin et al. 2010). Earlier studies based on these data have discussed OH emission (Goicoechea et al. 2011) as well as CH+ and SH+ lines (Nagy et al. 2013). The WADI and HEXOS programmes gathered spectroscopic data using HIFI and PACS instruments. In this article we have also used spectroscopic data from SPIRE that were obtained on the same objects as part of the SAG 4 GTKP (Habart et al. 2010; Köhler et al. 2014). Data from the literature and archives were also collected – particularly for H2. The full data sets used in this study are reported in Tables 1 and 2, for NGC 7023 and the Orion Bar respectively.

Table 1.

Observed data for NGC 7023, and dilution factor Ω.

Table 2.

Observed data for the Orion Bar and dilution factor Ω.

3.1. Herschel PACS observations

The PACS range spectroscopy observations of NGC 7023 and the Orion Bar were reduced using HIPE version 10. We used the standard pipeline to extract the spectrum from the central spaxel for the blue and red channels, including defringing. We applied the correction for point sources, considering that the high-J CO emission originates from dense structures of arcsec size and therefore smaller than the spaxel size, which is 3.2 × 3.2 arcsec2 (cf. Sect. 4.1). This is still an approximation, since the observed filamentary interface extends over several spaxels but there is no correction available for a semi-extended source. We checked that the intensities of the lines of interest remain very similar when performing the reduction with a more recent version of HIPE (e.g. v 13).

Line flux was calculated by fitting the lines observed on the central spaxel with a Gaussian function. Error bars were calculated by quadratically summing the different sources of errors: calibration error and error from the Gaussian fit associated with spectral rms. We used a 20% error for the flux calibration, which is intermediate between the values of 10–30% used in previous studies (Bernard-Salas et al. 2012, 2015; Okada et al. 2013).

3.2. Herschel HIFI observations

The HIFI observations were obtained using the Wide-Band Spectrometer, that provides a spectral resolution of 1.1 MHz (~0.6 km s−1 at 550 GHz). The half power beam width (HPBW) varies between 9″ at high frequency (1900 GHz) and 39″ at low frequency (550 GHz).

For NGC 7023, the observations were focused on specific frequency ranges, and cover a number of CO lines, CH+ (J = 1–0 and 2–1), HCO+ (J = 6–5) and the [CII] lines. The data reduction for NGC 7023 was straightforward and consisted in the subtraction of a linear baseline for each scan and then in averaging all the scans, including both the H and V polarisations. The data reduction procedure for the [CII] line is described in Ossenkopf et al. (2013). In the case of the Orion Bar, our observations consist of a spectral survey performed at the CO+ peak presented in Nagy et al. (2013) and further discussed in Nagy et al. (2017).

Line intensities for both sources were calculated by using the beam efficiencies of Roelfsema et al. (2012) and by performing a Gaussian fit to the line profiles. Ossenkopf et al. (2013) argued that TA is more appropriate for extended emission, and Tmb for point sources. Emission in some lines is extended (e.g. [CII]), whereas it is not in others (e.g. high-J CO). In this work, we use the mean value of TA and Tmb because the bright interface is in between point-source and extended. In addition, it is impractical to have a different treatment for every individual line. Finally, the fit error being negligible (<1%), our observational uncertainty is defined as the quadratic sum of the spectral rms and the flux calibration error, which is reported in Roelfsema et al. (2012). A sample of the observed 12CO and 13CO lines observed with Herschel-HIFI for both sources is shown in Fig. 2.

thumbnail Fig. 2.

12CO and 13CO lines observed with HIFI towards NGC 7023 (left) and the Orion Bar (right). The vertical red line indicates the systemic velocity of the source.

Open with DEXTER

3.3. Herschel SPIRE observations

We have used the SPIRE FTS fully sampled maps to extract complementary data at the observed positions. The data reduction is presented in details in Köhler et al. (2014) for NGC 7023, and in Parikka et al. (2017) for the Orion Bar. One important aspect is the use of the super-resolution method SUPREME (Ayasso et al., in prep.) to achieve higher spatial resolution than standard SPIRE observations, more specifically 11.9″ at 200 μm, 19.0″ at 400 μm and 24.1″ at 600 μm. The total error on the integrated line intensities was estimated to be 30%, which includes the calibration uncertainties and the line fitting errors (see Köhler et al. 2014, for more information).

3.4. Additional data

We used additional observations to further constrain our models. For NGC 7023, the pure rotational lines of H2 were observed with the low spectral resolution modules of Spitzer-IRS (Werner et al. 2004; Houck et al. 2004) and with ISO-SWS (de Graauw et al. 1996). For the Spitzer data (Werner et al. 2004), we used the CUBISM reduction tool (Smith et al. 2007) including the slit-loss correction function for extended sources, considering that the filamentary interfaces in Orion Bar and NGC 7023 are narrow but extended sources; their width is indeed marginally resolved at the Spitzer spatial resolution. The line fluxes were then extracted by fitting Gaussian profiles to the spectra extracted towards the H2 peak position over four pixels. Uncertainties take into account fit and calibration errors. The same lines were observed with ISO, and we used the intensity values reported in Fuente et al. (1999). The intensities for the vibrationally-excited lines of H2, observed at the CFHT and the Perkins Telescope, are taken from Lemaire et al. (1996, 1999), and Martini et al. (1999) respectively. Finally, we used the HCO+J = 1 – 0 observations obtained at the Plateau de Bure Interferometer (Fuente et al. 1996a), integrating line intensities between 1 and 5 km s−1.

For Orion Bar, we used pure rotational H2 lines from ISO-SWS data (Bertoldi priv. comm.; Habart et al. 2004) and ground-based TEXES data (Allers et al. 2005). The latter intensities are slightly lower compared to the ISO measurements. This can be due to the fact that the two instruments observed at different positions. The ISO-SWS data provides a position closer to the CO+ peak but lower spatial resolution. For vibrationally excited H2, we used data obtained with the BEAR instrument at the CFHT (Joblin, Maillard, Noel, unpublished) and the observations from van der Werf et al. (1996).

4. Observation results

4.1. Combination of observational data

For the purpose of comparing the observations with PDR models (cf. Sect. 5), we needed to combine the results obtained with different instruments. This is particularly challenging because of different calibrations, beam sizes and observational techniques. In this section, we discuss the procedure we used to deal with these differences.

Beam dilution factors. The data used in this work gather line intensities that were measured using different beam sizes. This raises the question of the morphology of the emitting structures and to which extend these structures fill the different beams. In NGC 7023, the bright NW interface is composed of a filamentary structure with a spatial characteristic width of a few arcsecs as clearly seen in vibrationally excited H2 emission (Lemaire et al. 1996). From these maps one can derive a relevant width of 2″ for the bright PDR interface that is centred towards the H2 peak. In the Orion Bar, the sharp molecular edge has been observed in vibrational and rotational H2 transitions (Tielens et al. 1993; van der Werf et al. 1996; Walmsley et al. 2000; Allers et al. 2005). These observations suggest a narrow (few arcsec) and patchy interface that is also seen in ALMA maps obtained at a spatial resolution between 1 and 5″ (Goicoechea et al. 2016, 2017). The ALMA data show that emission lines of some ions such as HCO+ 4–3 and SH+ 1–0 are located in the narrow layer given by vibrationally excited H2. Although the morphology is more complex than a single structure, we assume in the following that the observed emission from warm molecular tracers arises from a 2″ filamentary interface, consistently with the case of NGC 7023. To derive dilution factors, we therefore assumed for both PDRs that the emitting structure is a filament that follows the interface with infinite length and a 2″ thickness. For each observation, we calculated the fractional coverage of the beam by this filament. This is a simplistic procedure but the best we can do to overcome the lack of spatial information. The values of the derived dilution factors Ω are reported in Tables 1 and 2. The observed intensities were then divided by this factor Ω in order to be compared with the models.

Cross-calibration factors In addition to the dilution factors we needed to apply scaling factors to some of the data sets. The 12CO ladder for NGC 7023 (see the data in Fig. 3) reveals discrepancies between the different instruments that cannot be compensated by our dilution factors. In the following, we have considered that HIFI fluxes are the references and scale the intensities from SPIRE and PACS. As common lines are observed by SPIRE and HIFI, we simply searched for the scaling factor giving the least square error, and divided all line intensities observed by SPIRE (including species other than CO) by a factor of 0.54 as a result. PACS observations can not be directly compared, but the rotational diagram of 12CO reveals an unphysical jump between the PACS and SPIRE/HIFI lines. We determined the factor giving the best linear alignment of the PACS observations with the SPIRE/HIFI lines on this rotational diagram, and divided all line intensities observed by PACS by a factor of 1.3 as a result. In the case of the H2 observations in NGC 7023 it is not possible to satisfactorily merge the Spitzer and ISO data. This is because H2 emission is quite extended and the ISO beam contains emission from regions other than just the filament of interest. Spitzer values are taken as references and we divided all ISO observations by a factor of 2.54 to get the best agreement between the two data sets. In the case of the Orion Bar, we find no obvious reason for such adjustments. In particular no systematic discrepancy is visible in the CO ladder (cf. Fig. 3). Thus, no such adjustment was applied to the observations of the Orion Bar.

thumbnail Fig. 3.

Observed intensities of 12CO (left) and 13CO (right) in the Orion Bar (top panel) and NGC 7023 NW (bottom panel).

Open with DEXTER

Finally, after correcting for beam dilution and cross-calibration factors, we combined the line intensities from the different instruments by simply taking an average of the different observations of a given line. We computed error bars on this mean value as the interval between the minimum and maximum values in the error ranges of the different instruments.

The original (uncorrected) data of the different instruments are presented in Tables 1 and 2 and shown on Fig. 3 for 12CO and 13CO. Tables 4 and 5 present the data after correction of dilution, cross calibration and averaging of the different instruments. These corrected and averaged data are used in all figures except Fig. 3.

Table 3.

Input parameters used in the Meudon PDR code.

Table 4.

Comparison of models and observations (combined) for NGC 7023.

Table 5.

Comparison of models and observations (combined) for the Orion Bar.

4.2. CO rotational diagrams

Figure 4 presents the rotational diagrams and local thermodynamical equilibrium (LTE) fits of the 12CO and 13CO observations. The rotational diagrams present column density estimates in the upper rotational levels, Nu, without corrections for opacity effects with where I is the observed intensity, Aul and vul are the Einstein coefficient and frequency of the transition from the upper to lower levels and gu is the degeneracy of the upper level. As such, the derived column densities are lower limits when the lines are optically thick. In the case of 12CO, the fit was restricted to the highest excitation lines, starting at Jup = 15, in order to obtain information on the warmest gas. For the Orion Bar, we derived an excitation temperature of 147 ± 9 K and a column density of N(12CO) = (9.0 ± 3.9) × 1017 cm−2. In NGC 7023, the excitation temperature is 112 ± 6 K with N(12CO) = (1.7 ± 0.7) × 1017 cm−2. We note that the column density of warm CO is higher (at least a factor of 4) in the Orion Bar as compared to NGC 7023. The 12CO rotational temperature is also somewhat higher in the Orion Bar. Both facts favour the detection of higher-J CO transitions in Orion Bar relative to NGC 7023. We can compare these results with previous studies. Köhler et al. (2014) studied the CO lines measured with SPIRE in NGC 7023. Using the non-LTE radiative transfer code RADEX (van der Tak et al. 2007), they derived a kinetic temperature in the range 65–130 K and a column density of N(12CO) = 2–3×1018 cm−2. They concluded that the emitting structure has a dilution factor of 0.1. All these results are compatible with ours if we keep in mind our strategy to gain contrast on the hottest molecular interface, including taking into account higher-J CO lines. Our derived rotational temperature of 141 K for the Bar is also consistent with the kinetic temperature in the range 100–150 K that was reported by Nagy et al. (2017) in the analysis of the Herschel/HIFI spectral line survey of the Orion Bar mentioned above. It is also compatible with the lower limit of the kinetic temperature derived by Goicoechea et al. (2016) at the dissociation front.

thumbnail Fig. 4.

Rotational diagram of 12CO (left) and 13CO (right) lines observed in the Orion Bar (top panel) and NGC 7023 (bottom panel) PDRs. For 12CO, the excitation temperature and total column density are computed in the range J up = [15, 23] (Orion Bar), and J up = [15, 19] (NGC 7023). For 13CO, the lines used were J up = [5, 16] (Orion Bar), and J up = [5, 10] (NGC 7023).

Open with DEXTER

The 13CO diagrams reveal gas at an average rotational temperature of ~80 K, which is cooler than the 12CO temperature discussed above. In addition, Fig. 5 reports the value of the 12CO over 13CO line intensity ratio across the velocity profile of the 12CO 10–9 line. This ratio is found to be weak (~4) near the line centre, whereas it reaches the isotopic 12C to 13C ratio value of ~50 in the wings. This is characteristic of strong opacity effects. Köhler et al. (2014) calculated the optical depth of the 12CO lines in NGC 7023 using the RADEX code and concluded that the lines are optically thick up to Jup ~ 13–14. This is our justification for selecting only lines from Jup = 15 and higher in our fit of the 12CO rotational diagrams (cf. Fig. 4).

thumbnail Fig. 5.

Line profile for the J = 10–9 transition of 12CO (dashed line) and 13CO (dotted line) in the Orion Bar (top panel) and NGC 7023 NW (bottom panel). The relative intensity of the 12CO to 13CO lines across the velocity profiles is also shown (plain line).

Open with DEXTER

The above results strongly suggest that CO emission in the beam stems from a two-component medium: an extended cool or warm component and the hot and sharp interface. In 13CO, the former component is well seen, whereas it is partly hidden due to optical depth effects in 12CO emission. This results in 12CO emission seeming to stem mainly from the hot component. On the other hand, the hot component is difficult to observe in the high-J lines of 13CO due to the low expected signal. As an example, in the Orion Bar an intensity ratio of 5.6 can be derived from the observed J = 16–15 and J = 20–19 12CO lines (see Table 2). By applying the same factor to the 13CO lines, we can predict an intensity of 8 × 10−10 W m−2 sr−1 for the J = 20–19 13CO line, which cannot be detected as it is a factor of three weaker than the error bar of the J = 16–15 13CO line. In conclusion, only high-J 12CO lines can be used to characterise the physical conditions at the warm and bright interface of PDRs.

5. Models

5.1. Meudon PDR model

We compared our observations to PDR models using an updated version of the Meudon PDR code (Le Petit et al. 2006)1. This 1D PDR code simulates the physical and chemical processes at stationary state in a plane-parallel slab of gas and dust. At each position in the cloud, the code computes the temperatures of gas and grains, the chemical densities and, for the most important species, the non-LTE level populations. Then a post-treatment gives access to column densities and line intensities.

We updated several atomic and molecular data. In particular, it is worthwhile to notice for the present study that this version of the PDR code implements the CO–H2 collision rates of Yang et al. (2010) and the CO–H collision rates of Balakrishnan et al. (2002). The chemical network includes 213 species linked by 5067 gas-phase chemical reactions (except H2 formation on grains). Table 3 summarises the elemental abundances used in this paper. The formation reaction of H2 on grains plays a critical role in the computation of the chemical structure of PDRs. Several prescriptions are available in the Meudon PDR code: a mean formation rate of based on Copernicus and FUSE H2 observations in diffuse clouds (Jura 1974; Gry et al. 2002), Langmuir-Hinshelwood (LH) and Eley-Rideal (ER) mechanisms treated with a rate equation formalism (Le Bourlot et al. 2012) and a new stochastic approach that considers the impact of grain temperature fluctuations on H2 formation (both LH and ER) in a master equation formalism (Bron et al. 2014, 2016). In this paper, we have used the Le Bourlot et al. (2012) formalism2.

The code considers two external sources of energy, the radiation field and cosmic rays. The external radiation field can be the combination of a beamed stellar radiation field that simulates neighbouring stars and an isotropic ambient radiation field. This isotropic component is composed of the Mathis et al. (1983) field for the far UV to IR part, scaled by a factor, and black bodies that simulate dust IR emission and the cosmic microwave background. The Meudon PDR code allows us to use different radiation fields on each side of the cloud. We call the side illuminated by the stars responsible for the PDR the front side, and the other the back side. For NGC 7023, we used a beamed stellar radiation field built from the Kurucz stellar spectra (Kurucz 1993) as explained in Pilleri et al. (2012). This leads to a value of the UV intensity at the edge of the NW PDR of G0 = 2600 in Habing units for a distance of d ≃ 0.143 pc between the star and the dissociation front. In addition, the scattered light in the surrounding region was determined to represent about 4% of the direct stellar light using the scattered light measured at 475 nm by the Hubble Space Telescope (Witt et al. 2006). We thus assumed that the back side of the PDR is illuminated by an isotropic radiation field with G0 = 100, making the simplifying assumption that there is no dependence of the scattering with wavelength. This back side illumination was found to contribute marginaly to the line intensities. For the Orion Bar, several estimations of the UV flux have been proposed in the literature. As discussed by Allers et al. (2005), the estimated UV flux intensity depends on several parameters such as the inclination of the Bar. Based on Tielens & Hollenbach (1985) and Marconi et al. (1998), we fixed the radiation field impinging on the PDR so that, at the edge of the PDR, G0 = 2 × 104 in Habing units. We also assumed an ambient UV field illuminating the back side of the cloud with 10% of the front side illumination. Finally, the cosmic ray flux is introduced as a cosmic-ray ionisation rate ζ of H2 molecules. As we lack information about this flux in NGC 7023 and Orion Bar PDRs, we use an intermediate value of ζ = 5 × 10−17 s−1 that lies between estimations in diffuse and dense gas (Indriolo et al. 2015; Padovani et al. 2013; Le Petit et al. 2004; McCall et al. 1999). We checked that an increase of this ionisation rate by a factor of ten has no impact on CO line intensities.

The determination of photo-reaction rates requires to compute the specific intensity in the UV at each position in the cloud. The Meudon PDR code solves the radiative transfer equation on a wavelength grid (from the UV to the radio domain) at each position, including absorption and scattering by dust as well as absorption in lines or in the continuum for some chemical species. For instance, continuum absorption by C photoionisation turns out to have a non-negligible impact on the PDR structure (Rollins & Rawlings 2012) and is considered, as well as H line absorptions. For other species (as H2 and CO) line self-shielding is computed using the Federman et al. (1979) approximation. Absorption by grains is implemented using parametrised extinction curves (Fitzpatrick & Massa 1986). In PDRs, extinction curves are usually flatter in the far UV than the mean Galactic extinction curve. This is characteristic of larger-than-standard grain sizes. The adopted extinction curve is the one of HD 38087 in Fitzpatrick & Massa (1990) and RV = 5.62 as measured in NGC 7023 by Witt et al. 2006, close to the value determined for Orion Bar, 5.5 (Marconi et al. 1998). We also assume a column density to reddening ratio NH/E(B − V) = 1.05 × 1022 cm−2 mag−1, an intermediate value between the standard value 5.8 × 1021 and the value determined for ρ Oph, 15.4 × 1021 cm−2 mag−1 (Bohlin et al. 1978). Grains are simulated as a mixture of spherical amorphous carbonaceous and silicate grains following a MRN size distribution (Mathis et al. 1977), with minimum and maximum radii 3×10−7 and 3 × 10−5 cm, and with a dust-to-gas mass ratio of 0.01. The dust scattering properties are from Laor & Draine (1993).

At each position, the code computes the equilibrium gas temperature from the balance of total heating and cooling rates. The heating mechanisms considered are the photo-electric effect on grains, cosmic rays heating, and exothermic chemical reactions. Gas-grain collisions and H2 vibrational de-excitation can heat or cool the gas depending on the local physical conditions. For the photo-electric effect, we used the Bakes & Tielens (1994) prescription. Cooling rates were obtained by computing the radiative emission of the main coolants (15 species included), among which C+, O, C, CO, and its isotopologues. For these species, non-LTE level populations are computed taking into account collisional excitations and de-excitations, spontaneous emission, non-local radiative pumping by line and continuum photons, and chemical formation and destruction in specific levels as described in Gonzalez Garcia et al. (2008).

Finally, several prescriptions can be used to define the gas density profile as a function of depth. In this work, we tried two prescriptions: constant density models and constant pressure models. Parameters used in the models of NGC 7023 and Orion Bar are summarised in Table 3. The modelling strategy for the two objects is described in the following sections.

5.2. Fitting strategy

Since we focussed our analysis on the physical conditions at the PDR edge, we investigated the models that best account for the emission of tracers specific to this region, such as high-J 12CO and H2 lines. We investigated different scenarios. First, we ran PDR models at constant density but they proved to be incompatible with the observations. In particular, they were unable to reach the observed excitation temperature of high-J CO levels. We find that a density change across the PDR was necessary, with the density increasing with PDR depth. As earlier suggested by Marconi et al. (1998), constant pressure models are found to provide a satisfying density gradient due to the temperature drop when going towards the inside of the cloud.

Both PDRs appear as bright narrow interfaces (~2″) in vibrational H2 emission (cf. Fig. 1), which could be the result of an overdense surface layer seen roughly edge-on. We thus adopted for our models a geometry in which the PDR is observed with a high viewing angle of 60°; this angle being defined with 0° being face-on and 90° edge-on. The value of 60° gives an approximation of a nearly edge-on PDR and is the maximum inclination that can be used to derive line intensities in the 1D PDR Meudon code. The uncertainty on this angle could lead to an additional scaling factor on all line intensities. We thus allowed for a free global scaling factor f on the model intensities when fitting the model. A value of this factor larger than one would indicate a more edge-on configuration. In addition, this factor can correct for systematic errors we made on the assumed geometry, in particular viewing angle and the dilution factors.

We searched for the best fitting model in a grid of isobaric models for varying values of the thermal pressure P th and global scaling factor f. In both cases, we used the observed high-J 12CO lines from Jup = 11, the rotational H2 lines (S(0) to S(5)) and CH+ (1–0 to 6–5 in Orion Bar and 1–0 to 3–2 in NGC 7023) rotational lines. We thus have 17 (NGC 7023) or 23 (Orion Bar) observational constraints that the best model must simultaneously reproduce with two free-parameters, Pth and f. As these tracers only constrain the warm molecular layer of the PDR, we fixed a total AV value of 10 (NH ~ 2 × 1022 cm−2 for a face-on geometry and a factor of two higher assuming 60° inclination) in our models. We then compared the results of the best model with other available observations, i.e. lines from 13CO, C+, O, C, vibrational H2, HD, and in addition HCO+ in NGC 7023, and OH in the Orion Bar.

To provide an idea on how much the observations constrain the thermal pressure, we also report in the following the results obtained by using a pressure that differs by a factor of 1.5 (lower and higher) from the pressure obtained in the best fit model. For these cases, the scaling factor was adjusted to provide the best fit.

5.3. Model results

5.3.1. NGC 7023

The best fit was found for a model with Pth = 108 K cm−3 and with a global scaling factor of f = 0.7. Comparison of line intensities computed by this model to observed values is presented in Fig. 6. This model shows excellent agreement with the high-J 12CO lines (above Jup = 10), with the H2 pure rotational lines (except the S(3) line), and with the CH+ lines (except the (1–0) line). We also note that the model is able to fit both the ortho- and para-H2 lines and therefore can account for the observed non-equilibrium ortho-to-para ratio (Fuente et al. 1999, see Sect. 6.1 for further discussion). The discrepancy obtained for the H2 S(3) line (a factor of two brighter in the model than in the observations) likely indicates that the actual dust extinction in the line of sight towards NGC 7023 has a stronger silicate feature than assumed in our model (cf. Table 3).

thumbnail Fig. 6.

Excitation diagram of the different tracers observed in NGC 7023 after dilution correction (red squares) and the best fit model (black, Pth = 108 K cm−3, = 10, global scaling factor = 0.7). The best model has been chosen to optimise the fitting of the 12CO (high-J lines), H2 and CH+ lines only. In the last panel, the intensity values are normalised by the mean observed value for each line. The grey lines show the obtained variability when the thermal pressure is divided (dashed lines) or multiplied (plain lines) by a factor of 1.5. The best value for the scaling factor was found to be 1.2 for the model at Pth/1.5 and 0.47 for the model at Pth × 1.5.

Open with DEXTER

Further evidence of the adequacy of the model is given by its success in reproducing additional lines that were not used in the adjustment of the model (see also Fig. 6). The 12CO lines below Jup = 11 show an excellent agreement and the H2 vibrational lines (1–0) S(1) and (1–0) S(2) are reproduced within a factor <2. For HCO+, the model reproduces the J = 1–0 line within a factor of 2, but underestimates the J = 6–5 line by a factor of ~4. This line is very sensitive to the density due to the high dipole moment of HCO+. For instance for a gas kinetic temperature of T = 50 K a change of a factor of two in the density n(H2) leads to an increase of the J = 6–5 line by a factor of two while the J = 1–0 remains practically unchanged. This deviation suggests that the density provided by the isobaric model is good within a factor of a few. The [OI] 145 μm, [CI] 609 μm and HD J = 1–0 lines are well reproduced. Other lines on this same figure exhibit stronger differences of at most a factor of ~6: the [CII] 158 μm line is underestimated and the [OI] 63 μm line is overestimated. Herschel observations have shown that the [CII] emission comes not only from the H2 filaments but also from the surrounding more diffuse gas strongly emitting in the PAH features (Joblin et al. 2010). The adopted beam dilution factor is not justified in this case, and this correction could account for most of the discrepancy between the observed and calculated fluxes. The [OI] 63 μm line is overestimated in the model but this line is known to be optically thick and affected by self-absorption in most cases. Its study therefore requires a detailed analysis in velocity components. Finally, the model underestimates the 13CO lines by a factor of ~3. This is consistent with the arguments raised in Sect. 4.2, in particular that the 13CO emission arises from other regions inside the beam than the bright sharp interface.

Overall, we can conclude that the model is able to reproduce 17 observational constraints with only two free parameters, which is indicative of the adequacy of the model and of the physics it contains to describe the warm molecular edge region of the PDR. The grey lines in Fig. 6 demonstrate that the gas thermal pressure is best constrained by the high-J 12CO lines.

5.3.2. Orion Bar

The best fit was found for Pth = 2.8 × 108 K cm−3, and a global scaling factor of f = 1.3. The results of this model compared to the observations are presented in Table 5 and Fig. 7. The model provides a satisfactory fit of the high-J 12CO lines, CH+ lines, pure rotational lines of H2 (with maximal differences by a factor of two except for the H2 S(0) line that is underestimated by a factor of ~4). We thus obtain a reasonable agreement, within a factor of a few, for all the tracers used in the fitting procedure.

thumbnail Fig. 7.

Excitation diagram of the different tracers observed in the Orion Bar (red squares) and the best model (black, P th = 2.8 × 108 K cm−3, = 10 and global scaling factor f = 1.3). The best model has been chosen to optimise the fitting of the 12CO (high-J lines), rotational H2 and CH+ lines only. In the last panel, the intensity values are normalised by the mean observed value for each line. The grey lines show the obtained variability when the thermal pressure is divided (dashed lines) or multiplied (plain lines) by a factor of 1.5. The best value for the scaling factor was found to be 2.0 for the model at Pth/1.5 and 0.9 for the model at Pth × 1.5.

Open with DEXTER

Most of the other tracers (not used in the fitting procedure) show similar differences. The low lying 12CO lines (Jup < 11) are underestimated by up to a factor of ~3. The OH, HD, [OI] and [CI] lines are reproduced within a factor of 4. Some tracers, however, show stronger discrepancies such as the 13CO lines, which are underestimated by up to a factor of ~ 8 except for the highest line, and the [CII] 158 μm line that is underestimated by a factor of 16. For these lines at least part of the discrepancy is related to the correction by the beam dilution factor. We have discussed previously that some of the emission from 13CO is expected to arise from the surrounding molecular cloud (cf. Sect. 4.2). Emission in the [CII] line is also known to be extended in this region (see Goicoechea et al. 2015b). In the case of the vibrational H2 lines, the (1–0) S(1) line is strongly overestimated (factor of ~6) while the (2–1) S(1) line is underestimated by a factor of ~2 only.

To summarise, the agreement of the model with observations is not as good for Orion Bar as for NGC 7023 and leaves some areas of concern. On average, from the spatial structure that is reported in the next section we can derive that the emission lines coming from the most external layers are overpredicted by the model (CH+, H2 lines except low H2 rotational lines), whereas those arising from more internal layers (low- and mid-J CO lines, low H2 rotational lines) are underpredicted. This suggests that a single pressure is not sufficient to describe the evolution of gas density and temperature across the hot to warm irradiated interface (see also the grey curves in Fig. 7).

5.3.3. Spatial structure

In Figs. 8 and 9, we show the spatial structure of the calculated PDRs in terms of the gas temperature and H nuclei number density (upper panels), of the abundances of the different species (lower panels), and of the emission regions of the different lines (middle panels). Our calculations were performed on a slab of gas of AV = 10. However, we provide here information for the region up to AV ~ 5 in which the tracers we selected emit. Although the thickness and absolute position of the transitions depend on the characteristics of each PDR, the general trends can be summarised as follows: the atomic part has a density of a few 104 cm−3 (close to 105 cm−3 in Orion Bar) followed by a hot to warm (T = 1000 − 100 K) and relatively dense (few 105 to few 106 cm−3) molecular part. In the Orion Bar model, the atomic region is found to be significantly more extended (~6 × 10−3 pc) than in NGC 7023 (~0.5 × 10−3 pc). It has been truncated on Fig. 9 for better readability. In both PDRs, the hot to warm molecular region is found to start significantly before the H/H2 transition. Indeed a sharp increase of the H2 abundance is seen at 0.6 × 10−3 pc whereas the H/H2 transition occurs around 1.5 × 10−3 pc, in NGC 7023. For the Orion Bar, these values are 5.7 × 10−3 pc and 6.5 × 10−3 pc, respectively. We note that this increase in H2 abundance impacts the penetration of UV photons, which leads to a significant decrease of the heating rate and therefore of the temperature in the PDR. The size of the hot to warm molecular region is ~1.5 × 10−3 pc in Orion Bar, which is a factor of two lower than in NGC 7023. Taking into account the distances to the studied objects this would correspond to a thickness of ~2″ for NGC 7023 and ~0.8″ for the Orion Bar, which agrees well with the observations of the filaments in NGC 7023 and the dilution procedure we adopted to correct the observed intensities.

thumbnail Fig. 8.

NGC 7023 PDR model (Pth = 108 K cm−3, global scaling factor = 0.7). Top panel: evolution of the H nuclei number density and gas temperature with AV or distance (in 10−3 pc). Centre panel: spatial profile of the local emissivities of the main tracers. The emissivities have been scaled so that their maximum in the cloud is 1. Bottom panel: spatial profiles of the abundances of the species of interest in the model.

Open with DEXTER

thumbnail Fig. 9.

Orion Bar PDR model (Pth = 2.8 × 108 K cm−3, global scaling factor = 1.3). Top panel: evolution of the H nuclei number density and gas temperature with AV or distance (in 10−3 pc). Centre panel: spatial profile of the local emissivities of the main tracers. The emissivities have been scaled so that their maximum in the cloud is 1. Bottom panel: spatial profiles of the abundances of the species of interest in the model.

Open with DEXTER

Emission in the excited CO lines peaks at the back end of the warm molecular region, close to the C+/C/CO transition, with higher-J lines peaking closer to the surface. The emission profiles of the highest lines (e.g. J = 19–18 or 21–20) exhibit a wide left tail extending towards the H/H2 transition and accounting for roughly half of the total emission. This fraction of the emission comes from a low fraction (fractional abundance ~10−6 of CO existing before the C+/C/CO transition. In this region, we find CO formation to be initiated by CH+ formation as explained in Sect. 6.1. Emission in the H2 rotational lines spans the whole warm molecular region, with S(5) peaking close to the H/H2 transition and S(0) peaking closer to the high-J CO emission peaks. Vibrational H2 emission (e.g (1–0) S(1) line) peaks at the very edge of the hot molecular region, where the H2 abundance starts to increase steeply due to self-shielding processes. CH+ emission coincides with the higher H2 rotational lines (such as S(5)) and the vibrational line (1–0) S(1). Finally, the fine structure lines [CII] (158 μm ) and [OI] (63 and 145 μm ) arise in the whole hot and warm molecular region. The [CII] 158 μm and [OI] 63 μm lines also show significant emission in the atomic region in our model.

6. Discussion

6.1. Formation and excitation of CO at the PDR bright edges

The main objective of this paper is to analyse the emission in high-J CO lines in PDRs. This emission, as observed in the two template PDRs, NGC 7023 and Orion Bar, points either to an enhanced formation of CO in the PDR warm external layers or to an increased temperature in the layers containing CO. As will be discussed in more details below, this question is intimately related to another important question, which is the formation rate of H2 in PDRs. Several authors have previously studied H2 rotational emission in PDRs and found that a significant fraction of H2 has an excitation temperature in the range 400–700 K, much higher than expected from models (Parmar et al. 1991; Draine & Bertoldi 1999; Allers et al. 2005; Habart et al. 2011). Those authors suggested that this is due either to an increase of the photoelectric heating efficiency or a larger H2 formation rate on grains at high temperatures or to dynamical effects such as advection. The effect of an advancing photodissociation front was also invoked in NGC 7023 NW to account for the non-equilibrium values of the ortho-to-para ratio (Lemaire et al. 1999; Fuente et al. 1999; Fleming et al. 2010; Le et al. 2017). Since our models of NGC 7023 and Orion Bar successfully account for both CO high-J and H2 line intensities (as well as the ortho-to-para ratio), it is worthwhile investigating the underlying reason.

In our models, we used the prescription of Le Bourlot et al. (2012) that considers the formation on grains via the LH and the ER mechanisms. The resulting formation rate at the edge of the NGC 7023 and Orion Bar models are found to be three to four times higher than the rate determined for the diffuse gas of ~3 × 10−17 cm3 s−1 (Jura 1974; Gry et al. 2002). The ER mechanism is found to dominate at the edge of the PDR (Röllig et al. 2013) because the grains are warm. We note that Bron et al. (2014) have shown that the LH mechanism can also be efficient because of the temperature fluctuations of small grains due to photon absorption, although this mechanism has not been taken into account in our models. Thanks to the high H2 formation rate, the atomic to molecular transition is shifted towards the edge of the cloud, in other words into a hotter gas where H2 can be more excited. In our models, we also used a fully efficient conversion between ortho- and para-H2 on dust grains, meaning that all H2 molecules adsorbing on grains are converted before desorbing, which happens to be a good prescription to account for the observed relative H2 line intensities and therefore the observed ortho-to-para ratio. This efficient conversion on grains is justified by the temperature fluctuations of the grains as shown in Bron et al. (2016).

The presence of H2 in the hot to warm gas impacts the chemistry including the formation of CO since several routes initiated by the formation of CH+ become efficient. The C+ + H2 reaction has an activation threshold of ~4500 K that can be overcome, first by the kinetic energy of the reactants in hot gas, and second by the internal energy of FUV-pumped H2. In our models, we compute this reaction rate with the prescription by Agúndez et al. (2010) that takes into account H2 excitation. H2 ro-vibrational level populations of the ground electronic state are computed explicitly taking into account collisional excitation and de-excitation, radiative emission and pumping, UV pumping in electronically excited states followed by radiative de-excitation in the ground electronic state, excitation at formation following Sizun et al. (2010) and chemical destruction. This detailed treatment of H2 excitation was already found to provide a good agreement between predicted CH+ and SH+ line intensities and observations in the Orion Bar (Agúndez et al. 2010; Nagy et al. 2013; Goicoechea et al. 2017). Once CH+ is present in the gas, efficient ion-neutral reactions take place to produce species such as the CH+, , chain. Then reacts with O to give HCO+ followed by electronic recombination that gives CO. So the formation of CO is found to occur via this hot chemistry at the edge of PDRs as soon as H2 starts to self-shield (see also Goicoechea et al. 2016).

As can be seen in Figs. 8 and 9, CO emits in high rotational transitions as soon as H2 is formed efficiently in the PDR. In the frame of our stationary PDR models, ~50% of the J = 19–18 line intensity originates behind the H/H2 transition and before the C+/C/CO one. The remaining 50% are produced at the C+/C/CO transition. In our models, the gas temperature between the H/H2 and C+/C/CO transitions is a few hundred Kelvin, enough to excite collisionally CO in high rotational states. The main heating mechanism of the gas is the photo-electric effect on grains. Allers et al. (2005) modified both the UV dust extinction and the photoelectric efficiency in order to increase the gas temperature and therefore account for the H2 emission in the Orion Bar. This is in line with the fact that the abundance of polycyclic aromatic hydrocarbons (PAH) is found to increase at the border of PDRs (see Pilleri et al. 2012, and references therein). Since PAHs are major contributors to the photoelectric heating, we implemented in the code several ad-hoc profiles for the photoelectric heating that describe at best the observed abundance variation. However we can not conclude that these modified models were necessary to account for the observations since the new implementation of H2 formation was already leading to a sufficient warm molecular layer. In fact the analysis of mid-IR observations shows that the increased abundance of PAHs is found at AV ≤ 1 whereas our model predicts that excited CO emission is found at a different depth of AV ~ 2–3 (Figs. 8 and 9). We therefore did not include a profile in the photo-electric heating and used our standard implementation based on Bakes & Tielens (1994) in the final analysis presented here.

Finally, approximations in the computation of the self and mutual shielding of pre-dissociating UV lines can affect the calculated CO line intensities. In our models, we used the Federman et al. (1979) prescription to estimate the self-shielding. This formalism considers only the self-shielding of lines by themselves. UV absorption lines of CO and its isotopologues are only shifted by hundredths of Angstrom. Shielding of CO lines by H2 lines, and mutual line shielding of CO and its isotopologues can reduce the photo-destruction rates of CO and its isotopologues. The Meudon PDR code allows to compute exactly such shieldings but at the price of a very significant computing time (Goicoechea & Le Bourlot 2007). We performed some test calculations since a complete grid exploration was out of reach. We find that in the NGC 7023 and Orion Bar parameter range, this effect contributes to a maximum of a factor of three on high-J 13CO lines, such as the 15–14 transition, and is not sufficient to explain the mismatch between the modelled and observed 13CO line intensities. This supports the notion that most of this mismatch arises from the mixing of components in the beam, as discussed in Sect. 4.2.

6.2. Structure of the brightest PDR interfaces

As presented in Sect. 5.3 our models manage to account for most of the observed lines that trace the bright PDR interface within a factor of a few (even better in the case of NGC 7023). We can then conclude that, to the first order, it is possible to explain the line intensities emitted at the edge of bright PDRs (G0 ~ 103–104) with a single interface modelled by a stationary isobaric PDR model at high thermal pressure (~108 K cm−3). A PDR model at Pth = 108 K cm−3 has been used for the Orion Bar in studies related to the one presented here but only for a limited data set and without exploring the role of the parameters (Nagy et al. 2013; Cuadrado et al. 2015). This idea of using such isobaric models has been raised previously by several authors (Marconi et al. 1998; Lemaire et al. 1999; Allers et al. 2005) but, most of the times, the alternative clump-interclump scenario has been used or invoked (Stutzki & Guesten 1990; Parmar et al. 1991; Meixner & Tielens 1993; Tauber et al. 1994; Hogerheijde et al. 1995; van der Werf et al. 1996; Usuda et al. 1996; Andree-Labsch et al. 2017). In this scenario, clumps with density of ~106–107 cm−3 that include about 10% of the gas, are embedded in a more diffuse gas ~104–105 cm−3, meaning that the difference of density between the two components is a factor ~100. In our isobaric model at high pressure (~few 108 K cm−3) a comparable gradient in density naturally arises from the hypothesis of constant pressure. It is clear that the assumption of constant pressure equilibrium for our two PDRs is an approximation but the strength of our isobaric model is, however, that it can describe the structure by a single parameter while the clumpy models have to select or fit a very large and partially arbitrary parameter space. Some lines are found to be more responsive to the pressure than others as shown in Fig. 10. This is the case of the 0–0 S(5) H2 line and the high-J CO lines. On the opposite, the [CII] line is rather constant at 5 × 106 < Pth <108 K cm−3 and its intensity then decreases significantly at higher pressures.

thumbnail Fig. 10.

Variation of the intensity of different lines of interest with the gas thermal pressure for the illumination conditions of NGC 7023 NW, as computed by the PDR model.

Open with DEXTER

Nevertheless, a stationary model is not expected to capture all the complexity of a strongly irradiated PDR, such as the Orion Bar, where out-of-equilibrium effects take place (Bertoldi & Draine 1996). However there is no point in refining our model considering the poor spatial resolution in most tracers. The pressure we derived at the Orion Bar edge, is likely an effective pressure providing the best compromise between fitting the emission of surface tracers (excited H2 lines and CH+ lines) and that of CO (mid- to high-J lines mainly emitted at AV ~ 2–4). The computed value of the H2 1–0 S(1) to 2–1 S(1) ratio, which is significantly too high compared to the observations (a factor of approximately 10) suggests that the gas density is too high at the surface of the PDR, considering that the 1–0 S(1) line intensity increases with collisional relaxation of higher UV-pumped vibrational levels. The existence of a pressure gradient is in line with the work of Goicoechea et al. (2016) who showed that the gas kinematics derived from the ALMA maps suggests gas flowing from the high-pressure molecular layers (Pth ~ 2 × 108 K cm−3) to the atomic layers (Pth ~ 5 × 107 K cm−3). The advection of molecular gas through the PDR edge would impact both the chemical and thermal structures and therefore the calculated line intensities (Bertoldi & Draine 1996). In this case one can expect a lower value for the H2 1–0 S(1) to 2–1 S(1) ratio since more emission will arise from UV-pumped levels. Störzer & Hollenbach (1998) calculated that on average the H2 lines are affected by a factor of three.

6.3. A Pth − G0 relation in PDRs

Our work shows that the emission of the warm molecular gas at the PDR edge can be attributed to a slab at a high thermal pressure of 108 and 3 × 108K cm−3 in NGC 7023 NW and Orion Bar, respectively. There is a trend that this pressure increases with the G0 value. To further explore this relation, we compiled data from the literature. While doing so, one needs to be sure that the collected values are consistent with our study, which means that this pressure was derived from relevant tracers. In particular we avoided studies that were only based on the analysis of the rotational lines of H2.

The data presented in Fig. 11 was obtained from the following studies. Habart et al. (2005) studied in detail the Horsehead PDR and have shown that there is a density gradient at the interface with the HII region, which can be modelled by a thermal pressure equilibrium at Pth ~ 4 × 106 K cm−3. The value of Pth = 5 × 106 K cm−3 for NGC 7023 E corresponds to the maximum value derived by Köhler et al. (2014) from their combined analysis of dust and excited CO lines. Pérez-Beaupuits et al. (2010) included some mid-J CO lines in their analysis of the warm gas in M17 SW and concluded that the high-density gas (nH = 5 × 105cm−3) has a maximum temperature of 230 K. We therefore used these values to derive a pressure of Pth = 1.2 × 108 K cm−3. Finally, we included the massive star forming region W49A in which Nagy et al. (2012) studied the warm molecular gas and derived a kinetic temperature map using H2, CO excited lines and a volume density map using HCN excited lines. The maximum thermal pressure can be derived at the centre of the map with a value of Pth = 5.4 × 108 K cm−3.

thumbnail Fig. 11.

Relation between the thermal pressure in the dense structures of PDRs and the UV intensity G 0, see text for references. The dashed lines show the range of values obtained by Bron et al. (2018) in their photoevaporating PDR models.

Open with DEXTER

Figure 11 shows that Pth indeed increases with G0. Considering the very uncertain error bars, we do not provide here a fitting of the reported points that would provide a more quantitative scaling of Pth with G0. Nevertheless, a visual inspection of Fig. 11 leads to Pth/G0 = 1–4 × 104 K cm−3 except for W49A that falls below this range. We note that this graph can help to rationalise the results presented in Stock et al. (2015). In the two PDRs studied, S 106 and IRAS 23133+6050, the CO SLEDs are close to those measured in Orion Bar, which could be explained by a UV radiation field, G0, of a few 104. We can also comment on the results obtained by Indriolo et al. (2017) on the CO SLEDs in prototypical massive star-forming regions in which both a high-UV field and shocks are thought to excite the gas. Our study suggests that excitation by UV photons also plays a major role in the case of Orion S and W49N. Indeed both objects have similar CO SLEDs and W49N/A is found to follow to some extent the PthG0 trend shown in Fig 11. In these objects, shocks are however also involved; they are likely to be the major driver for emission in CO lines with Jup > 25 and are revealed when line profiles can be resolved. For instance, Tahani et al. (2016) showed that for the CO J = 16–15 line towards Orion S, there is a narrow (4 km s−1) component associated with the PDR and a broad (15 km s−1) component associated with shock excitation, both having similar integrated intensities. More detailed modelling would be necessary to disentangle the contribution of both excitations on the CO SLEDs.

The obtained PthG0 relation can also give us further insights into the (unclear) origin of the density structures that are found at the edge of H II regions (case of Orion Bar) or of atomic regions (case of NGC 7023 in which no H II region is present). It suggests that the UV radiation field plays a major role in the compression of the PDR. As the pressures found in the PDRs are significantly higher than the pressures found in the H II regions (e.g. Pth = 6 × 107 K cm−3 in the Orion Bar, Goicoechea et al. 2016), pressurisation by the thermal pressure of the H II region (e.g. in an expanding H II region) is not sufficient to explain the trend. Photoevaporation, in which photoheated gas at the ionisation and dissociation fronts expands into the central cavity and exerts by reaction a force on the neutral and/or molecular part of the cloud (Bertoldi 1989; Bertoldi & Draine 1996), could induce compression of the molecular part of the PDR and explain the pressure difference with the central ionised or atomic cavity. In addition, the tight correlation with G0 independently of the presence of an ionisation front (case of NGC 7023) close to the PDR seems to indicate that non-ionising (FUV) photons can be at least as efficient as ionising photons for this photoevaporation process. These considerations have found theoretical support in a recent study by Bron et al. (2018). The authors find that photoevaporation of the illuminated edge of the molecular cloud can indeed lead to high pressures and account for the Pth-G0 trend. To illustrate this result, we show in Fig. 11 the range of values (dashed lines representing Pth/G0 = 5 × 103 and 8 × 104 K cm−3, respectively) obtained by the authors using their time-dependent hydrodynamical PDR code. The agreement with the observations is striking and opens new perspectives to study dynamical evolution of the strongly illuminated edges of molecular clouds in massive star-forming regions.

7. Conclusion

Thanks to Herschel, we have measured the CO SLEDs in two prototypical PDRs: NGC 7023 NW (observed 12CO lines from Jup = 4–19) and the Orion Bar (observed 12CO lines from Jup = 4–23). The excitation temperature deduced for Jup≥15 from the rotational diagrams are 112 and 147 K, respectively, showing the presence of warm CO gas at the irradiated PDR edge. We have used the Meudon PDR code and more specifically stationary isobaric PDR models to account for high-J 12CO lines as well as for H2 and CH+ lines. The thermal pressure value and a global scaling factor alone were used as free parameters. The best models were obtained for a gas thermal pressure of Pth ~ 108 K cm−3 and provide a good agreement with the observed values. The prediction made by these models for other lines from HCO+, O, C, C+, HD, and OH has also been found to be compatible with the observed values.

Compared to previous works, we found that by simulating in detail the H2 formation process on grains, its level excitation, and considering state-to-state chemistry for key reactions, it is possible to explain line emission of molecules at the edge of PDRs with a stationary model without the introduction of adhoc hypothesis as clumps or shocks. One of the key mechanisms to account for warm CO in PDRs is the high efficiency of H2 formation at the PDR edge, which brings the H/H2 transition closer to the interface. The gas temperature ranges from 1000 K where H2 starts to self-shield to 100 K at the C+/C/CO transition. CO starts to form as soon as H2 appears in the PDR. This is due to the fact that the warm temperature, high density and presence of FUV-pumped H2 allow the formation of CH+ via the H2 + C+ reaction, opening a hot chemistry channel that leads to the formation of CO. As a consequence, in the framework of these stationary models, a significant fraction (~50% for the models of the two PDRs presented here) of high-J CO emission is produced before the C+/C/CO transition. In any case, the separation between the H/H2 and C+/C/CO transitions is predicted to be very small (less than one arcsec. at the distance of Orion) and this is in line with ALMA observations (Goicoechea et al. 2016).

Our results impact our view of the irradiated edge of star-forming regions and the feedback of star formation on its parental cloud. In the two prototypical PDRs, NGC 7023 NW and Orion Bar, we found that the FUV photons from nearby massive stars have enough energy to explain CO excitation in mid- and high-J levels. No additional energy source as mechanical heating is required. A comparison of NGC 7023 NW and Orion Bar with other typical PDRs shows a correlation between the thermal pressure at the edge of PDRs and the intensity of the UV radiation field. A similar correlation was recently reported by Wu et al. (2018) in their spatial study of the Carina nebula. This seems to indicate that the UV radiation field is the main driver of the high pressure at the edge of PDRs. This high-pressure edge is very sharp of the order of a few 10−3 pc (one arcsec at the distance of Orion) and is expected to evolve in the surrounding gas at lower pressure. The presence of pressure gradients and advection processes is invoked but could not be further exploited due to the limited spatial resolution of the observations.

Implications for extragalactic studies have been discussed by Indriolo et al. (2017) and are supported by our study. As the high-pressure edges are thin structures, their observation suffers from beam dilution effects. The availability of high-spatial resolution observations is therefore crucial to refine PDR models and is only feasible while studying Galactic PDRs. Perspectives in this topic include the coming James Webb Space Telescope that will soon give us access to subarcsecond-resolution observations for a large number of H2 rotational and ro-vibrational lines. Observations of high-J CO lines with high spatial resolution are however currently out of reach and this could be the case for quite some time. The use of alternative species, which could be observed by ALMA, such as HCO+ or HCN, should be explored. From the modelling side, further improvements should include a better description of the grain populations and their properties since the penetration of the UV field plays a central role in the energetics and chemistry of PDRs. An additional big step would consist of coupling the chemistry and energetics with dynamical processes that can describe the evolution of these high-pressure edges (Gorti & Hollenbach 2002).

Acknowledgments

This work was in part supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP, and co-funded by CEA and CNES. It was supported by the German Deutsche Forschungsgemeinschaft, DFG project number SFB 956, C1. J.R.G and J.C. thank Spanish MINECO for funding supports under grants AYA2012-32032 and AYA2017-85111-P. E. B. thanks the ERC grant ERC-2013-Syg-610256-NANOCOSMOS for support.


2

We have not used the most sophisticated treatment at our disposal, Bron et al. (2014), for computing-time reasons.

References

All Tables

Table 1.

Observed data for NGC 7023, and dilution factor Ω.

Table 2.

Observed data for the Orion Bar and dilution factor Ω.

Table 3.

Input parameters used in the Meudon PDR code.

Table 4.

Comparison of models and observations (combined) for NGC 7023.

Table 5.

Comparison of models and observations (combined) for the Orion Bar.

All Figures

thumbnail Fig. 1.

Composite images of the Orion Bar (top panel) and NGC 7023 (bottom panel). Red indicates the 8 μm emission observed with Spitzer. Green shows the vibrationally excited emission of H2 from Lemaire et al. (1996) and van der Werf et al. (1996) for NGC 7023 and Orion Bar, respectively. Blue shows the 12CO emission, J = 6 − 5 for the Orion Bar (Lis & Schilke 2003) and J = 1−0 for NGC 7023 (Gerin et al. 1998). The circles represent the HPBW of Herschel at 550 GHz and 1900 GHz. The square indicates the position of the central spaxel of the PACS observations.

Open with DEXTER
In the text
thumbnail Fig. 2.

12CO and 13CO lines observed with HIFI towards NGC 7023 (left) and the Orion Bar (right). The vertical red line indicates the systemic velocity of the source.

Open with DEXTER
In the text
thumbnail Fig. 3.

Observed intensities of 12CO (left) and 13CO (right) in the Orion Bar (top panel) and NGC 7023 NW (bottom panel).

Open with DEXTER
In the text
thumbnail Fig. 4.

Rotational diagram of 12CO (left) and 13CO (right) lines observed in the Orion Bar (top panel) and NGC 7023 (bottom panel) PDRs. For 12CO, the excitation temperature and total column density are computed in the range J up = [15, 23] (Orion Bar), and J up = [15, 19] (NGC 7023). For 13CO, the lines used were J up = [5, 16] (Orion Bar), and J up = [5, 10] (NGC 7023).

Open with DEXTER
In the text
thumbnail Fig. 5.

Line profile for the J = 10–9 transition of 12CO (dashed line) and 13CO (dotted line) in the Orion Bar (top panel) and NGC 7023 NW (bottom panel). The relative intensity of the 12CO to 13CO lines across the velocity profiles is also shown (plain line).

Open with DEXTER
In the text
thumbnail Fig. 6.

Excitation diagram of the different tracers observed in NGC 7023 after dilution correction (red squares) and the best fit model (black, Pth = 108 K cm−3, = 10, global scaling factor = 0.7). The best model has been chosen to optimise the fitting of the 12CO (high-J lines), H2 and CH+ lines only. In the last panel, the intensity values are normalised by the mean observed value for each line. The grey lines show the obtained variability when the thermal pressure is divided (dashed lines) or multiplied (plain lines) by a factor of 1.5. The best value for the scaling factor was found to be 1.2 for the model at Pth/1.5 and 0.47 for the model at Pth × 1.5.

Open with DEXTER
In the text
thumbnail Fig. 7.

Excitation diagram of the different tracers observed in the Orion Bar (red squares) and the best model (black, P th = 2.8 × 108 K cm−3, = 10 and global scaling factor f = 1.3). The best model has been chosen to optimise the fitting of the 12CO (high-J lines), rotational H2 and CH+ lines only. In the last panel, the intensity values are normalised by the mean observed value for each line. The grey lines show the obtained variability when the thermal pressure is divided (dashed lines) or multiplied (plain lines) by a factor of 1.5. The best value for the scaling factor was found to be 2.0 for the model at Pth/1.5 and 0.9 for the model at Pth × 1.5.

Open with DEXTER
In the text
thumbnail Fig. 8.

NGC 7023 PDR model (Pth = 108 K cm−3, global scaling factor = 0.7). Top panel: evolution of the H nuclei number density and gas temperature with AV or distance (in 10−3 pc). Centre panel: spatial profile of the local emissivities of the main tracers. The emissivities have been scaled so that their maximum in the cloud is 1. Bottom panel: spatial profiles of the abundances of the species of interest in the model.

Open with DEXTER
In the text
thumbnail Fig. 9.

Orion Bar PDR model (Pth = 2.8 × 108 K cm−3, global scaling factor = 1.3). Top panel: evolution of the H nuclei number density and gas temperature with AV or distance (in 10−3 pc). Centre panel: spatial profile of the local emissivities of the main tracers. The emissivities have been scaled so that their maximum in the cloud is 1. Bottom panel: spatial profiles of the abundances of the species of interest in the model.

Open with DEXTER
In the text
thumbnail Fig. 10.

Variation of the intensity of different lines of interest with the gas thermal pressure for the illumination conditions of NGC 7023 NW, as computed by the PDR model.

Open with DEXTER
In the text
thumbnail Fig. 11.

Relation between the thermal pressure in the dense structures of PDRs and the UV intensity G 0, see text for references. The dashed lines show the range of values obtained by Bron et al. (2018) in their photoevaporating PDR models.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.