© V. Maillard et al. 2021

1 Introduction

In the atomic envelop surrounding a molecular cloud (Wannier et al. 1983; Andersson et al. 1991), the transition where molecular hydrogen turns into atomic hydrogen is called the atomic-to-molecular hydrogen (H/H₂) transition, or the photodissociation front. Knowledge about the location of the H/H₂ transition gives us information on the total fraction of molecular gas in a cloud. It is indeed in the cold and dense molecular medium that stars are formed through gravitational collapse of the gas. The H₂ mass fraction is linked to star-formation thresholds by the Schmidt–Kennicutt relations (Schmidt 1959; Kennicutt 1998, see also Tacconi et al. 2020 for a review on molecular gas mass in galaxies and its link to star formation rates and the evolution of galaxies). Secondly, the intensity of the lines of the molecular tracers depends on the depth of the H/H₂ transition due to the specific chemistry occurring in the warm H₂ gas.

Estimating the total H I column density at the surface of neutral interstellar clouds is a key problem for tackling several important astrophysical tasks: to deduce the H₂ mass fraction intervening in the Schmidt–Kennicutt law of star-formation (Villanueva et al. 2017), to deduce the CO-to-H₂ conversion factor (Schruba et al. 2017; Dessauges-Zavadsky et al. 2017; Pineda et al. 2017), or even to study the dynamics of a cloud (Valdivia et al. 2016; Beuther et al. 2020). Some direct estimates of the total H I column density were also performed through 21 cm observations (Kim et al. 1998, 2003; Stanimirović et al. 1999, 2014; Staveley-Smith et al. 2003; Stil et al. 2006; Peek et al. 2011; Bihr et al. 2015; Beuther et al. 2016).

A first analytical model of the H/H₂ transition was presented by Sternberg (1988) who derived an analytic formula for the total H I column density for a cloud with planar geometry illuminated by a beamed radiation field. He provided results as a function of the far-ultraviolet (FUV) radiation intensity, the cloud gas density, and the metallicity. In two papers, Krumholz et al. (2008, 2009) presented new models for the H/H₂ transition taking into account the metallicity dependence of the molecular mass fraction in galaxy disks and multidimensional geometry of finite-sized clouds embedded in ambient isotropic fields. More recently, Sternberg et al. (2014) and Bialy & Sternberg (2016) updated the analytic theory for one-dimensional (1D) planar slabs and spheres for beamed and isotropic radiation fields, aiming to apply them to global galaxy evolution studies. Their results were compared and validated by Meudon PDR Code (Le Petit et al. 2006) computations. These models are widely used to predict the column density of the atomic surface layer (Wong et al. 2009; Bolatto et al. 2011; Lee et al. 2012, 2015; Roman-Duval et al. 2014; Bialy et al. 2015, 2017a; Schruba et al. 2018; Noterdaeme et al. 2019; Klimenko et al. 2020). They have been extended to other problems such as the determination of the C⁺ /H₂ ratio (Nordon & Sternberg 2016) or the HD/H₂ ratio in the diffuse interstellar gas (Balashev & Kosenko 2020). Also, Bialy et al. (2017b) explored the effects of turbulence on the H/H₂ transition using analytic and numerical methods with an application to the Perseus molecular cloud. Recently, Sternberg et al. (2021) extended the theory to the H/H₂ transition inthe dust-free primordial interstellar gas.

In the analytic theories presented above, the medium was assumed to be in a static and stationary state. However, clues about dynamical effects were recently found. Herschel observations of excited lines in strongly UV-illuminated galactic photodissociation regions (PDRs) (Joblin et al. 2018; Wu et al. 2018) revealed a thin compressed surface layer (10⁻³ pc) with a thermal pressure well above its environment at the edge of the PDR (P_th ~ 10⁸ K cm^-3). The gas thermal pressure was found to be strongly correlated to the UV field intensity. Moreover, Atacama Large Millimeter Array (ALMA) observations (Goicoechea et al. 2016, 2017) revealed the spatial structure of the Orion Bar PDR with an excellent resolution, showing that hot molecular tracers were emitting from a thin layer at high values of pressure and temperature. It was proposed that photoevaporation at the ionization front could explain the high pressures found in the neutral PDR and their correlation between the thermal pressure and the UV field intensity. Numerical models of photoevaporating PDRs with the time-dependent and dynamical Hydra PDR Code (Bron et al. 2018) reproduced the pressure structure and the correlation of the pressure with the UV radiation field.

When a star illuminates a neutral cloud, its surface is heated by the radiation field. If the pressure in the H II region is not sufficient to contain the heated gas (e.g., dense globules embedded in the H II region, Bertoldi 1989 and Bertoldi & Draine 1996, blister H II regions where the ionized gas is free to escape to the surrounding diffuse medium, Israel 1978 and Tenorio-Tagle 1979), a photoevaporation flow of ionized gas streams from the neutral cloud into the H II region. The theory of ionization fronts was first derived by Kahn (1954) following the suggestion of Oort & Spitzer (1955), and is presented in textbooks (e.g., Spitzer 1978, Chapter 12.1, Draine 2011b, Chapter 37). This theory shows that steady ionization fronts cannot exist for a specific range of ionizing photon flux to neutral gas density ratios, F_EUV∕n_H, which thus separates two types of allowed solutions, called R-type and D-type. A R-type ionization front propagates supersonically relative to the neutral gas, and can occur either during the initial phase of development of an H II region or in a later phase when encountering a steep decreasing density ramp (Franco et al. 1989, 1990). A D-type front propagates subsonically relative to the neutral gas, the gas is accelerated when passing through the ionization front and a pressure jump is maintained between the ionized and the neutral gas. When conditions fall in the forbidden range between R-type and D-type, a (usually low-velocity) shock propagates ahead of the ionization front in the neutral gas, raising the neutral gas density to reach D-type conditions. After an initial phase, this shock front overtakes the H/H₂ transition. Nearly critical D-type fronts are thought to be omnipresent in evolved, blister-type (Israel 1978, or “champagne flow", Tenorio-Tagle 1979) H II regions found in many star forming regions and where archetypal dense PDRs are found (e.g., the Orion Bar in the Orion Nebula, O’dell 2001). Advection of the neutral gas from the neutral PDR region through these D-type fronts can thus maintain the H/H₂ chemistry out of equilibrium.

This photoevaporation mechanism, proposed to explain the dynamics of bright PDRs illuminated by young stars at the edges of molecular clouds bordering blister H II regions, has been studied in the past to describe the outer layer of small neutral molecular globules embedded in compact H II regions by Bertoldi & Draine (1996) and Störzer & Hollenbach (1998). The upcoming James Webb Space Telescope (JWST) will give us access to high angular resolution observations of many H₂ rotational and ro-vibrational lines that will help probe the high pressure layer of PDRs and constrain the photoevaporation mechanism.

Bertoldi & Draine (1996) studied the PDR structure at the surface of a clump embedded in a H II region by taking into account the expansion of newly ionized gas in the H II region separating the cloud from one or more O stars. Thanks to the curvature of the modeled clump, the gas can expand freely, so that the ionizing radiation erodes the gas. This translates into a propagating ionization front, which is the transition between ionized hydrogen H⁺ and atomic hydrogen. They linked the Extreme-UV (EUV, >13.6 eV) and Far-UV (FUV, $\in \left[6-13.6\right]$ eV) radiation flux arriving at the respective fronts (ionization and photodissociation) with their velocities. When the ionization front is slower than the dissociation front, an atomic region is developed, yet the velocity of the dissociation front is too high for stationary models of PDRs to still be applicable. In the other case, where the dissociation front is slower than the ionization front, the fronts are not separated and the atomic region does not exist. For a wide range of parameters, the atomic region is controlled by nonequilibrium effects, leading either to its nonexistence or to a smaller value of total H I column density than predicted with stationary models.

A following paper, Störzer & Hollenbach (1998) presented a model of the thermal and chemical structure of the PDR including a propagating ionization front. They solved the full time-dependent structure with a constant ionization front velocity and a uniform density that does not vary with time. They obtained similar conclusion to those of Bertoldi & Draine (1996) about the merging of the fronts, and found that the H₂ rotational and vibrational line emission can be enhanced by a factor of 3.

As of today, with Herschel, ALMA or upcoming JWST observations, dynamical effects can be highlighted in PDRs, leading to a need for an extension of Sternberg et al. (2014) to include the effects of photoevaporation. In this paper we extend the H/H₂ transition model to take into account the dynamics induced by the photoevaporation mechanism in a new semi-analytical theory. As in Störzer & Hollenbach (1998), we model the propagation of the ionization front with a constant velocity, but find the structure of the PDR by solving the stationary equation of the chemistry including advection. This equation depends only on the physicalconditions that are the velocity of the ionization front, the FUV field intensity, the gas density, and metallicity.

Our goal is to compute the impact of the dynamical effects on the position and width of the H/H₂ transition, and its consequences on total atomic hydrogen column density. In Sect. 2, we discuss our semi-analytical approach to solve the equation driving the H₂ abundance. In Sect. 3, we present the results of the dynamical models, first for a standard metallicity and then as a function of the metallicity, and propose a new formula for the total atomic hydrogen column density. In Sect. 4, we discuss our semi-analytical results and compare them to observations and to Hydra PDR Code models.

2 Analytic overview

In this section, we present an analytic model including the most important processes controlling the H/H₂ transition in interstellar clouds exposed to FUV radiation fields and considering the propagation of the ionization front.

2.1 Model geometry

We consider a 1D plane-parallel semi-infinite slab of neutral gas irradiated by stellar UV photons. The ionization front, from where we start our model, is assumed punctual similarly to other analytical models: Sternberg et al. (2014); Krumholz et al. (2008, 2009) or to usual PDR codes: Meudon PDR Code, Le Petit et al. (2006), so that the EUV component, responsible for the ionization of hydrogen up to the ionization front, is assumed to be fully extinguished, and the radiation field penetrating the slab is limited to FUV photons, which will create a photodissociation front of H₂.

The photoevaporation dynamics of the ionized gas allows the ionization front to advance into the neutral gas with a velocity v_IF. In the reference frame of the ionization front, this is equivalent to an advection of the neutral gas through the PDR and the ionization front with an equal velocity in the opposite direction (toward the ionization front). To account for the impact of this dynamic, we include a uniform advection velocity toward the ionization front, in a similar fashion to the study of Störzer & Hollenbach (1998). We assume that the shock front preceding the ionization front has already separated enough from the ionization front so that the H/H₂ transition is included in the shocked layer. For the purpose of our model describing the H/H₂ transition, we can thus assume a constant advection velocity in the neutral region. This velocity is taken as a free parameter, and the possible values and mechanisms controlling this velocity are discussed in Sect. 4.

We also assume the neutral region to have a uniform density n_H, irradiated by a unidirectional flux of dissociating photons from the 912–1108 Å Lyman–Werner (LW) band. In this study, we assume stationary state, since the timescale to reach a stationary PDR structure $t_{\text{eq}}\approx 500 \text{yr} \left[\left({10}^6 {\text{cm}}^{-3}\right)/n_{\text{H}}\right]$ (Hollenbach & Natta 1995) is short compared to a cloud lifetime. This timescale is derived from a static view, but gets shorter when the propagation of the ionization front is taken into account as seen in Sect. 2.2. We can therefore use the stationary state hypothesis.

2.2 Physical model

To consider the impact of photoevaporation on the H/H₂ transition, it is thus necessary to expand on previous studies that considered the H/H₂ transition ina static slab of gas. The propagation of the ionization front takes the form of an advection velocity of the neutral gas. We do not solve the dynamics of the gas through the ionization front, but only account for the impact of this dynamics by assuming a steady flow with a constant advection velocity. We can then solve for the stationary but out-of-equilibrium state of the H/H₂ chemistry inthe presence of this steady flow.

Aside from advection, the processes taken into account for the formation and destruction of H₂ are its formation on dust and photodissociation by FUV photons. Similarly to Sternberg et al. (2014), we neglect here the impact of cosmic rays,which would maintain a very small fraction of atomic hydrogen deep inside the cloud. In the following, we use most of Sternberg et al. (2014) notation conventions.

We take the rate of the H₂ formation on dust per hydrogen nucleon as the average value for diffuse gas of Jura (1974), that we assume to be linearly proportional to the metallicity: $$R=3\times {10}^{-17}{Z}^{\prime } {\text{cm}}^3{\text{s}}^{-1},$$(1)

where Z′ is the metallicity defined as the ratio between heavy elements abundances in the cloud and in the solar photosphere.

Let N be the hydrogen nuclei column density so that N = N₁ + 2N₂ with N₁ the column density of atomic hydrogen and N₂ the column density of H₂. We also introduce n₁ and n₂ the number densities of H and H₂ respectively. For the photodissociation of H₂, we assume that the dissociating FUV radiation is only absorbed by grains and H₂ (we neglect absorption by gas phase species such as atomic hydrogen, carbon or CO, refer to Appendix B for more details) so that the H₂ photodissociation rate at a column density N is: $${\mathcal{P}}_{{\text{H}}_2}(N)={\mathcal{P}}_0{G}_0{f}_{\text{shield}}(N_2){\text{e}}^{-{\sigma }_{\text{g}}N},$$(2)

with $${\mathcal{P}}_0=3.3\times {10}^{-11} {\text{s}}^{-1},$$(3)

and $${\sigma }_{\text{g}}=1.9\times {10}^{-21}{Z}^{\prime } {\text{cm}}^2,$$(4)

where G₀ is the energy density of the radiation field integrated between 6 and 13.6 eV and expressed in units of the corresponding integral for the ISRF of Habing (1968)¹, u_Habing = 2.65 × 10⁻¹⁴ erg cm⁻³ (Draine 2011b). We use here the G₀ at the cloud surface (i.e. seeing the incoming radiation field from 2π sr). f_shield is the H₂ self-shielding function. The parameter ${\mathcal{P}}_0$ is only an approximate value of the H₂ dissociation rate at the cloud surface for G₀ = 1 as the true value depends on H₂ ro-vibrational level populations. The value computed by Sternberg et al. (2014) with respect to an isotropic Draine ISRF (Draine 1978) in free space is D₀ = 5.8 × 10⁻¹¹ s⁻¹. Our ${\mathcal{P}}_0$ value is the corresponding value for G₀ in Habing units and taken at the cloud surface. Moreover, σ_g defined in Eq. (4) is the dust effective LW-photon absorption cross section per hydrogen nucleus. We take the same value as Sternberg et al. (2014), derived from a standard interstellar extinction curve with a total-to-selective extinction ratio R_V ≡ A_V∕E(B − V) = 3.1. We take a linear dependence of σ_g to metallicity following the latest prescriptions seen in Draine et al. (2007) and Wiseman et al. (2017) and first given in Franco & Cox (1986) with a value at solar abundances very close to the one used here. Sternberg et al. (2014) also used σ_g proportionalto metallicity. A constant dust-to-metal ratio for a wide range of metallicities supports this linear dependence, as found by Zafar & Watson (2013), and more recently by Chiang et al. (2021) for nearby galaxies with a ratio of about 0.5. Mattsson et al. (2014) indicate that the dust-to-metal ratio is kept constant by stellar dust production and dust growth, counteracting dust destruction, except perhaps for low metallicities. Finally, our expression for the H₂ self-shielding function is taken from Draine & Bertoldi (1996): $$\begin{array}{lll}{f}_{\text{shield}}(N_2)\hfill & =\hfill & \frac{0.965}{{(1+y/b_5)}^2}\hfill \\ \hfill & \hfill & +\frac{0.035}{{(1+y)}^{0.5}}\times \exp \left[-8.5\times {10}^{-4}{(1+y)}^{0.5}\right],\hfill \end{array}$$(5)

with y = N₂∕5 × 10¹⁴ cm⁻², b₅ = b∕1 km s⁻¹, and b = 2 km s⁻¹ being the Doppler broadening parameter.

The expression of the equation of the density of H₂ is: $$\frac{\text{d}{n}_2}{\text{d}t}=R{n}_{\text{H}}{n}_1-G_0{\mathcal{P}}_0{\text{e}}^{-{\sigma }_{\text{g}}N}{f}_{\text{shield}}(N_2)n_2+v_{\text{IF}}\frac{\text{d}{n}_2}{\text{d}x}.$$(6)

where v_IF > 0 is the velocity of the cloud in the ionization front frame and n_H is the uniform and stationary density of the gas. On the right-hand side, the first term represents H₂ formation on dust, the second term is its photodissociation by FUV, taking into account H₂ self-shielding and dust absorption, and the third term is the advection in our 1D geometry. The advection is a positive term in this equation as it acts as a source term for H₂. This can be interpreted as a stationary flow of H₂ coming from the infinite reservoir that is the molecular region.

The timescale to reach a stationary state is inversely proportional to formation, so that the usual static approach (v_IF = 0) gives $t_{\text{eq}}=1/(2R{n}_{\text{H}})\approx 500 \text{yr} \left[\left({10}^6 {\text{cm}}^{-3}\right)/n_{\text{H}}\right]$ as discussed earlier in Sect. 2.1. But in our advection case, the advection acts like another formation term, so that t_eq is reduced. The time to reach a stationary state gets shorter as the velocity of the ionization front increases.

We introduce τ_g = σ_gN the dust optical depth in the cloud and τ_1,2 = σ_gN_1,2 the optical depths due to the dust associated with the atomic H and the H₂ gas, respectively. We also introduce x_1,2 = n_1,2∕n_H the fractional abundances of atomic H and the H₂. We can now rewrite Eq. (6) in dimensionless form and consider the stationary state $\left(\partial /\partial t=0\right)$: $$0=1-\frac{\text{d}{\tau }_2}{\text{d}{\tau }_{\text{g}}}\left(2+\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{f}_{\text{shield}}\left(\frac{{\tau }_2}{{\sigma }_{\text{g}}}\right){\text{e}}^{-{\tau }_{\text{g}}}\right)+\frac{v_{\text{IF}}{\sigma }_{\text{g}}}{R}\frac{{\text{d}}^2{\tau }_2}{\text{d}{\tau }_{\text{g}}^2}.$$(7)

In Eq. (7), we observe the appearance of two dimensionless quantities representing the physical conditions. The first is ${\mathcal{P}}_0{G}_0/R{n}_{\text{H}}$, which gives the ratio of the H₂ formation timescale to the H₂ dissociation timescale. This parameter (multiplied by the H₂ self-shielding function) matches the αG parameter from Sternberg (1988) and Sternberg et al. (2014). As explained in these papers, it determines the cloud FUV optical depth due to the dust associated with the H I gas. We find a new dimensionless parameter v_IF σ_g∕R, which corresponds to the ratio of the H₂ formation timescale to the timescale necessary to advect a dust column density equivalent to τ_g = 1. The physical conditions of the cloud G₀ and n_H only appear as the G₀ ∕n_H ratio in the first dimensionless parameter, while v_IF appears separately in the second dimensionless parameter. We then present our results as functions of G₀ ∕n_H and of v_IF.

We now explore both the case with a null ionization front velocity, which we call “static case” in the following (Sect. 2.3), and the case with a nonzero ionization front velocity, which we call “propagating ionization front case” in the following abbreviated as “advection case” (Sect. 2.4).

2.3 Static case

We first rederive the static solution found in Sternberg et al. (2014). In this situation, v_IF = 0. Eq. (7) then becomes: $$\frac{\text{d}{\tau }_2}{\text{d}{\tau }_{\text{g}}}=\frac{1}{2+\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\exp \left(-{\tau }_{\text{g}}\right)f_{\text{shield}}\left({\tau }_2/{\sigma }_{\text{g}}\right)}.$$(8)

It is possible to solve analytically Eq. (8) by using the following simplified H₂ self-shielding function presented in Draine & Bertoldi (1996): $$f_{\text{shield}}(N_2)=\{\begin{array}{cc}1& \text{if} N_2<{10}^{14} {\text{cm}}^{-2},\\ {\left(\frac{N_2}{{10}^{14} {\text{cm}}^{-2}}\right)}^{-3/4}& \text{else}.\end{array}$$(9)

The proof of the solution to Eq. (8) with Eq. (9) as the H₂ self-shielding function is presented in Appendix A. We obtained that τ₂ can be written as the following function of τ₁ by using that τ_g = τ₁ + 2τ₂: $${\tau }_2({\tau }_1)=\frac{1}{2}{\Gamma }^{-1}\left[\frac{1}{4},\Gamma \left(\frac{1}{4},{\sigma }_{\text{g}}\times \left({10}^{14} {\text{cm}}^{-2}\right)\right)+(1-{\text{e}}^{{\tau }_1})/\xi \right],$$(10)

with Γ the incomplete Gamma function and ξ a dimensionless parameter defined as: $$\xi =\frac{{\left({\sigma }_{\text{g}}\times {10}^{14} {\text{cm}}^{-2}\right)}^{3/4}}{\sqrt[4]{2}}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}.$$(11)

By rearranging Eq. (10) to have ${\tau }_1\left({\tau }_2\right)$ and taking the limit τ_g → +∞ (in which case τ₂ → +∞), we can then obtain the expression for the total H I column density N_1,tot, also presented in Appendix A. It gives: $$N_{1,\text{tot},\text{static}}=\frac{1}{{\sigma }_{\text{g}}}\ln \left[1+2.97\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{\left({\sigma }_{\text{g}}\times {10}^{14} {\text{cm}}^{-2}\right)}^{3/4}\right].$$(12)

To finally obtain τ₁ at the position of the transition we need to solve the system given by Eq. (A.2) (valid only at the position of the transition) and Eq. (A.13) (valid everywhere). We solved that system numerically to obtain the results shown in Fig. 1. The black solid line represents the total atomic hydrogen column density N_1,tot. We introduce two other quantities, which are the column densities of atomic hydrogen in the atomic and molecular regions. The first one is calculated by integrating the atomic hydrogen density only from the beginning of the cloud (i.e., the ionization front) to the H/H₂ transition defined as the position where n₁ = n₂. The second is taken from the transition to the end of the cloud. The sum of these two quantities gives N_1,tot. This separation allows to compare the contribution to N_1,tot of the atomic fraction remaining after the H/H₂ transition. The black dashed line shows the column density of atomic hydrogen integrated in the atomic region. The black dotted line is the column density of atomic hydrogen integrated over the molecular region. A direct comparison of these three curves shows that for the weak field limit² with G₀∕n_H < 10⁻² cm³, the column density of atomic hydrogen is mainly built up in the molecular region. For the strong field limit with G₀ ∕n_H > 10⁻² cm³, N_1,tot is mainly fromthe atomic region. This feature is very well explained in Sternberg et al. (2014) as a consequence of the shape of the H/H₂ transition. Indeed, on Fig. 2, one can find the abundances of H and H₂ as functions of A_V for two different G₀ ∕n_H ratios. One can see that in the weak-field regime (low G₀∕n_H ratio, black solid lines), the atomic hydrogen is present into the molecular medium in significant proportions. In fact, the integral in both the atomic and molecular regions shows that they roughly contribute to the same amount of atomic hydrogen column density. In this weak field regime, the FUV absorption is dominated by H₂ self-shielding. In the strong-field regime in dashed lines, the atomic medium has a far higher contribution to the atomic hydrogen column density than the molecular region. In such conditions, the dominant source of far-UV absorption is dust extinction.

We also show in Fig. 1 the solution of Sternberg et al. (2014). The small difference between the black line from our study and the red one from Sternberg et al. (2014) comes from the value of ${\mathcal{P}}_0$ used and discussed above.

To test the impact of the second and more accurate expression of the shielding function from Draine & Bertoldi (1996) (Eq. (5)), we solve numerically Eq. (8) and compare the solutions obtained with the two expressions of the self-shielding function in Fig. 3. We see on Fig. 3 that the more accurate function of H₂ self-shielding as the red curves leads to slight differences, mainly concerning the distribution of the atomic hydrogen before and after the transition. Apart from a little factor in the weak-field regime, the total atomic hydrogen column density is very similar with the two H₂ self-shielding functions. We now only use Eq. (5) as the H₂ self-shielding function in the following.

Fig. 1 N1 integrated over the atomic region (black dashed line), N1 integrated over the molecular region (black dotted line), N1,tot (black solid line) as functions of G0∕nH for a static and stationary approach, calculated by solving Eq. (8) analytically with Eq. (9). The red line is N1,tot from Sternberg et al. (2014) calculated with a similar approach.

Fig. 2 Spatial profiles of H and H2 abundances obtained by integrating Eq. (8). The solid lines represent a weak field static solution and the dashed linesrepresent a strong field static solution.

Fig. 3 N1 integrated over the atomic region (dashed line), N1 integrated over the molecular region (dotted line), N1,tot (solid line) as functions of G0∕nH for a static and stationary approach. Black lines are for Eq. (9) of H2 self-shielding, and red lines are for Eq. (5), both from Draine & Bertoldi (1996).

2.4 Propagating ionization front case

To study the case where the ionization front propagates into the cloud (from now on, v_IF ≠0, referred to as the “advection case” in the following), Eq. (7) is solved numerically. As we saw in Eq. (7) the solution does not depend on G₀ and n_H independently but on the G₀ ∕n_H ratio. So, the free parameters that are explored are G₀∕n_H, Z′ (appearing in both R and σ_g expressions) and v_IF. The resolution of this equation gives the spatial profile of τ₂ from which we derive τ₁, and then x₁ and x₂.

We take v_IF as a free parameter, but as explained in Bertoldi & Draine (1996), this quantity is directly linked to the flux of ionizing photons and so depends on the stellar type of the illuminating star. We discuss this point more thoroughly in Sect. 4.1. The results of our numerical computations are presented in Sect. 3.

3 Results

3.1 Advection effect on the transition

The profile of the abundances of H and H₂ have been numerically computed by a direct integration of Eq. (7). Until Sect. 3.4, we assume solar metallicity Z′ = 1. We can see an example of the shift between the static and advection profiles discussed above on Fig. 4. It compares the spatial profiles of atomic hydrogen and molecular hydrogen with and without advection for physical conditions typical of a bright PDR, with a G₀ ∕n_H ratio of 0.5 cm³ and ionization front velocity v_IF = 1 km s⁻¹. In this case, we clearly see that the H/H₂ transition is closer to the ionization front on the left than the static model. The higher the ionization front velocity is, the closer the H/H₂ transition is to the left edge of the cloud. The atomic region is therefore smaller as well as the accumulated atomic hydrogen column density. For nonzero ionization front velocity, the dissociation front differs from a weaker radiation field because the transition issharper. To illustrate that, we can compare the solid lines of Fig. 4 representing the H/H₂ transition for a nonzero velocity to the dashed line in Fig. 2 representing the H/H₂ transition without advection. These two cases have the H/H₂ transition roughly at the same position, but in the advection case the transition is much sharper.

Fig. 4 Spatial profiles of H and H2 abundances obtained by integration of Eq. (7) with typical PDR conditions, i.e., G0 ∕nH = 0.5 cm3 and vIF= 1 km s−1. The dashed lines are the static solution and the solid lines are the advection one.

3.2 Merging fronts

With a high-enough velocity, the atomic region disappears and we observe in our models a merging of the dissociation front with the ionization front. This merging was first predicted by Bertoldi & Draine (1996). When the ionization front propagates too fast in the medium, the dissociation front cannot separate from the ionization front. In these conditions, they form a merged structure, and the atomic region does not exist.

We solved Eq. (7) for a grid of different v_IF and differentG₀∕n_H ratio and determined the critical value of v_IF above which front merging occurs for each value of G₀∕n_H. These computed critical values are presented as red circles in Fig. 5. To the right of the boundary defined by these circles, the fronts are merged. In the following we compare the fronts merging criteria from previous papers to our results and derive a new expression for the boundary.

Bertoldi & Draine (1996) derive their criterion by comparing the velocity of the dissociation front to the velocity of the ionization front assuming they are initially merged. The two fronts can only separate if the initial dissociation front velocity (which depends on G₀∕n_H) is larger than the ionization front velocity. In their derivation of the criterion (seen in Fig. 5 as the dotted line) they neglected H₂ self-shielding, which is why their criterion is a proportionality relation.

Another fronts merging criterion was provided in Störzer & Hollenbach (1998). To derive their criterion, they compared the dynamical timescale of the flow of H₂, representing how long it takes for H₂ to flow across the atomic region with the H₂ photodissociation timescale. If the dynamical timescale is shorter than the photodissociation timescale, the advected H₂ could flow across the PDR without being dissociated and would reach the ionization front, the fronts would be merged. This approach was made without taking into account H₂ formation on grains and dust absorption. They also assumed a fixed H₂ column density value when taking H₂ self-shielding into account (using Eq. (9)). In Fig. 5, the black dashed line represents Störzer & Hollenbach front merging criterion.

As we can see, both of the previously discussed fronts merging criteria are close in term of magnitude, but neither of them can explain the shape of the red circles with the smooth transition between two straight lines, corresponding to two regimes. By not taking into account the H₂ self-shielding, Bertoldi & Draine (1996) are closer to the regime at high values of v_IF, whereas Störzer & Hollenbach (1998), by taking a fixed value of the H₂ self-shielding, are closer to the regime at low-values of v_IF. One can guess that the two regimes and the transition between the two come from the treatment of H₂ self-shielding. Indeed, for low values of the G₀∕n_H ratio, Sternberg et al. (2014) explained that H₂ self-shielding is the dominant source of FUV absorption, so that advection has a stronger effect in this case. A lower value of v_IF is needed to merge the fronts. We found a closer criterion (the black solid line) by finding an approximation of the slope of x₂ at the position of the transition and at the merging of the fronts. We can rewrite Eq. (7) at the position of the transition where x₁ = x₂ = 1∕3 and when the fronts are merged, namely when τ_g = τ₂ = 0. With these assumptions, H₂ and dust shieldings are equal to 1, leading to: $$v_{\text{IF}}=\left(\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}-1\right)\times \frac{1}{3}\frac{R}{{\sigma }_{\text{g}}}{\left(\frac{\text{d}{x}_2}{\text{d}{\tau }_{\text{g}}}\right)}^{-1}.$$(13)

From here we need to find the x₂ derivative by isolating dx₂∕dτ_g in the previous equation and by taking the merging conditions found in the semi-analytical models. The fit gives us: $$\frac{\text{d}{x}_2}{\text{d}{\tau }_{\text{g}}}=\frac{1}{3}\frac{1}{4.5}\frac{\left(\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}-1\right)}{\exp \left({\tau }_{\text{g}}^0\right)-1},$$(14)

with ${\tau }_{\text{g}}^0$ the depth of the H/H₂ transition in the static case for this G₀∕n_H ratio. Then we can deduce that: $$v_{\text{IF}}=\frac{4.5R}{{\sigma }_{\text{g}}}\left[\exp \left({\tau }_{\text{g}}^0\right)-1\right],$$(15)

for which we need to numerically compute ${\tau }_{\text{g}}^0$ first. As it then depends on computing a semi-analytical model, we do not have an analytical expression for the merging criteria.

To get a more convenient expression for our criterion ${\left(G_0/n_{\text{H}}\right)}_{\text{crit}}$ as a function of v_IF, we fit the merging criterion from Eq. (15). The fit is presented as the blue dashed line of Fig. 5, and is given by: $$\begin{array}{lll}{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}(v_{\text{IF}})\hfill & =\hfill & 3.0\times {10}^{-2}\times \left(\frac{v_{\text{IF}}}{1.0 \text{km}{\text{s}}^{-1}}\right)\hfill \\ \hfill & \hfill & +4.5\times {10}^{-4}\times \ln \left(1+\frac{v_{\text{IF}}}{1.0\times {10}^{-3} \text{km}{\text{s}}^{-1}}\right).\hfill \end{array}$$(16)

Fig. 5 Criterion describing, for a given ionization front velocity, the G0∕nH ratio leading to a merge of the ionization and dissociation fronts. The ionization front and the dissociation front are merged to the right of the criteria. The black dotted line is the criterion derived by Bertoldi & Draine (1996), the black dashed line is the one of Störzer & Hollenbach (1998). The red circles are the effective drop from Fig. 6. The black line is our new semi-analytical criterion from Eq. (15) and the blue dashed line is a fit of the criterion from Eq. (16).

3.3 Total atomic hydrogen column density

In Fig. 6, we can see the total atomic hydrogen column density N_1,tot in solid lines, the atomic region contribution to N_1,tot in dashed lines and the molecular region contribution in dotted lines for different velocities ranging from v_IF = 0 to 1 km s⁻¹ as functions of the G₀∕n_H ratio. We can make a few observations. First and as expected, shifting the transition even slightly toward the ionization front makes N_1,tot decrease from the static case. A smaller shift from a smaller velocity leads to a fainter deviation from the static case. Secondly, we see that for the advection cases, as the G₀ ∕n_H ratio decreases, the atomic region contribution to N_1,tot (in dashed lines) decreases, but the drop is much more abrupt in the advection cases than in the static one. This is the translation of the merging of the fronts discussed above (Sect. 3.2). As the H/H₂ transition approaches and merges with the ionization front, N₁ integrated in the atomic region drops to 0. For G₀∕n_H ratio under the cutoff, the atomic region does not exist. However, N_1,tot is not 0, because a fraction of atomic hydrogen can still be found in the molecular region, although not being the dominant specie.

We can notice that we retrieve the two regimes discussed in Sect. 2.3 coming from the competition of the two components of N_1,tot. The first one is for high G₀∕n_H ratio, where N_1,tot is dominated by the atomic hydrogen accumulated in the atomic region. The other is for low G₀ ∕n_H, where N_1,tot converges toward a power law representing the atomic hydrogen integrated in the molecular region. The switch between the two regimes follows the merging criteria discussed earlier, so it highly depends on the velocity. As seen on Fig. 5, higher velocity translates in a switch happening at higher G₀ ∕n_H ratio.

In the top panel of Fig. 7, we see ΔN_1,tot = N_1,tot,static − N_{1,tot,dynamical} as the black solid line for v_IF = 0.1 km s⁻¹. The column density N_1,tot,static is in black dashed line for comparison. As we can see, the two regimes discussed above are clearly visible. The change in the total atomic hydrogen column density ΔN_1,tot is proportional to the G₀ ∕n_H ratio as long as the fronts are merged, and when they are dissociated, ΔN_1,tot is constant. We highlight this behavior with the superposition of a constant function (red dashed line) to the dissociated fronts regime and a linear function (blue dashed line) to the merged fronts regime.

On the middle panel of Fig. 7, we show as red circles the value of ΔN_1,tot in the dissociated fronts regime as a function of the ionization front velocity v_IF. The models were computed for a very high G₀∕n_H ratio (100 cm³). The black line is a formula only depending on v_IF. On the bottom panel of Fig. 7, the blue circles give the ΔN_1,tot∕N_1,tot,static ratio in the merged fronts regime (where ΔN_1,tot∕N_1,tot,static is roughly constant) as a function of v_IF. Here these values were computed with a very low G₀∕n_H (10⁻⁵ cm³). The black line is a fit made on the blue circles as well as low-metallicity results as introduced and explained later in Sect. 3.4. Looking at Fig. 5, we see that for such a low G₀ ∕n_H ratio, every velocities investigated lead to merged fronts. But if v_IF is taken small enough, the fronts start to dissociate so the fit does not apply in that case. That is what we start to see with the point where v_IF = 10⁻³ km s⁻¹: our fit does not account very well for very low velocities because the fit is built to replicate strongly merged fronts conditions. Finally, we present the value of ΔN_1,tot in each of the two regimes (separated fronts and merged fronts) in the following Eq. (17): $$\begin{array}{lll}\hfill & \hfill & \Delta N_{1,\text{tot}}\left[{\text{cm}}^{-2}\right]=N_{1,\text{tot},\text{static}}-N_{1,\text{tot},\text{dynamical}}\hfill \\ \hfill & \hfill & =\{\begin{array}{l}\frac{1}{1.5{\sigma }_{\text{g}}}\ln \left(1+\frac{v_{\text{IF}}{\sigma }_{\text{g}}}{2R}\right) \text{for} \frac{G_0}{n_{\text{H}}}>{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}\left(v_{\text{IF}}\right),\hfill \\ \hfill \\ \frac{N_{1,\text{tot},\text{static}}}{\pi /2}\arctan \left(\sqrt{\frac{v_{\text{IF}}}{0.05 \text{km}{\text{s}}^{-1}}}\right) \text{for} \frac{G_0}{n_{\text{H}}}<{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}\left(v_{\text{IF}}\right).\hfill \end{array}\hfill \end{array}$$(17)

Fig. 6 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines) and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio. The black lines are the static result. The blue and red lines are the advection results respectively for vIF = 0.1 and 1.0 km s−1.

Fig. 7 Top panel: ΔN1,tot (black solid line) as the difference between static (black dashed line for comparison) and advection cases as a function of the G0 ∕nH ratio for vIF = 0.1 km s−1. The blue and red dashed lines are fits of ΔN1,tot in the weak and strong-field limit, respectively. Middle panel: ΔN1,tot between static and advection cases as a function of vIF. The points were obtained for a strong field of G0∕nH = 102. The black lines are the formula presented in Eq. (17). Bottom panel: same as middle panel but for a weak field of G0∕nH = 10−5. In this limit, Eq. (17) is a fit made on all the points, including metallicity, discussed in Sect. 3.4.

Fig. 8 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines) and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio for, from top to bottom, respectively, Z′ = 0.01 in blue, 0.1 in red and 1 in black, obtained for vIF = 0.1 km s−1. The dotted-dashed lines are the static results.

3.4 Metallicity

From Eqs. (1) and (4), we can see that by decreasing the metallicity Z′, R and σ_g are also decreased because of their proportionality to Z′. When both of these equations are injected in Eq. (7), we first see that a smaller R increases our dimensionless parameter ${\mathcal{P}}_0{G}_0/R{n}_{\text{H}}$, it is then equivalent to a larger G₀∕n_H ratio by a factor 1∕Z′. Metallicity does not affect the advection term because the two Z′ dependencies cancel each other in the R∕σ_g ratio. Because the ${\mathcal{P}}_0{G}_0/R{n}_{\text{H}}$ ratio depends on Z′ but not the advection, the advection needs to be stronger to merge the fronts, causing the merging criterion from Fig. 5 to be shifted to the right for a lower metallicity.

Figure 8, represents, for v_IF = 0.1 km s⁻¹, the atomic hydrogen column density integrated in the atomic region (dashed lines), the molecular region (dotted lines), and the total (solid lines), along with the static models (dotted-dashed lines), for different metallicities Z′ = 0.01, 0.1, 1, respectively in blue, red, and black. The column density N_1,tot is larger as Z′ decreases. Indeed, the parameter Z′ acts like a scaling factor because σ_g intervenes in the column density N₁ = σ_gτ₁.

Finally, a decreased metallicity has an impact on H₂ self-shielding, as the self-shielding term is f_shield(τ₂∕σ_g) and σ_g decreases as Z′. As a consequence, H₂ self-shielding is more important relative to dust-shielding. This effect changes the shape and position of the curvature seen for Z′ = 1 in the merging criterion from Fig. 5.

We computed the merging criteria for the three metallicities, Z′ = 0.01, 0.1, 1 and presented them in Fig. 9. The merging values obtained with our semi-analytical model are the red circles. We see that the effect is not only a scaling factor but is a bit more complex, yet is very well reproduced by our analytical formula from Eq. (15) as the black lines. We can see that merging the fronts requires a stronger ionization front velocity as metallicity decreases. The fit presented in Eq. (16) does not take into account the metallicity, but can still be used by adding a dependence to metallicity as follows: $$\begin{array}{lll}{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}(v_{\text{IF}},Z^{\prime })\hfill & =\hfill & 1.5\times {10}^{-2}\left(1+Z^{\prime }\right)\left(\frac{v_{\text{IF}}}{1.0 \text{km}{\text{s}}^{-1}}\right)\hfill \\ \hfill & \hfill & +4.5\times {10}^{-4}{Z^{\prime }}^{2/3}\ln \left(1+\frac{v_{\text{IF}}/\sqrt{Z^{\prime }}}{1.0\times {10}^{-3} \text{km}{\text{s}}^{-1}}\right).\hfill \end{array}$$(18)

This fit can be seen in Fig. 9 as the blue dashed line.

The two regimes dichotomy discussed previously for the Z′ = 1 case is still applicable to the other cases, as seen on the top panel of Fig. 10, with the merged or dissociated fronts regimes. We clearly observe the two regimes with a constant value for dissociated fronts and a power law for merged fronts. On the middle panel, we see ΔN_1,tot for a few metallicities. Our formula from Eq. (17) fits very well the shift in this dissociated fronts regime, reproducing the scaling factor. On the bottom panel, the semi-analytical results were fit to obtain the black curve. We present hereafter the final formula for both of the regimes and including metallicity: $$\begin{array}{lll}\hfill & \hfill & \Delta N_{1,\text{tot}}\left[{\text{cm}}^{-2}\right]=N_{1,\text{tot},\text{static}}-N_{1,\text{tot},\text{dynamical}}\hfill \\ \hfill & \hfill & =\{\begin{array}{l}\frac{1}{1.5{\sigma }_{\text{g}}}\ln \left(1+\frac{v_{\text{IF}}{\sigma }_{\text{g}}}{2R}\right)\text{for}\frac{G_0}{n_{\text{H}}}>{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}\left(v_{\text{IF}},Z^{\prime }\right),\hfill \\ \hfill \\ \frac{N_{1,\text{tot},\text{static}}}{\pi /2}\arctan \left(\sqrt{\frac{v_{\text{IF}}/\sqrt{Z^{\prime }}}{0.05 \text{km}{\text{s}}^{-1}}}\right)\text{for}\frac{G_0}{n_{\text{H}}}<{\left(\frac{G_0}{n_{\text{H}}}\right)}_{\text{crit}}\left(v_{\text{IF}},Z^{\prime }\right).\hfill \end{array}\hfill \end{array}$$(19)

In the dissociated fronts regime, where $G_0/n_{\text{H}}>{\left(G_0/n_{\text{H}}\right)}_{\text{crit}}\left(v_{\text{IF}},Z^{\prime }\right)$ and where ΔN_1,tot is constant with respect to G₀∕n_H, the effect of metallicity is just a scaling factor that affects our fit through the dependence of σ_g on Z′. In the merged fronts regime ($G_0/n_{\text{H}}<{\left(G_0/n_{\text{H}}\right)}_{\text{crit}}\left(v_{\text{IF}},Z^{\prime }\right)$), the effect or metallicity can be described as a shift to the left of the fit as Z′ decreases, modeled by the $\sqrt{Z^{\prime }}$ factor in the fit.

Fig. 9 Criterion describing, for a given ionization front velocity, the G0∕nH ratio leading to a merge of the ionization and dissociation fronts. The red circles are the empirical merging values. The black lines are our new semi-analytical criterion from Eq. (15) and the blue dashed lines are the fits from Eq. (18). The ionization front and the dissociation front are merged to the right of the criteria.

Fig. 10 Top panel: ΔN1,tot (solid lines) as a function of the G0∕nH ratio for vIF = 0.1 km s−1 and for Z′ = 0.01 (blue), Z′ = 0.1 (red), and Z′ = 1 (black). The static N1,tot are the dashed lines. Middle panel: ΔN1,tot between static and advection cases as a function of vIF. The points were obtained for a strong field of G0∕nH = 102 cm3. The black lines are the analytical formula presented in Eq. (19). Bottom panel: Same as middle panel but for a weak field of G0∕nH = 10−5 cm3. In this limit, Eq. (19) is a fit made on all the points.

4 Discussion

4.1 The link between G₀/n_H and v_IF

As explained in Sect. 2.4, we treated v_IF as a free parameter independent of G₀. We investigate in this section how v_IF is linked to our G₀∕n_H ratio. First, mass conservation across the ionization front (Eq. (37.1), Draine 2011b) yields: $$v_{\text{IF}}=\frac{F_{\text{EUV}}({\text{H}}^{+}/\text{H})}{n_{\text{H}}},$$(20)

with F_EUV(H⁺∕H) the flux of ionizing photons arriving at the ionization front. This flux depends on the absorption of EUV photons in the ionized region, on the distance between the star and the PDR and on the spectral type of the star. The G₀ parameter similarly depends on the distance of the star, the spectral type of the star and the absorption of FUV photons in the ionized region.

We note that G₀∕n_H and v_IF are determined by taking the EUV and FUV fluxes at the ionization front. The star emission rate of EUV and FUV photons are respectively S_EUV and S_FUV (photons s⁻¹) and the distance between the star and the ionization front is d_star. Following Bertoldi & Draine (1996), we call q_EUV the factor representing EUV absorption by dust and hydrogen recombination in the ionized gas, and q_FUV the factor representing FUV dust absorption in the ionized gas. The ionization front velocity is then expressed as (Bertoldi & Draine 1996): $$v_{\text{IF}}=\frac{S_{\text{EUV}}}{4\pi d_{star}^2{q}_{\text{EUV}}{n}_{\text{H}}},$$(21)

and, for a radiation field coming from a point source (star), the G₀ ∕n_H ratio is³ : $$\frac{G_0\times F_{\text{ISRF}}}{n_{\text{H}}}=\frac{S_{\text{FUV}}}{4\pi d_{star}^2{q}_{\text{FUV}}{n}_{\text{H}}},$$(22)

with F_ISRF = 1.27 × 10⁷cm⁻² s⁻¹. Dividing the two expressions and using Eq. (12) from Bertoldi & Draine (1996): $$\frac{v_{\text{IF}}}{G_0/n_{\text{H}}}=\frac{S_{\text{EUV}}}{S_{\text{FUV}}}\times \frac{\exp \left(({\sigma }_{\text{g}}-{\sigma }_{\text{EUV}})N_{\text{HII}}\right)}{1+\frac{5.07\times {10}^{-20}\times N_{\text{HII}}}{1+{\sigma }_{\text{EUV}}{N}_{\text{HII}}/3}}\times F_{\text{ISRF}},$$(23)

where N_HII is the total gas column density integrated in the ionized region and σ_EUV the EUV dust extinction cross section. In this expression, the link between the ionization front velocity and G₀ ∕n_H is quite clear, depending only on the properties of the star and the absorption by dust and gas before the ionization front.

As stated in Bertoldi & Draine (1996), the result depends a lot on the dust absorption properties. In the following, we use σ_g from Eq. (4). For σ_EUV, the EUV dust extinction, we take the value of Bertoldi & Draine (1996) for R_V = 3.1: $${\sigma }_{\text{EUV}}=1.64\times {10}^{-21} {\text{cm}}^2.$$(24)

As we see in Sect. 4.4, for a various set of objects, τ_EUV = σ_EUVN_HII is close to 1 in a range going from about 0.1 to 10. We decided to take a look at the link between v_IF and the G₀∕n_H ratio for this range of EUV dust extinction. As the link also depends on the ratio between the EUV and FUV emission rates, we first assume S_EUV∕S_FUV = 1, corresponding to an early-type O star.

Figure 11 represents the link between v_IF and the G₀ ∕n_H ratio for different τ_EUV, 0.1 in dotted line, 1 in dashed line, and 10 in dotted-dashed line. This range is compared to the merging criterion in solid line. We see that for a low value of τ_EUV the conditions are always very close to front merging. Indeed, a low τ_EUV means that EUV absorption is low, resulting in a fast propagation of the ionization front into the cloud. For higher values of τ_EUV, the fronts are merged for weak fields, with G₀∕n_H ratio under 10⁻² cm³. For stronger fields, the link between v_IF and the G₀∕n_H ratio does not move far away from the merging criterion, with less than a decade of difference. In such conditions, fronts are clearly dissociated, yet dynamical effects are still visible on N_1,tot.

That is what is presented on Fig. 12 for the same range of τ_EUV. We see that ΔN_1,tot is very closeto the static solution, meaning that dynamical effects are very strong. They are higher for low τ_EUV values as we are closer to the front merging conditions, but are still very high as τ_EUV increases. We also observe thatthe effects are more important in the strong field limit, because for such conditions the ionization front velocity is higher, meaning that the H/H₂ transition is sharper, translating into a lower atomic fraction in the molecular region.

As we can see in Eq. (20), the velocity of a steady ionization front is inversely proportional to the gas density immediately behind the IF. In the classical scenario of a D-type front preceded by a shock front, this density is thus not the initial density of the molecular cloud, but the density after shock compression. In addition, density gradients preexisting in the molecular cloud lead either to an accelerating ionization front (propagating in the direction of a decreasing density gradient, as e.g., in the formation of a champagne flow H II region) or a decelerating ionization front (propagating in the direction of an increasing density gradient). As shown by Franco et al. (1989, 1990), respectively for spherically symmetric and disk-like clouds, an ionization front propagating from the inside of the structure is accelerated along the decreasing density ramp, and if the ramp is steeper than r^−3∕2, the D-type ionization front overtakes the shock front and becomes a R-type front. This acceleration originating from the density structure leads to the aforementioned blisters and champagne flows. The stationary hypothesis for the PDR structure only remains valid if the timescale of variation of n_H and v_IF due to the density gradient and the acceleration of the ionization front is larger than the timescale to reach a stationary PDR structure. As n_H and v_IF are linked to each other by Eq. (20), the variation of one of them causes the variation of the other, so that their timescales of variation are equal. The timescale of variation of the density n_H is: $${\tau }_{n_{\text{H}}}=\frac{n_{\text{H}}}{\text{d}{n}_{\text{H}}/\text{d}t}=\frac{n_{\text{H}}}{v_{\text{IF}}\times \text{d}{n}_{\text{H}}/\text{d}x}.$$(25)

By comparing it to the timescale to reach a stationary PDR structure, which is ${(2R{n}_{\text{H}})}^{-1}$ (Hollenbach & Natta 1995), one finds that to maintain a stationary state, the logarithmic derivative of the density has to verify: $$\frac{1}{n_{\text{H}}}\frac{\text{d}{n}_{\text{H}}}{\text{d}x}<\frac{2R{n}_{\text{H}}}{v_{\text{IF}}}.$$(26)

This criterion represents the inverse of the characteristic length L of variation of density and gives us the condition on the logarithmic derivative of the density to maintain a stationary state. It means that the criterion is equivalent to: $$L>5.5\times {10}^{-3}\left(\frac{{10}^5 {\text{cm}}^{-3}}{n_{\text{H}}}\right)\left(\frac{v_{\text{IF}}}{1 \text{km}{\text{s}}^{-1}}\right) \text{pc}.$$(27)

Thus, for dense gas, preexisting density variations on scales of a few mpc would be needed to invalidate our stationary hypothesis. In comparison, dense filaments have been argued to have a typical scale of 0.1 pc (Arzoumanian et al. 2011), so that in a wide range of environments, the acceleration or deceleration of the ionization front is slow enough for the stationary hypothesis to remain reasonable. However, smaller scale preexisting substructures might exist in some environments, in which the acceleration or deceleration of the ionization front would not permit the establishment of a stationary H/H₂ structure. Potential observational signatures of such nonstationary PDRs remain to be explored.

Finally, our study adopts a 1D geometry. In 2D and 3D geometries, other structures can however emerge from instabilities induced by the propagation of the front itself. These instabilities have been proposed to explain the bright rims or elephant trunks structures. The linear analysis of thin-shell instabilities by Giuliani (1979, 1980) were confirmed by numerical simulations by Garcia-Segura & Franco (1996). They also showed the importance of cooling in the neutral gas in the appearance of dynamic instabilities that lead to a rapid shell fragmentation. With the merging of the fragments, dense and massive clumps are created with a variety of shapes. Williams (1999) showed that the clump cast a shadow that creates a tail of nonionized gas behind the clump, leading to further instabilities until the clump vanishes. In Whalen & Norman (2008), 3D radiation hydrodynamics simulations are presented to model the propagation of ionization fronts. Theyalso concluded that the cooling efficiency in the shocked neutral medium plays a key role in the development of instabilities. Strong cooling leads to long clumpy structures whereas weak cooling leads to turbulent flows.

Fig. 11 Comparison between the merging criterion from Fig. 5 and the link between vIF and the G0∕nH ratio from Eq. (23).

Fig. 12 ΔN1,tot range for different τEUV compared with N1,tot from the static solution, using the link between vIF and the G0∕nH ratio from Eq. (23). The closer ΔN1,tot is to N1,tot, static, the strongerthe dynamical effects.

4.2 D-type ionization fronts and critical IF velocity

While Eq. (20) is necessarily true across a steady ionization front, the classical theory of ionization fronts (Kahn 1954, Draine 2011b, chap. 37) shows that a steady (plane-parallel) ionization front can only exist if the resulting velocity is either lower than v_D or higher than v_R, with $$v_{\text{R}}=c_{\text{II}}+\sqrt{c_{\text{II}}^2-c_{\text{PDR}}^2-\frac{v_{\text{A}}^2}{2}},$$(28)

and $$v_{\text{D}}=c_{\text{II}}-\sqrt{c_{\text{II}}^2-c_{\text{PDR}}^2-\frac{v_{\text{A}}^2}{2}},$$(29)

where $c_{\text{PDR}}=\sqrt{\gamma T/\mu }$ is the isothermal sound speed in the atomic region right after the ionization front with γ the adiabatic index, T the PDR temperature, μ the mean mass per particle, c_II the isothermal sound speed of the H II region and v_A is the Alfvén speed in the atomic region when a magnetic field is present. The Alfvén speed is simply taken as a free parameter here as we do not solve the magneto-hydrodynamics of the ionization front propagation.

When the conditions (F_EUV and n_H) are such that v_IF ≥ v_R, the ionization front is called R-type. A R-type ionization front propagates supersonically relative to the neutral gas, and the gas mainly remains frozen in place as the IF propagates. This corresponds to very early phase of development of an H II region, or to a later stage of acceleration of the front along a decreasing density ramp, as mentioned in Sect. 4.1. On the other hand, when the conditions are such that v_IF ≤ v_D, the ionization front is called D-type. Two types of D-type solutions are possible. Strong D-type fronts correspond to near pressure equilibrium between the ionized and neutral regions as can happen in the late stage of the H II region expansion, approaching the static stage of a pressure-confined H II region. In weak D-type fronts, the ionized gas is accelerated to roughly the ionized gas sound speed, with a corresponding pressure jump from the ionized to the neutral region. The atomic region thus has a higher pressure than the ionized region. When the conditions are such that we should have v_D < v_IF < v_R, no steady ionization front is possible. A shock front develops and precedes the ionization front, increasing the density of the neutral gas immediately before the ionization front so that v_IF = F_EUV∕n_H ≤ v_D.

The values of v_IF that we can consider in our model are thus necessarily lower than v_D or higher than v_R. We consider a typical PDR temperature of 1000 K, giving c_PDR = 3.7 km s⁻¹. The Alfvén speed in the PDR is taken in a wide range of 0 to 5 km s⁻¹, corresponding to a transverse magnetic field of 0 to a few 100 μG. With c_II = 10 km s⁻¹ we deduce that our upper boundary for v_IF in D-type fronts ranges between 0.7 to 1.4 km s⁻¹, while the lower boundary for v_IF in R-type fronts ranges between 18.6 and 19.3 km s⁻¹. Due to their high velocities relative to the neutral gas, R-type fronts can only exist in very short-lived phases and are thus most likely not relevant for our goal of characterizing the H I column density at the surface of dense molecular clouds in star forming regions where photoevaporation can be present. In situations where photoevaporation flows (also called photoabalation flows) are present, as in blister-type H II regions, the fronts are expected to be close to D-critical, with v_IF ≈ v_D (Henney et al. 2005; Henney 2007).

Equations (28) and (29) assume a transverse magnetic field. Ionization fronts with oblique magnetic fields have been studied in Williams et al. (2000) and Williams & Dyson (2001). By rewriting the jump conditions for obliquely magnetized ionization fronts, they determined that a high transverse field shrinks the forbidden range of front velocities. Solutions of R and D-type ionization fronts and transitions between the regimes were investigated both from a mathematical perspective and physical relevance. They concluded that these obliquely magnetized ionization fronts may lead to transverse velocities and modulation of magnetic field, which we do not consider in this study.

4.3 Radiation pressure

Radiation pressure, caused by momentum transferred from the photons that are absorbed and scattered to the gas and dust, can be nonnegligible in H II regions. Krumholz & Matzner (2009) studied the impact of the radiation pressure on the expansion of both embedded and blister H II regions and concluded that radiation pressure could only become important for H II regions that formed around very massive star clusters. However, their study assumed that all of the radiation momentum was given to the neutral shell directly, without considering the impact of the radiation pressure on the ionized gas. This impact on the ionized region has been exhaustively studied for static, pressure-confined H II regions by Draine (2011a), showing that it can cause the apparition of a central cavity with most of the ionized gas compressed into an ionized shell close to the ionization front. The impact of the radiation pressure on the expansion phase of spherical, embedded H II regions has been investigated by Kim et al. (2016), assuming the static stratification found by Draine (2011a) to remain valid during the expansion. They showed that the compression of the ionized gas in an outer ionized shell caused by radiation pressure can significantly affect the expansion of the H II region due to the raised pressure in this ionized shell. Numerical simulations of the expansion of an embedded, spherical H II region including radiative pressure on dust, gas-dust drift, and dust photodestruction (for the conditons of RCW 120) where presented by Akimkin et al. (2017). They found that radiation pressure hasa significant impact on the dynamics of expansion of the H II region in the conditions of their model, and that dust drift can lead to the expulsion of a significant fraction of the dust from the H II region (see also Gustafsson 2018).

While radiation pressure in embedded, spherical H II region has thus been studied significantly in recent years, the impact of radiation pressure on photoevaporation flows (such as found in blister-type H II regions) has however not been investigated in detail. Krumholz & Matzner (2009) consider the case of blister H II region, but as they neglect the impact of radiation pressure on the ionized gas, they implicitly assume the photoevaporation flows to be unaffected by radiation pressure, and radiation pressure and photoevaporation to act independently on the neutral shell. Henney (2007) argue that radiation pressure is not important for typical galactic blister H II regions based on an order of magnitude comparison of the radiative pressure-induced acceleration to the acceleration in a photoevaporation flow. We assume here radiative pressure to be negligible for typical galactic dense PDRsfound at the borders of blister H II regions. In more extreme conditions, radiative pressure could be expected to quench photoevaporation by containing the flow of ionized gas. This could result in a decreased v_IF, and thus a reduced impact on the H/H₂ transition as described by the results of our semi-analytical model (Sect. 2 and 3).

4.4 Comparison with observations

Most observations of PDRs illuminated by B stars show that the H II region is very close to the star and far from the dense PDR, so that photoevaporation is not relevant in that case. Furthermore Diaz-Miller et al. (1998) indicated that in B stars or lower mass cases, the H II region is expected to be far closer to the star than the dissociation front. As the ionization and dissociation fronts are independent in this configuration, the semi-analytical model developed in this paper is no longer relevant. Still, if we do apply our semi-analytical model to the parameters of well-known B-star PDRs (NGC7023,...), it indeed yields a very small ionization front velocity, on the order of less than 10⁻⁴ km s⁻¹ and a G₀ ∕n_H ratio around 10⁻³–10⁻² cm³, leading to negligible effects of the dynamics. In the following comparison, we limit ourselves to PDRs illuminated by O stars.

We present in Table 1 the parameters relevant for our study to confront our conclusions to actual PDRs illuminated by O stars. For these PDRs, we deduced the column density between the star and the ionization front by multiplying the density estimated in the H II region n_HII by the distance between the star and the PDR or the width of the region if mentioned. The ratio between the star emission rates of EUV and FUV photons are determined by making a log-linear interpolation on Table 1 from Diaz-Miller et al. (1998). For Carina, there are a lot of stars illuminating the PDR, with a maximum temperature of 43 900 K (Wu et al. 2018). We took the approximation of a ratio equal to 1 as expected for such a population of stars. Knowing τ_EUV and S_EUV∕S_FUV for these objects, we use Eq. (23) to calculate v_IF. The result of the computation is given as the last column from Table 1 and on Fig. 13 with error bars. We can see that for all these objects, the dynamic effects have to be important, as they are close to the merging criterion. For the Horsehead and California, the error bars even spread into the merged fronts space.

The PDRs seems to form two groups in the v_IF versus G₀∕n_H space, although whether this distinction is an artifact of our small sample or a physical effect remains to be confirmed. For the high-excitation PDRs (high G₀∕n_H) in our sample (Orion Bar, Carina, S 106, and IRAS 23133+6050), our estimated values of v_IF seem to show a clear saturation of the ionization front velocity at the expected value for D-critical fronts (indicated by the dashed vertical line), as predicted by ionization front theory (see Sect. 4.2). Such D-critical fronts are indeed expected at the dense walls of blister H II regions (Henney et al. 2005; Henney 2007). This group of PDRs clearly falls in the nonmerged fronts regime, which is compatible with ALMA observations of the Orion Bar showing a relatively extended atomic region separating the ionization front and the H/H₂ transition (Goicoechea et al. 2016). As found in Sect. 3.1, we nevertheless predict moderate but non negligible effects of the advection dynamics on the H/H₂ transition. A second group at low G₀∕n_H ratio is composed of the Horsehead nebula, California, and M17. The ionization front velocity estimates for these PDRs seems to be close to the merging boundary, although the large error bars do not allow to conclude whether they fall in the merged fronts regime or inthe nonmerged front regime. The velocities found correspond to noncritical weak D-type ionization fronts. If the line-up of these PDRs on the merging boundary is confirmed by more accurate estimations, it might indicate a specific dynamics of merged fronts causing a saturation at the critical merging velocity. Based on these results, we thus expect strong effects of the advection dynamics with merged, or close to merged, ionization front and dissociation front in these low excitation PDRs illuminated by O stars. Strong effect on the emission lines of H₂ could in particular be expected in these low G₀∕n_H PDRs, in contrast to the conclusions of Störzer & Hollenbach (1998) who focused on high G₀ ∕n_H conditions representative of Orion Bar-like PDRs.

Fig. 13 Comparison between the merging criterion from Fig. 5 and deduced parameters G0 ∕nH and vIF for observations. The vertical bar corresponds to the typical value vD = 1 km s−1 deduced from Eq. (29).

4.5 Comparison to Hydra PDR Code simulations

In order to compare our N_1,tot semi-analytical predictions to Hydra PDR Code simulations, we implemented an advection of the cloud through a punctual ionization front. We plot in Fig. 14 the predictions of N_1,tot and N₁ integrated over the atomic and molecular regions for a number of simulations as circles. The agreement is very good with a maximum deviation of a factor 1.2 for every point but one. The latter is for N₁ integrated over the atomic region for v_IF = 1 km s⁻¹ and G₀∕n_H = 0.04. It has a deviation of 2, which can easily be explained by the lack of resolution of the Hydra PDR Code simulation. Because we examine a case very close to front merging, the resolution between the ionization and photodissociation fronts may be insufficient for further precision, yet it is in excellent agreement with the semi-analytical prediction.

5 Conclusions

To sum up, we built a semi-analytical model of the H/H₂ transition ina 1D plane-parallel PDR illuminated by FUV radiation with advancing ionization front resulting from the photoevaporation mechanism. We investigated the stationary equilibrium state in the ionization front reference frame where its propagation translates into advection of the neutral gas through the front. We did not solve the hydrodynamics of the gas through the front, but included the advection term to account for the impact of the steady flow. The advection positively adds up to H₂ formation on dust and its photodissociation by FUV photons. Our radiative transfer takes into account not only the extinction by dust but also the H₂ self-shielding.

We provide a simple set of analytical expressions allowing to compute the total H I column density N_1,tot in the presence of photoevaporation: Knowing G₀, n_H, the metallicity Z′, the stellar type of the illuminating star and N_HII the total proton column density in the H II region, one can first estimate the velocity of the ionization front v_IF from Eq. (23). One then has to find whether the fronts are merged or not by comparing the G₀ ∕n_H ratio to the critical value found in Eq. (18). Equation (12) then gives the static prediction value of N_1,tot, and Eq. (19) finally gives the desired dynamical prediction of N_1,tot.

We present a summary of our results:

• The advection of the gas through the ionization front induces a change in the location of the H/H₂ transition, which is closer to the ionization front compared to the static case. This effect reduces the width of the atomic region;

• If the ionization front velocity exceeds a certain value, the atomic region disappears. The photodissociation front and the ionization front then form a merged structure. The merging of the fronts has been studied in previous papers (Bertoldi & Draine 1996; Störzer & Hollenbach 1998). We provide a new, more accurate criterion to discriminate whether or not the fronts are merged;

• The total atomic hydrogen column density N_1,tot integrated in the PDR is decreased in the presence of a propagating ionization front. When the fronts are not merged, the absolute difference in N_1,tot between the static case and the case of a propagating ionization front depends only on the velocity of the ionization front. When the fronts are merged, it is the relative difference that depends only on the velocity of the ionization front. We provide analytical approximations in both cases;

• The metallicity plays a key role in this study. A lower metallicity changes the merging criterion. Indeed, in such a case, a higher ionization front velocity is needed to merge the fronts. It also changes the magnitude of the total atomic hydrogen column density. The two-regimes dichotomy is still valid with different metallicities;

• Observations of PDRs illuminated by O-stars are close to the merging criteria, where the dynamical effects are strong. For this kind of objects, we expect N_1,tot to be much lower in the advection approach than with a static one, especially for low excitation PDRs.

By taking into account the photoevaporation, we found that the H/H₂ transition is closer to the ionization front than modeled by previous static studies. Some H₂ is then found in a hotter region, richer in FUV photons. One can then expect H₂ ro-vibrational lines to be more intense as collisions and UV-pumping are more efficient. As a lot of H₂ ro-vibrational lines will be observed by the upcoming JWST, taking into account the dynamical effects highlighted in this paper will be essential for the modeling of observations. The effects are expected to cause larger differences for low excitation PDRs, such as the Horsehead Nebulae.

Fig. 14 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines), and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio. The different colors (black, blue, and red) represent vIF = 0, 0.1 and 1 km s−1. The different points are obtained with the Hydra PDR Code.

Acknowledgements

We thank the anonymous referee for helpful comments and suggestions. This work was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Solving Eq. (8) with Eq. (9) as H₂ self-shielding function

Let us find the position of the H/H₂ transition defined as the po- sition where x₁ = x₂, and, as a reminder, x_1,2 = d τ_1,2∕dτ_g.

Writing $x_1^0$, $x_2^0$, ${\tau }_1^0$, ${\tau }_2^0$, ${\tau }_g^0$ the values of abundances and optical depths at the position of the transition and by using the definition of the transition and $1=x_1^0+2x_2^0$, we can deduce that $x_2^0=1/3$. We can then rewrite Eq. (8) at the position of the transition as follows: $$\frac{1}{3}=\frac{1}{2+\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\exp \left(-{\tau }_{\text{g}}^0\right)f_{\text{shield}}({\tau }_2^0/{\sigma }_{\text{g}})}.$$(A.1) $${\tau }_{\text{g}}^0={\tau }_1^0+2{\tau }_2^0=\ln \left[\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{f}_{\text{shield}}({\tau }_2^0/{\sigma }_{\text{g}})\right].$$(A.2)

We also rewrite Eq. (8) at any position (and not just at the transition): $$\frac{d{\tau }_{\text{g}}}{d{\tau }_2}=2+\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\exp \left(-{\tau }_{\text{g}}\right)f_{\text{shield}}({\tau }_2/{\sigma }_{\text{g}}).$$(A.3)

Using that: $$\frac{d{\tau }_{\text{g}}}{d{\tau }_2}=\frac{d\left({\tau }_1+2{\tau }_2\right)}{d{\tau }_2}=\frac{d{\tau }_1}{d{\tau }_2}+2,$$(A.4)

we deduce that : $$\frac{d{\tau }_1}{d{\tau }_2}=\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\exp \left(-{\tau }_{\text{g}}\right)f_{\text{shield}}({\tau }_2/{\sigma }_{\text{g}}).$$(A.5) $$\frac{d{\tau }_2}{d{\tau }_1}={\left(\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{f}_{\text{shield}}({\tau }_2/{\sigma }_{\text{g}}){\text{e}}^{-({\tau }_1+2{\tau }_2)}\right)}^{-1}.$$(A.6)

To solve this, we use Eq. (9) as the H₂ self-shielding func- tion. This function is composed of two sections. One for low values of H₂ column densities with N₂ < 10¹⁴ cm⁻², and one for N₂ > 10¹⁴ cm⁻². Solving Eq. (A.6) consequently requires to first solve it for the first case, giving us a solution from τ₂ = 0 to ${\tau }_2={\tau }_2^{\text{shield}}={\sigma }_{\text{g}}\times \left({10}^{14} {\text{cm}}^{-2}\right)=1.9\times {10}^{-7}$. After that, we solve it for the second case, from ${\tau }_2={\tau }_2^{\text{shield}}$ to τ₂ → +∞. For the first case where the H₂ self-shielding is equal to 1, Eq. (A.6) can be rewritten as: $$\frac{d{\tau }_1}{d{\tau }_2}=\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{\text{e}}^{-({\tau }_1+2{\tau }_2)},$$(A.7)

with the following solution giving τ₁ as a function of τ₂: $${\tau }_1({\tau }_2)=\ln \left(1+\frac{1}{2}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{\text{e}}^{-2{\tau }_2}\right).$$(A.8)

The boundary where ${\tau }_2={\tau }_2^{\text{shield}}$ gives: $${\tau }_1\left({\tau }_2^{\text{shield}}\right)=\ln \left(1+\frac{1}{2}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{\text{e}}^{-2{\tau }_2^{\text{shield}}}\right).$$(A.9)

In the other case, where the self-shielding of H₂ is a decreasing power-law, Eq. (A.6) is rewritten as: $$\frac{d{\tau }_2}{d{\tau }_1}={\left(\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}{\left(\frac{{\tau }_2}{{\tau }_2^{\text{shield}}}\right)}^{-3/4}{\text{e}}^{-({\tau }_1+2{\tau }_2)}\right)}^{-1},$$(A.10)

with the following solution: $$\begin{array}{lll}{\tau }_2\left({\tau }_1\right)\hfill & =\hfill & \frac{1}{2}{\Gamma }^{-1}[\frac{1}{4},\Gamma \left(\frac{1}{4},{\tau }_2^{\text{shield}}\right)\hfill \\ \hfill & \hfill & +\left(1+\frac{1}{2}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\exp \left(2{\tau }_2^{\text{shield}}\right)-{\text{e}}^{{\tau }_1}\right)/\xi ],\hfill \end{array}$$(A.11)

with Γ the incomplete Gamma function so that $\frac{\partial \Gamma (a,x)}{\partial x}=-x^{a-1}\times {\text{e}}^{-x}$, Γ⁻¹ the inverse incomplete Gamma function and: $$\xi =\frac{{\left({\tau }_2^{\text{shield}}\right)}^{3/4}}{\sqrt[4]{2}}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}.$$(A.12)

The ${\tau }_2^{\text{shield}}$ parameter is very small compared to typical values of τ₁ or τ₂ at the H/H₂ transition. We can thus neglect the exponential term involving ${\tau }_2^{\text{shield}}$ in Eq. A.11, so that: $${\tau }_2\left({\tau }_1\right)=\frac{1}{2}{\Gamma }^{-1}\left[\frac{1}{4},\Gamma \left(\frac{1}{4},{\tau }_2^{\text{shield}}\right)+(1-{\text{e}}^{{\tau }_1})/\xi \right].$$(A.13)

To finally obtain τ₁ at the position of the transition we need to solve the system given by Eq. (A.2) (valid only at the position of the transition) and Eq. (A.13) (valid everywhere). We solved that system numerically to obtain the results shown in Figure 1.

Moreover, as we are interested in the total atomic hydrogen column density N_1,tot, we derive here Eq. (12). Expressing Eq. (A.13) as τ₁ (τ₂): $${\tau }_1({\tau }_2)=\ln \left(1+\xi \left(\Gamma \left(\frac{1}{4},{\tau }_2^{\text{shield}}\right)-\Gamma \left(\frac{1}{4},2{\tau }_2\right)\right)\right),$$(A.14)

from which we can derive N_1,tot by taking the limit τ₂ → +∞. Indeed, the fraction of atomic hydrogen is decreasing as τ_g → +∞, so at infinity the gas is strictly in molecular form, τ₂ or τ_g is the same in this limit. In the end, we obtain: $$\begin{array}{lll}{N}_{1,\text{tot}}\hfill & =\hfill & \frac{1}{{\sigma }_{\text{g}}}\ln \left(1+\xi \Gamma \left(\frac{1}{4},{\tau }_2^{\text{shield}}\right)\right)\hfill \\ \hfill & \hfill & =\frac{1}{{\sigma }_{\text{g}}}\ln \left(1+2.97{\left({\sigma }_{\text{g}}\times {10}^{14} {\text{cm}}^{-2}\right)}^{3/4}\frac{{\mathcal{P}}_0}{R}\frac{G_0}{n_{\text{H}}}\right).\hfill \end{array}$$(A.15)

Appendix B C I extinction

Far-UV photons are susceptible of being absorbed by atomic carbon C I in the continuum, thus possibly contributing to the shielding of H₂. The resulting opacity was neglected in this paper, and we show in this appendix that while extinction by gas phase carbon can be important for the global structure of a PDR, its impact on the H/H₂ transition is negligible.

In the FUV band, the atomic carbon absorbs from 912 Å to 1101 Å with an almost constant absorption cross-section σ_C ≃ 1.7 × 10⁻¹⁷ cm² (Nahar & Pradhan 1997). In this band and at a given position, it diminishes the flux by a factor of $\exp \left(-{\sigma }_{\text{C}}{N}_{\text{C}}\right)$, with N_C the column density of C I integrated from the edge of the cloud to the given position. Adding its effect to the semi-analytical model of this paper would thus require taking into account a more extensive chemical network in order to compute the column density of atomic carbon. The hypothesis maximizing this extinction would be to assume that all the carbon is in atomic form. With an abundance around x(C) = 10⁻⁴, one can rewrite the extinction term as $\exp \left(-{\sigma }_{\text{C}}{N}_{\text{C}}\right)=\exp \left(-{\sigma }_{\text{C}}x(\text{C})N\right)$ with σ_C x(C) ≃ 1.7 × 10⁻²¹ cm², which is very close to dust ex- tinction cross-section σ_g = 1.9 × 10⁻²¹ cm². So, at best, atomic carbon extinction would be equal to dust extinction. Sternberg et al. (2014) showed that the conditions for which dust extinction dominates over H₂ self-shielding are met in the strong field limit. In the weak field limit, where G₀ ∕n_H ≤ 10⁻² cm³, the H/H₂ transition is driven by H₂ self-shielding and not by dust, and atomic carbon extinction can thus be expected to have no impact on the H/H₂ transition. In the strong field limit (G₀ ∕n_H ≥ 10⁻² cm³), its impact depends on the abundance of atomic carbon before the H/H₂ transition.

To investigate this impact, we used the Meudon PDR Code. The Meudon PDR Code (Le Petit et al. 2006) simulates a stationary, plane-parallel PDR, computing self-consistently the radiative transfer both in the continuum and in lines (Goicoechea & Le Bourlot 2007), the chemistry with an exhaustive chemical network and the thermal balance of the gas. We obtained the re- sults presented in Fig. B.1 and B.2 for different G₀ ∕n_H ratio, from 1 to 10⁻⁴ cm³. On Fig. B.1, we have the abundances profiles of H, H₂, C⁺, C, and CO as functions of A_V. Fig. B.2 presents the extinction terms of dust, C, H₂, and the combination of the three as functions of A_V. For dust and carbon, the extinction terms are exp(−σ_g N) and exp(−σ_C x(C) N) respectively and, for H₂, it is equal to the shielding function f_shield(H₂). These extinction terms are not extracted directly from the PDR models because the PDR code computes the sum of extinctions at each wavelength and it is not convenient to separate the various contributions. Instead, since here we just need simple estimations of these extinction terms, we computed them with the above expressions, and for the carbon extinction term, we use the abundance of C computed at each position by the code.

Fig. B.1 Abundance profiles of H in cyan, H2 in blue, C+ in orange, C in red, and CO in violet as functions of AV for models with G0∕nH = 1, 10−2, 10−3, and 10−4 cm3 obtained using Meudon PDR Code.

Fig. B.2 Extinction factor of dust in blue, of C in red, of H2 self-shielding in green, and of all of them combined in black, as functions of AV. C I and H2 column densities are taken from Meudon PDR Code simulations for models with G0 ∕nH = 1, 10−2, 10−3, and 10−4 cm3.

In the high G₀∕n_H regime (see e.g., G₀ ∕n_H = 1 cm³, top panel), atomic carbon is not very abundant at the H/H₂ transition position with x(C) ≃ 10⁻⁷ for G₀∕n_H = 1 cm³, so that extinction due to C I is smaller than dust extinction. When we compare the actual extinction factors in the top panel of Fig. B.2, we see that indeed, C I extinction is negligible until around A_V ≃ 2.3, deeper than the transition, found around A_V = 1. As G₀∕n_H is decreased, the abundance of atomic carbon in the atomic region increases. In the weak field case, dust and atomic carbon extinctions are very close because carbon is almost entirely in neutral form. But as expected for the weak field regime (Sternberg et al. 2014), one can see that from the edge of the cloud to the H/H₂ transition, the column density is far too small for the extinction factor (both by dust and by atomic carbon) to significantly differ from 1. FUV extinction is mainly driven by H₂ self-shielding, with negligible contributions of carbon or dust extinction. In in- termediate cases, like panels 2 and 3, dust slowly takes a leading importance, but extinction due to C I always remains negligible up to the H/H₂ transition. Deeper in the cloud, atomic carbon continuum extinction of course becomes important, and in particular, we see that it cannot beneglected for the transition to CO, but this region is not considered in this paper.

In conclusion, we find atomic carbon extinction to be negligible in all cases for the purpose of the study of the H/H₂ transition. In the case of high G₀∕n_H, the very low fraction of atomic carbon in the atomic region makes it negligible compared to dust ex- tinction. In the case of low G₀∕n_H, where the fraction of atomic carbon in the atomic region increases to significant levels, the H/H₂ transition is dominated by H₂ self-shielding and appears at a very low column density, so that extinction by both dust and atomic carbon is negligible.

References

All Tables

All Figures

	Fig. 1 N1 integrated over the atomic region (black dashed line), N1 integrated over the molecular region (black dotted line), N1,tot (black solid line) as functions of G0∕nH for a static and stationary approach, calculated by solving Eq. (8) analytically with Eq. (9). The red line is N1,tot from Sternberg et al. (2014) calculated with a similar approach.
In the text

	Fig. 2 Spatial profiles of H and H2 abundances obtained by integrating Eq. (8). The solid lines represent a weak field static solution and the dashed linesrepresent a strong field static solution.
In the text

	Fig. 3 N1 integrated over the atomic region (dashed line), N1 integrated over the molecular region (dotted line), N1,tot (solid line) as functions of G0∕nH for a static and stationary approach. Black lines are for Eq. (9) of H2 self-shielding, and red lines are for Eq. (5), both from Draine & Bertoldi (1996).
In the text

	Fig. 4 Spatial profiles of H and H2 abundances obtained by integration of Eq. (7) with typical PDR conditions, i.e., G0 ∕nH = 0.5 cm3 and vIF= 1 km s−1. The dashed lines are the static solution and the solid lines are the advection one.
In the text

Fig. 5 Criterion describing, for a given ionization front velocity, the G0∕nH ratio leading to a merge of the ionization and dissociation fronts. The ionization front and the dissociation front are merged to the right of the criteria. The black dotted line is the criterion derived by Bertoldi & Draine (1996), the black dashed line is the one of Störzer & Hollenbach (1998). The red circles are the effective drop from Fig. 6. The black line is our new semi-analytical criterion from Eq. (15) and the blue dashed line is a fit of the criterion from Eq. (16).

In the text

	Fig. 6 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines) and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio. The black lines are the static result. The blue and red lines are the advection results respectively for vIF = 0.1 and 1.0 km s−1.
In the text

Fig. 7 Top panel: ΔN1,tot (black solid line) as the difference between static (black dashed line for comparison) and advection cases as a function of the G0 ∕nH ratio for vIF = 0.1 km s−1. The blue and red dashed lines are fits of ΔN1,tot in the weak and strong-field limit, respectively. Middle panel: ΔN1,tot between static and advection cases as a function of vIF. The points were obtained for a strong field of G0∕nH = 102. The black lines are the formula presented in Eq. (17). Bottom panel: same as middle panel but for a weak field of G0∕nH = 10−5. In this limit, Eq. (17) is a fit made on all the points, including metallicity, discussed in Sect. 3.4.

In the text

	Fig. 8 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines) and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio for, from top to bottom, respectively, Z′ = 0.01 in blue, 0.1 in red and 1 in black, obtained for vIF = 0.1 km s−1. The dotted-dashed lines are the static results.
In the text

Fig. 9 Criterion describing, for a given ionization front velocity, the G0∕nH ratio leading to a merge of the ionization and dissociation fronts. The red circles are the empirical merging values. The black lines are our new semi-analytical criterion from Eq. (15) and the blue dashed lines are the fits from Eq. (18). The ionization front and the dissociation front are merged to the right of the criteria.

In the text

Fig. 10 Top panel: ΔN1,tot (solid lines) as a function of the G0∕nH ratio for vIF = 0.1 km s−1 and for Z′ = 0.01 (blue), Z′ = 0.1 (red), and Z′ = 1 (black). The static N1,tot are the dashed lines. Middle panel: ΔN1,tot between static and advection cases as a function of vIF. The points were obtained for a strong field of G0∕nH = 102 cm3. The black lines are the analytical formula presented in Eq. (19). Bottom panel: Same as middle panel but for a weak field of G0∕nH = 10−5 cm3. In this limit, Eq. (19) is a fit made on all the points.

In the text

	Fig. 11 Comparison between the merging criterion from Fig. 5 and the link between vIF and the G0∕nH ratio from Eq. (23).
In the text

	Fig. 12 ΔN1,tot range for different τEUV compared with N1,tot from the static solution, using the link between vIF and the G0∕nH ratio from Eq. (23). The closer ΔN1,tot is to N1,tot, static, the strongerthe dynamical effects.
In the text

	Fig. 13 Comparison between the merging criterion from Fig. 5 and deduced parameters G0 ∕nH and vIF for observations. The vertical bar corresponds to the typical value vD = 1 km s−1 deduced from Eq. (29).
In the text

	Fig. 14 N1 integrated over the atomic region (dashed lines), over the molecular region (dotted lines), and over the whole PDR, N1,tot (solid lines) as functions of the G0∕nH ratio. The different colors (black, blue, and red) represent vIF = 0, 0.1 and 1 km s−1. The different points are obtained with the Hydra PDR Code.
In the text

	Fig. B.1 Abundance profiles of H in cyan, H2 in blue, C+ in orange, C in red, and CO in violet as functions of AV for models with G0∕nH = 1, 10−2, 10−3, and 10−4 cm3 obtained using Meudon PDR Code.
In the text

	Fig. B.2 Extinction factor of dust in blue, of C in red, of H2 self-shielding in green, and of all of them combined in black, as functions of AV. C I and H2 column densities are taken from Meudon PDR Code simulations for models with G0 ∕nH = 1, 10−2, 10−3, and 10−4 cm3.
In the text