Merged H/H₂ and C⁺/C/CO transitions in the Orion Bar Free

ABSTRACT

High-resolution ALMA images towards the Orion Bar show no discernible offset between the peak of H₂ emission in the photodissociation region (PDR) and the ¹³CO(3–2) and HCO⁺(4–3) emission in the molecular region. This implies that positions of H₂ and CO dissociation fronts are indistinguishable in the limit of ALMA resolution. We use the chemo-dynamical model MARION to show that the ALMA view of the Orion Bar, namely, no appreciable offset between the ¹³CO(3–2) and HCO⁺(4–3) peaks, merged H₂ and CO dissociation fronts, and high intensity of HCO⁺(4–3) emission, can only be explained by the ongoing propagation of the dissociation fronts through the molecular cloud, coupled to the dust motion driven by the stellar radiation pressure, and are not reproduced in the model where the dissociation fronts are assumed to be stationary. Modelling line intensities, we demonstrate that after the fronts have merged, the angular separation of the ¹³CO(3–2) and HCO⁺(4–3) peaks is indeed unresolvable with the ALMA observations. Our model predictions are consistent with the results of the ALMA observations about the relation of the bright HCO⁺(4–3) emission to the compressed gas at the border of the PDR in the sense that the theoretical HCO⁺(4–3) peak does correspond to the gas density enhancement, which naturally appears in the dynamical simulation, and is located near the H₂ dissociation front at the illuminated side of the CO dissociation front.

1 INTRODUCTION

Young massive stars vigorously interact with surrounding molecular clouds. Expansion of ionized gas around a young massive star(s) is accompanied by a shock wave propagation into the molecular cloud, as has been shown both analytically (Spitzer 1978) and with numerical simulations (e.g. Hosokawa & Inutsuka 2005). But observational evidence of this process in real H ii regions is scarce and ambiguous, (Zavagno et al. 2007; Kirsanova et al. 2008; Gordon & Sorochenko 2009; Anderson et al. 2015), as it is not easy to prove that a compression wave, accompanied by ionization and dissociation fronts, does move into the molecular cloud.

The Orion Bar is an ionized border of a molecular cloud in the South-Eastern part of the Orion Nebula (e.g. O’Dell 2001), highly irradiated by far-ultraviolet (FUV, 6 eV < hν < 13.6 eV) photons. Regions of this kind are usually referred to as photodissociation regions (PDRs). The primary photoionization source of the Orion Nebula is a massive star θ¹C Ori of the spectral type O7 (Sota et al. 2011) or O6.5V (Tsivilev, Parfenov & Sobolev 2014). Along with the primary star, other Trapezium stars provide significant contribution to the FUV irradiation of the Orion Bar PDR (Ferland et al. 2012). The Orion Bar is seen almost edge-on (see e.g. O’Dell 2001), so that we can investigate how irradiation from the nearby stars penetrates the molecular cloud and changes its chemical and physical structure (see e.g. studies by Hogerheijde, Jansen & van Dishoeck 1995; Jansen et al. 1995; Andree-Labsch, Ossenkopf-Okada & Röllig 2017).

Sellgren, Tokunaga & Nakada (1990) have found spatial stratification, distinguishing an almost parallel ionization front, 3.3  |$\mu$|m PAH emission, and 2 |$\mu$|m H₂ emission. Tielens et al. (1993), using an equilibrium model of a plane-parallel irradiated cloud, have reproduced the ∼10 arcsec offset between the PAH and H₂ emission, and also between the H₂ and CO emission. The study of Tielens et al. (1993) is based on previous equilibrium model calculations of Tielens & Hollenbach (1985a,b), which are now considered as a classic model of a dense PDR.

However, recent high-resolution ALMA millimetre-band images presented by Goicoechea et al. (2016) draw a different picture, in which there is no appreciable offset between a peak of H₂ vibrational emission and an edge of ¹³CO(3–2) and HCO⁺(4–3) emission towards the Orion Bar. Locations of the ¹³CO(3–2) and HCO⁺(4–3) emission peaks coincide with the region of brightest H₂ emission that appears close to the H₂ dissociation front (see Extended Data figure 2 in Goicoechea et al. 2016). The location of the ¹³CO(3–2) emission edge supposedly coincides with a location of a CO dissociation front. Goicoechea et al. (2016) were not able to explain the ALMA data with a stationary PDR model, so they suggested that time-dependent PDR modelling could explain both the coincidence of H₂, ¹³CO(3–2), and HCO⁺(4–3) emission peaks and high brightness of HCO⁺(4–3) emission. Bright HCO⁺(4–3) emission requires gas number density of about 10⁶ cm⁻³ [see e.g. discussion of the critical density for HCO⁺ in Shirley 2015, based on calculations of Flower (1999)], and relative HCO⁺ abundance greater than 10⁻⁹.

The goal of the present study is to verify whether time-dependent astrochemical calculations can indeed reproduce the observed features of the Orion Bar, namely, no appreciable offset between the ¹³CO(3–2) and HCO⁺(4–3) emission peaks, merged H₂ and CO dissociation fronts, and bright HCO⁺(4–3) emission. To study the Orion Bar PDR, we utilize a chemo-dynamical model of an expanding H ii region, implemented as a computer code MARION. Its development is described in Kirsanova, Wiebe & Sobolev (2009), Pavlyuchenkov, Kirsanova & Wiebe (2013), Akimkin et al. (2015), and Akimkin et al. (2017). Results of the chemo-dynamical computation are then passed to the RADEX model (van der Tak et al. 2007) to calculate intensities of ¹³CO(3–2) and HCO⁺(4–3) lines and to compare them to the results of the ALMA observations presented by Goicoechea et al. (2016).

2 MODEL CALCULATION

MARION is a 1D spherically symmetric model of the dynamical and chemical evolution of gas and dust in the vicinity of a massive star. It needs to be adapted for the study of the Orion Bar PDR. We simulate the plane-parallel geometry of the molecular cloud edge by adopting a very large radius for a first inner cell, which is supposed to represent the cloud surface. The cloud thickness is 0.225 pc, and it is divided into 1500 cells, so the spatial resolution is about 30 au. We take the cloud thickness larger than e.g. 0.1 pc as found by Salgado et al. (2016) in order to provide a sufficient model time for a shock wave to reach a quasi-stationary solution prior to breaking out of the computational domain.

The radiation field is assumed to be generated by massive stars from the Orion Nebula Cluster. We constructed radiation field using parameters of five massive stars that were defined as ‘M42 inner’ by Ferland et al. (2012). The resultant radiation field contains no photons with λ ≤ 91.2 nm to save computational time. In other words, we only consider the PDR, but not the ionized region. We calculated the ‘M42 inner’ radiation field at a distance of 0.22 pc, which corresponds to the projected distance 111 arcsec from the main ionizing star to the Bar according to Pellegrini et al. (2009) and Andree-Labsch et al. (2017), and then passed the radiation field to the described above plane-parallel modification of the MARION code. The intensity of the radiation field is χ = 4.4 × 10⁴ in units of the Draine field (Draine 1978) in the wavelength range of 91.2 nm < λ ≤ 200 nm. This χ value is consistent with the FUV radiation flux near the Orion Bar PDR (Tielens & Hollenbach 1985b).

We note that at λ ≥ 200 nm this procedure produces radiation field that is orders of magnitude lower than the ‘universal’ Draine field (Draine 1978) with a long wavelength modification introduced by van Dishoeck & Black (1982), scaled by χ = 4.4 × 10⁴. Even if we supplement the radiation field generated by the five massive stars with radiation field of less massive stars from the Trapezium Cluster (Muench et al. 2002), we still do not obtain the same intensity enhancement in the λ ≥ 200 nm range as the scaled sum of the Draine field and the long wavelength modification. The usage of the radiation field constructed from the real stellar population is important for our calculations due to radiation pressure term in the dynamical equations.

After the upgrade by Akimkin et al. (2015, 2017), various dust size distributions can be accommodated in the MARION model. Here we use 48 dust size bins and WD16 (Weingartner & Draine 2001) grain size distribution for R_V = 5.5 (Bohlin & Savage 1981). The PAH mass is assumed to account for 4 per cent of the total dust mass. With these assumptions, A_V depends linearly on the hydrogen number density with the ratio A_V/N_H = 5.3 × 10⁻²² cm². As in the MARION code dust can move with respect to gas, the value of A_V/N_H varies with distance and time. To calculate photoreaction rates, we use |$\tau (5000\, \mathring{\rm A})$| computed for the current dust spatial distribution, when needed [see how |$\tau _{\rm V} \equiv \tau (5000\, \mathring{\rm A})$| changes with time in Section 5]. Heating and cooling processes, which are included in the calculation, are listed in Table 1. Radiation pressure on dust grains and momentum transfer from grains to gas are treated as in Akimkin et al. (2017).

Parameters of the model cloud and initial chemical abundances are summarized in Table 2. Gas number density n_gas at the outset of the calculation corresponds to the value estimated by Goicoechea et al. (2016). We initially put all carbon and oxygen into CO and O₂ and apply ‘H ii’ elemental abundances from the cloudy model (Ferland et al. 2013). Chemical reactions on dust surfaces along with accretion and desorption processes are not considered to save computational time. Gas-phase chemical network is taken from a website of PDR-benchmarking workshop.¹ File rate99_edited_incl_crp.dat (Röllig et al. 2007) with cross-sections for most photoreactions is taken from the Leiden photo site (Heays, Bosman & van Dishoeck 2017). The network is based on the UMIST99 ratefile (Le Teuff, Millar & Markwick 2000). A full list of species includes H, H⁺, H₂, |$\mathrm{H}_2^+$|⁠, |$\mathrm{H}_3^+$|⁠, O, O⁺, OH⁺, OH, O₂, |$\mathrm{O}_2^+$|⁠, H₂O, H₂O⁺, H₃O⁺, C, C⁺, CH, CH⁺, CH₂, |$\mathrm{CH}_2^+$|⁠, CH₃, |$\mathrm{CH}_3^+$|⁠, CH₄, |$\mathrm{CH}_4^+$|⁠, |$\mathrm{CH}_5^+$|⁠, CO, CO⁺, HCO⁺, He, He⁺, e⁻. Cosmic ray (CR) ionization rate is taken to be 5 × 10⁻¹⁷ s⁻¹, but we checked higher values of the ionization rate up to 5 × 10⁻¹⁶ s⁻¹ as mentioned in Section 5.

Röllig et al. (2007) provide eight test cases for the benchmarking of PDR models. Prior to the study of the Orion Bar, we perform all these tests with the corresponding values for element abundances and using MARION in a static mode. Results for the test cases with variable gas and dust temperatures are presented in Appendix A.

We use the RADEX model (van der Tak et al. 2007) to calculate ¹³CO(3–2) and HCO⁺(4–3) emission intensities at each time step and in each computational cell. Main parameters of the RADEX calculation are the following. Background radiation temperature is 2.73 K in each computational cell. The line width for calculation of HCO⁺(4–3) integrated intensity is fixed to 2 km s⁻¹ according to the results by Goicoechea et al. (2016). We use LAMDA database (Schöier et al. 2005) in the RADEX calculation and account for collisions of CO and HCO⁺ with H₂ molecules. Collisions of CO with atomic hydrogen (Yang et al. 2013) were also included into the calculations with RADEX. RADEX method number 2 is utilized. To convert spatial offsets to angular offsets and to produce a synthetic observational picture we adopt the distance to the Orion Bar of 414 pc (Menten et al. 2007) and the line-of-sight extent of the Bar of 0.12 pc (Werner et al. 1976). Model distributions of ¹³CO(3–2) and HCO⁺(4–3) emission are convolved with 1 arcsec ALMA beam at 350 GHz.

3 H/H₂ AND C⁺/C/CO TRANSITIONS IN A STATIONARY CALCULATION

We begin with calculating a stationary model with parameters described above to compare the computed locations of H₂ and CO dissociation fronts to those from the study by Tielens et al. (1993). For that purpose, we switch off the dynamical components of the MARION model to calculate a stationary PDR. The run continued for 1.4 Myr to assure stationary conditions (they are actually established much earlier). Results of the calculation are shown in Fig. 1. The r value in Figs 1, 2, 5, and 6 represents a distance from the inner boundary of the computational domain. In the used set-up of the MARION model, ionizing photons are not taken into account, so we cannot pinpoint the location of the ionization front. The dissociation front is defined as a location, where an abundance of either H₂ or CO drops by a factor of two relative to its initial value given in Table 2.

In Fig. 1, the H₂ dissociation front is located right next to the border of the computational domain (on this and subsequent figures it is marked with a vertical green line), while the CO dissociation front (marked with a vertical red line) is shifted deeper into the cloud. The offset between the H₂ and CO dissociation fronts is about 12 arcsec. This result is in agreement with previous equilibrium calculations by Tielens et al. (1993) and Tauber et al. (1994). The location of the CO dissociation front marks the C⁺/C/CO transition. Thus, the separation between the H/H₂ and C⁺/C/CO transitions in the stationary model would have been definitely resolved with 1 arcsec resolution.

Right-hand panel of Fig. 1 shows results of the RADEX calculation. Goicoechea et al. (2016) mention that the peak value of the ¹³CO(3–2) line intensity is a measure of the gas temperature in molecular clouds. We plot the gas temperature on the right panel of Fig. 1 and do find it to be very close to the ¹³CO(3–2) intensity (within 5 K) everywhere beyond the CO dissociation front. Brightness of the ¹³CO(3–2) emission peak is 2.5 times smaller in our results than the value obtained by Goicoechea et al. (2016), and the HCO⁺(4–3) integrated intensity is about 40 times smaller than the value found in the Orion Bar. Thus, the stationary model reproduces neither the merging of the H₂, ¹³CO(3–2) and HCO⁺(4–3) emission peaks, seen in the ALMA observations, nor the observed line intensities.

4 TIME-DEPENDENT LOCATION OF H/H₂ AND C⁺/C/CO TRANSITIONS

In this section we present results of a time-dependent PDR model for a cloud that is initially completely molecular (see Table 2) and then instantly becomes irradiated by strong FUV field. Radiation heats up gas and drives a compression wave into the cloud. At the same time, radiation pressure pushes dust grains into the cloud, and the momentum transfer from dust to gas significantly affects the gas dynamics. FUV quanta dissociate molecules, and both H₂ and CO dissociation fronts also move through the cloud. We stop the calculation when the H₂ dissociation front crosses the outer border of the computational domain that occurs after about 13 kyr. Fig. 2 shows physical, chemical, and emission properties of the model cloud for three representative time moments, which correspond to 2.3, 5.0, and 8.3 kyr since the beginning of the calculation and are further referred to as early, intermediate, and late time moments. The moments show three varieties of mutual locations of H₂ and CO dissociation fronts, mutual locations of ¹³CO(3–2) and HCO⁺(4–3) emission peaks, and their intensities.

At the early time moment gas density distribution in the considered region is still nearly flat, but a density enhancement starts to appear in the leftmost part of the region. This is where the H₂ dissociation front is located. The CO dissociation front initially is located somewhat deeper into the cloud; the separation between the H₂ and CO dissociation fronts is ≈18 arcsec. Transition from H to H₂ occurs closer to the border of the cloud because of effective self-shielding of H₂ molecules against photo-dissociation. While some radiation capable of dissociating CO molecules is absorbed by H₂ molecules (mutual shielding), photons dissociating CO can penetrate deeper into the cloud with the nearly flat density distribution. This is why the CO dissociation front forms far ahead of the H₂ dissociation front at this time (see van Dishoeck & Black 1988). It is interesting that abundance of HCO⁺ is slightly larger than CO abundance in the leftmost part of the computational domain, between the HCO⁺(4–3) peak and the H₂ dissociation front. The theoretical HCO⁺(4–3) emission is much fainter at this moment than the observed emission.

As the structure of the transition region between the cloud and the PDR evolves, HCO⁺(4–3) emission brightens in the direction of the H₂ dissociation front, as shown in the middle row of Fig. 2. The HCO⁺(4–3) emission peak nearly coincides with the gas density enhancement that continues to grow ahead of the H₂ dissociation front. The separation between the H₂ and CO dissociation fronts is about 11 arcsec, which is still larger than ALMA resolution at 350 GHz. Maximum HCO⁺(4–3) integrated intensity is ≈70 K km s⁻¹, which is in agreement with the average value observed in the Orion Bar. ¹³CO(3–2) peak emission is seen towards the CO dissociation front in the gas that is heated by radiation but is still mostly undisturbed. There is also a secondary ¹³CO(3–2) emission peak at r = 17 arcsec that is nearly coincident with the gas density enhancement and the HCO⁺(4–3) emission peak. Later this minor ¹³CO(3–2) emission peak will merge with the bulk of the CO-emitting gas.

After about 7–8 kyr of evolution (the late time moment), the structure of the transition region is finally established. Its main feature is that the H₂ and CO dissociation fronts are nearly merged with each other (see the bottom rows in Figs 2 and 3). The distance between the fronts is only 1.0–1.5 arcsec. Since this moment, the merged H/H₂ and C⁺/C/CO transitions move together away from the source of FUV radiation.

The separation between the HCO⁺(4–3) and ¹³CO(3–2) emission peaks is also about 1 arcsec after the structure is established. This value is close to ALMA resolution at 350 GHz and is, thus, consistent with the ALMA observations. The n_gas enhancement exceeds an order of magnitude at this moment. The density and temperature values in the compressed layer are in agreement with estimates obtained by Goicoechea et al. (2016) using the equilibrium PDR model and HCO⁺(4–3) and ¹³CO(3–2) emission analysis (in Fig. 2 these estimates are shown with hatched rectangles). Similar gas number density and temperature in the PDR were obtained by Sorochenko & Tsivilev (2000) after the analysis of carbon radio recombination lines and infrared fine-structure lines of C ii and O i (n_gas = (1.3 ± 0.4) · 10⁵ cm⁻³, T_gas = 215 ± 32 K). Both the peak intensity of the ¹³CO(3–2) line and the maximum integrated intensity of the HCO⁺(4–3) line are consistent with the corresponding observed values averaged through the front in the Orion Bar. This consistency is preserved thereafter within a factor of three as illustrated in Fig. 4.

The theoretical ¹³CO(3–2) peak intensity does not differ from the observed value by more than a factor of two over the entire considered time span. This is because the peak intensity of an optically thick transition is close to the gas temperature in the molecular region, and this temperature does not vary much at the CO dissociation front and beyond. The same result is also obtained in the stationary run. Thus, the ¹³CO(3–2) emission is indeed a reliable gas temperature measure in the Orion Bar, as suggested by Goicoechea et al. (2016).

The situation with the peak of HCO⁺(4–3) emission is less trivial. At the intermediate and late time moments, it nearly coincides with the peak of the gas number density. This is consistent with the suggestion made by Goicoechea et al. (2016) on the correspondence of the HCO⁺(4–3) emission peak to the gas density peak. The integrated intensity of the HCO⁺(4–3) line persistently grows through the early and the intermediate times, reaching the value of about 100 K at the late time moment and then decreases down to about 30 K. There is a significant scatter in the computed parameters of the HCO⁺ line, because they are very sensitive to the environmental conditions (a smoothed curve is shown in Fig. 4). One can see in the bottom row of Fig. 2 that the peak of the HCO⁺ emission is located at the illuminated side of the CO dissociation front. Evidently, HCO⁺ molecules bind some carbon atoms freed from photo-dissociated CO molecules. The density peak is located in the region right next to the CO dissociation front, where HCO⁺ abundance sharply decreases with distance. Even minor variations in the combination of HCO⁺ abundance and the total gas density may cause large scatter in the emission parameters. We note that the chemical network we use for this study is reduced and adopted to reproduce the locations of CO and H₂ dissociation fronts in the equilibrium calculations (Röllig et al. 2007). But it appears to be good enough to reproduce HCO⁺ abundances as well. We note that the HCO⁺(4–3) emission has optical depth greater than one.

In Fig. 5 we show column densities of the species which contain most carbon atoms after the merging of the fronts. Various components of the carbon hydrogenation chain becomes abundant, with the most abundant CH₂, CH, CH|$_3^+$|⁠, and CH⁺. In Fig. 5 we also compare our theoretical results with observed abundances of several species in the Orion Bar measured by Nagy et al. (2013, 2017) toward the ‘CO⁺ peak’ with the Herschel telescope. Spatial resolution of the Herschel data is worse than the ALMA data discussed here. So, the abundances from Herschel are averaged over much larger area than the horizontal axis range we show in Fig. 5. It is seen that we reproduce abundances of atomic carbon, C⁺, CH⁺, HCO⁺ within an order of magnitude at least between the H₂ and CO dissociation fronts. We also reproduce abundance of CH, which has an uncertainty range much smaller (about a quarter of magnitude) than the other mentioned molecules, at least around the CO dissociation front.

We also found that the ratio of C/CO relative abundances in the PDR are in general agreement with 0.05–0.1 values reported by White & Sandell (1995). The ratio rises from 0.01 at the H₂ dissociation front to 0.1 at the CO dissociation front. We consider the abundance of CO⁺ between the H₂ and the CO dissociation fronts as being in agreement with Stoerzer, Stutzki & Sternberg (1995) and Fuente et al. (2003), who determined the abundance using a single-dish telescope.

The used chemical network does not contain hydrocarbon molecules with more than one carbon atom as well as species containing both carbon and nitrogen. Such a restriction of the network does not liberate the carbon atoms from the chemical chains leading to formation of CO, HCO⁺, and CH₄. This is probably why we obtain high abundances of CH₂ and CH|$_3^+$| which are intermediaries of the chemical chain for the CH₄ formation. Richer chemical network is needed to model the detailed PDR chemistry which is beyond the scope of the present study.

The conclusion of this section is the following. In a dynamically evolving PDR, the H₂ and CO dissociation fronts rapidly merge and move together thereafter. By the end of the computation run a thin compressed gas layer forms, which is no more than 5 arcsec wide (0.01 pc), and the thickness of the topmost density peak is about 2 arcsec (0.005 pc). CO and HCO⁺ emission radial profiles look like sharp peaks. Their widths and respective intensities are similar to what is observed in the Orion Bar. The observed bright HCO⁺(4–3) emission is only matched by the theoretical one after the dissociation fronts merge and propagate into the surrounding molecular gas.

To explore the importance of the dynamics we made an additional test. The physical structure (gas and dust number densities) were taken as the initial conditions for the stationary calculations similar to the one presented in Section 3. Thermal structure for gas and dust distributions was calculated from the pressure equilibrium condition. We found that no bright HCO⁺(4–3) emission appears in the shocked gas at such conditions because of low T_gas in the region with high n_gas.

5 DISCUSSION

Our goal is to verify whether the observed closeness of the H₂ and CO dissociation fronts in the Orion Bar is related to the dynamical evolution of the object. Indeed, in the previous sections we showed that the near coincidence of ¹³CO(3–2) and HCO⁺(4–3) emission peaks observed in the Orion Bar is not reproduced in the stationary PDR model and appears naturally in the dynamical model. Merging of the emission peaks and the corresponding dissociation fronts should not be considered as a certain temporary feature of the presented simulation. While we do perform time-dependent calculations, the model time cannot be interpreted as an ‘age’ of the Orion Bar. First, we do not know when the PDR has been formed. Secondly, the Orion Nebula is a blister-type H ii region, so that the 1D MARION model cannot be applied to this object ‘as is’. Our results should rather be interpreted as a ‘stationary dynamic situation’, which can only be realized when we take into account the action of pressure on to the molecular gas.

This result is consistent with conclusions that have been reached in some previous studies. For example, Hosokawa & Inutsuka (2005) presented hydrodynamic simulations of expanding H ii regions that included calculation of the H₂ and CO dissociation. They showed that the dissociation fronts of H₂ and CO become closer to each other during the expansion and nearly merge in 400 kyr after the beginning of the expansion for physical conditions pertaining to H ii region Sh2-104. The calculation of Hosokawa & Inutsuka (2005) is the first chemo-dynamical simulation of an expanding H ii region with a real hydrodynamic model. We should also mention a study of Störzer & Hollenbach (1998) who used analytic approach to show that in non-equilibrium calculations the C⁺/CO transition layer becomes closer to the H₂ dissociation front compared to the equilibrium model. They also concluded that non-equilibrium effects in the Orion Bar are probably small. This conclusion was in agreement with contemporary observational results, however, now we may argue that dynamical effects are important in highly FUV-irradiated gas.

Studies of PDRs commonly predict that the HI/H₂ transition occurs at A_V ≈ 2 mag, and the C⁺/C/CO transition occurs at A_V ≈ 4 mag (Tielens & Hollenbach 1985a; Tauber et al. 1994). Our results are generally consistent with this picture as seen in the left panel of Fig. 6, where A_V = τ_V/1.086, and τ_V is computed for the current dust distribution. The diagram shows locations of the H₂ and CO dissociation fronts in terms of A_V. It is seen that after the stationary dynamic situation has been established on a time-scale of about 7–8 kyr, which is significantly smaller than the typical age of an expanding H ii region around an O-type star (e.g. Spitzer 1978), the fronts moves slowly from A_V  = 3 mag and 4 mag for H₂ and CO, respectively, to A_V  = 4 mag (H₂) and 5.5 mag (CO). The front convergence is caused by the appearance of the density enhancement. Because of the sharp density increase in the enhancement the transition region between the illuminated medium to the shaded medium becomes quite narrow, and all the major chemical transformations occur over a very limited distance range. This result is naturally absent in stationary models that assume constant density across the PDR. One may argue that the nearly coincident dissociation fronts can be naturally reproduced by a constant pressure model with a non-uniform density distribution. However, our tests indicate that it is impossible in this case to reproduce bright HCO⁺(4–3) emission.

Joblin et al. (2018) claim the existence of thin (few 10⁻³ pc corresponding to 0.5−1 arcsec) emission layers at high thermal pressures about 10⁸ K cm⁻³ in the Orion Bar PDR. Our computed physical conditions in the compressed gas layer for the late time moment presented in the Fig. 2 agree with their suggestion. It is noteworthy that the similar value of the thermal pressure in the compressed layer was obtained by Bron et al. (2018) using photoevaporation flow model. Edge-on orientation of the Orion Bar PDR does not allow testing gas kinematics, although studies of this kind could possibly give us a direct answer about the dynamical model, appropriate for the region. The Orion Bar is irradiated not only by O7 θ¹C Ori star from north-west but also by O9.5 θ²A Ori (Sota et al. 2011) from south-east. This other star may also influence the chemical structure of the Orion Bar. While θ¹C Ori dominate ionization in the inner part of the Orion Nebula (O’Dell, Kollatschny & Ferland 2017), we still believe that influence of θ²A Ori on the molecular abundances towards the south-east edge of the Bar should be further investigated.

It is worthwhile to relate our results to those presented by Pellegrini et al. (2009) and Shaw et al. (2009). These authors have explored the Orion Bar observations, presented by Tielens et al. (1993) and Tauber et al. (1994) for a different position, where H₂ and CO(1–0) emission peaks appear to be separated by an angular distance of about 10 arcsec. Pellegrini et al. (2009) and Shaw et al. (2009) have been able to reproduce both the location and the intensity of the emission peaks for the [S ii], H α, H₂, and CO(1–0) lines, using a stationary model with the enhanced magnetic field and CR ionization rate. The value of the CR ionization rate had to be increased by nearly four orders of magnitude in order to reproduce observed CO and H₂ intensities in simulations of Pellegrini et al. (2009) and Shaw et al. (2009).

The ALMA observations and our simulations draw a somewhat different picture, in which H₂, ¹³CO(3–2), and HCO⁺(4–3) emission peaks nearly coincide, without enhancing the cosmic ray ionization rate above the standard ISM value. Even more so, raising the ionization rate by just one order of magnitude may lead to some serious chemical consequences. We tested higher values of the CR ionization rate up to 5 × 10⁻¹⁶ s⁻¹, considering recent results of Sorochenko & Smirnov (2010), Indriolo & McCall (2012), and Neufeld & Wolfire (2017). The peak values of the ¹³CO(3–2) and HCO⁺(4–3) emission intensity in the PDR remain the same as in our basic model. The difference, significant for the present study, is the high HCO⁺ abundance in the undisturbed molecular region, where the abundance is defined by reactions with H|$_3^+$| ion.

These different results may be somehow related to the non-uniform matter distribution within the Orion Bar. Goicoechea et al. (2016) and Andree-Labsch et al. (2017) demonstrate a clumpiness of the Orion Bar. So, some disagreement between our results and the observed picture can be caused by the fine structure of the PDR unresolved in our model. Also, our restricted simulations do not include modelling of the Orion Nebula’s ionization front so we cannot relate mutual positions of the ionization and dissociation fronts. In conclusion we emphasize that the model suggested by Pellegrini et al. (2009) and Shaw et al. (2009) successfully describes spatially resolved observations in optics and infrared while our model reproduces the spatially resolved observations in millimeter wavelengths. In this sense, these two models compliment each other. A more advanced modelling that includes both dynamics and chemistry, considers the ionized, atomic and molecular gas, and takes into account more molecular line data (e.g. Parikka et al. 2018; Goicoechea et al. 2019), will definitely help to resolve the issue.

Fig. 6 (left) shows that a value of A_V corresponding to the location of the CO dissociation front drops by one magnitude at the late time moment when the H₂ and CO dissociation fronts merge at the border of the PDR. This drop reflects a fact that dust grains and gas can move relative to each other in the used version of the MARION model. This is further illustrated in the right panel of Fig. 6, which shows how relation between the hydrogen column density and visual extinction changes in the PDR and the molecular cloud. Due to dust motion through the gas and variable dust-to-gas mass ratio, we see different N_H/A_V ratios. Within the molecular cloud this ratio is the same as in the original dust model. At the density peak and within the PDR, dust is depleted due to the drift driven by radiation pressure. The important difference is that at the early and intermediate time moments the CO dissociation front is located farther into the cloud, far ahead of the dust depletion zone, while at the late time moment and afterwards the CO dissociation front is located within the dust depletion zone. This is why A_V is lower at the front location. In the absence of the abundant dust, the CO protection from photodissociation is mostly provided by H₂ shielding. Thus, CO survival is directly related to H₂ survival. This effect further narrows the gap between the two dissociation fronts. The enhanced N_H/A_V ratio found in the present study is in agreement with previous findings by Abel et al. (2004, 2006, 2016) in the Orion Veil, which is physically associated with the Trapezium stars and the Orion Nebula.

To strengthen this conclusion, we repeated the calculation assuming that dust is frozen to gas. In this case, the evolution of the PDR including the front convergence proceeds somewhat slower. The separation between the fronts becomes smaller than the ALMA resolution after about 13 kyr, which is still far less than other relevant time-scales of the PDR evolution.

Salgado et al. (2016) report the decrease of dust opacity in the Orion Bar PDR. Specifically, they find a factor of four and five decrease of infrared and ultraviolet dust opacities, respectively, in the ionized gas compared to the value for diffuse interstellar gas. Factor of five decrease of the infrared dust opacity is found towards the PDR itself. These results support our conclusion on the enhanced N_H/A_V ratio in the Orion Bar PDR. Salgado et al. (2016) obtain the increase of the dust optical depth at 19.7 |$\, \mu$|m from 0.025 in the ionized gas to 0.2 in the PDR. These values are also in agreement with the variation of τ_19 μm between the H₂ and CO dissociation fronts obtained with the presented MARION model.

Some important remarks should be made before we conclude. First, we tried various initial conditions for the calculations. Specifically, we tested initial densities in the range of 5 × 10⁴−10⁵ cm⁻³ and χ in the range of 10⁴−10⁵. We also considered the Orion Veil A_V/N_H = 3.3 × 10⁻²² cm² ratio for the calculations of the Orion Bar PDR (Abel et al. 2004, 2006, 2016) and found no difference in terms of the ¹³CO(3–2) and HCO⁺(4–3) emission intensity in the PDR with our basic model. The mentioned test calculations bring us to a conclusion that radial profiles of various molecular lines, similar to those presented by Goicoechea et al. (2016) for CO and HCO⁺, and the abundance profiles are highly desirable to testify the predictions of PDR models and confine the physical parameters. Ever for such an extensively studied region as the Orion Nebula and its surroundings we did not find the abundance profiles in the literature to compare with the predictions from the Fig. 5.

6 CONCLUSION

We argue that dynamical effects in the FUV-irradiated gas, namely, the formation of the compressed layer at the border of the PDR and dust drift driven by the radiation pressure, are responsible for the observational appearance of the Orion Bar on the ¹³CO(3–2) and HCO⁺(4–3) maps at 350 GHz produced by ALMA. The interaction between the heated PDR and the cold molecular cloud leads to the formation of the density enhancement that slowly moves into the cloud. As the mutual location of the transition regions depends on A_V, and its value grows on a short distance scale within the compressed layer, the regions become close to each other, so that their angular separation is too small to be resolved with the ALMA observations. Another important factor that speeds up the formation of the observed configuration is the dust expulsion from the PDR by the stellar radiation pressure. The resultant chemical and thermal structure leads to the correspondingly coincident locations of the H₂, ¹³CO(3–2) and HCO⁺(4–3) emission peaks. Finally, high gas number density and high relative abundance of HCO⁺ cause bright HCO⁺(4–3) line emission, which appears on the illuminated side of the CO dissociation front. These observed features are not reproduced by the stationary model and appear naturally in dynamical model of the PDR. We also confirm suggestions of Goicoechea et al. (2016) that HCO⁺(4–3) emission corresponds the gas density peak in the compressed region and that ¹³CO(3–2) line emission is a reliable measure of the gas temperature in the Orion Bar.

SUPPORTING INFORMATION

Figure A1. Comparison of MARION calculations with the results of the PDR benchmark meeting.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We are thankful to Ya. N. Pavlyuchenkov for fruitful discussions of the radiation transfer with RADEX. The anonymous referee is gratefully acknowledged for his/her comments and suggestions which allows us to improve the paper.

The study is supported by Russian Foundation for Basic Research, grant number 16-02-00834-a.

Footnotes

REFERENCES

APPENDIX: COMPARISON WITH OTHER PDR CODES

We performed tests of the MARION model using results of the PDR-benchmarking workshop held in 2004 (Röllig et al. 2007). All F1–F4 (fixed gas and dust temperature) and V1–V4 (variable gas and dust temperature) models are tested. We show here comparison of the MARION results with Röllig et al. (2007) for V1–V4 models. In Fig. A1 results of MARION are shown by a solid line and results of Meudon, HTBKW, KOSMA-τ, Leiden, Meijerink, and Sternberg are shown by dashed lines. All the references to the codes and initial conditions for the test calculations are given in Röllig et al. (2007).