Structure of photodissociation fronts in star-forming regions revealed by observations of high-J CO emission lines with Herschel

In bright photodissociation regions (PDRs) associated to massive star formation, the presence of dense"clumps"that are immersed in a less dense interclump medium is often proposed to explain the difficulty of models to account for the observed gas emission in high-excitation lines. We aim at presenting a comprehensive view of the modeling of the CO rotational ladder in PDRs, including the high-J lines that trace warm molecular gas at PDR interfaces. We observed the 12CO and 13CO ladders in two prototypical PDRs, the Orion Bar and NGC 7023 NW using the instruments onboard Herschel. We also considered line emission from key species in the gas cooling of PDRs (C+, O, H2) and other tracers of PDR edges such as OH and CH+. All the intensities are collected from Herschel observations, the literature and the Spitzer archive and are analyzed using the Meudon PDR code. A grid of models was run to explore the parameter space of only two parameters: thermal gas pressure and a global scaling factor that corrects for approximations in the assumed geometry. We conclude that the emission in the high-J CO lines, which were observed up to Jup=23 in the Orion Bar (Jup=19 in NGC7023), can only originate from small structures of typical thickness of a few 1e-3 pc and at high thermal pressures (Pth~1e8 K cm-3). Compiling data from the literature, we found that the gas thermal pressure increases with the intensity of the UV radiation field given by G0, following a trend in line with recent simulations of the photoevaporation of illuminated edges of molecular clouds. This relation can help rationalising the analysis of high-J CO emission in massive star formation and provides an observational constraint for models that study stellar feedback on molecular clouds.


Introduction
Photodissociation regions (PDRs) are key regions in the study of the interstellar medium. They 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 Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. 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, C, as well as H 2 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 pho-Article number, page 1 of 20 arXiv:1801.03893v2 [astro-ph.GA] 15 Apr 2018 tons from nearby massive stars involving the photo-electric effect on grains (Bakes & Tielens 1994;Weingartner & Draine 2001) as well as the collisional de-excitation of H 2 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 that simulate radiative transfer, chemistry, 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 as the Meudon PDR code, can simulate very detailed micro-physical processes 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. 2014Bron et al. , 2016. The analysis of line emission in PDRs is intimately connected to considerations on the morphology of these regions. Earlier observations of a few mid-/high-J CO lines from groundbased 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 (n H ∼ 10 6 -10 7 cm −3 ) are embedded in a less dense interclump medium (n H ∼ 10 3 -10 5 cm −3 ). In following studies on the bright PDRs M17 SW, Orion Bar and Carina, clumpy PDR models were used to analyse quite successfully the observations (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 J up = 4 to J up = 50. This allows us to build full CO spectral line energy distributions (SLEDs) 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 on the Galactic Center with observations towards Sgr A* , as well as in the Sgr B2 cores, B2(M) and B2(N) . In external galaxies, CO SLEDs have been obtained for Seyfert galaxies, starburst galaxies, (ultra)luminous infrared galaxies ((U)LIRGs) and AGNs (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. 2012Kazandjian et al. , 2015. In PDRs, UV photons are expected to be the major actor and these CO lines offer the possibility to further constrain 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 (n H ∼ 10 6 cm −3 ) immersed in a strong far-UV radiation field (G 0 ∼ 10 5 in units of the Habing field; Habing 1968) and an interclump medium at lower density and irradiation field (n H ∼ 10 4 cm −3 , G 0 ∼ 10 4 ). However, the high value derived for the UV field on the clump (a factor of 10 higher compared to the interclump medium) is striking.
In this work, we try to bring some new insights into 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 . NGC 7023 NW is well known to exhibit bright H 2 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 take benefit of the Herschel HIFI, SPIRE and PACS data. In order to include further tracers of the warm/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-/high-J CO, H 2 , CH + , [CII], [OI], [CI], HD, OH and HCO + lines. We analyse these observations using the latest version of the Meudon PDR code (Le Petit et al. 2006). Our goal is to get 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 modeling 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-/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.  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 G 0 = 2600 and a density of n H ∼ 4 × 10 3 cm −3 . Later observations have shown that the star formation process has shaped a cavity inside the molecular cloud with walls consisting of dense PDRs, at the north-west (NW), at the south and at the east (Fuente et al. 1992;Rogers et al. 1995;Fuente et al. 1998;Gerin et al. 1998 Lemaire et al. (1996) and van der Werf et al. (1996) for NGC 7023 and Orion Bar, respectively. Blue shows the 12 CO 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.

The prototypical PDRs: NGC 7023 and Orion Bar
bright thin filaments were revealed by high angular resolution images in the Extended Red Emission (ERE), vibrationally excited H 2 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) (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, n H ∼ 10 5 − 10 6 cm −3 (Martini et al. 1997;Lemaire et al. 1996;Fuente et al. 1996a), which are embedded in a more diffuse gas with n H ∼ 10 3 − 10 4 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. Whereas 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 turned out to be one of the best sites to study the physical and chemical processes taking place in a PDR.

The 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 H ii region that penetrates into the parent molecular cloud. The UV intensity impinging on the PDR has been estimated to be G 0 = 1 − 4 × 10 4 in Habing units (Tielens & Hollenbach 1985;Marconi et al. 1998). The distance of the Orion nebula has been measured with great precision by Menten et al. (2007) from trigonometric parallax, yielding a value of 414 ± 7 pc. An angular distance of 1 arsecond 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 that changes from face-on to almost edge-on and then face-on. 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 × 10 4 cm −3 in order to account for the observed offsets of ionization front and molecular lines (see also Wyrowski et al. 1997). Most models of the molecular emission in the Orion Bar used high-density clumps (n H ∼ 10 6 − 10 7 cm −3 ) embedded in a lower density gas (n H = 5×10 4 −10 5 cm −3 ).
Clumps were introduced first to explain the emission of excited lines of CO and warm H 2 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 H ii 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 H ii 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 ionization front and in front of the molecular interface traced by vibrational H 2 . This is also well illustrated in Fig. 1 that combines the IRAC 8 µm map dominated by PAH emission with the H 2 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 (n H ∼ 10 6 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 H 2 vibrational emission. The H/H 2 and C + /C/CO transitions should therefore be very close at this interface, which is incompatible with stationary PDR models at n H ∼ 5 × 10 4 cm −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 centered on the bright H 2 filaments, a position referred to as the H 2 peak (Joblin et al. 2010(Joblin et al. , α 2000 =21 h 01 m 32.4 s , δ 2000 =+68 • 10 25.0 ). The Orion Bar was observed at the CO + peak position (Stoerzer et al. 1995(Stoerzer et al. , α 2000 (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 programs gathered spectroscopic data using HIFI and PACS instruments. We also use in this article 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 litterature and archives were also collected in particular for H 2 . The full data sets, which are used in this study, are reported in Tables 1 and 2 for NGC 7023 and the Orion Bar, respectively.

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 arcsec 2 (cf. Sec.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. 2012Okada et al. 2013).

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 in 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 T A is more appropriate for extended emission, and T mb 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 T A and T mb because the bright interface is in between point-source and extended. In addition, it is not practical 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 12 CO and 13 CO lines observed with Herschel-HIFI for both sources is shown in Fig. 2.

Herschel SPIRE observations
We have used the SPIRE FTS fully sampled maps to extract complementary data at the CO + peak. 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).

Additional data
We used additional observations to further constrain our models. For NGC 7023, the pure rotational lines of H 2 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 Article  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 toward the H 2 peak position over 4 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 H 2 , observed at the CFHT and the Perkins Telescope, are taken from Lemaire et al. (1996), Lemaire et al. (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 H 2 lines from ISO-SWS data (Bertoldi, private communication;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 H 2 , 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).

Combination of observational data
For the purpose of comparing the observations with PDR models (cf. Sect. 5), we need 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 H 2 emission (Lemaire et al. 1996). From these maps one can derive a relevant width of 2 for the bright PDR interface that is centered toward the H 2 peak. In Orion Bar, the sharp molecular edge has been observed in vibrational and rotational H 2 transitions 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. 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 H 2 . 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 in both PDRs 576.268 GHz 2.0 ± 0.6 (−9) ---0.07 J=6-5 691.473 GHz 5.3 ± 1.6 (−9) 7.5 ± 0.7 (−9) --0.08 J=7-6 806.652 GHz Fuente et al. (1996) Lemaire et al. 1996-d Martini et al. 1999 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 have to apply scaling factors on some of the data sets. The 12 CO 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 consider 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 12 CO 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 H 2 observations in NGC 7023 it is not possible to satisfactorily merge the Spitzer and ISO data. This is due to the fact that H 2 emission is quite extended and the ISO beam contains emission from other regions 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.
Finally, after correcting for beam dilution and crosscalibration 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   T ex = 147 ± 9 K N ( 12 CO) = (9.0 ± 3.9) × 10 17 cm −2    Figure 4 presents the rotational diagrams and LTE fits of the 12 CO and 13 CO observations. The rotational diagrams present column density estimates in the upper rotational levels, N u , without corrections for opacity effects with N u /g u = 4 π I A ul hν ul g u where I is the observed intensity, A ul and ν ul are the Einstein coefficient and frequency of the transition from the upper to lower levels and g u is the degeneracy of the upper level. As such, the derived column denities are lower limits when the lines are optically thick. In the case of 12 CO, the fit was restricted to the highest excitation lines, starting at J up =15, in order to obtain information on the warmest gas. In the Orion Bar, we derived an excitation temperature of 147 ± 9 K and a column density of N( 12 CO)=(9.0 ± 3.9) × 10 17 cm −2 . In NGC 7023, the excitation temperature is 112 ± 6 K with N( 12 CO)=(1.7 ± 0.7) × 10 17 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 12 CO 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( 12 CO)=2 − 3 × 10 18 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 that was mentioned earlier. It is also compatible with the lower limit of the kinetic temperature derived by Goicoechea et al. (2016) at the dissociation front.

The CO rotational diagrams
The 13 CO diagrams reveal gas at an average rotational temperature of ∼80 K, which is cooler than the 12 CO temperature discussed above. In addition, Figure 5 reports the value of the 12 CO over 13 CO line intensity ratio across the velocity profile of the 12 CO 10-9 line. This ratio is found to be weak (∼ 4) near the line center, whereas it reaches the isotopic 12 C/ 13 C 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 12 CO lines in NGC 7023 using the RADEX code and concluded that the lines are optically thick up to J up ∼ 13−14. This justifies why we selected only lines from J up =15 and higher in our fit of the 12 CO rotational diagrams (cf. Fig. 4).
The above results strongly suggest that CO emission in the beam stems from a two-component medium: an extended cool/warm component and the hot and sharp interface. In 13 CO, the former component is well seen, whereas it is partly hidden due to optical depth effects in 12 CO emission. This results in 12 CO emission looking like stemming mainly from the hot component. On the opposite, the hot component is difficult to observe in the high-J lines of 13 CO 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 12 CO lines (see Tab. 2). By applying the same factor to the 13 CO lines, we can predict an intensity of 8 10 −10 W m −2 sr −1 for the J=20-19 13 CO line, which cannot be detected considering that it is a factor of 3 weaker relative to the error bar of the J=16-15 13 CO line. In conclusion, only high-J 12 CO lines can be used to characterise the physical conditions at the warm and bright PDR interface.

The Meudon PDR model
We compared our observations to PDR models using an updated version of the Meudon PDR code (Le Petit et al. 2006, https://ism.obspm.fr). 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.
Several atomic and molecular data have been updated. In particular, it is worthwhile to notice for the present study that this version of the PDR code implements the CO-H 2 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 H 2 formation on grains). Table 3 summarizes the elemental abundances used in this paper. The formation reaction of H 2 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 3 × 10 −17 √ T/100 cm 3 s −1 based on Copernicus and FUSE H 2 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 H 2 formation (both LH and ER) in a master equation formalism (Bron et al. 2014(Bron et al. , 2016. In this paper, we use the Le Bourlot et al. (2012) formalism 1 .
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 neighboring 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 front side of the PDR the one illuminated by the stars responsible for the PDR and back side the other one. 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 G 0 = 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 assume that the back side of the PDR is illuminated by an isotropic radiation field with G 0 = 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, G 0 = 2 × 10 4 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 ionization rate ζ of H 2 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 ionization rate by a factor 10 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 photoionization 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 H 2 and CO) line selfshielding is computed using the Federman et al. (1979) approximation. Absorption by grains is implemented using parametrized extinction curves (Fitzpatrick & Massa 1986). In PDRs, extinction curves are usually flatter in the far UV than the mean Galactic extinction curve. That is characteristic of larger than standard grain sizes. The adopted extinction curve is the one of HD 38087 in Fitzpatrick & Massa (1990) and R V = 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 N H /E(B-V) = 1.05×10 22 cm −2 mag −1 , an intermediate value between the standard value 5.8 × 10 21 and the value determined for ρ Oph, 15.4 × 10 21 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 H 2 vibrational de-excitation can heat or cool the gas depending on the local physical conditions. For the photo-electric effect, we use the Bakes & Tielens (1994) prescription. Cooling rates are 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 summarized in Tab. 3. The modeling strategy for the two objects is described in the following sections.

Fitting strategy
Since we focus our analysis on the physical conditions at the PDR edge, we investigate the models that best account for the emission of tracers specific to this region, such as high-J 12 CO and H 2 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 found 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 were 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 H 2 emission (cf. Fig. 1), which could be the result of References.
(1) For the NGC 7023 PDR model, we adopt a beamed stellar radiation field as described in the text. If converted in a Habing scaling factor, this corresponds to G 0 = 2600, (2)  an overdense surface layer seen roughly edge-on. We thus adopt 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 1 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 12 CO lines from J up =11, the rotational H 2 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, P th and f . As these tracers only constrain the warm molecular layer of the PDR, we fixed a total A V value of 10 (N H ∼ 2 × 10 22 cm −2 for a face-on geometry and a factor 2 higher assuming 60 • inclination) in our models. We then compare the results of the best model with other available Table 4. Comparison of models and observations (combined) for NGC 7023. Values in parentheses are powers of ten.

NGC 7023
The best fit was found for a model with P th = 10 8 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 12 CO lines (above J up = 10), with the H 2 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-H 2 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 H 2 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. Tab. 3). 10 20

Other tracers
Best model P th, best × 1.5 P th, best /1.5 Observations Fig. 6. Excitation diagram of the different tracers observed in NGC 7023 after dilution correction (red squares) and the best fit model (black, P th = 10 8 K cm −3 , A tot V = 10, global scaling factor = 0.7). The best model has been chosen to optimise the fitting of the 12 CO (high-J lines), H 2 and CH + lines only. In the last panel, the intensity values are normalized 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 P th /1.5 and 0.47 for the model at P th ×1.5.
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 12 CO lines below J up = 11 show an excellent agreement and the H 2 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 ∼ 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 2 in the density n(H 2 ) leads to an increase of the J=6-5 line by a factor 2 while the J=1-0 remains practically unchanged. This deviation suggests that the density provided by

Other tracers
Model P th, best × 1.5 P th, best / 1.5 Observations Fig. 7. Excitation diagram of the different tracers observed in the Orion Bar (red squares) and the best model (black, P th = 2.8 × 10 8 K cm −3 , A tot V = 10 and global scaling factor f = 1.3). The best model has been chosen to optimise the fitting of the 12 CO (high-J lines), rotational H 2 and CH + lines only. In the last panel, the intensity values are normalized 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 P th /1.5 and 0.9 for the model at P th ×1.5. Local emissivity (normalized) [CO] [CH + ] Fig. 8. NGC 7023 PDR model (P th = 10 8 K cm −3 , global scaling factor = 0.7). Top: Evolution of the H nuclei number density and gas temperature with A V or distance (in 10 −3 pc). Center: 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: Spatial profiles of the abundances of the species of interest in the model. and affected by self-absorption in most cases. Its study therefore requires a detailed analysis in velocity components. Finally, the model underestimates the 13 CO lines by a factor of ∼ 3. This is consistent with the arguments raised in Sec. 4.2, in particular that the 13 CO 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 12 CO lines.

Orion Bar
The best fit was found for P th = 2.8 × 10 8 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 12 CO lines, CH + lines, pure rotational lines of H 2 (with maximal differences by a factor of 2 except for the H 2 S(0) line that is underestimated by a factor of ∼ 4). We thus obtained a reasonable agreement, within a factor of a few, for all the tracers used in the fitting procedure.
Most of the other tracers (not used in the fitting procedure) show similar differences. The low lying 12 CO lines (J up < 11) are underestimated by up to a factor of ∼ 3. The OH, HD, [OI] and [CI] lines are reproduced within a factor of 3-4. Some tracers show however stronger discrepancies such as the 13 CO 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 13 CO is expected to arise from the surrounding molecular cloud (cf. Sec. 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 H 2 lines, the (1-0) S(1) line is strongly overestimated (factor of ∼ 6 − 7) while the (2-1) S(1) line is underestimated by a factor of ∼2 only.
To summarize, 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 + , H 2 lines except low H 2 rotational lines), whereas those arising from more internal layers (low-and mid-J CO lines, low H 2 rotational lines) are underpredicted. This suggests that a [CO] [CH + ] Fig. 9. Orion Bar PDR model (P th = 2.8 × 10 8 K cm −3 , global scaling factor = 1.3). Top: Evolution of the H nuclei number density and gas temperature with A V or distance (in 10 −3 pc). Middle: 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: Spatial profiles of the abundances of the species of interest in the model. single pressure is not sufficient to describe the evolution of gas density and temperature across the hot/warm irradiated interface (see also the grey curves in Fig. 7).

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 pannels), 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 A V =10. However, we provide here information for the region up to A V ∼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 summarized as follows: the atomic part has a density of a few 10 4 cm −3 (close to 10 5 cm −3 in Orion Bar) followed by a hot/warm (T = 1000 − 100 K) and relatively dense (few 10 5 to few 10 6 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/warm molecular region is found to start significantly before the H/H 2 transition. Indeed a sharp increase of the H 2 abun-dance is seen at 0.6 × 10 −3 pc whereas the H/H 2 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. Note that this increase in H 2 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/warm molecular region is ∼ 1.5 × 10 −3 pc in Orion Bar, which is a factor of 2 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.
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/H 2 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 Sec. 6.1. Emission in the H 2 rotational lines spans the whole warm molecular region, with S(5) peaking close to the H/H 2 transition and S(0) peaking closer to the high-J CO emission peaks. Vibrational H 2 emission (e.g (1-0) S(1) line) peaks at the very edge of the hot/warm molecular region, where the H 2 abundance starts to increase steeply due to self-shielding processes. CH + emission coincides with the higher H 2 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/warm molecular region. The [CII] 158 µm and [OI] 63 µm lines also show significant emission in the atomic region in our model.

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 H 2 in PDRs. Several authors have earlier studied H 2 rotational emission in PDRs and found that a significant fraction of H 2 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). The authors suggested that this is due either to an increase of the photoelectric heating efficiency or a larger H 2 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 are successful to account for both CO high-J and H 2 line intensities (as well as the ortho-to-para ratio), it is worth to investigate the underlying reason for that.
In our models, we use the prescription of Le Bourlot et al. (2012) that considers the formation on grains via the Langmuir-Hinshelwood (LH) and the Eley Rideal (ER) mechanisms. The resulting formation rate at the edge of the NGC 7023 and Orion Bar models are found to be 3-4 times higher than the rate determined for the diffuse gas of ∼ 3 × 10 −17 cm 3 s −1 (Jura 1974;Gry et al. 2002). The ER mechanism is found to dominate at the edge of the PDR ) 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 H 2 formation rate, the atomic/molecular transition is shifted towards the edge of the cloud, i.e. in a hotter gas where H 2 can be more excited. In our models, we also used a fully efficient conversion between ortho-and para-H 2 on dust grains, i.e. all H 2 molecules adsorbing on grains are converted before desorbing, which happens to be a good prescription to account for the observed relative H 2 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 H 2 in the hot/warm gas impacts the chemistry including the formation of CO since several routes initiated by the formation of CH + become efficient. The C + + H 2 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 H 2 . In our models, we compute this reaction rate with the prescription by Agúndez et al. (2010) that takes into account H 2 excitation. H 2 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 H 2 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 + , CH + 2 , CH + 3 chain. Then CH + 3 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 H 2 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 H 2 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/H 2 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 two transition layers (H/H 2 and C + /C/CO) is a few hundreds 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 H 2 emission in the Orion Bar. This is in line with the fact that the abundance of polycyclic aromatic hydrocarbons (PAHs) 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 could not conclude that these modified models were necessary to account for the observations since the new implementation of H 2 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 A V ≤1 whereas our model predicts that excited CO emission is found at a different depth of A V ∼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 use 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 H 2 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 found that in the NGC 7023 and Orion Bar parameter range, this effect contributes to a maximum of a factor 3 on high-J 13 CO lines, such as the 15-14 transition, and is not sufficient to explain the mismatch between the modelled and observed 13 CO line intensities. This supports that most of this mismatch arises from the mixing of components in the beam, as discussed in Sec. 4.2.

Structure of the brightest PDR interfaces
As presented in Sec. 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 (G 0 ∼ 10 3 − 10 4 ) with a single interface modelled by a stationary isobaric PDR model at high thermal pressure (∼ 10 8 K cm −3 ). A PDR model at P th =10 8 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;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 ∼ 10 6 − 10 7 cm −3 that include about 10% of the gas, are embedded in a more diffuse gas ∼ 10 4 − 10 5 cm −3 , i.e. the difference of density between the two components is a factor ∼ 100. In our isobaric model at high pressure (∼ few 10 8 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) H 2 line and the high-J CO lines. On the opposite, the [CII] line is rather constant at 5 10 6 < P th < 10 8 K cm −3 and its intensity then decreases significantly at higher pressures.
Still, 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 H 2 lines, CH + lines) and that of CO (mid-/high-J lines mainly emitted at A V ∼2-4). The computed value of the H 2 1-0 S(1)/2-1 S(1) ratio, which is significantly too high compared to the observations (a factor ∼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 (P th ∼ 2 × 10 8 K cm −3 ) to the atomic layers (P th ∼ 5 × 10 7 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 H 2 1-0S(1)/2-1S(1) ratio since more emission will arise from UV-pumped levels. Störzer & Hollenbach (1998) calculated that on average the H 2 lines are affected by a factor of 3.

A P th -G 0 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 10 8 and 3 × 10 8 K cm −3 in NGC7023 NW and Orion Bar, respectively. There is a trend that this pressure increases with the G 0 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 H 2 .
The data presented in Fig.11 is coming from the following studies. Habart et al. (2005) have studied in details the Horsehead PDR and have shown that there is a density gradient at the interface with the H ii region, which can be modelled by a thermal pressure equilibrium at P th ∼ 4 × 10 6 K cm −3 . The value of P th = 5 × 10 6 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 (n H = 5 × 10 5 cm −3 ) has a maximum temperature of 230 K. We therefore used these values to derive a pressure of P th = 1.2 × 10 8 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 H 2 , CO excited lines and a volume density map using HCN excited lines. The maximum thermal pressure can be derived at the center of the map with a value of P th = 5.4 × 10 8 K cm −3 . Figure 11 shows that P th indeed increases with G 0 . 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 P th with G 0 . Nevertheless, a visual inspection of Fig. 11 leads to P th /G 0 =1-4 10 4 K cm −3 except for W49A that falls below this range. It is interesting to note that this graph can help rationalizing 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, G 0 , of a few 10 4 . 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 P th -G 0 trend shown in Fig 11. In these objects, shocks are however also involved; they are likely the major driver for emission in CO lines with J up >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 P th -G 0 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 Hii regions (e.g. P th = 6 × 10 7 K cm −3 in the Orion Bar, Goicoechea et al. 2016), pressurization by the thermal pressure of the Hii region (e.g. in an expanding Hii region) is not sufficient to explain the trend. Photoevaporation, in which photoheated gas at the ionization 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 ionized/atomic cavity. In addition, the tight correlation with G 0 independently of the presence of an ionization front (case of NGC 7023) close to the PDR seems to indicate that non-ionizing (FUV) photons can be at least as efficient as ionizing photons for this photoevaporation process. These considerations have found theoretical support in a recent study by Bron et al. (2018). The authors found that photoevaporation of the illuminated edge of the molecular cloud can indeed lead to high pressures and account for the P th -G 0 trend. To illustrate this result, we show in Fig. 11 the range of values (dashed lines representing P th /G 0 = 5 × 10 3 and 8 × 10 4 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.

Conclusion
Thanks to Herschel, we have measured the CO SLEDs in two prototypical PDRs: NGC 7023 NW (observed 12 CO lines from J up =4 to 19) and the Orion Bar (observed 12 CO lines from J up =4 to 23). The excitation temperature deduced for J up ≥ 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 12 CO lines as well as for H 2 and CH + lines. Only the thermal pressure value and a global scaling factor were used as free parameters. The best models were obtained for a gas thermal pressure of P th ∼ 10 8 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 simulating in detail the H 2 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 ad-hoc hypothesis as clumps or shocks. One of the key mechanisms to account for warm CO in PDRs is the high efficiency of H 2 formation at the PDR edge, which brings the H/H 2 transition closer to the interface. The gas temperature ranges from 1000 K where H 2 starts to self-shield to 100 K at the C + /C/CO transition. CO starts to form as soon as H 2 appears in the PDR. This is due to the fact that the warm temperature, high density and presence of FUV-pumped H 2 allow the formation of CH + via the H 2 + C + reaction, opening a hot chemistry channel that leads to the formation of CO. As a consequence, in the frame of these stationary models, a significant fraction (∼50 % in 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/H 2 and C + /C/CO transition layers is predicted to be very small (less than one arsec. 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 starforming 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 midand 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. Implication for extragalactic studies has been discussed by Indriolo et al. (2017) and is supported by our study. The highpressure edges being 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 H 2 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, has to be explored. From the modelling side, further improvement 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 in coupling the chemistry and energetics with dynamical processes that can describe the evolution of these high-pressure edges (Gorti & Hollenbach 2002).