Dynamical effects of the radiative stellar feedback on the H I-to-H2 transition

The atomic-to-molecular hydrogen (H/H2) transition has been extensively studied as it controls the fraction of gas in a molecular state in an interstellar cloud. This fraction is linked to star-formation by the Schmidt-Kennicutt law. While theoretical estimates of the column density of the H I layer have been proposed for static photodissociation regions (PDRs), Herschel and well-resolved ALMA (Atacama Large Millimeter Array) observations have revealed dynamical effects in star forming regions, caused by the process of photoevaporation. We extend the analytic study of the H/H2 transition to include the effects of the propagation of the ionization front, in particular in the presence of photoevaporation at the walls of blister H II regions, and we find its consequences on the total atomic hydrogen column density at the surface of clouds in the presence of an ultraviolet field, and on the properties of the H/H2 transition. We solved semi-analytically the differential equation giving the H2 column density profile by taking into account H2 formation on grains, H2 photodissociation, and the ionization front propagation dynamics modeled as advection of the gas through the ionization front. Taking this advection into account reduces the width of the atomic region compared to static models. The atomic region may disappear if the ionization front velocity exceeds a certain value, leading the H/H2 transition and the ionization front to merge. For both dissociated and merged configurations, we provide analytical expressions to determine the total H I column density. Our results take the metallicity into account. Finally, we compared our results to observations of PDRs illuminated by O-stars, for which we conclude that the dynamical effects are strong, especially for low-excitation PDRs.


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-tomolecular hydrogen (H/H 2 ) transition, or the photodissociation front. Knowledge about the location of the H/H 2 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 2 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 2 transition due to the specific chemistry occurring in the warm H 2 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 2 mass fraction intervening in the Schmidt-Kennicutt law of star-formation (Villanueva et al. 2017), to deduce the CO-to-H 2 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. 1998Stanimirović et al. 1999Stanimirović et al. , 2014Staveley-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 2 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 presented new models for the H/H 2 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. 2012Lee et al. , 2015Roman-Duval et al. 2014;Bialy et al. 2015Bialy et al. , 2017aSchruba 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 2 ratio (Nordon & Sternberg 2016) or the HD/H 2 ratio in the diffuse interstellar gas (Balashev & Kosenko 2020). Also, A&A 656, A65 (2021) Bialy et al. (2017b) explored the effects of turbulence on the H/H 2 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 2 transition in the 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 −3 pc) with a thermal pressure well above its environment at the edge of the PDR (P th ∼ 10 8 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(Goicoechea et al. , 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  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 andBertoldi &, blister H II regions where the ionized gas is free to escape to the surrounding diffuse medium, Israel 1978 andTenorio-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. 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(Franco et al. , 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 2 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 2 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  and Störzer & Hollenbach (1998). The upcoming James Webb Space Telescope (JWST) will give us access to high angular resolution observations of many H 2 rotational and rovibrational lines that will help probe the high pressure layer of PDRs and constrain the photoevaporation mechanism.  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, ∈ [6 − 13.6] 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  about the merging of the fronts, and found that the H 2 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 2 transition model to take into account the dynamics induced by the photoevaporation mechanism in a new semianalytical 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 physical conditions 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 2 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 2 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.

Analytic overview
In this section, we present an analytic model including the most important processes controlling the H/H 2 transition in interstellar clouds exposed to FUV radiation fields and considering the propagation of the ionization front. A65, page 2 of 18 V. Maillard et al.: Dynamical effects of the radiative stellar feedback on the H I-to-H 2 transition

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 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 2 .
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 2 transition is included in the shocked layer. For the purpose of our model describing the H/H 2 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 eq ≈ 500 yr 10 6 cm −3 /n H (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.

Physical model
To consider the impact of photoevaporation on the H/H 2 transition, it is thus necessary to expand on previous studies that considered the H/H 2 transition in a 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 2 chemistry in the presence of this steady flow.
Aside from advection, the processes taken into account for the formation and destruction of H 2 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 2 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: 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 1 + 2N 2 with N 1 the column density of atomic hydrogen and N 2 the column density of H 2 . We also introduce n 1 and n 2 the number densities of H and H 2 respectively. For the photodissociation of H 2 , we assume that the dissociating FUV radiation is only absorbed by grains and H 2 (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 2 photodissociation rate at a column density N is: with and where G 0 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) 1 , u Habing = 2.65 × 10 −14 erg cm −3 (Draine 2011b). We use here the G 0 at the cloud surface (i.e. seeing the incoming radiation field from 2π sr). f shield is the H 2 self-shielding function. The parameter P 0 is only an approximate value of the H 2 dissociation rate at the cloud surface for G 0 = 1 as the true value depends on H 2 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 0 = 5.8 × 10 −11 s −1 . Our P 0 value is the corresponding value for G 0 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 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 proportional to 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 2 self-shielding function is taken from : with y = N 2 /5 × 10 14 cm −2 , b 5 = b/1 km s −1 , and b = 2 km s −1 being the Doppler broadening parameter. A&A 656, A65 (2021) The expression of the equation of the density of H 2 is: 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 2 formation on dust, the second term is its photodissociation by FUV, taking into account H 2 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 2 . This can be interpreted as a stationary flow of H 2 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 eq = 1/(2Rn H ) ≈ 500 yr 10 6 cm −3 /n H 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 = σ g N the dust optical depth in the cloud and τ 1,2 = σ g N 1,2 the optical depths due to the dust associated with the atomic H and the H 2 gas, respectively. We also introduce x 1,2 = n 1,2 /n H the fractional abundances of atomic H and the H 2 . We can now rewrite Eq. (6) in dimensionless form and consider the stationary state (∂/∂t = 0): In Eq. (7), we observe the appearance of two dimensionless quantities representing the physical conditions. The first is P 0 G 0 /R n H , which gives the ratio of the H 2 formation timescale to the H 2 dissociation timescale. This parameter (multiplied by the H 2 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 2 formation timescale to the timescale necessary to advect a dust column density equivalent to τ g = 1. The physical conditions of the cloud G 0 and n H only appear as the G 0 /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 0 /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).

Static case
We first rederive the static solution found in Sternberg et al. (2014). In this situation, v IF = 0. Eq. (7) then becomes: It is possible to solve analytically Eq. (8)  Bertoldi (1996): The proof of the solution to Eq. (8) with Eq. (9) as the H 2 selfshielding function is presented in Appendix A. We obtained that τ 2 can be written as the following function of τ 1 by using that τ g = τ 1 + 2τ 2 : with Γ the incomplete Gamma function and ξ a dimensionless parameter defined as: By rearranging Eq. (10) to have τ 1 (τ 2 ) and taking the limit τ g → +∞ (in which case τ 2 → +∞), we can then obtain the expression for the total H I column density N 1,tot , also presented in Appendix A. It gives: To finally obtain τ 1 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 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 2 transition defined as the position where n 1 = n 2 . 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 2 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 2 with G 0 /n H < 10 −2 cm 3 , the column density of atomic hydrogen is mainly built up in the molecular region. For the strong field limit with G 0 /n H > 10 −2 cm 3 , N 1,tot is mainly from the atomic region. This feature is very well explained in Sternberg et al. (2014) as a consequence of the shape of the H/H 2 transition. Indeed, on Fig. 2, one can find the abundances of H and H 2 as functions of A V for two different G 0 /n H ratios. One can see that in the weak-field regime (low G 0 /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 2 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 P 0 used and discussed above.
To test the impact of the second and more accurate expression of the shielding function from  2 We define the weak and strong field limits as done in Sternberg et al. (2014), where αG < 1 for the weak field limit, corresponding to roughly G 0 /n H < 10 −2 cm 3 . For the strong field limit, αG > 1, corresponding to G 0 /n H > 10 −2 cm 3 . (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 2 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 2 self-shielding functions. We now only use Eq. (5) as the H 2 self-shielding function in the following.

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 0 and n H independently but on the G 0 /n H ratio. So, the free parameters that are explored are G 0 /n H , Z (appearing in both R and σ g expressions) and v IF . The resolution of this equation gives the spatial profile of τ 2 from which we derive τ 1 , and then x 1 and x 2 .
We take v IF as a free parameter, but as explained in , 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.

Advection effect on the transition
The profile of the abundances of H and H 2 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 0 /n H ratio of 0.5 cm 3 and ionization front velocity v IF = 1 km s −1 . In this case, we clearly see that the H/H 2 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 2 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 is sharper. To illustrate that, we can compare the solid lines of Fig. 4 representing the H/H 2 transition for a nonzero velocity to the dashed line in Fig. 2 representing the H/H 2 transition without advection. These two cases have the H/H 2 transition roughly at the same position, but in the advection case the transition is much sharper.

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 . 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 different G 0 /n H ratio and determined the critical value of v IF above which front merging occurs for each value of G 0 /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.  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 0 /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 2 self-shielding, which is why their criterion is a proportionality relation. Storzer & Hollenbach, 1998Bertholdi & Draine, 1996 Fit of the criterion New semi-analytical criterion Modeled merging values Fig. 5. Criterion describing, for a given ionization front velocity, the G 0 /n H 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 , 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).
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 2 , representing how long it takes for H 2 to flow across the atomic region with the H 2 photodissociation timescale. If the dynamical timescale is shorter than the photodissociation timescale, the advected H 2 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 2 formation on grains and dust absorption. They also assumed a fixed H 2 column density value when taking H 2 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 2 self-shielding,  are closer to the regime at high values of v IF , whereas Störzer & Hollenbach (1998), by taking a fixed value of the H 2 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 2 self-shielding. Indeed, for low values of the G 0 /n H ratio, Sternberg et al. (2014) explained that H 2 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 2 at the position of the transition and at the merging of the fronts. We can rewrite Eq. (7)  transition where x 1 = x 2 = 1/3 and when the fronts are merged, namely when τ g = τ 2 = 0. With these assumptions, H 2 and dust shieldings are equal to 1, leading to: From here we need to find the x 2 derivative by isolating dx 2 /dτ g in the previous equation and by taking the merging conditions found in the semi-analytical models. The fit gives us: with τ 0 g the depth of the H/H 2 transition in the static case for this G 0 /n H ratio. Then we can deduce that: for which we need to numerically compute τ 0 g 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 (G 0 /n H ) 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:

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 −1 as functions of the G 0 /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 0 /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 2 transition approaches and merges with the ionization front, N 1 integrated in the atomic region drops to 0. For G 0 /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 0 /n H ratio, where N 1,tot is dominated by the atomic hydrogen accumulated in the atomic region. The other is for low G 0 /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 0 /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 −1 . 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 0 /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 0 /n H ratio (100 cm 3 ). 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 0 /n H (10 −5 cm 3 ). 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 0 /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 −3 km s −1 : 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):

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 P 0 G 0 /Rn H , it is then equivalent to a larger G 0 /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 P 0 G 0 /Rn 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 −1 , 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 1 = σ g τ 1 .
Finally, a decreased metallicity has an impact on H 2 selfshielding, as the self-shielding term is f shield (τ 2 /σ g ) and σ g decreases as Z . As a consequence, H 2 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: 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 semianalytical results were fit to obtain the black curve. We present hereafter the final formula for both of the regimes and including metallicity: In the dissociated fronts regime, where G 0 /n H > (G 0 /n H ) crit (v IF , Z ) and where ∆N 1,tot is constant with respect to G 0 /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 H < (G 0 /n H ) crit (v IF , Z )), the effect or metallicity can be described as a shift to the left of the fit as Z decreases, modeled by the √ Z factor in the fit.

The link between G 0 /n H and v IF
As explained in Sect. 2.4, we treated v IF as a free parameter independent of G 0 . We investigate in this section how v IF is linked to our G 0 /n H ratio. First, mass conservation across the ionization front (Eq. (37.1), Draine 2011b) yields: 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 0 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 0 /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 −1 ) and the distance between the star and the ionization front is d star . Following , 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 : and, for a radiation field coming from a point source (star), the G 0 /n H ratio is 3 : with F ISRF = 1.27 × 10 7 cm −2 s −1 . Dividing the two expressions and using Eq. (12) from : 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 0 /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 , 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  for R V = 3.1: As we see in Sect. 4.4, for a various set of objects, τ EUV = σ EUV N 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 0 /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.
3 See footnote 1 on the FUV range relevant for H 2 photodissociation.  Figure 11 represents the link between v IF and the G 0 /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 0 /n H ratio under 10 −2 cm 3 . For stronger fields, the link between v IF and the G 0 /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 That is what is presented on Fig. 12 for the same range of τ EUV . We see that ∆N 1,tot is very close to 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 that the effects are more important in the strong field limit, because for such conditions the ionization front velocity is higher, meaning that the H/H 2 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. (1989Franco et al. ( , 1990, respectively for spherically symmetric and disk-like clouds, an ionization front propagating from the inside of the structure is accelerated along A65, page 10 of 18 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: By comparing it to the timescale to reach a stationary PDR structure, which is (2 R n H ) −1 (Hollenbach & Natta 1995), one finds that to maintain a stationary state, the logarithmic derivative of the density has to verify: 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 × 10 −3 10 5 cm −3 n H v IF 1 km s −1 pc.
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 2 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 (1979Giuliani ( , 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. They also 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.

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 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 and where c PDR = γT/µ 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 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 −1 . The Alfvén speed in the PDR is taken in a wide range of 0 to 5 km s −1 , corresponding to a transverse magnetic field of 0 to a few 100 µG. With c II = 10 km s −1 we deduce that our upper boundary for v IF in D-type fronts ranges between 0.7 to 1.4 km s −1 , while the lower boundary for v IF in R-type fronts ranges between 18.6 and 19.3 km s −1 . 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 blistertype 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 Dtype 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.

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, pressureconfined 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, gasdust drift, and dust photodestruction (for the conditons of RCW 120) where presented by Akimkin et al. (2017). They found that radiation pressure has a 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 PDRs found 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 2 transition as described by the results of our semi-analytical model (Sect. 2 and 3).

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 Bstar PDRs (NGC7023,...), it indeed yields a very small ionization front velocity, on the order of less than 10 −4 km s −1 and a G 0 /n H ratio around 10 −3 -10 −2 cm 3 , 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.
As found in Sect. 3.1, we nevertheless predict moderate but non negligible effects of the advection dynamics on the H/H 2 transition. A second group at low G 0 /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 in the 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 2 could in particular be expected in these low G 0 /n H PDRs, in contrast to the conclusions of Störzer & Hollenbach (1998) who focused on high G 0 /n H conditions representative of Orion Bar-like PDRs.

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 1 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 1 integrated over the atomic region for v IF = 1 km s −1 and G 0 /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.

Conclusions
To sum up, we built a semi-analytical model of the H/H 2 transition in a 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 2 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 2 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 0 , 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 0 /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 2 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 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 2 transition is closer to the ionization front than modeled by previous static studies. Some H 2 is then found in a hotter region, richer in FUV photons. One can then expect H 2 ro-vibrational lines to be more intense as collisions and UV-pumping are more efficient. As a lot of H 2 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.
with σ C x(C) 1.7 × 10 −21 cm 2 , which is very close to dust extinction cross-section σ g = 1.9 × 10 −21 cm 2 . 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 2 self-shielding are met in the strong field limit. In the weak field limit, where G 0 /n H ≤ 10 −2 cm 3 , the H/H 2 transition is driven by H 2 self-shielding and not by dust, and atomic carbon extinction can thus be expected to have no impact on the H/H 2 transition. In the strong field limit (G 0 /n H ≥ 10 −2 cm 3 ), its impact depends on the abundance of atomic carbon before the H/H 2 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 results presented in Fig. B.1 and B.2 for different G 0 /n H ratio, from 1 to 10 −4 cm 3 . On Fig. B.1, we have the abundances profiles of H, H 2 , C + , C, and CO as functions of A V . Fig. B.2 presents the extinction terms of dust, C, H 2 , 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 2 , it is equal to the shielding function f shield (H 2 ). 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.
In the high G 0 /n H regime (see e.g., G 0 /n H = 1 cm 3 , top panel), atomic carbon is not very abundant at the H/H 2 transition position with x(C) 10 −7 for G 0 /n H = 1 cm 3 , 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 0 /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 2 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 2 self-shielding, with negligible contributions of carbon or dust extinction. In intermediate 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 2 transition. Deeper in the cloud, atomic carbon continuum extinction of course becomes important, and in particular, we see that it cannot be neglected 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 2 transition. In the case of high G 0 /n H , the very low fraction of atomic carbon in the atomic region makes it negligible compared to dust extinction. In the case of low G 0 /n H , where the fraction of atomic carbon in the atomic region increases to significant levels, the H/H 2 transition is dominated by H 2 self-shielding and appears at a very low column density, so that extinction by both dust and atomic carbon is negligible. Extinction factor of dust in blue, of C in red, of H 2 selfshielding in green, and of all of them combined in black, as functions of A V . C I and H 2 column densities are taken from Meudon PDR Code simulations for models with G 0 /n H = 1, 10 −2 , 10 −3 , and 10 −4 cm 3 .