ULTRAVIOLET ESCAPE FRACTIONS FROM GIANT MOLECULAR CLOUDS DURING EARLY CLUSTER FORMATION

The escape of UV photons from massive stars in young star clusters within molecular clouds drives many critical processes in the interstellar medium (ISM) and intergalactic medium (IGM). The radiation released by stars contributes to the interstellar radiation field (ISRF), which has the highest energy densities at optical and UV wavelengths (Draine 2011), the strength of which was first estimated by Habing (1968) to be ∼4 × 10⁻¹⁴ erg cm⁻³ for 12.4 eV photons. Later authors have further characterized the strength of the UV portion of the ISRF by including wavelength dependence (Draine 1978; Mathis et al. 1983).

The ISRF, as well as its interactions with gas and dust, is responsible for determining the chemical, thermal, and ionization state of the ISM via photoionization, photodissociation, photoelectric heating, and absorption and re-emission by dust grains (Draine 2011). Since most UV photons are generated by massive stars in the range 10–100 M_⊙, they contribute significantly to the strength of the ISFR and significantly alter the state of the ISM in their vicinity, even when considering their short lifetimes.

It has also become clear in recent years that UV ionizing photons from galaxies hosting active galactic nuclei (AGNs) are not sufficient to completely reionize the IGM by z = 6 (Fan et al. 2006; Robertson et al. 2013). Instead, fainter dwarf galaxies, with masses as low as ∼10⁸ M_⊙, are needed to provide the remaining UV photons via their stellar content (Wise et al. 2014; Xu et al. 2015). These low-mass galaxies may contribute up to ∼40% of the total ionizing photons required for reionization (Wise et al. 2014).

In order to contribute to reionization, ionizing photons produced in these galaxies must escape into the IGM (Robertson et al. 2010). The exact fraction of photons, f_esc, that escape their host galaxies, however, is a debated topic. For bright, high-redshift galaxies, measured via the Lyman continuum, f_esc ∼ 7% (Siana et al. 2015), but this number can be as high as ∼30% for fainter Lyα-emitting galaxies (Nestor et al. 2013). Estimates of f_esc from the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC) based on H ii region mapping suggest global escape fractions of 4% and 11%, respectively (Pellegrini et al. 2012).

Simulations that attempt to quantify f_esc for both high- and low-mass galaxies have been performed, but these results often vary by orders of magnitude. For example, Paardekooper et al. (2011) found f_esc < 1% for high-redshift dwarf galaxies, while later numerical works have found f_esc > 10% (Razoumov & Sommer-Larsen 2010; Ferrara & Loeb 2013; Paardekooper et al. 2015). Moreover, f_esc can vary by orders of magnitude over the lifetime of the galaxy (Paardekooper et al. 2011).

As illustrated by the numerical simulations in Paardekooper et al. (2011), the distribution of dense gas in star-forming regions is one of the main constraints on f_esc from a galaxy. This suggests that detailed modeling of f_esc from dense regions within galaxies is required to fully understand the trends observed in more global simulations. Giant molecular clouds (GMCs) are the densest regions of the galactic ISM, and they are the sites where all known star formation takes place. Studying the escape of UV photons from GMCs is therefore also important for a better understanding of cosmological reionization.

The GMC environment is complex, consiting of filaments produced by supersonic turbulence out of which stars and clusters ultimately form (Bertoldi & McKee 1992; Lada & Lada 2003; Mac Low & Klessen 2004; McKee & Ostriker 2007; André et al. 2014; p. 27, Klessen & Glover 2016). Stars that form in this environment can then alter their surroundings via the emission of radiation, producing H ii regions. The complexity of this problem necessitates the use of numerical simulations. While simulations of GMCs that include star formation and radiative transfer have been completed (Dale et al. 2005; Krumholz et al. 2010; Murray et al. 2010; Peters et al. 2010a; Bate 2012; Klassen et al. 2012b; Walch et al. 2013), these studies do not examine the fraction of photons that escape the cloud.

In this paper, we address the critical question of UV escape fractions from turbulent molecular clouds by computing f_esc from 10⁶ M_⊙ GMCs. We employ our suite of simulations that simulated star cluster formation and radiative feedback within young, 10⁶ M_⊙ GMCs that have varying initial virial parameters (Howard et al. 2016). We model the early evolution of star clusters, defined here as less than 5 Myr, since the effects of supernovae are not included. We found that, despite producing large H ii regions, the inclusion of radiative feedback only suppressed the formation of clusters by a few percent. In comparison, varying the initial virial parameter from 0.5 to 5 (i.e., bound to unbound) reduced the efficiency of cluster formation by ∼34%. The high final star formation efficiencies, which range from 18% to 34%, suggest that radiative feedback alone is not responsible for limiting star formation but that initially unbound clouds better reproduce locally observed GMCs.

Given that we have computed the structure and dynamics of cluster-forming clouds undergoing radiative feedback, we can now address the question of what fraction of the UV photons produced by the massive stars in clusters escapes the molecular cloud.

We present maps of the ionizing photon flux escaping the cloud to demontrate its highly nonuniform nature in space. We also present f_esc (used hereafter to represent the escape fraction from a GMC) during the first 4 Myr of the GMC's evolution, which is shown to be highly variable in time and peaks at ∼35% with a long-term average value of ∼15%. The variable nature of f_esc is attributed to H ii regions that dramatically vary in both shape and size due to the dynamical nature of the gas and embedded clusters.

Below, we provide a brief description of our numerical methods and subgrid model for star cluster formation.

We have simulated a 10⁶ M_⊙ GMC using the adaptive mesh refinement (AMR) code FLASH (Fryxell et al. 2000), which includes self-gravity, radiative transfer, star cluster formation, and cooling processes (see Howard et al. 2016 for more detail). This cloud mass was chosen in particular because high-mass GMCs contain most of the molecular mass in the Milky Way and are host to the most massive stellar clusters (Mac Low & Klessen 2004; McKee & Ostriker 2007; Klessen & Glover 2016).

The cloud is initially overlaid with a turbulent velocity field that is composed of a mixture of solenoidal and compressive turbulence with a Burgers spectrum (as in Girichidis et al. 2011). We selected a configuration with an initial virial paramater of 3, corresponding to an initial Mach number of 73. We chose this simulation in particular out of the suite presented in Howard et al. (2016) because we found that initially unbound clouds best reproduce the properties of massive GMCs in the Milky Way. The radius of the cloud is 33.8 pc. The initial average density of the GMC is n = 100 cm⁻³, with a density profile that is uniform in the inner half of the cloud and decreases as r^−3/2 in the outer half.

The package PARAMESH is used for the adaptive mesh portion of FLASH (Fryxell et al. 2000). The grid is refined at locations with sharp density or temperature contrasts to improve the resolution near filaments and H ii regions. The minimum cell size in our simulation is 0.13 pc.

To model gas cooling, we employ the method from Banerjee et al. (2006), which treats cooling via molecular line emission, gas–dust interactions, H₂ dissociation, and radiative diffusion in the optically thick limit. The cooling rates from Neufeld et al. (1995) are used to treat molecular line emission, while the treatment in Goldsmith (2001) cools the gas via gas–dust transfer.

Radiative transfer is treated via a hybrid-characteristics ray tracer developed by Rijkhorst et al. (2006) and adapted for astrophysical use by Peters et al. (2010a). This scheme treats both ionizing and nonionizing radiation and makes use of the DORIC package (Mellema & Lundqvist 2002) to solve the ionization equations. While the DORIC package is capable of treating a large number of species, we consider hydrogen to be the only gas component for simplicity. The flux of ionizing photons, F_*, from an individual source is given by

$$\begin{eqnarray}&&{F}_{* }=\displaystyle \frac{{S}_{* }}{4\pi {r}^{2}}{e}^{-\tau },\end{eqnarray} \tag{ 1 }$$

where S_* is the cluster's ionizing photon rate in s⁻¹, r is the distance between the source and cell of interest, and τ is the intervening optical depth. The opacity to nonionizing radiation is represented by the Planck mean opacities from Pollack et al. (1994), which are used because the ray tracer has no frequency dependence. We adopt a single UV opacity in neutral gas of κ = 775 cm² s⁻¹ from Li & Draine (2001). This opacity is scaled by the neutral fraction of the gas, so completely ionized regions have an opacity of zero.

We make use of sink particles (Federrath et al. 2010) to model star cluster formation with a custom subgrid model to represent star formation within the clusters (Howard et al. 2014). We adopt a threshold density for formation of 10⁴ cm⁻³, which is based on observations of star-forming clumps (Lada & Lada 2003). Our subgrid model within cluster sink particles (henceforth referred to as clusters) divides the cluster mass into two types: stars, and the remaining gas mass (denoted as the reservoir). We convert the reservoir to stars by randomly distributing the mass into main-sequence stars via a Chabrier (2005) initial mass function (IMF) with an efficiency of 20% per free-fall time, where the free-fall time is taken to be 0.36 Myr. The IMF is sampled every 1/10 of a free-fall time to allow cluster propeties to evolve smoothly over time. Newly accreted gas is added to the reservoir (i.e., gas that is available for star formation during the next IMF sampling step). The masses of all stars formed in the cluster are recorded, and analytical fits provided by Tout et al. (1996) are used to determine each star's total and ionizing luminosity. The cluster's luminosity is then the sum of its consituents, which is then used by the ray tracer.

In order to reduce the computational time, we apply a mass threshold of 1000 M_⊙ in stars (which typically have ∼1 O star), below which clusters do not radiate. Clusters below the threshold continue to accrete gas and form new stars, but they are not included in the radiative transfer calculation.

To study the spatial distribution of the escaping UV flux from the cloud, we produce maps of the ionizing flux across a spherical surface, which are presented in Figure 1. The radius of this sphere corresponds to the intial cloud radius of 33.8 pc, and all clusters are contained within the surface. The maps were made using a Hammer projection, which was chosen because it is an equal-area projection. We also include the locations of the 10 most luminous clusters (accounting for 93% of the final ionizing luminosity), projected to the closest location on the sphere, in white circles. Note that the clusters are not actually located on this spherical surface, but are contained within its volume.

Zoom In Zoom Out Reset image size

The first panel, plotted at 1.5 Myr, shows the ionizing flux shortly after the first clusters begin to radiate. A large fraction of the surface is not receiving any ionizing photons, shown by the white patches. This is because at this time, the clusters have only recently formed (meaning that their total ionizing luminosity is low compared to their final values).

In the same panel, the regions that are receiving ionizing photons are concentrated in the upper right quadrant. Note that the most luminous clusters appear in a grouping toward the right side as well, suggesting that these clusters are responsible for much of the emission observed outside the cloud. There is also some flux associated with the cluster in the bottom left quadrant of this panel.

At 2.5 Myr, we see that the entire surface is now being traversed by UV photons from the clusters. The flux of photons, however, is not spatially uniform. Since the flux on the sphere's surface depends on the intervening column density, the presence of dense clumps and filaments manifests itself as regions with lower flux. We note that the simulation has virialized (α = 1) at 2.5 Myr, so any further turbulence is driven by gravitational collapse of the gas (see Howard et al. (2016) for details and Klessen & Hennebelle (2010) for a more general discussion of accretion-driven turbulence).

As the total ionizing luminosity increases and the total mass in gas decreases, the presence of these dark filaments becomes less pronounced. At 3.18 Myr, only the left side of Figure 1 shows regions with low flux. The grouping of clusters on the right of this figure is likely responsible for the higher flux in that region. From 3.75 Myr onward, the flux is more spatially uniform due to increased cluster luminosities and lower total gas mass.

The above visualizations show that the ionizing flux can vary significantly over both space and time within a GMC.

In Figure 2, we focus on the evolution of f_esc from the cloud. We define f_esc as the total number of photons crossing the spherical surface previously discussed in Figure 1, divided by the summed total of all photons being generated by the clusters. Its time evolution is shown in the top left panel of Figure 2. Note that we only include the clusters that are above the mass threshold discussed in Section 2, since these are the clusters that are used by the radiative transfer scheme.

Zoom In Zoom Out Reset image size

The total escape fraction remains low at approximately 3% between 1.5 and 2.5 Myr. After 2.5 Myr, f_esc rises to a peak of 30% at 3.25 Myr, followed by a sudden drop. The escape fraction begins to rise again, reaching a peak of 37% at 3.8 Myr. The average f_esc from the first rise at 2.5 Myr to the end of the simulation, shown by the horizontal line, is 15%.

The rising f_esc and subsequent rapid drops are not due to changes in the ionizing photon output from the clusters, which is shown in the top right panel of Figure 2. These clusters are accreting new gas vigourously from their surroundings and building new, massive stars as time progresses. The increase in the ionizing photon output is steady and shows no distinct features that correspond to the features seen in f_esc.

Rather, the ionization structure of the gas is responsible for the variable f_esc. In the bottom left panel of Figure 2, we plot the fraction of gas mass that has an ionization fraction of greater than 95%. This figure clearly mirrors the features seen in f_esc, with an increasing f_esc corresponding to an increase in the mass fraction of ionized gas. Recall that UV opacity in ionized regions is significantly lower than in neutral regions.

While the change in the ionized mass fraction is low, peaking at ∼3%, the H ii regions can spatially occupy a significant fraction of the simulation volume, typically filling large voids that are interspersed between dense filaments. Since clusters are the source of radiation and therefore tend to exist in H ii regions, photons can travel large distances due to the reduced opacity, resulting in higher photon fluxes near the boundary of the simulation volume.

However, the size and shape of H ii regions are not constant. Both observations (De Pree et al. 2014, 2015) and simulations (Peters et al. 2010a, 2010b; Galván-Madrid et al. 2011; Klassen et al. 2012a) show that the size of an H ii region can fluctuate on short timescales, a phenomenon described as "flickering." The dynamic and anisotropic nature of the gas, in combination with dynamic clusters, can result in H ii regions becoming shielded to radiation due to changes in density between the source and the ionized regions. The formerly irradiated gas then recombines, causing the H ii to flicker.

A visual demonstration of this flickering is displayed in Figure 3, which shows three-dimensional (3D) images of density (in green) and H ii regions (in red). The green density contours show gas at ∼100 cm⁻³, which is the typical density of the filaments out of which the clusters form. The entire simulation volume is shown, and the box side length is 83 pc.

Zoom In Zoom Out Reset image size

The left panel of Figure 3 shows the state of the simulation at 3.18 Myr, corresponding to the first pronounced peak in f_esc. A large H ii region has developed on one side of the cloud that extends away from the dense, central gas to the boundary of the simulation volume. The middle panel, shown at 3.31 Myr, shows the decrease in the size of the H ii region that is responsible for the deep trough in f_esc at 3.25 Myr. The right panel of Figure 3, plotted at the second peak of f_esc at 3.75 Myr, shows that the H ii region has expanded again to a similar size to that seen in the first panel.

To investigate the cause of the variable H ii region size, we focus our analysis on one luminous cluster that is associated with the H ii region. This cluster is the second most luminous in the simulation with a final ionizing luminosity of 1.40 × 10⁵¹ s⁻¹. The most luminous cluster was not chosen because it is deeply embedded in the dense, central gas and therefore its associated H ii region is small in comparison to the one that extends to the boundary of the simulation volume, as seen in Figure 3.

We drew lines of sight that originate at the cluster's position and extend a distance of 20 pc through the large H ii region. This was done at two times, one just before the H ii region collapses for the first time (at 3.25 Myr) and one immediately after the collapse (∼35,000 yr after the first image). We can then examine how the density and the recombination rate differ before and after the H ii region collapses along these lines of sight.

We find that the radiative recombination rate along the lines of sight increases significantly immediately after the H ii region collapses, increasing from ∼5 × 10⁻⁸ to 1 × 10⁻⁶ cm⁻³ s⁻¹. The radiative recombination rate is given by αn², assuming an ionization degree of 100%, where α is the radiative recombination coefficient and n² is the square of the number density. The recombination coefficient varies with temperature as

$$\begin{eqnarray}&&\alpha =2.59\times {10}^{-13}{\left(\displaystyle \frac{T}{{10}^{4}{\rm{K}}}\right)}^{-0.7},\end{eqnarray} \tag{ 2 }$$

where T is the gas temperature in kelvin.

The radiative recombination rate increases after the collapse for two reasons. First, the density immediately surrounding the cluster increases, likely due to the turbulent nature of the surrounding gas. Second, as the region cools, the recombination coefficient increases.

We also examined the quantity αx²n², where x is the ionization fraction of the gas, which removes the assumption of a 100% ionization fraction. In this case, we see the opposite trend and the recombination rate drops from 5 × 10⁻⁸ to 2 × 10⁻¹³ cm⁻³ s⁻¹. Despite the increase in density and the recombination coefficient, the recombination rates decrease significantly due to the low ionization fraction of the gas after the H ii region collapses.

We find that n² increases by a factor of ∼1.5–4 along the lines of sight within a 1 pc radius of the cluster's location. This increased density limits the amount of radiation that propagates to larger radii. The gas can then recombine and cool from ∼10⁴ K, typical of H ii regions, to ∼10 K, which is the temperature floor adopted in the simulation.

This can be visualized by examining the neutral column density from the luminous cluster through the H ii region, as shown in Figure 4. The top panels show a slice of density (left) and ionization fraction (right) centered on the luminous cluster associated with the H ii region that collapses at 3.25 Myr. These images are plotted before the H ii region collapses. At this time, the cluster is no longer deeply embedded in the massive cold filament out of which it formed in the first place. The bottom panels of Figure 4 show a Hammer projection of the neutral gas column density across a spherical surface of radius 20 pc centered on the same cluster before the H ii region collapse (left) and after the collapse (right). A 20 pc radius circle is shown in the density and ionization fraction slices for reference. The column density projections clearly show that the region that was previously ionized has increased in column density after the H ii region collapses.

Zoom In Zoom Out Reset image size

The increase in density surrounding the massive cluster can be understood through turbulent shocks in the surrounding ionized gas. We measured the local gas velocity dispersion at the location of the cluster immediately before the H ii region collapse to be 10.1 km s⁻¹ (corresponding to a Mach number of 1.14). Thus, a density fluctuation can cross the cluster's radius of 0.78 pc in ∼76,000 yr. This is comparable to the ∼35,000 yr it takes for the H ii region to collapse. A passing shock could therefore lead to a local density enchancement, causing the H ii region to collapse in the observed time. As the gas recombines, it cools and shields regions further along the line of sight.

This further analysis supports our claim that the dynamic nature of the gas, which at this point in time causes the density to increase within the H ii region, is responsible for the strong fluctuations in the size of the H ii regions being produced by luminous clusters, and hence in the value of f_esc.

The above discussion has focused on f_esc from the entire molecular cloud, which is a useful quantity when studying the buildup of the ISRF or estimating the global escape fraction from a galaxy as a whole. We may instead investigate the escape fraction from smaller regions surrounding luminous clusters to follow the evolution of f_esc as a function of distance from the cluster. We are limited in this regard due to the fact that the ray tracer used to compute the radiative transfer only tracks the total flux in each grid cell with no directional information about the incoming rays. This means that if we calculate the flux across a small spherical surface centered on a luminous cluster, it will likely include contributions from sources outside the sphere.

Still, we can calculate f_esc across a surface if the majority of the total ionizing luminosity is being generated within its volume. This minimizes the contribution to the total flux from outside sources. We find that the 10 most luminous clusters generate 90% of the total ionizing luminosity and are located a maximum of 24.3 pc from the simulation center. We therefore repeat the f_esc calculation for a sphere of radius 25 pc instead of 33.8 pc, which was presented in Figure 2.

For reference, we show a two-dimensional projection of the position of all clusters in the left panel of Figure 5. The 10 most luminous clusters are shown by the stars, and all clusters are colored by their Z-position in the cloud. We only show the positions at one time, 3.18 Myr, which corresponds to the first peak of f_esc in Figure 2, to illustrate that the clusters are not strongly grouped together but instead cover the cloud's entire extent.

Zoom In Zoom Out Reset image size

The right panel of Figure 5 shows f_esc across a spherical surface of radius 25 pc. Comparing to the top left panel of Figure 2, we see that f_esc across the smaller surface well within the cloud has a similar temporal evolution with pronounced peaks at 3.25 and 3.75 Myr, but the values are ∼2 times larger. The average f_esc from 2.5 to 4.2 Myr is 41%, compared to 15% for ionizing radiation that escapes from the surface of the cloud. The early evolution (i.e., less than 2.5 Myr) of f_esc is also significantly enhanced, likely due to the generation of small H ii regions surrounding the luminous clusters that are not large enough to extend to the edge of the simulation volume.

Overall, these results suggest that f_esc decreases with cloud radius as one moves through the cloud and out its surface. This trend has been found observationally by Pellegrini et al. (2012), who noted that the global f_esc from the LMC and SMC is estimated to be 4% and 11%, respectively, while f_esc from individual star-forming regions can be as high as ∼60%.

The values presented in Figure 2 are particularly important for the global ISRF and ISM structure since they represent the escape fraction from the surface of an entire 10⁶ M_⊙ GMC. Clouds of this mass are host to the most massive stellar clusters, which dominate the stellar feedback and overall luminosity of a galaxy (Harris & Pudritz 1994; Mac Low & Klessen 2004; McKee & Ostriker 2007; Klessen & Glover 2016).

We computed the UV photon escape fraction from a turbulent, 10⁶ M_⊙ GMC using the astrophysical code FLASH. The cloud, taken from Howard et al. (2016), is initially unbound with a virial parameter of 3, and sink particles are used to model the formation of star clusters. Our simulations end just before supernova explosions could disrupt the star clusters and remove surrounding gas.

Our analysis indicates that the flux is both highly anisotropic, due to the filamentary and clumpy nature of the intervening turbulent gas, and highly variable in time. As time progresses, the flux naturally increases since the clusters contain more massive stars.

The integrated escape fraction, defined as the total number of photons leaving the cloud divided by the total number of photons being produced by all clusters, also varies significantly over time. For the first 2.5 Myr of evolution, f_esc remains low at ∼5%. There are two distinct peaks in the escape fraction at 3.25 and 3.8 Myr, with a maximum escape fraction of 30% and 37% at the two peaks. The average f_esc from the onset of large H ii regions at 2.5 Myr to the end of the simulation is 15%. The average f_esc increases to 41% if we instead consider a smaller surface well inside the cloud at a radius of 25 pc, as compared to the GMC's surface, which is at a radius of 33.8 pc.

The peaks of f_esc and subsequent troughs are tied to the local gas density structure surrounding the luminous clusters, which determines the size of the H ii region they produce. At a peak, the H ii is large and extends toward the boundary of the simulation volume. At a trough, the H ii region is only a small fraction of its previous size. The collapse of the H ii region is tied to the turbulent nature of the gas surrounding the luminous clusters, which causes the density to vary with time. As it increases, so do the recombination rates, and the H ii regions shrink.

We argue that calculations of the photon escape fraction on galactic scales require knowledge of f_esc for individual GMCs and for the star-forming complexes in their interior. For many applications (e.g., cosmic reionization), this is computationally prohibitive, and we suggest to use an average value of f_esc = 15% for a 10⁶ M_⊙ GMC with fluctuations of a factor of two superimposed on timescales of about 1 Myr.

We thank Bill Harris and Eric Pelligrini for interesting discussions.

C.S.H. and R.E.P. thank the the Zentrum für Astronomie der Universität Heidelberg (ZAH), Institut für Theoretische Astrophysik (ITA), and the Max-Planck-Institut für Astronomie (MPIA), for their generous support during R.E.P.'s sabbatical leave (2015/16) and C.S.H.'s extended visit (2015 October–November).

R.E.P. is supported by Discovery Grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada. C.S.H. acknowledges financial support provided by the Natural Sciences and Engineering Research Council (NSERC) through a postgraduate scholarship. The FLASH code was in part developed by the DOE-supported Alliances Center for Astrophysical Thermonuclear Flashes (ASCI) at the University of Chicago. This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada.

R.S.K. acknowledges support from the European Research Council under the European Communitys Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant "STARLIGHT: Formation of the First Stars" (project no. 339177). R.S.K. further acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) in the Collaborative Research Center SFB 881 The Milky Way System (subprojects B1, B2, and B8) and in the Priority Program SPP 1573 Physics of the Interstellar Medium (grant nos. KL 1358/18.1 and KL 1358/19.2).