Photoevaporation of Jeans-unstable molecular clumps

TL;DR

This study uses 3D simulations to analyze how FUV radiation affects the evolution and star-forming potential of Jeans-unstable molecular clumps, revealing that radiation influences H$_2$ content but does not prevent collapse.

Contribution

It provides the first detailed 3D hydrodynamical simulations including non-equilibrium chemistry to study FUV photoevaporation of Jeans-unstable molecular clumps.

Findings

FUV radiation reduces H$_2$ mass by 40-90% initially.

Clumps remain gravitationally unstable despite mass loss.

FUV intensity regulates the residual H$_2$ available for star formation.

Abstract

We study the photoevaporation of Jeans-unstable molecular clumps by isotropic FUV (6 eV $< {\rm h}\nu$ < 13.6 eV) radiation, through 3D radiative transfer hydrodynamical simulations implementing a non-equilibrium chemical network that includes the formation and dissociation of H$_2$. We run a set of simulations considering different clump masses ($M=10-200$ M$_\odot$) and impinging fluxes ($G_0=2\times 10^3-8\times 10^4$ in Habing units). In the initial phase, the radiation sweeps the clump as an R-type dissociation front, reducing the H$_2$ mass by a factor 40-90%. Then, a weak ($\mathcal{M}\simeq 2$) shock develops and travels towards the centre of the clump, which collapses while loosing mass from its surface. All considered clumps remain gravitationally unstable even if radiation rips off most of the clump mass, showing that external FUV radiation is not able to stop clump collapse. However, the FUV intensity regulates the final H$_2$ mass available for star formation: for example, for $G_0 < 10^4$ more than 10% of the initial clump mass survives. Finally, for massive clumps ($\sim 100$ M$_\odot$) the H$_2$ mass increases by 25-50% during the collapse, mostly because of the rapid density growth that implies a more efficient H$_2$ self-shielding.

Figures

Tables

Table 1. Table 1 : List of photochemical reactions included in our chemical network.

Photochemical reactions
$H+\gamma \to H^{+}+e$	$H_2^{+}+\gamma \to H^{+}+H$
$\mathrm{He}+\gamma \to {\mathrm{He}}^{+}+e$	$H_2^{+}+\gamma \to H^{+}+H^{+}+e$
${\mathrm{He}}^{+}+\gamma \to {\mathrm{He}}^{++}+e$	$H_2^{+}+\gamma \to H+H$ (direct)
$H^{-}+\gamma \to H+e$	$H_2^{+}+\gamma \to H+H$ (Solomon)
$H_2+\gamma \to H_2^{+}+e$

Table 2. Table 2 : Summary of the 3D simulation run in this work. Given a mass M 𝑀 M , the corresponding radius R 𝑅 R and density (average number density ⟨ n ⟩ delimited-⟨⟩ 𝑛 \langle n\rangle and central number density n c subscript 𝑛 𝑐 n_{c} ) are determined, as detailed in Sec. 3.1 . Simulations of clumps with the same mass differ for the intensity of the external source of FUV radiation G 0 subscript 𝐺 0 G_{0} , that is calculated at the clump surface (see Sec. 3.2 ). The free-fall time is also reported for reference.

	$M\left[M_{\odot }\right]$	$R\left[\mathrm{pc}\right]$	$⟨n⟩\left[{\mathrm{cm}}^{-3}\right]$	$n_c\left[{\mathrm{cm}}^{-3}\right]$	$t_{\mathrm{ff}}\left[\mathrm{Myr}\right]$	$G_0$
clump_M50_noRad	50	0.5	$3.9\times {10}^3$	$6.6\times {10}^3$	0.81	0
clump_M50_G2e3	50	0.5	$3.9\times {10}^3$	$6.6\times {10}^3$	0.81	$2\times {10}^3$
clump_M50_G3e4	50	0.5	$3.9\times {10}^3$	$6.6\times {10}^3$	0.81	$3\times {10}^4$
clump_M50_G8e4	50	0.5	$3.9\times {10}^3$	$6.6\times {10}^3$	0.81	$8\times {10}^4$
clump_M10_G3e4	10	0.2	$1.5\times {10}^4$	$2.1\times {10}^4$	0.41	$3\times {10}^4$
clump_M100_G3e4	100	0.7	$2.5\times {10}^3$	$4.4\times {10}^3$	1.01	$3\times {10}^4$
clump_M200_G3e4	200	1.0	$1.6\times {10}^3$	$3.3\times {10}^3$	1.26	$3\times {10}^4$

Equations

$${\Gamma }_{{\mathrm{H}}_2}=7.5\times 10^{-11}\frac{n_{\gamma }(12.24-13.51\mathrm{eV})}{2.256\times 10^4\mathrm{c}{\mathrm{m}}^{-3}}{\mathrm{s}}^{-1},$$

$${\tau }_i=\sum _j{n}_j\Delta x_{\mathrm{cell}}{\sigma }_{ij},$$

$$g(s)=\frac{1}{(2\pi {\sigma }^2{)}^{1/2}}\exp \left[-\frac{1}{2}\left(\frac{s-s_0}{\sigma }\right)\right],$$

$${\sigma }^2=\ln \left(1+b^2{\mathcal{M}}^2\frac{\beta }{\beta +1}\right),$$

$$R=\frac{1}{2}{\lambda }_{\text{textsc}j}=\frac{1}{2}\frac{\pi {\sigma }^2+\sqrt{36\pi c_s^{2}G{L}^2{m}_p\mu n+{\pi }^2{\sigma }^4}}{6GL{m}_p\mu n}$$

$$n(r)=\begin{cases}n_{c}&\quad r<R/2\\ n_{c}(2r/R)^{-1.5}&\quad R/2\leq r<R\\ \end{cases}.\$$

$$v_{\mathrm{front}}=\frac{\varphi }{n}\simeq \frac{1.6\times 10^{-3}{G}_0}{⟨n⟩⟨h\nu {⟩}_{\text{textsc}fuv}}\simeq 10^9\mathrm{cm}{\mathrm{s}}^{-1}>c_{\mathrm{red}}$$

$$N_{\mathrm{trans}}=0.7\ln \left[{\left(\frac{\alpha G}{2}\right)}^{1/0.7}+1\right]{\left(1.9\times 10^{-21}Z\right)}^{-1}{\mathrm{cm}}^{-2}$$

$$\alpha G=0.35G_0\left(\frac{100\mathrm{c}{\mathrm{m}}^{-3}}{n_{\text{textsc}h}}\right){\left(\frac{9.9}{1+8.9Z/{\mathrm{Z}}_{\odot }}\right)}^{0.37}.$$

$$\log (f_{M_{{\mathrm{H}}_2}})=-0.85{\log }^2{G}_0+0.22\log G_0-0.38$$

$$x(t)={\left(\frac{1}{n{\alpha }_{\text{textsc}b}(T)}{\int }_{13.6eV}^{\infty }{F}_{\nu }{a}_{\nu }\mathrm{d}\nu \right)}^{1/2}$$

Full text

††pagerange: Photoevaporation of Jeans-unstable molecular clumps–C††pubyear: 2019

1 Introduction

Stars are known to form in clusters inside giant molecular clouds (GMCs), as a consequence of the gravitational collapse of overdense clumps and filaments (Bergin et al., 1996; Wong et al., 2008; Takahashi et al., 2013; Schneider et al., 2015; Sawada et al., 2018). The brightest (e.g. OB) stars have a strong impact on the surrounding interstellar medium (ISM), since their hard radiation field ionizes and heats the gas around them, increasing its thermal pressure. As a result, the structure of the GMC can be severely altered due to the feedback of newly formed stars residing inside the cloud, with the subsequent dispersal of low-density regions. Collapse can then occur only in dense regions able to self-shield from impinging radiation (Dale et al., 2005, 2012a, 2012b; Walch et al., 2012).

The ISM within the Strömgren sphere around a star-forming region is completely ionized by the extreme ultra-violet (EUV) radiation, with energy above the ionization potential of hydrogen ($h\nu>13.6$ eV). The typical average densities of H \scriptstyle\rm II\regions are $\langle n\rangle\simeq 100\,\,{\rm{cm^{-3}}}$: this results in ionization fractions $x_{\textsc{hii}}<10^{-4}$ and a final gas temperature $T>10^{4-5}$ K. Far-ultraviolet (FUV) radiation (photon energy $6\,\mathrm{eV}<h\nu<13.6\,\mathrm{eV}$) penetrates beyond the H \scriptstyle\rm II\region, thus affecting the physical and chemical properties of the ISM up to several parsecs. This region is usually referred to as the Photo-Dissociation Region (PDR; Tielens & Hollenbach, 1985; Kaufman et al., 1999; Le Petit et al., 2006; Bron et al., 2018). Typical fluxes in the FUV band due to OB associations may have values as high as $G_{0}=10^{4-5}$ (Marconi et al., 1998), in units of the Habing flux111The Habing flux ($1.6\times 10^{-3}\,{\rm erg}\,{\rm s}^{-1}{\rm cm}^{-2}$) is the average interstellar radiation field of our Galaxy in the range [6 eV, 13.6 eV] (Habing, 1968). A PDR is characterized by a layer with neutral atomic hydrogen, photo-dissociated by Lyman-Werner photons ($11.2\,\mathrm{eV}<h\nu<13.6\,\mathrm{eV}$), and a deeper layer where gas self-shielding allows hydrogen to survive in molecular form. There are many observational pieces of evidence that the structure of PDRs are not homogeneous, with gas densities spanning many orders of magnitude from $10^{2}\,\,{\rm{cm^{-3}}}$ to $10^{6}\,\,{\rm{cm^{-3}}}$. In particular, small isolated cores of few solar masses and sizes $\apprle$ 0.1 pc are commonly observed (Reipurth, 1983; Hester et al., 1996; Huggins et al., 2002; Mäkelä & Haikala, 2013). Radiative feedback by FUV radiation could explain their formation via shock-induced compression (Lefloch & Lazareff, 1994).

The effect of FUV radiation on clumps is twofold: (1) FUV radiation dissociates the molecular gas, which then escapes from the clump surface at high velocity (photoevaporation); (2) radiation drives a shock which induces the clump collapse (radiation-driven implosion, RDI Sandford et al., 1982). The first effect reduces the clump molecular mass, hence decreasing the mass budget for star formation within the clump. Instead, the latter effect may promote star formation by triggering the clump collapse. Hence, the net effect of radiative feedback on dense clumps is not trivial and deserves a careful analysis.

The flow of gas from clumps immersed in a radiation field has been studied by early works both theoretically (Dyson, 1968; Mendis, 1968; Kahn, 1969; Dyson, 1973) and numerically (Tenorio-Tagle, 1977; Bedijn & Tenorio-Tagle, 1984). Bertoldi (1989) and Bertoldi & McKee (1990) developed semi-analytical models to describe the photoevaporation of atomic and molecular clouds induced by ionizing radiation. In their models they also include the effects of magnetic fields and self-gravity. They find that clumps settle in a stationary cometary phase after the RDI, with clump self-gravity being negligible when the magnetic pressure dominates with respect to the thermal pressure (i.e. $B>6\,\mu\mathrm{G}$), or when the clump mass is much smaller than a characteristic mass $m_{\mathrm{ch}}\simeq 50\,{\rm M}_{\odot}$. They focused on gravitationally stable clumps, thus their results are not directly relevant for star formation.

Later, the problem was tackled by Lefloch & Lazareff (1994), who performed numerical simulations which however only included the effect of thermal pressure on clump dynamics. Gravity was then added for the first time by Kessel-Deynet & Burkert (2003). For an initially gravitationally stable clump of 40 ${\rm M}_{\odot}$, they find that the collapse can be triggered by the RDI; nevertheless, they notice that the collapse does not take place if a sufficient amount of turbulence is injected ($v_{\rm rms}\simeq 0.1\,{\rm km/s}$). Bisbas et al. (2011) also ran simulations of photoevaporating clumps, with the specific goal of probing triggered star formation. They find that star formation occurs only when the intensity of the impinging flux is within a specific range ($10^{9}\,{\rm cm}^{-2}{\rm s}^{-1}<\Phi_{\textsc{euv}}<3\times 10^{11}\,{\rm cm}^{-2}{\rm s}^{-1}$ for a $5\,{\rm M}_{\odot}$ initially stable clump 222Assuming for example the spectrum of a $10^{4}\,{\rm L}_{\odot}$, this EUV flux corresponds roughly to a flux $G_{0}=10-3\times 10^{4}$ in the FUV band. ). All these works include the effect of ionizing radiation only, while FUV radiation feedback is instead relevant for clumps located outside the H \scriptstyle\rm II\region of a star (cluster).

In a previous work (Decataldo et al., 2017, hereafter D17), we have constructed a 1D numerical procedure to study the evolution of a molecular clump, under the effects of both FUV and EUV radiation. We have followed the time evolution of the structure of the iPDR (ionization-photodissociation region) and we have computed the photoevaporation time for a range of initial clump masses and intensity of impinging fluxes. However, since D17 did not account for gravity, ${\mathrm{H}_{2}}$ dissociation is unphysically accelerated during the expansion phase following RDI. Those results have been compared with the analytical prescriptions by Gorti & Hollenbach (2002), finding photoevaporation times in agreement within a factor 2, although different simplifying assumptions where made in modelling the clump dynamics.

The same setup by D17 has been used by Nakatani & Yoshida (2018) to run 3D simulations with on-the-fly radiative transfer (RT) and a chemical network including $\mathrm{H}^{+}$, ${\mathrm{H}_{2}}$, $\mathrm{H}^{+}$, $\mathrm{O}$, $\mathrm{CO}$ and $\mathrm{e}^{-}$. Without the inclusion of gravity in their simulations, they find that the clump is confined in a stable cometary phase after the RDI, which lasts until all the gas is dissociated and flows away from the clump surface. Nevertheless, they point out that self-gravity may affect the clump evolution when photoevaporation is driven by a FUV-only flux, while the EUV radiation produces very strong photoevaporative flows which cannot be suppressed by gravity.

In this paper, we attempt to draw a realistic picture of clump photoevaporation by running 3D hydrodynamical simulations with gravity, a non-equilibrium chemical network including formation and photo-dissociation of ${\mathrm{H}_{2}}$, and an accurate RT scheme for the propagation of FUV photons. We focus on the effect of radiation on Jeans-unstable clumps, in order to understand whether their collapse is favoured or suppressed by the presence of nearby stars emitting in the FUV range.

The paper is organized as follows. In Sec. 2, we describe the numerical scheme used for the simulations, and, in Sec. 3, the initial conditions for the gas and radiation. We analyse the evolution of the clump for different radiative fluxes and masses is Sec. 4 and Sec. 5, respectively. A cohesive picture of the photoevaporation process is given in Sec. 6. Our conclusions are finally summarized in Sec. 7.

2 Numerical scheme

Our simulations are carried out with ramses-rt333https://bitbucket.org/rteyssie/ramses, an adaptive mesh refinement (AMR) code featuring on-the-fly RT (Teyssier, 2002; Rosdahl et al., 2013). RT is performed with a momentum-based approach, using a first-order Godunov solver and the M1 closure relation for the Eddington tensor. The basic version of ramses-rt thermochemistry module accounts only for the photoionization of hydrogen and (first and second) photoionization of helium. We have used the chemistry package krome 444https://bitbucket.org/tgrassi/krome (Grassi et al., 2014) to implement a complete network of H and He reactions, including neutral and ionized states of He, H and ${\mathrm{H}_{2}}$.

We track the time evolution of the following 9 species: H, $\mathrm{H}^{+}$, $\mathrm{H}^{-}$, $\mathrm{H}_{2}$, $\mathrm{H}_{2}^{+}$, He, $\mathrm{He}^{+}$, $\mathrm{He}^{++}$ and free electrons. Our chemical network includes 46 reactions in total555The included reactions and the respective rates are taken from Bovino et al. (2016): reactions 1 to 31, 53, 54 and from 58 to 61 in Tab. B.1 and B.2, photoreactions P1 to P9 in Tab. 2., and features neutral-neutral reactions, charge-exchange reactions, collisional dissociation and ionization, radiative association reactions and cosmic ray-induced reactions (we consider a cosmic ray ionization rate $\zeta_{\mathrm{H}}=3\times 10^{-17}\,{\rm s}^{-1}$, the reference value in the Milky Way (Webber, 1998). We follow ${\mathrm{H}_{2}}$ formation on dust, adopting the Jura rate at solar metallicity $R_{f}({\mathrm{H}_{2}})=3.5\times 10^{-17}n_{\mathrm{H}}n_{\mathrm{tot}}$ (Jura, 1975). There are 9 reactions involving photons, listed in Tab. 1: photoionization of H, He, $\mathrm{He}^{+}$, $\mathrm{H}^{-}$ and $\mathrm{H}_{2}$ to $\mathrm{H}_{2}^{+}$, direct photodissociation of $\mathrm{H}_{2}^{+}$ and $\mathrm{H}_{2}$ and the two-step Solomon process (Draine & Bertoldi, 1996). The Solomon process rate is usually taken to be proportional to the total flux at 12.87 eV (Glover & Jappsen, 2007; Bovino et al., 2016), but this is correct only if the flux is approximately constant in the Lyman-Werner band (11.2 - 13.6 eV), as pointed out by Richings et al. (2014a). These authors find that the most general way to parametrize the dissociation rate is

[TABLE]

where $n_{\gamma}(\Delta E_{\mathrm{bin}})$ is the photon density in the energy interval $\Delta E_{\mathrm{bin}}$. We work in the on-the-spot approximation, hence photons emitted by recombination processes are neglected.

Given the chemical network and the included reactions, the code can in principle be used with an arbitrary number of photon energy bins. In the particular context of photoevaporating clumps, we decided to make only use of two bins with energies in the FUV domain, i.e. [6.0 eV, 11.2 eV] and [11.2 eV, 13.6 eV]. As we consider molecular clumps located outside stellar H \scriptstyle\rm II\regions, we expect that EUV radiation does not reach the surface of the clump. On the other hand, we neglect photons with energies $<6.0$ eV since they do not take part in any chemical reactions of interest in our case.

Fig. 1 summarizes the approach we used to couple RT module in ramses with the non-equilibrium network adopted via krome (as also done in Pallottini et al., 2019). At each time-steps ($\Delta t$), photons are first propagated from each cell to the nearest ones by ramses-rt (A). Then, the gas-radiation interaction step (B) is executed, sub-cycling in absorption (B1) and chemical evolution (B2) steps with a time-steps $\Delta t_{i}<\Delta t$, such that the flux is not reduced by more than $20\%$ at each substep ($\Delta t_{i}$).

In step B1, we account for (1) photons that take part in chemical reactions, (2) $\mathrm{H}_{2}$ self-shielding and (3) dust absorption. The optical depth of a cell in the radiation bin $i$ (excluding the Solomon process) is computed by summing over all photo-reactions:

[TABLE]

where $n_{j}$ is the number density of the photo-ionized/dissociated species in the reaction $j$, $\Delta x_{\rm cell}$ is the size of the cell, and $\sigma_{ij}$ is the average cross section of the reaction $j$ in the bin $i$. For the Solomon process, the self-shielding factor $S_{\mathrm{self}}^{\mathrm{H}_{2}}$ is taken from Richings et al. (2014b), and it is related to the optical depth by $\tau_{\mathrm{self}}^{\mathrm{H}_{2}}=-\log(S_{\mathrm{self}}^{\mathrm{H}_{2}})/\Delta x_{\mathrm{cell}}$. Absorption from dust is included, with opacities taken from Weingartner & Draine (2001). We have used the Milky Way size distribution for visual extinction-to-reddening ratio $R_{V}=3.1$, with carbon abundance (per H nucleus) $b_{C}=60$ ppm in the log-normal populations666https://www.astro.princeton.edu/~draine/dust/dustmix.html.

After every absorption substep, krome is called in each cell (step B2): photon densities in each bin are passed as an input, together with the current chemical abundances and the gas temperature in the cell, and krome computes the new abundances after a time-steps accordingly.

In App. A, we show two successful tests that we performed to validate our scheme for the coupling between ramses-rt and krome:

• A.1

An ionized region, comparing the results with the analytical solution;

• A.2

The structure of ${\rm{H_{2}}}\,\,$in a PDR, compared with the standard benchmarks of Röllig et al. (2007).

3 Simulation setup

3.1 Gas

The computational box is filled with molecular gas777The gas has helium relative mass abundance $X_{\rm He}=25\%$. of density $n=100\,\,{\rm{cm^{-3}}}$ and metallicity $Z={\rm Z}_{\odot}$. A dense clump ($n=10^{3}-10^{4}\,\,{\rm{cm^{-3}}}$) is then located at the centre of the domain, with the same initial composition of the surrounding gas (Fig. 2).

Clumps in GMCs are self-gravitating overdensities. Observations of GMCs of different sizes and masses (Hobson, 1992; Howe et al., 2000; Lis & Schilke, 2003; Minamidani et al., 2011; Parsons et al., 2012; Barnes et al., 2018; Liu et al., 2018) show that clumps have a wide range of physical properties: radii range from $0.1-10$ pc, densities can be $10-10^{4}$ times the average density of the GMC; typical masses range from few solar masses to few hundreds ${\rm M}_{\odot}$.

The density distribution of clumps in GMCs has been studied with numerical simulations of supersonic magnetohydrodynamic turbulence (Padoan & Nordlund, 2002; Krumholz & McKee, 2005; Padoan & Nordlund, 2011; Federrath & Klessen, 2013), yielding a log-normal PDF (Probability Distribution Function):

[TABLE]

where $b\simeq 0.3-1$ is a factor taking into account the kinetic energy injection mechanism which is driving the turbulence (Molina et al., 2012). Including self-gravity, the high-density end of the PDF is modified with a power-law tail $g(n)\sim n^{-\kappa}$, with $\kappa\sim 1.5-–2.5$ (Krumholz & McKee, 2005; Padoan & Nordlund, 2011; Federrath & Klessen, 2013; Schneider et al., 2015).

If a value $n$ of the density is drawn from the PDF $g(s)$, the corresponding radius of the clump can be estimated with the turbulent Jeans length (Federrath & Klessen, 2012):

[TABLE]

where $c_{s}$ is the isothermal sound speed, $L$ is the size of the GMC, $m_{p}$ the proton mass and $\mu$ the mean molecular weight. The corresponding clump mass is estimated by assuming a spherical shape and uniform density.

In Fig. 3 the clump radius is plotted as a function of the number density, for different Mach numbers $\mathcal{M}$, GMC size fixed at $L=25$ pc and $c_{s}$ computed with a temperature $T=10$ K and molecular gas ($\mu=2.5$). The coloured points correspond to the position of clumps with different masses $M=10-200\,{\rm M}_{\odot}$ in the $R$-$n$ diagram. In GMCs with the same Mach number, clumps with larger mass are less dense: indeed, as a rough approximation, we have that $R\propto n^{-1}$ and $M\sim m_{p}\mu nR^{3}\propto n^{-2}$. It is also interesting to check how the properties of clumps with the same mass vary for different parent GMCs. Considering the 10 ${\rm M}_{\odot}$ clump (red point), we notice that decreasing $\mathcal{M}$, its position in the $R$-$n$ diagram shifts towards lower density and larger radius. Thus, in a $\mathcal{M}=1$ cloud, clumps with mass higher than 10 ${\rm M}_{\odot}$ would have a density as low as the average cloud density, implying that clumps more massive than 10 ${\rm M}_{\odot}$ do not exist at all in such a cloud.

Here we consider clumps residing in a GMC with size $L=25$ pc, average temperature $T=10$ K and Mach number $\mathcal{M}=15$. In particular, we explore the range of masses represented in Fig. 3, i.e. $M=10,50,100,200\,{\rm M}_{\odot}$. Masses, radii and average densities of these clumps are summarized in Tab. 2. Clumps are modelled as spheres located at the centre of the computational box; their initial density profile is constant up to half of the radius, and then falls as a power law:

[TABLE]

where for each clump the core density $n_{c}$ is chosen such that the total mass is the selected one (Tab. 2). The profile in Eq. 5 has been used in simulations of molecular clouds (e.g. Krumholz et al., 2011) and it is physically motivated by observations of star-forming clumps (Beuther et al., 2002; Mueller et al., 2002). The choice of this profile implies that there is a discontinuity in the gas density (and so the pressure) at the interface with the external ISM. In the simulations, this will produce a slight expansion of the clump just before the stellar radiation hits the clump surface.

A turbulent velocity field is added to the clump in the initial condition. We generate an isotropic random Gaussian velocity field with power spectrum $P(k)\propto k^{-4}$ in Fourier space, normalising the velocity perturbation so that the virial parameter $\alpha=5\,v_{\rm rms}^{2}\,R_{c}/G\,M_{c}$ is equal to 0.1, as measured for some clumps with mass around $10^{2}\,{\rm M}_{\odot}$ (e.g. Parsons et al., 2012). In three dimensions, the chosen power spectrum gives a velocity dispersion that varies as $\ell^{1/2}$ with $\ell$ the length scale, which is in agreement with Larson scaling relations (Larson, 1981).

3.2 Radiation sources

We set-up a roughly homogeneous radiation field around the clump by placing 50 identical point sources (i.e. stars), randomly distributed on the surface of a sphere centred on the clump and with radius larger than the clump radius ($R_{\mathrm{sources}}=R_{\mathrm{clump}}+1\,\mathrm{pc}$)888The choice of 1 pc as a distance of the sources from the clump is arbitrary. In fact the aim is to get a specific $G_{0}$ at the clump surface, which can be obtained either varying the source luminosity or the source distance. The number of sources is also not relevant, provided that it is large enough to ensure a nearly isotropic flux on the clump surface., as depicted in Fig. 2. Each star has a black-body spectrum and a bolometric luminosity in the range $L_{\star}=10^{4}-10^{6}\,{\rm L}_{\odot}$, according to the desired FUV flux at the clump surface (Tab. 2).

In our simulations we use the GLF scheme999To solve numerically the propagation of photons, different function for the intercell photon flux can be used. ramses implements the Harten-Lax-van Lee (HLL) function (Harten et al., 1983; González et al., 2007) and the Global Lax Friedrich (GLF), obtained by setting the HLL eigenvalues to the speed of light. for the propagation of photons, since it is more suitable for isotropic sources (while the HLL scheme@footnotemark introduces asymmetries, see Rosdahl et al. 2013). Our configuration of sources is prone to the problem of “opposite colliding beams”: the photon density is higher than expected in the cells where two or more fluxes come from opposite directions. This problem is due to the M1 closure relation (see González et al. 2007; Aubert & Teyssier 2008 and Rosdahl et al. 2013, in particular their Fig. 1) and it is detailed in App. B. To circumvent this issue, we compute the average flux on the clump surface at the beginning of the simulation, which could be higher than that expected from an analytical calculation. In Tab. 2 we list the resulting FUV flux at the clump surface, for the different set-ups.

3.3 Resolution

The coarse grid has a resolution of $128^{3}$ cells, which implies a cell size $\Delta x_{\mathrm{cell}}\simeq 0.047$ pc. We include one AMR level according to a refinement criterion based on the ${\mathrm{H}_{2}}$ abundance gradients: a cell is refined if the ${\mathrm{H}_{2}}$ abundance gradient with neighbouring cells is higher than $10\%$. Thus the effective resolutions is increased up to $256^{3}$ cells with size $\Delta x_{\mathrm{cell}}\simeq 0.023$ pc. For the control run without radiation, the resolution is increased by 4 additional levels of refinement in a central region of radius 0.05 pc.

The expected time-scale of photoevaporation is of the order of 1 Myr (Gorti & Hollenbach, 2002; Decataldo et al., 2017). Simulations are carried out with a reduced speed of light $c_{\mathrm{red}}=10^{-3}\,c$, where $c$ is the actual speed of light, in order to prevent exceedingly small time-stepss, which would result in a prohibitively long computational time. Indeed, time-stepss are settled by the light-crossing time of cells in the finest grid, hence in our simulations $\Delta t\simeq 75$ yr with reduced speed of light. The reduced speed of light affects the results of simulations when $c_{\mathrm{red}}$ is lower than the speed of ionization/dissociation fronts (Deparis et al., 2019; Ocvirk et al., 2019), which is given by $v_{\mathrm{front}}=\phi/n$, where $\phi$ is the photons flux in cm*-2s-1*. In our set of simulations, clump_M50_G8e4 has the highest $G_{0}/\langle n\rangle$ ratio, yielding

[TABLE]

This points out that the propagation of the dissociation front (DF) is not treated accurately in our simulations. Nevertheless, the DF propagates at $v_{\mathrm{front}}$ only for a time around few kyr, after which the front stalls and the photoevaporation proceeds for about 0.1-1 Myr. Hence, most of the simulation is not influenced by the reduced speed of light approximation, and the error concerns only the speed of the dissociation front (and not the thermochemical properties of the photo-dissociated gas).

We have checked that the resolution of our simulation is sufficient to describe the effect of radiation on dense gas. In fact, FUV radiation dissociates and heats a shell of gas at the clump surface, hence the simulation is physically accurate only if the thickness of this layer is resolved. The thickness of the dissociated shell corresponds to the depth $\delta_{\textsc{hi}}$ of the H \scriptstyle\rm I\/${\mathrm{H}_{2}}$ transition in a PDR. An analytic approximation for the column density of the transition ($N_{\mathrm{trans}}=\delta_{\textsc{hi}}n_{\textsc{h}}$) is given by the expressions (Bialy & Sternberg, 2016)

[TABLE]

[TABLE]

This expression is obtained by considering the ${\mathrm{H}_{2}}$ formation-dissociation balance in a slab of gas with constant density, accounting for ${\mathrm{H}_{2}}$ self-shielding (Sternberg et al., 2014) and dust absorption, without the effect of cosmic rays.

In Fig. 4, we have plotted $\delta_{\textsc{hi}}$ as a function of the FUV flux and the gas density. The SED of the impinging radiation is that of a $10^{4}\,{\rm L}_{\odot}$ star, with a distance scaled to obtain the flux $G_{0}$ shown in the $x$ axis of the plot. The red line marks the contour corresponding to the maximum resolution of our simulations, i.e. $\Delta x=0.023$ pc. Hence, referring to Tab. 2 for the values of the mean density, $\langle n\rangle$, we can see that the effect of radiation on clumps with mass larger than 50 ${\rm M}_{\odot}$ is well resolved for every value of $G_{0}$ of interest. Simulations with $M=10,50{\rm M}_{\odot}$ are close to the resolution limit if $G_{0}<10^{4}$. Keep in mind, though, that due to the imposed profile Eq. 5, the density in the outer regions is smaller than $\langle n\rangle$.

3.4 Set of simulations

We run in total seven 3D simulations of photoevaporating dense molecular clumps (Tab. 2). The first simulation of a 50 ${\rm M}_{\odot}$ clump does not include radiation and it is used for comparison. Then we run a set of simulations of clumps with the same mass ($M=50\,{\rm M}_{\odot}$) and different intensities of the radiation field ($G_{0}=2\times 10^{3}-8\times 10^{4}$). Finally, in the last set of simulations, we consider clumps with masses varying in the range $M=10-200\,M_{\odot}$ and constant impinging radiation ($G_{0}\simeq 3\times 10^{4}$).

4 Fiducial clumps (50 ${\rm M}_{\odot}$)

We study first the evolution of a fiducial clump with initial mass of $50\,M_{\odot}$. The clump has an initial radius $R=0.5$ pc and a profile described by Eq. 5 with a core density $n_{c}\simeq 6.6\times 10^{3}\,\,{\rm{cm^{-3}}}$. We run a simulation with no external radiation source, to be used as a control case. We then perform three simulations introducing a nearly isotropic radiation field with different intensities $G_{0}$.

The clump radius is defined as the distance from the centre where the 99.7% of the molecular gas mass is enclosed. In the run with no radiation, all the box is fully molecular, so we had to adopt a different definition of radius, namely where the ${\mathrm{H}_{2}}$ density reaches 10% of its maximum value at the centre.

4.1 Run with no radiation

In Fig. 5, the clump without radiation sources is shown with grey lines: the plot tracks the evolution of the radius (left-hand panel) and the central density (right-hand panel). The radius decreases by a factor of 100 during the collapse, leading to an increase of the maximum core density by 3 orders of magnitude. The clump reaches the minimum radius slightly before the free-fall time $t_{\mathrm{ff}}$ (shown by the vertical dotted line in the right-hand panel), and the maximum compression state is indicated by the grey circle. As opposed to works on protostellar accretion (Larson, 1969; Shu, 1977), we do not follow the clump evolution after the collapse, since our aim here is only to compare the implosion phase with that of clumps exposed to radiation.

4.2 Runs with different impinging flux

We now want to compare the gravity-only simulation with the simulations where radiation is impinging on the clump surface. For the three simulations clump_M50_G2e3, clump_M50_G3e4 and clump_M50_G8e4, the Habing flux $G_{0}$ through the clump surface is $\simeq 2\times 10^{3}$, $3\times 10^{4}$, and $8\times 10^{4}$, respectively.

Slices of number density, thermal pressure and radial velocity for the $G_{0}=2\times 10^{3}$ run are shown in Fig. 6, at times $t=0.1$ Myr and $t=0.5$ Myr. In the first snapshot ($t=0.1$ Myr), it is possible to see that a thick layer of gas ($n\sim 500\,\,{\rm{cm^{-3}}}$) is pushed away at high velocity (Mach number $\mathcal{M}=1.5$), due to the underneath high-pressure gas ($P\sim 3\times 10^{5}\,\mathrm{K}\,\,{\rm{cm^{-3}}}$). This is due to radiation heating the clump surface, so that the pressure is increased and drives an expansion of an outer shell of gas. On the other hand, the overpressure (seen as the highest pressure layer in the first snapshot) also drives a shock inward. The shocked gas is flowing towards the centre at a velocity that is higher than the rest of the clump gas, which is undergoing gravitational collapse. The qualitative behaviour of the gas is the same found in D17, albeit that work considered the effect of ionizing radiation but no gravity. Nevertheless, the effect of the FUV radiation is qualitatively the same: a surface layer of the clump is heated, so that it expands and drives a shock towards the centre of the clump. In the second snapshot, the clump is collapsed to a dense small core, while the photoevaporating flow has reached a higher distance from the centre.

4.3 Radius and density evolution

In Fig. 5, the evolution of clump radius and clump central density are shown for all the four runs. The radius decreases suddenly from the initial value, because of the photodissociation front moving into the cloud as an R-type front: the molecular hydrogen is thus photodissociated and the gas remains almost unperturbed, until the flux is attenuated (both because of photodissociation and dust absorption) and the front stalls. This phase lasts few kyrs (see Eq. 6). After that, the front drives a shock front compressing the gas ahead of it (D-type shock front, Spitzer, 1998), so that the clump shrinks further because of the radiation-driven shock wave. Even if the flux is different in the 3 runs with radiation, the temperature at which the clump surface is heated does not vary much (200-250 K), hence the clump contraction proceeds almost at the same speed. The simulations are then stopped when the clump radius reaches a minimum of \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}2\Delta x_{\mathrm{min}}, i.e. when the clump collapses to a size below the adopted resolution.

The right-hand panel of Fig. 5 shows the clumps central density at different times. For the coloured lines, a circle marks the moment in which the clump has reached the minimum radius (i.e. about the size of a cell), and it corresponds to the higher compression state of the clump. The density reached in the implosion phase is lower when the radiation is stronger. This happens because for the simulations with strong flux the dissociated shell in the R-type phase is thicker, hence the remaining molecular collapsing core is less massive.

4.4 Mass evolution

In order to compare the efficiency of the photoevaporation process with different radiation intensities, we plot in Fig. 7 the total clump mass $M_{\mathrm{tot}}$ (solid line) and the ${\mathrm{H}_{2}}$ mass $M_{\mathrm{H}_{2}}$ (dotted line) as a function of time.

The lines start from the time $t_{\textsc{r}}$ which marks the end of the R-type dissociation front propagation. Subsequently, photodissociation continues at the surface of the clump, generating a neutral flow such that both $M_{\mathrm{tot}}$ and $M_{\mathrm{H}_{2}}$ decrease with time. Notice that the ratio $M_{\mathrm{tot}}/M_{\mathrm{H}_{2}}$ is not constant, since the clump molecular fraction $x_{{\mathrm{H}_{2}}}$ changes with time, even in the interior of the clump.

In particular, $x_{{\mathrm{H}_{2}}}$ is lower just after the R-type phase, when clumps have lower density. Indeed, both the bins [6.0, 11.2] eV (bin a) and [11.2, 13.6] eV (bin b, corresponding to Lyman-Werner band) play a role in the dissociation of ${\mathrm{H}_{2}}$, as we have verified in some test simulations with only one radiation bin: the bin b dissociates the ${\mathrm{H}_{2}}$ through the Solomon process, while the bin a causes an increase of the gas temperature, hence increasing the collisional dissociation of ${\mathrm{H}_{2}}$. While the absorption of radiation in the bin b is very strong because of ${\mathrm{H}_{2}}$ self-shielding, radiation in the bin a is basically absorbed by dust only101010With the absorption cross sections adopted, the optical depth in the bin a is $\tau_{a}=1$ when $N_{\textsc{h}}=7\times 10^{20}\,{\rm cm}^{-2}$, while for LW radiation $\tau_{b}=1$ when $N_{\textsc{h}}=10^{15}\,{\rm cm}^{-2}$ and can penetrate even in the clump interior. At later times, clumps become denser and the gas cools more efficiently. As a result, ${\mathrm{H}_{2}}$ formation is promoted, yielding the maximum value $x_{{\mathrm{H}_{2}}}\simeq 0.76$ within the clump.

4.5 Stability of the molecular core

We have explicitly verified that the molecular core of the clumps is Jeans unstable when the simulations are stopped, implying that the clump could eventually collapse. This conclusion does not depend on the resolution of the simulation, as we have verified by running the simulation clump_M50_G8e4 with $2\times$ and $4\times$ the standard resolution (see App. C). Hence, the final clump mass is an upper limit to the final stellar mass $M_{\star}$. For the 50 ${\rm M}_{\odot}$ clumps, we find $M_{\star}\simeq 15\,{\rm M}_{\odot}$ in the lowest flux case ($G_{0}=2\times 10^{3}$) and about $M_{\star}\simeq 3\,{\rm M}_{\odot}$ in the highest flux case ($G_{0}=8\times 10^{4}$). Following the standard Shu classification (Shu et al., 1987; André, 1994; André & Motte, 2000), star formation proceeds through an early accretion phase from the surrounding envelope ($t\simeq 1$ kyr), a late accretion phase via a disk ($t\simeq 0.2-1$ Myr) and a protostellar phase ($t\simeq 10$ Myr). We cannot investigate the effect of photoevaporation during these phases with the resolution and physics considered in our simulation suite, but we expect that the mass that goes into stars will be in general lower than $M_{\star}$ due to photoevaporating gas from the protostellar system.

5 Clumps with different masses

After focusing on clumps with mass $50\,{\rm M}_{\odot}$, we also analyse the effect of photoevaporation on clumps with lower (10 ${\rm M}_{\odot}$) and higher (100 and 200 ${\rm M}_{\odot}$) masses, with an initially impinging flux of $G_{0}\simeq 3\times 10^{4}$ (see Tab. 2).

In the left-hand panel of Fig. 8 the clump radius is plotted as a function of time. Clumps with larger masses take more time to collapse, as they have a lower initial density and hence longer free fall time. The maximum density $n_{\mathrm{max}}$ as a function of time is shown in the right-hand panel, with higher values reached by more massive clumps.

In Fig. 9 (left-hand panel), the total mass and the ${\mathrm{H}_{2}}$ mass in the clumps are plotted. The $10\,{\rm M}_{\odot}$ and $50\,{\rm M}_{\odot}$ have the same trends seen in Sec. 4, with both total and ${\mathrm{H}_{2}}$ mass decreasing during the implosion phase. Instead, the behaviour is different for the $100\,{\rm M}_{\odot}$ and $200\,{\rm M}_{\odot}$ clumps:

• •

$100\,{\rm M}_{\odot}$: $M_{{\mathrm{H}_{2}}}$ increases slightly (by 25%) during the implosion, while $M_{\mathrm{tot}}$ decreases slowly.

• •

$200\,{\rm M}_{\odot}$: $M_{{\mathrm{H}_{2}}}$ increases substantially (by 50%) during the RDI, while $M_{\mathrm{tot}}$ does not decrease considerably and shows a slight raise after 0.4 Myr.

The different behaviour of $M_{{\mathrm{H}_{2}}}$ in the two more massive clumps, before the gravitational collapse, is due to the effect of radiation in the bin [6-11.2] eV. Photons in this band make their way to the clump interior, which is less dense in the centre with respect to smaller clumps, and increase the gas temperature so that the ${\mathrm{H}_{2}}$ decreases in the centre. Nevertheless, when the clump collapses, the gas self-shields from this radiation and ${\mathrm{H}_{2}}$ abundance increases again in the centre. This boost of the ${\mathrm{H}_{2}}$ abundance compensates for the photoevaporative loss.

In the final part of the RDI, the two massive clumps behave differently, with $M_{\rm tot}$ raising for the 200 $M_{\odot}$ clump. This is due to the fact that the escape velocity from the $200\,{\rm M}_{\odot}$ becomes higher than the typical velocity of the outflowing gas. Fig. 9 (right-hand panel) shows the time evolution of the escape velocity $v_{\mathrm{esc}}$ from the surface of the clump, for the four simulations with similar impinging flux. We have marked with a horizontal line the isothermal sound speed $c_{s}$ at $T\simeq 250$ K, which is a typical temperature of the heated atomic surface of the clumps. The circles mark the time when clump collapse. Only the 200 ${\rm M}_{\odot}$ clump has $v_{\mathrm{esc}}>c_{s}$ at the end of the implosion phase, clarifying why this clump can accrete further in spite of the radiation impinging on its surface.

The same considerations of Sec. 4.5 hold: clumps are Jeans unstable at the end of the RDI, hence the remaining mass in the clump is an estimate of the mass going to form stars. Considering that photoevaporation can also remove mass from the protostellar system, such final mass should be regarded as an upper limit to the stellar mass. While the clump mass is generally reduced by photoevaporation during the RDI, due to its self-gravity the most massive clump in our set of simulations ($200\,{\rm M}_{\odot}$) manages to retain its core mass after the R-type propagation of the DF. Hence in this case photoevaporation is completely inefficient in limiting star formation.

6 General picture

The 3D simulations that we have run show that photoevaporating clumps undergo three main evolutionary phases, summarized in Fig. 10:

R-type DF propagation: FUV radiation penetrates the clump as an R-type front, with the clump density structure unaltered;

• -

RDI: the heated atomic shell drives a shock inward, so that the clump implodes;

• -

gravitational collapse: the molecular core is Jeans unstable, so it undergoes a gravitational collapse, with photoevaporation regulating the mass going into stars.

In the following, we analyse the details of these phases, trying to generalize the results of the simulations to a range of clump masses and impinging flux.

6.1 R-type dissociation front propagation

The first phase of the photoevaporative phenomenon is the dissociation of clump molecules by the propagation of the DF as an R-type front. During this phase, a clump shell is converted in atomic form during the DF propagation, without any dynamical effect on the gas. This phase is very short (less than 0.02 Myr), and it is responsible for the sudden decrease of $M_{\mathrm{H}_{2}}$ with respect to its initial value, as we already pointed out in Fig. 7 and Fig. 9.

This phase determines substantially the fate of a clump. In fact, if the remaining ${\mathrm{H}_{2}}$ mass is very small, the clump can be quickly eroded in the following photoevaporative phase. To make a prediction of the ${\mathrm{H}_{2}}$ mass in the clump after the R-type DF propagation, we have run a set of 1600 1D simulations, with clump masses varying between 6 and $10^{3}\,{\rm M}_{\odot}$ and FUV fluxes in the range $G_{0}=10^{2-6}$. As initial setup of the 1D simulations, we have considered a 1D stencil passing through the centre of a three-dimensional clump. Hence, radiation is injected from both sides of this 1D box, with the prescribed flux. The results are shown in Fig. 11 (leftmost panel), where the colours mark the ratio of molecular mass after the R-type DF propagation to the initial molecular mass ($f_{M_{{\mathrm{H}_{2}}}}$). With respect to the results of the 3D simulations, the 1D simulations slightly underestimate $f_{M_{{\mathrm{H}_{2}}}}$, because of the different geometry. The plot shows that the fraction $f_{M_{{\mathrm{H}_{2}}}}$ depends only weakly on the initial clump mass (apart from the low mass clumps, $M<30{\rm M}_{\odot}$), so that a relation between $f_{M_{{\mathrm{H}_{2}}}}$ and $G_{0}$ can be derived by fitting the data:

[TABLE]

We also notice that for the range of masses and fluxes that we have chosen, no clump is suddenly dissociated by the FUV radiation, hence all the analysed clumps will eventually undergo the implosion phase.

We investigate the dependence of gas metallicity on the efficiency of photoevaporation. We run two additional sets of 1D simulations ($2\times 1600$) with lower values of $Z$, i.e. $0.5Z_{\odot}$, $0.2Z_{\odot}$. The results for $f_{M_{{\mathrm{H}_{2}}}}$ are shown in the central and right-hand panel of Fig. 11. Since the dust abundance is assumed to be proportional to $Z$, the low-$Z$ gas is more transparent to radiation, and molecular gas is dissociated more efficiently. Clumps are fully photodissociated by the R-type DF for $Z=0.5Z_{\odot}$ ($Z=0.2Z_{\odot}$) when $G_{0}>3\times 10^{4}$ ($G_{0}>3\times 10^{3}$). This finding agrees with results from Vallini et al. (2017) and Nakatani & Yoshida (2018), both pointing out that photoevaporation is more rapid in metal-poor clouds, i.e. the ones expect in high redshift galaxies (e.g. Pallottini et al., 2017; Pallottini et al., 2019).

6.2 Radiation-driven implosion

Clumps that are not completely dissociated by the propagation of the DF, will then attain a configuration with a molecular core surrounded by an atomic heated shell (this is the case for all clumps considered in this work). Since the latter has a higher pressure with respect to the molecular core, a shock propagates inward compressing the clump. This phase is generally called radiation-driven implosion.

In D17 we have developed an analytical solution for the propagation of the shock towards the centre of the clump. In that work, the shock parameters (as Mach number and compression factor) are computed as a result of the discontinuity between the cold molecular core and the heated atomic shell. The backup pressure from the atomic shell is kept constant during the evolution, which is equivalent to assume a constant heating of the shell by radiation. Hence, the shock velocity $v_{\mathrm{shock}}$ changes only because of spherical convergence, causing an increase with radius as $v_{\mathrm{shock}}\sim r^{-0.394}$. Nevertheless, this work shows that radiation is absorbed in the photoevaporative flow, hence the heating of the atomic gas is reduced and its pressure decreases accordingly. This effect acts in the opposite direction than the spherical convergence, reducing the shock speed.

In Fig. 12, the shock radial distance from the centre is plotted as a function of time. The simulation clump_M10_G3e4 has been excluded because the resolution is too low to resolve the shock position properly. The initial position of the shock is determined by the clump radius at the end of the DF propagation (cfr. left-hand panel of Figs 5 and 8). The plot shows that the shock speed is almost constant in time, implying that there is a balance between absorption of radiation which is backing up the shock, and spherical convergence of the shock. The simulations show that the shock Mach number (with respect to the gas ahead of the shock) is $\mathcal{M}=2.0\pm 0.3$, with a shock speed about half of the speed of sound in the heated clump shell ($c_{\textsc{pdr}}$). This is in contrast to Gorti & Hollenbach (2002) assumption, where the shock speed is approximated to be exactly $c_{\textsc{pdr}}$.

Finally, we notice that the RDI lasts more than the time needed for the radiation-driven shock to reach the centre of the clump, in contrast with Gorti & Hollenbach (2002) assumption that the clump stops contracting when the shock has reached the centre. This is evident especially for massive and larger clumps, and it is due to the fact that the shock moves towards the centre with a higher speed than the collapsing clump surface. Thus, when the shock has reached the centre, the surface is still moving inward.

During the RDI, two effects changing the ${\mathrm{H}_{2}}$ mass are competing: (1) the clumps loses mass from its surface, (2) the central density increases, raising the ${\mathrm{H}_{2}}$ abundance. The second effect dominates for clumps massive enough ($M\geq 100\,{\rm M}_{\odot}$). Furthermore, self-gravity inhibits the photoevaporative flow of even more massive clumps ($M\geq 200\,{\rm M}_{\odot}$), showing that in this case photoevaporation is not effective in reducing the total clump mass during the RDI.

6.3 Gravitational collapse

After the RDI, the molecular core is still Jeans unstable, thus we expect it to undergo a gravitational collapse with possible star formation (if the molecular mass of the core is sufficient). We do not include star formation routines in our simulations, as done for instance via seeding of a protostellar object and by following its accretion (Dale et al., 2007; Peters et al., 2010).

Other works on photoevaporation (Gorti & Hollenbach, 2002; Decataldo et al., 2017) do not include clump self-gravity in their analysis. As a result, in those works the clump reaches a minimum radius where thermal pressure balances the pressure of the shell heated by radiation. In this scenario, gas continues to photoevaporate from its surface, until all the molecular gas is dissociated. This is not seen in our simulations, since gravity leads to the clump self-collapse after the RDI.

The clump mass at the end of the R-type DF propagation (Fig. 11) is in general an upper limit to the mass $M_{\star}$ going into stars. In fact during the RDI, a fraction of the mass flows away from the clump. However, this does not happen for clumps with sufficiently large masses in which self-gravity prevents the gas from escaping the clump.

7 Conclusions

We have studied the photoevaporation of Jeans unstable clumps by Far-Ultraviolet (FUV) radiation, by running 3D RT simulations including a full chemical network which tracks the formation and dissociation of ${\mathrm{H}_{2}}$.

The simulations have been run with the adaptive mesh refinement code ramses (Teyssier, 2002) by using the ramses-rt module (Rosdahl et al., 2013), in order to perform momentum-based on-the-fly RT. The RT module has been coupled with the non-equilibrium chemical network generated with krome (Grassi et al., 2014), in order to consider photo-chemical reactions, as the dissociation of ${\mathrm{H}_{2}}$ via the two-step Solomon process. We have run seven simulations of dense clumps, embedded in a low density medium ($n=100\,\,{\rm{cm^{-3}}}$), with different clump masses ($M=10-200\,{\rm M}_{\odot}$) and different impinging FUV radiation fields ($G_{0}=2\times 10^{3}-8\times 10^{4}$). These clumps have central number densities $n_{c}\simeq 6\times 10^{3}-2\times 10^{4}$ and total column density $N_{\textsc{h}_{2}}\simeq 5\times 10^{21}{\rm cm}^{-2}$.

In all the cases, we find that the evolution a clump follows three phases:

1. 1)

R-type DF propagation: the density profile remains unaltered while most of the ${\mathrm{H}_{2}}$ mass is dissociated ($40-90\%$ of the ${\mathrm{H}_{2}}$ mass, depending on $G_{0}$). 2. 2)

RDI: the heated shell drives a shock inward ($\mathcal{M}\simeq 2$) promoting the clump implosion; at the same time, the heated gas at the surface evaporates with typical speed $1.5-2\,\mathrm{km}\,\mathrm{s}^{-1}$. 3. 3)

Gravitational collapse of the core: the clump collapses if the remaining ${\mathrm{H}_{2}}$ core is Jeans unstable after the RDI; if $M_{\mathrm{H}_{2}}$ is significantly higher than 1 ${\rm M}_{\odot}$, than we expect it to form stars.

During the RDI, both the molecular mass $M_{\textsc{h}_{2}}$ and the total mass $M$ decrease for the 10 and $50\,{\rm M}_{\odot}$ clumps. However, we find that $M_{\textsc{h}_{2}}$ tends to increase by $\sim 25-50\%$ for more massive clumps, due to the fact that previously dissociated ${\mathrm{H}_{2}}$ recombines when the clump collapses and the density increases. For the most massive clump only ($200\,{\rm M}_{\odot}$), photoevaporation is inefficient even in reducing the total mass $M$, since during the RDI the escape velocity becomes larger than the outflowing gas speed (comparable to the $\mathrm{HI}$ sound speed).

All the ${\mathrm{H}_{2}}$ cores are still Jeans unstable after the RDI. This shows that FUV radiative feedback is not able to prevent the gravitational collapse, although it regulates the remaining molecular gas mass. All the simulated clumps manage to retain a mass $M>2.5\,{\rm M}_{\odot}$, hence suggesting that star formation may indeed take place. The evolution of low-mass clumps follows what expected from analytical works (Bertoldi, 1989; Gorti & Hollenbach, 2002; Decataldo et al., 2017). However, our analysis clarifies that self-gravity has a non negligible effect for massive clumps ($\apprge 100{\rm M}_{\odot}$), limiting the mass loss by photoevaporation.

The dynamics of photoevaporating clumps can also have important consequences for their Far-Infrared (FIR) line emission (Vallini et al., 2017), [C $\scriptstyle\rm II$] in particular. In fact, a strong $G_{0}$ increases the maximum [C $\scriptstyle\rm II$] luminosity, as FUV radiation ionizes carbon in PDRs, albeit short-lived clumps may contribute less significantly to the parent GMC emission. To understand the effect of photoevaporation on line luminosity, simulations accounting for the internal structure of GMCs are required. This would allow to track the contribution of many clumps with different masses and subject to different radiation fields. We will address this study in a forthcoming paper.

Finally, we point out that photoevaporation is also a crucial effect in regulating the molecular mass in ultra-fast outflows launched from active galaxies. Indeed, dense molecular gas ($n\sim 10^{4-6}\,\,{\rm{cm^{-3}}}$) is observed up to kpc scale (Cicone et al., 2014; Bischetti et al., 2018; Fluetsch et al., 2019), and its origin and fate are currently under investigation (Ferrara & Scannapieco, 2016; Decataldo et al., 2017; Scannapieco, 2017; Richings & Faucher-Giguère, 2018). In addition to FUV radiation, these clumps are also subject to a strong EUV field (not considered in this work), causing a fast photoevaporative flow of ionized gas and a stronger RDI. Since radiation is coming from the nuclear region of the active galaxy, the radiation field seen by the clump is non-isotropic, affecting only the side of the clump facing the source. Moreover, clump masses are generally larger ($10^{3-4}\,{\rm M}_{\odot}$, see Zubovas & King 2014; Decataldo et al. 2017) and high turbulence is reasonably expected within the outflow. Nevertheless, clumps will still follow a similar evolutionary path, with R-type dissociation/ionization, RDI, and gravitational collapse if the imploded core is Jeans unstable. This can be relevant to explain the observed molecular outflow sizes, which depends on the lifetime of clumps undergoing photoevaporation. Moreover, in this scenario the recent hints of star formation inside the outflow (Maiolino et al., 2017; Gallagher et al., 2019; Rodríguez del Pino et al., 2019) seem plausible, in view of our finding that most clumps manage to retain a sufficient amount of mass to collapse and form stars.

Acknowledgments

This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the Deutsche Forschungsgemeinschaft (DFG) cluster of excellence “Origin and Structure of the Universe”. We thank A. Lupi, J. Rosdahl, and the participants of the “The Interstellar Medium of High Redshift Galaxies” MIAPP conference for fruitful discussions. The simulations have been run on the UniCredit Research & Development facilities; in particular we thank M. Paris for his technical support. DD and AF acknowledge support from the European Research Council (ERC) Advanced Grant INTERSTELLAR H2020/740120. LV acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant agreement No. 746119. We acknowledge use of the python programming language, astropy (Astropy Collaboration et al., 2013), matplotlib (Hunter, 2007), numpy (Van der Walt et al., 2011), pymses (Labadens et al., 2012).

Appendix A Benchmarks

A.1 H \scriptstyle\rm II\region

As a first test of our version of ramses-rt coupled with krome, we run a simulation of 1D slab of atomic gas invested by ionizing radiation. The test is aimed to check the coupling between radiative transfer and chemistry, hence we run in static mode, i.e. forcing the gas velocity in every cell to zero. The 1D computational domain has a size $L=0.2$ pc, resolved with 1024 cells, and it is filled with neutral hydrogen, at a density of 100 $\,{\rm{cm^{-3}}}$.

The radiation field comes from the left side of the slab, with the spectrum of a black body source at a distance of 1 pc, discretised in 10 bins (delimited by 13.6, 18.76, 23.9, 29.1, 34.2, 39.4, 44.5, 49.8, 54.8, 60.0, 1000.0 eV). The geometrical flux reduction is not considered, while dust cross section is set to $\sigma_{d}=2.94\times 10^{-22}\,\mathrm{cm}^{2}$. We note that dust is not relevant since the column density in the whole slab is lower than $10^{20}\,\mathrm{cm}^{-2}$. The simulations are stopped at 1 kyr, i.e. longer than the recombination time $t_{\mathrm{rec}}\simeq n/\alpha_{\textsc{b}}$.

The result is compared with a semi-analytical solution, obtained by dividing the slab of gas in 1024 cells and proceeding from the first cell to the $n$-th cell as follows:

compute the temperature $T$ in the $n$-th cell, using heating and cooling function provided by Gnedin & Hollon (2012); 2. 2.

compute the ionization fraction $x_{\textsc{}}(T)$ in the $n$-th cell, according to

[TABLE]

where $n$ is the H nuclei number density, $\alpha_{\textsc{b}}(T)$ is the case B recombination coefficient (taken from Abel et al., 1997), $F_{\nu}$ is the specific flux at the cell and $a_{\nu}$ is the photoionization cross section (values in Verner & Ferland, 1996). 3. 3.

compute the total H \scriptstyle\rm I\column density $N_{\textsc{hi}}$ from cell $1$ to $n$ 4. 4.

reduce radiation according to absorption due to $N_{\textsc{hi}}$

The comparison between the simulation and the semi-analytical model is shown in Fig. 13, showing the neutral fraction $x_{\textsc{hi}}$ as a function of the depth into the gas slab, for black bodies sources with different bolometric luminosities ($10^{3}$, $10^{4}$, $10^{5}$ ${\rm L}_{\odot}$). The solid lines shows the results of the simulations, while the dotted lines are obtained with the semi-analytical model. The curves for the $10^{3}$ and $10^{4}\,{\rm L}_{\odot}$ are in very good agreement, while there is a small difference for the the $10^{5}\,{\rm L}_{\odot}$ source near the Strömgren depth $\delta_{\rm Str}$. We have verified that this is due to slightly different gas temperatures around $\delta_{\rm Str}$, which in turn implies a difference in the recombination coefficient. Indeed, krome and Gnedin & Hollon (2012) do not use the same coefficients for computing heating and cooling functions.

A.2 PDR

In order to test the chemistry of molecular hydrogen, we have carried out a 1D simulation of a slab of molecular gas, with radiation coming from one side. The setup is the same of typical simulations carried out with 1D time-independent PDR codes, as kosma-$\tau$ (Stoerzer et al., 1996), cloudy (Ferland et al., 1998), ucl_pdr (Bell et al., 2005) and many others. A PDR-code comparison has been done by Röllig et al. (2007), that benchmarks the results obtained with the different codes starting from the same setup conditions.

In particular, here we compare the results of our code with the models V1 and V2 in the Röllig work. In the V1 and V2 models, the gas filling the computational domain has an H density $n_{\mathrm{H}}=10^{3}\,\,{\rm{cm^{-3}}}$, with He elemental abundance $A_{\mathrm{He}}=0.1$. The flux is coming from the left side, and it is given by $\chi$ times the Draine spectrum (Draine, 1978) in the FUV (6.0 - 13.6 eV); the V1 model has $\chi=10$, while V2 has $\chi=10^{5}$.

The dust cross section is $\sigma_{d}=1.75\times 10^{-21}\,\mathrm{cm}^{2}$ and does not depend on frequency, the cosmic ray H ionization rate is $\zeta_{\mathrm{H}}=5\times 10^{-17}\,\mathrm{s}^{-1}$ and the ${\mathrm{H}_{2}}$ formation and dissociation rates are set to $R_{\mathrm{H}_{2}}=3\times 10^{-18}T^{1/2}$ and $\Gamma_{\mathrm{H}_{2}}=5.10\times 10^{-11}\,\chi)$ respectively. Note that in Röllig et al. (2007), such prescription on the dissociation rate is adopted only for the codes that do not compute explicitly it by summing over all oscillator strengths.

We prepare our simulation to match the setup of the V1 and V2 models in the PDR benchmark, in particular by setting the same ${\mathrm{H}_{2}}$ formation and dissociation rate of molecular hydrogen. The main drawbacks with respect to PDR codes are that (1) we use only the two FUV bins [6.0-11.2] eV and [11.2-13.6] eV to sample the Draine spectrum, and (2) we do not track C, O and Si as separate species.

The results of the comparison are displayed in Fig. 14, for model V1 (upper panel) and model V2 (lower panel). The ${\mathrm{H}_{2}}$ number density and the $\mathrm{HI}$ number density are represented with a blue and a red line respectively; with solid lines we plot the result from our ramses-rt + krome simulation; with dashed lines we plot results111111Data from Röllig et al. (2007) is available at http://www.ph1.uni-koeln.de/pdr-comparison. for the cloudy test. The depth into the gas slab in the x axis is measured as a visual extinction $A_{V}$, which is proportional to the H column density, according to $A_{V}=6.289\times 10^{-22}\,N_{H}$. We obtain a good agreement with cloudy, matching the $\mathrm{HI}$ and ${\mathrm{H}_{2}}$ abundances at low/high $A_{V}$. The small difference in the residual ${\mathrm{H}_{2}}$ abundance in the V2 test is due to a difference in the temperature, affecting the ${\mathrm{H}_{2}}$ formation rate. The position of the $\mathrm{HI}-{\mathrm{H}_{2}}$ transition is well captured in the simulations, with a steeper transition profile due to the different ${\mathrm{H}_{2}}$ self-shielding recipe.

Appendix B Opposite colliding beams problem

Radiative transfer codes based on the M1 closure relations give unexpected results when radiation beams that travel in opposite directions collide (Rosdahl et al., 2013). The flux in the computational domain can be different to what expected by summing the flux of the single beams, showing an excess in the direction perpendicular to the beams.

With the purpose of testing the behaviour of the radiation field for the setup of our simulations, we have run a simulation with resolution $128^{3}$, with 50 stars with luminosity $10^{4}\,{\rm L}_{\odot}$ at a distance of 1.5 pc from the box centre and constant gas density inside the whole domain ($n=10^{-3}\,\,{\rm{cm^{-3}}}$). Dust is not included and gas absorption is negligible at this density, so the Habing flux can be easily calculated analytically for comparison: the flux in every cell is computed by summing the contribution to the flux from all sources, by scaling it with the square root of the distance.

The result after a time sufficient for radiation to cross the whole domain is showed Fig. 15, where flux in the FUV band (left-hand panel) is compared with the analytical solution for the same setup (middle panel). From the map of the relative difference between the two (right-hand panel), we observe that there is a factor of 2 at most. Nevertheless, in the region inside the spherical shell where stars reside, the flux is rather homogeneous. Hence, despite the opposite colliding beams problem, we can carry out our simulations accounting for the fact that the flux on the clump is higher than what expected from an analytical prescription.

Appendix C Convergence test

The standard resolution of our simulation suite is $\Delta x\simeq 0.047$ pc, refining up to $\Delta x\simeq 0.023$ pc according to a ${\mathrm{H}_{2}}$ density gradient criterion. Due to radiation-driven implosion and gravitational instability, all the clump mass collapses to one cell at some time $t_{c}$ which depends on the initial mass and the intensity of the radiation field. We stop our simulations at $t_{c}$, since the later time evolution of the clumps would be unresolved.

We argue that – despite the resolution limit – $t_{c}$ is a good estimate of radiation-driven implosion duration and, being the final clump Jeans unstable. Its mass at $t_{c}$ represents an upper limit to the gas mass that collapses and forms stars. To show this point, we perform a convergence test: we take the $50\,{\rm M}_{\odot}$ (with $G_{0}=8\times 10^{4}$) clump as a reference and we run the two additional simulations: (1) we increase the resolution by a factor 2 in a central region with radius 0.1 pc; (2) we add a further level of refinement in a region with radius 0.05 pc.

Fig. 16 shows the time evolution of the maximum density (upper panel) and the clump mass (lower panel) for the runs with different resolution. Simulations are stopped at the time $t_{c}$ when the clump reaches the minimum resolvable size (1-2 cells): $t_{c}$ is roughly the same in the three cases. In fact, by increasing the resolution, we can resolve the clump collapse to a later stage, but since the timescale for collapse becomes shorter, the clumps reach the minimum size in almost the same time. The lower panel confirms that the final clump mass is not affected by numerical resolution, implying that the estimate of the clump mass at the end of the RDI is reliable.