Lyman–Werner escape fractions from the first galaxies

Abstract

1 INTRODUCTION

In recent years, several highly luminous quasars have been observed, for example, a 1.2 × 10¹⁰ M_⊙ black hole (BH) at redshift z ≈ 6.3 (Wu et al. 2015) and a 9 × 10⁹ M_⊙ BH at redshift z ≈ 7.1 (Mortlock et al. 2011). These objects provide a stern challenge for current structure formation models as it is not clear yet how BHs can assemble so much mass during the first billion years of the Universe (see a review by Volonteri 2010).

In addition, the luminous Ly α emitter CR7 (Sobral et al. 2015), recently discovered at z = 6.6, is found to have a unique structure. A deep Hubble Space Telescope image shows that it consists of three clumps, separated by ≈ 5 kpc. The most luminous clump does not show any metal lines, but the two smaller, accompanying clumps do. CR7 is not only very luminous, but also shows a strong He ii 1640 Å line and therefore a hard spectrum. Sobral et al. (2015) interpreted this object as a chemically pristine protogalaxy forming a large cluster of Population III (Pop III) stars, but other authors have argued that it is more likely to be a massive accreting BH (Pallottini et al. 2015; Agarwal et al. 2016; Dijkstra, Gronke & Sobral 2016; Hartwig et al. 2016, but see Bowler et al. 2016).

Explaining these observations and understanding how supermassive black holes (SMBHs) form and grow at high redshift is a topic of active research. On one hand, SMBH seeds have been suggested to be remnants from Pop III stars (Fryer, Woosley & Heger 2001; Madau & Rees 2001; Inayoshi, Haiman & Ostriker 2016) with relatively low masses (10²–10³ M_⊙) that accrete close to the Eddington limit for large fractions of their lifetime. On the other hand, very massive seeds (10⁴–10⁶ M_⊙) could be created by the direct collapse of primordial gas (Loeb & Rasio 1994; Eisenstein & Loeb 1995; Koushiappas, Bullock & Dekel 2004; Begelman, Volonteri & Rees 2006; Lodato & Natarajan 2006; Latif et al. 2013, see Volonteri 2010 for a review).

A necessary condition for the direct collapse scenario is a very massive and hot (T > 10⁴ K) halo in which gas contracts quasi-isothermally without cooling and fragmenting. The central clump, to which the halo compresses, contains ≥10 per cent of the total gas mass (Bromm & Loeb 2003). The fate of this clump is not entirely clear: It may form a massive BH directly, or may form a supermassive star that subsequently collapses to form a BH. To keep the collapsed gas at high temperatures, the gas must be metal-free (Omukai, Schneider & Haiman 2008), and molecular hydrogen has to be prevented from forming or must be destroyed in that system, as it would otherwise cool a pristine halo down to temperatures of 150 K and encourage fragmentation (Bromm & Loeb 2003). For this, a critical level of external Lyman–Werner (LW) flux is crucial (e.g. Omukai 2001; Shang, Bryan & Haiman 2010; Wolcott-Green, Haiman & Bryan 2011; Agarwal et al. 2012).

Current work (Smidt, Whalen & Johnson, in preparation) favours the direct collapse scenario (but see Alexander & Natarajan 2014). Unlike BH seeds originating from Pop III stars, direct collapse BH (DCBH) seeds are not born starving (Whalen, Abel & Norman 2004), but instead retain their fuel supply (Alvarez, Wise & Abel 2009; Park & Ricotti 2011). Lower mass BH seeds are often kicked out of their host halo (Whalen & Fryer 2012). As the metal-free clump of CR7 is accompanied by two other clumps, the question arises whether this is an observed formation in place of such a DCBH.

With this study, we impose additional constraints to the DCBH seed formation scenario from simulations (see Dijkstra, Ferrara & Mesinger 2014, who show that the number density of DCBH seeds is highly dependent on the LW escape fraction). The questions we address in this work are as follows: How much LW radiation can actually escape an atomically cooling halo and contribute to the LW background? Is shielding by neutral hydrogen an important mechanism to prevent the build-up of the LW background?

We have recently performed a study on LW escape fractions from minihaloes (Schauer et al. 2015, hereafter S15; see also Kitayama et al. 2004), finding the escape fraction varying with halo and stellar mass from 0 to 95 per cent. We adopted two limiting cases: the near-field and the far-field. In the near-field case, we examine what would be seen by an observer in the rest frame of the host halo. This is a good approximation for when the relative velocity is smaller than the typical width of the LW lines. In the far-field case, we consider an observer who is significantly Doppler-shifted with respect to the halo, either due to a large peculiar velocity or simply due to the Hubble flow. This is the relevant case if we are interested in the global build-up of the LW background.

In this new study, we use our existing framework presented in S15 but focus on more massive haloes (10⁷–10⁸ M_⊙), typical hosts of the first galaxies. Here, instead of single stars, we focus on radiation from stellar clusters with time-varying spectral energy distributions (SEDs). We consider several different star formation efficiencies (SFEs) and stellar cluster masses, meant to be descriptive of the first galaxies at z ≈ 10.

Our paper is structured as follows. In Section 2, we briefly summarize the semi-analytic methods described in S15 and introduce our numerical simulations. Details of our parameter space, including the choice of haloes, SEDs and SFEs, are described in Section 3. We present our results in Section 4, where we show the dependence on the LW escape fractions for all halo, SED and SFE combinations over time in both the near- and far-field. In addition, we provide tabulated results with values averaged over the lifetime of stellar populations. Finally, we discuss our findings and conclude in Section 5.

2 METHODS

In this section, we describe the methods used to calculate the LW escape fractions in the near- and far-field limits. The influence of radiation from a central stellar cluster on the gas in a surrounding protogalaxy is modelled in 1D using the zeus-mp hydrodynamical code (Section 2.1). The simulation results are then post-processed to derive the LW escape fraction (Section 2.2).

2.1 zeus-mp

Our simulations are run with the radiation-hydrodynamics code zeus-mp (Whalen & Norman 2006). A full description of the code can be found in more detail in S15. Here, we only give a short overview of the most important features.

The version of zeus-mp used in this study couples photon-conserving ionizing UV transport self-consistently to non-equilibrium primordial chemistry and hydrodynamics to model the growth of cosmological ionization fronts (I-fronts; Whalen & Norman 2006). We include a network of nine chemical species (H, H⁺, He, He⁺, He^{2 +}, H⁻, H|$^{+}_{2}$|⁠, H₂ and e⁻). Heating and cooling processes like photoionization, collisional ionization, excitation and recombination, inverse Compton scattering from the cosmic microwave background, bremsstrahlung emission and H₂ cooling, mostly follow the implementation in Anninos et al. (1997). Details of the radiative reactions can be found in table 1 of Whalen & Norman (2008).

The spectra are tabulated as a function of time, in 120 energy bins with 40 uniform bins from 0.755 to 13.6 eV and 80 uniform logarithmically spaced bins from 13.6 to 90.0 eV. Since single stars were analysed in S15, we used SEDs that were constant in time. However, in this work, as we are focusing on stellar populations, we consider time-varying SEDs as described in Section 3.2.

For the simulation, we assume a spherical symmetry. The volume is divided into 1000 radial cells, reaching out to about twice the virial radius. We use a logarithmically spaced grid, allowing us to have high resolution within the inner regions of the halo without compromising the running time. For 8 out of a total of 30 simulations, we additionally refine a 20-pc region around the I-front with 20 000 cells to achieve numerical convergence. For all setups, we run comparison simulations with half the number of cells. In all cases, the results from the comparison simulations differ by less than 1 per cent from the original simulations, demonstrating that our results have converged numerically.

2.2 Semi-analytical post-processing

After the simulations have run, we calculate escape fractions of LW photons in two limiting cases, the near-field (Section 2.2.1) and the far-field (Section 2.2.2). For a detailed discussion of these two limits, we refer the reader to S15. We summarize the calculations in the sections below.

2.2.1 Near-field

The near-field limit applies in the vicinity of the halo, and is the escape fraction that would be measured by an observer who is at rest with respect to the halo, or moving with only a low velocity relative to it. In this limit, H₂ within the halo can potentially provide effective self-shielding of the LW lines.

Draine & Bertoldi (1996) and Wolcott-Green & Haiman (2011, hereafter WH11) have calculated the effects of H₂ self-shielding in the low-density and local thermodynamic equilibrium limits, respectively. They provide simple fitting functions describing the shielding as a function of H₂ column density and temperature. LW photons that do not dissociate H₂ molecules, or are not absorbed and re-emitted as lower energy photons, can escape the halo. For high column densities, H₂ can shield itself against LW radiation, or can be shielded by H. We can therefore equate the escape fraction to the shielding fraction of molecular, or molecular and neutral hydrogen.

For the optically thick case (N_H2 > 10¹³ cm⁻²), we use the shielding functions from WH11, which depend on the column density N_H2 and, via the Doppler width, on temperature and velocity dispersion of the gas in the halo. We consider two cases: one in which only H₂ contributes, and another where we additionally account for shielding from neutral hydrogen. In the H₂-only case, we use the shielding functions given in equation 8 of WH11. When we account for shielding from atomic hydrogen, we additionally use equations 11 and 12 from the same paper. For small column densities (N_H2 < 10¹³ cm⁻²), we set the escape fraction to 1.

2.2.2 Far-field

The far-field limit corresponds to the escape fraction that would be measured by an observer located far away from the halo, moving with a large peculiar velocity or Hubble flow velocity with respect to it. In this limit, the LW lines in the rest frame of the observer are Doppler-shifted and do not coincide with the LW lines in the halo.

To calculate the escape fraction in this limit, we choose the LW range of 11.2–13.6 eV (corresponding to a wavelength range of 912–1110 Å) and calculate what fraction of the light in this range will be absorbed by H₂ or H in the halo. If we denote this fraction as f_abs, then the far-field escape fraction is simply f_esc = 1 − f_abs.

To calculate f_esc, we first compute the dimensionless equivalent width of the full set of accessible LW lines. To do this, we follow Rodgers & Williams (1974) and calculate the equivalent width of every line, accounting for both Lorentz broadening and Doppler broadening.

We then sum the equivalent widths of the individual lines and add to them the dimensionless equivalent width of all Lyman series transitions of atomic hydrogen that lie within the energy range of interest, including transitions up to n = 10. Abgrall & Roueff (1989) and Abgrall et al. (1992) provide data on the total radiative de-excitation rates of the excited states, oscillator strengths and transition frequencies of molecular hydrogen. Wise & Fuhr (2009) give simulated data for atomic hydrogen. In a final step, we combine the dimensionless equivalent widths and follow the prescription in Draine & Bertoldi (1996) and divide by the total width of the LW band. This yields f_abs, from which f_esc follows trivially.

3 PARAMETER SPACE

In this section, we present our parameter space. We consider two haloes with masses of 5.6 × 10⁷ (Halo A) and 4 × 10⁸ M_⊙(Halo B), three SFEs and five SEDs that evolve over time for a given total stellar mass in the halo. The haloes in our study are taken from cosmological simulations by Latif & Volonteri (2015). We summarise our simulations in Table 1.

3.1 Haloes

We perform our simulations with data from cosmological simulations from Latif & Volonteri (2015). All 3D quantities are mapped on to 1D radial profiles. In Halo A, we smooth out a region at radius ≈20 pc, where a substructure is falling into the halo. In our 1D picture, this provides a more physical result. In Fig. 1, we show the temperature, density and velocity profiles as well as the initial H₂ abundance. Both haloes follow a similar density profile, but Halo A shows a higher abundance of molecular hydrogen and therefore lower temperatures in the centre. Both haloes are selected at redshift z = 10 and assumed to be metal-free. The virial radii of the objects are 1.3 and 4.3 physical kpc, respectively.

Figure 1.

Number density (upper left-hand panel), temperature (upper right-hand panel), H₂ abundance (lower left-hand panel) and velocity (lower right-hand panel) profiles at start-up. Halo A (black solid line) has a total mass of 5.6 × 10⁷ M_⊙; Halo B (red dashed line) has a total mass of 4 × 10⁸ M_⊙.

Open in new tabDownload slide

3.2 Spectral energy distribution

We study SEDs derived from direct sums of individual Pop III stars and from spectral synthesis models of Pop III and Pop III/Pop II stellar populations calculated with the yggdrasil and starburst99 codes. As the initial mass function (IMF) of first stars is still very uncertain (Clark, Glover & Klessen 2008; Stacy, Greif & Bromm 2010; Clark et al. 2011; Hirano et al. 2014; Stacy & Bromm 2014; Stacy, Bromm & Lee 2016), we consider different descriptions that cover the range currently discussed in the literature. A total of five IMFs for the stars are considered. We use different methods for creating the time-dependent SEDs; for a flat IMF (9–500 M_⊙), we use spectra from Pop III stars alone, for a Salpeter slope IMF (50–500 M_⊙) and a lognormal IMF (1–500 M_⊙), we make use of the code yggdrasil and for Kroupa IMFs (0.1–100 M_⊙), we use a combination of yggdrasil and starburst99; all details can be found in the following subsections. For all SEDs, we assume an instantaneous starburst at the beginning. When the runtime exceeds the lifetime of a star, that star does not explode in a supernova, but is simply removed from the spectrum. We discuss the validity of this assumption in Section 4.3. We show the SED for our fiducial halo-SFE case in Fig. 2 for the first and last moment of the simulation for all five SEDs in terms of dQ/dE, where dQ/dE is the number of photons emitted per second and per electron volt. Depending on the SED, the simulation runtimes differ. The SED with a Salpeter slope ranging from 50 to 500 M_⊙ lasts for 3.6 Myr, corresponding to the lifetime of the least massive star. For all other SEDs, we adopt a runtime given by the lifetime of a 9 M_⊙ star, i.e. 20.2 Myr (Schaerer 2002). A summary of the SEDs can be seen in Fig. 3, where we show the number of ionizing photons and LW photons produced as a function of time.

Figure 2.

SEDs for an SFE of 0.5 per cent in Halo A. Different colours and linestyles mark different SEDs: a Salpeter IMF (olive solid line), a flat IMF (dark blue dashed line), a lognormal IMF (green dot–dashed line), a Kroupa IMF without nebular emission (black dotted line) and a Kroupa IMF with nebular emission (light blue dot–dot-dot–dashed line). The upper panel shows the initial spectrum, and the lower panel shows the final spectrum at the cut-off of the simulation at 3.6 and 20.2 Myr, for the Salpeter and other IMFs, respectively.

Open in new tabDownload slide

Figure 3.

Ionizing (upper panel) and LW (lower panel) photon production rates for all SEDs for an SFE of 0.5 per cent in Halo A. The colour-coding is the same as in Fig. 2.

Open in new tabDownload slide

3.2.1 Flat slope IMF from Pop III stellar spectra

We create a stellar population from a flat slope IMF within the mass range of 9–500 M_⊙, using the stellar SEDs and same methodology as in Mas-Ribas, Dijkstra & Forero-Romero (2016); we refer the reader to that work for a detailed description of the calculations. These authors used parameters from Schaerer (2002) and Marigo et al. (2001), and assumed the stars to be metal-free and to reside on the main sequence. We account for the stellar evolution over time by simply removing from the calculation the stars that reach the end of their lives. This SED is considered a radiation source until the least massive star disappears, t_end = t_9 M⊙ = 20.2 Myr.

3.2.2 Salpeter and lognormal IMF created with yggdrasil

To generate SEDs for Pop III star clusters with the yggdrasil spectral synthesis code, we sum the spectra of its constituent stars (Zackrisson et al. 2011; Rydberg et al. 2015). The cluster is assumed to form instantaneously. In summing the spectra of individual stars to obtain an SED, we consider two top-heavy IMFs. One is a Salpeter-like extremely top-heavy IMF ranging from 50 to 500 M_⊙ with a typical mass of 100 M_⊙ (Schaerer 2002). The second is a lognormal distribution IMF with a typical mass of 10 M_⊙ and the standard deviation of the underlying normal distribution σ = 1.0, ranging from 1 to 500 M_⊙ (Raiter, Schaerer & Fosbury 2010, see Tumlinson 2006). yggdrasil provides spectra for star clusters for ages of up to 3.6 Myr for the extremely top-heavy IMF, which is the lifetime of the longest lived star in the cluster. The lognormal distribution is turned off after 20.2 Myr.

3.2.3 Kroupa IMF with and without nebular emission created with yggdrasil/starburst99

The purely stellar SEDs are based on synthetic spectra from Raiter et al. (2010) for absolute metallicities Z = 10⁻⁷ and 10⁻⁵ and starburst99 (Leitherer et al. 1999) with Geneva high mass-loss tracks in the case of Z = 0.001, 0.004, 0.008 and 0.020. Nebular emission has been added using the cloudy photoionization code (Ferland et al. 2013) with parameters as in Zackrisson et al. (2011), under the assumption of no Lyman continuum leakage and no dust. The stellar SEDs have been rescaled in all cases to be consistent with the Kroupa (2001) stellar IMF.

We point out here that there are some physical limitations of the Kroupa plus nebular emission case. In principle, the radiation in the halo would be a mixture of starlight and nebular emission. No data sets of mixed SED grids are available at this time, and we simply add the nebular to the stellar emission. Thus, we urge the reader to treat this case as an extreme limit where the nebular emission has been maximized. Excluding the Kroupa plus nebular emission case from our work does not qualitatively change our results.

3.3 Star formation efficiencies

For our lower limit SFE of 0.1 per cent in Halo A, we get a total of ≈10⁴ M_⊙ in stars. Mas-Ribas et al. (2016) show that for a total stellar mass of 10³ M_⊙, there is a factor of few difference in the LW flux for different stochastically sampled IMFs, but this is much reduced for ≈10⁴ M_⊙ in stars.

4 RESULTS

Consistent with S15, we find that the LW escape fraction depends on the evolution of the ionization front (I-front). In a completely ionized halo, the escape fraction is 100 per cent, but when the I-front retreats or only slowly expands, molecular hydrogen production is enhanced, allowing it to self-shield. The propagation of the I-front is therefore an important quantity for the calculation of the escape fraction. We identify three different ways the I-front can behave and discuss one prominent example of each type in detail in Section 4.1. In the following sections, we will tabulate the time-averaged escape fractions in both the near-field (Section 4.2) and the far-field (Section 4.4) as well as the escape fraction averaged over the LW flux (Section 4.3). We also present estimates for the ionizing escape fractions based on the position of the I-front in Section 4.5. The escape fractions we provide are lower limits: Our 1D simulations cannot capture asymmetric break-out along directions of minimum column density and we expect radiation to break out faster in some low-density regions.

4.1 Individual ionization front behaviour

4.1.1 Outbreaking ionization front

In the case of a strong radiation source, the ionization front quickly breaks out of the halo. Simultaneously, the escape fraction jumps up to 100 per cent because there is no molecular hydrogen within the ionized region. The lower limit case SFE (1 per cent) for our Salpeter SED in Halo B provides a prominent example. Fig. 4 displays the position of the I-front and the LW escape fraction in the near-field as a function of time.

Figure 4.

Position of the I-front (black solid line; left axis) and LW escape fraction in the near-field (red dot–dashed line; right axis) as a function of time for a 1.0 per cent SFE with a Salpter IMF in Halo B. The virial radius is shown by a black dotted line.

Open in new tabDownload slide

Figure 5.

Abundances of the primordial species as a function of radius for four output times for a 1.0 per cent SFE with a Salpter IMF in Halo B. The black solid line shows the number density of the gas against radius (corresponding to the label on the right-hand side of the plots).

Open in new tabDownload slide

4.1.2 Outbreaking and returning ionization front

Similar to the case just mentioned, here, the ionization front quickly advances through the halo, leading to a sudden increase in the LW near-field escape fraction from 0 to 100 per cent at the beginning of the simulation. However, after some time, the radiation source weakens and produces fewer ionizing photons. As a result, these haloes quickly recombine in an outside-in fashion, triggering the production of molecular hydrogen. The near-field escape fraction drops to 0 per cent when recombination sets in.

Figure 6.

Position of the I-front (black solid line) and LW escape fraction in the near-field (red dot–dashed line) as a function of time for a 1.0 per cent SFE with a flat IMF in Halo B.

Open in new tabDownload slide

Figure 7.

Abundances of the primordial species as a function of radius for six output times for a 1.0 per cent SFE with a flat IMF in Halo B.

Open in new tabDownload slide

4.1.3 Slowly moving ionization front

In the case of a weak radiation source, a D-type I-front advances slowly through the halo. At its border, the many free electrons trigger massive production of H₂, and the dense thin shell that builds up reaches optical depths that can shield the LW radiation (see also Ricotti, Gnedin & Shull 2001). Our calculations depend critically on the thickness and peak abundance of the H₂ shell. This requires specific consideration and we have to implement a region of high numerical resolution at the position of the I-front. We adopt a challenging technique where we resolve a 20-pc wide region that is bracketing and moving with the I-front with a resolution of 0.001 pc. In order to keep the computational time manageable, the first 5 per cent of the simulations are run with our standard resolution, introducing a maximum error of 5 per cent in the time-averaged LW escape fraction.

As an example, in Fig. 8, we show the I-front position and the escape fraction as a function of time for a 0.1 per cent SFE and Salpeter IMF in Halo A. The initial 5 per cent of the simulation had to be run without the refinement around the I-front and therefore shows an LW near-field escape fraction that is probably too low. In the upper left-hand panel of Fig. 9, we show the abundances of all species at 0.4 Myr. The I-front and the shock front move together. A thin shell of molecular hydrogen with peak abundances reaching above 10⁻⁶ forms at the position of the I-front that is able to shield some of the LW radiation. The thin H₂ shell continues moving with the I-front and shock front (upper right-hand panel and middle left-hand panel). From 2.3 Myr on (middle right-hand panel), the spectrum starts to weaken and the molecular hydrogen shell starts to get stronger, preventing more LW radiation from escaping the halo. At 3.0 Myr (lower left-hand panel), the He^{2 +} abundance decreases and the I-front starts to turn back. The shell of molecular hydrogen gets thicker, and at the end of the simulation at 3.6 Myr (lower right-hand panel), all LW radiation is prevented from escaping the halo (in the near-field approximation). Averaged over the runtime of the simulation, the LW escape fraction in this example is 63 per cent. For a slowly moving I-front, we cannot give a simple approximation, in contrast to the outbreaking or returning I-front cases.

Figure 8.

Position of the I-front (black solid line) and LW escape fraction in the near-field (red dot–dashed line) as a function of time for a 0.1 per cent SFE with Salpeter IMF in Halo A.

Open in new tabDownload slide

Figure 9.

Abundances of the primordial species as a function of radius for six output times for a 0.1 per cent SFE with a Salpeter IMF in Halo A.

Open in new tabDownload slide

4.2 Time-averaged escape fractions in the near-field

For all SFE, SED and halo combinations, we average the LW escape fraction over the runtime of the simulation, i.e. 20.2 Myr (or 3.6 Myr for the Salpeter IMF). We list the corresponding escape fractions in Table 3.

Shielding by neutral hydrogen makes only a slight difference: The LW escape fractions are at most 16 per cent higher if we only consider H₂ self-shielding compared to the case where we also account for shielding from H. In many cases, the difference is negligible. This result is in contrast to S15, where we found significantly lower LW escape fractions when including shielding by neutral hydrogen.

As expected, a higher SFE increases the LW escape fraction. This can be observed in both haloes and for all SEDs. Outbreaking I-fronts lead to LW escape fractions ≥ 95 per cent. For the slowly moving or returning I-fronts, the values can be as small as 3 per cent. This is mostly observed for the case of a flat SED ranging from 9 to 500 M_⊙, where many massive stars die at the beginning of the simulation, leading to a much weaker ionizing source of photons later on. As soon as the more massive stars have died, the production of molecular hydrogen sets in and the halo becomes optically thick for the majority of the runtime of the simulations. In all the other simulations, the change in spectrum is not as extreme, and ionizing photons are produced for a longer period, leading to higher LW escape fractions.

Our simulations assume that no stars explode as supernovae. Instead, at the end of their lives, we simply remove their contributions from the SED. We made this simplifying assumption for several reasons. First, we expect the dynamics of the supernova remnant to be highly sensitive to the 3D matter distribution in the protogalaxy (see e.g. Mac Low & Ferrara 1999). It is therefore questionable how well we can represent their effects in the simplified geometry of a 1D model. Secondly, considerable uncertainties still exist regarding the final masses and explosion energies of Pop III stars. In particular, the effects of rotation-driven mass-loss may be significant but not yet fully understood (Ekström, Meynet & Maeder 2008). Finally, the supernovae that explode may be unable to clear all of the gas from the protogalaxy, thus resulting in a system polluted by metals and dust. Our assumption of a primordial gas composition may then become highly inaccurate, depending on the final metallicity of the gas.

Nevertheless, once stars begin to explode as supernovae, it is possible that they will rapidly drive gas out of the protogalaxy, leading to escape fractions quickly reaching 100 per cent. To account for this possibility, we have computed LW near-field escape fractions that are averaged over only the first 2 Myr, rather than over the lifetime of the longest lived massive star (Table 4). 2 Myr corresponds to the lifetime of stars of a few hundred solar masses (Schaerer 2002) and we can safely assume that no supernovae exploded earlier.

Most LW escape fractions averaged over 2 Myr are high; for SFEs with a flat IMF in Halo B, this behaviour is most pronounced. The four outbreaking I-front cases show slightly reduced LW escape fractions, but only down to 97 per cent. This change occurs because the initial few time-steps, until an outbreak with escape fractions smaller than 100 per cent, now carry a slightly larger weight. For the slowly moving I-fronts, the LW escape fraction decreases (with the exception of a 0.1 SFE flat IMF in Halo A). Here, the LW escape fraction slowly increases over time as the I-front moves through the halo. These last results need to be taken as lower limits as we have to run the slowing moving I-front simulations with higher resolution, but can do so only after 1 Myr (0.18 Myr for the Salpeter IMF) due to the high computational costs of the required refinement technique. LW escape fractions of returning I-fronts decrease or increase depending on when the I-front first breaks out (in the case of a quiet late Kroupa IMF) and starts to return. We therefore cannot give a general trend.

4.3 LW flux-averaged escape fractions in the near-field

The LW flux of the SED is decreasing over time (compare Fig. 3) by factors of up to 4 × 10⁴ in the case of the flat IMF. We therefore provide the LW escape fraction weighted by the LW flux emitted by the stellar source rather than averaged over time. We list these values in Table 5 and show them in Fig. 10. Except for an SFE of 0.1 per cent, the escape fraction is always larger than 75 per cent. The difference to the time-averaged LW escape fraction from Section 3.2 is most significant for the flat IMF, as in that case, the LW escape fraction drops when the LW flux drops and the overall contribution is small. Additional shielding by atomic hydrogen does not decrease the LW escape fraction in this limit by more than 0.5 per cent. Therefore, we only provide the escape fraction obtained by considering both H₂ and H shielding together.

Figure 10.

Flux-averaged LW escape fractions in the near-field for our parameter space. Different SEDs are shown by different symbols, Halo A is colour-coded in light blue, and Halo B is colour-coded in black.

Open in new tabDownload slide

4.4 Time-averaged escape fractions in the far-field

In the far-field, the LW escape fractions are generally larger than in the near-field, varying from 59 per cent to 100 per cent. Our values are listed in Table 6. Including shielding by neutral hydrogen can reduce the LW escape fraction by up to 29 per cent (comparing 0.1 per cent SFE with a Kroupa SED in Halo A) and thus should not be forgotten in simulations. For the fiducial SFE, all LW escape fractions are higher than 75 per cent and we conclude that the host haloes of the first galaxies do not suppress the build-up of the LW background by a large factor.

4.5 Lower limits for ionizing escape fractions

Since we track the position of the I-front through the halo, we are able to calculate the escape fractions of ionizing photons as well. The values need to be taken as a first-order approximation only, since our 1D setup prevents us from seeing the outbreak of ionized cones like authors recently do with high-resolution 3D studies (Ritter et al. 2012; Paardekooper, Khochfar & Dalla Vecchia 2015; Stacy et al. 2016). The escape fraction of ionizing photons is determined by the position of the I-front alone. For each simulation output, we check if the I-front has crossed the virial radius, in which case we take f_{esc, ion}(t) = 1, otherwise, f_{esc, ion}(t) = 0. Averaging these values over the runtime of the simulation yields the ionizing escape fractions f_{esc, ion}, listed in Table 7.

The values of f_{esc, ion} are similar to f_{esc, LW} in the near-field for the case of an outbreaking or a returning I-front, since the LW escape fractions depend critically on the position of the I-front. For the case of a slowly moving I-front that never crosses the virial radius, we obtain ionizing escape fractions of zero while the LW escape fractions vary from 3 per cent to 85 per cent in the near-field. In these cases, the outer halo is optically thick to ionizing photons, but optically thin to LW photons.

5 DISCUSSION AND CONCLUSION

We studied a multidimensional grid of parameters to calculate the LW escape fractions in both the near-field and the far-field limit. We consider two haloes of 5.6 × 10⁷ and 4.0 × 10⁸ M_⊙, which are irradiated by a stellar population in their centre. Different SEDs, ranging from a Salpeter IMF between 50 and 500 M_⊙ to a Kroupa IMF between 0.1 and 100 M_⊙, combined with SFEs between 0.1 per cent and 5.0 per cent, are taken into account.

We find that the near-field LW escape fraction depends on the motion of the I-front and we therefore group our results into three different ways the I-front can behave:

• For an outbreaking I-front, the LW escape fraction in the near-field is |$f_\mathrm{esc}^{\mathrm{(o)}} \ge$| 95 per cent. This is obtained for SFEs greater than 0.5 per cent with a Salpeter slope IMF ranging from 50 to 500 M_⊙: The I-front overruns the halo and all gas is ionized for the rest of the simulation.

• In the case of a spectrum that weakens substantially over time, haloes recombine outside-in. Large numbers of free electrons still present in the outer parts of the halo lead to the aforementioned molecular hydrogen formation in the outer part of the halo. The halo becomes optically thick. The LW escape fraction can be approximated by |$f_\mathrm{esc}^{\mathrm{(r)}} \approx t_\mathrm{return} / t_\mathrm{total}$|⁠, where t_return is the time when the recombination sets in and t_total is the runtime of the simulation.

• For a slowly moving I-front, the LW escape fraction can be as small as 3 per cent and as large as 88 per cent, depending on a thin shell of molecular hydrogen that builds up at the position of the I-front. In this case, the stellar radiation acts as an agent of negative feedback that enhances the abundance of molecular hydrogen and suppresses the LW background (see also Ricotti et al. 2001, discussing H₂-shell formation in the IGM).

Out of 30 simulations, we find 4 overrunning, 17 returning and 9 slowly moving I-fronts.

In the far-field, all LW escape fractions are much larger than in the near-field, as we are shifting out of the line centres and average the LW escape fraction over the whole LW energy range. Here, shielding by neutral hydrogen starts to play a role since Lyman lines of H β and higher, which lie inside the LW range, become very broad (Haiman, Abel & Rees 2000). Neutral hydrogen is able to reduce the LW escape fraction in the far-field by up to 29 per cent. In general, all LW escape fractions in the far-field are higher than 75 per cent in the fiducial case and the halo only slightly hinders the build-up of the LW background.

We also tabulated values of the ionizing escape fractions for our parameter space, depending on the position of the I-front relative to the virial radius. For slowly moving I-fronts, this value equals zero, since the I-front does not reach out to the virial radius at any time in our simulation. For simulations where the I-front quickly breaks out of the halo, the ionizing escape fractions are larger than 95 per cent. If the I-front first breaks out and later returns due to outside-in recombination of the halo, the escape fraction ranges from 10 per cent to 85 per cent.

The biggest drawback of this study is our 1D approach. In 3D simulations, ionizing photons break cones into the haloes, following low-density regions. We expect our LW escape fractions to be close to 1 in these directions, leading to an overall higher f_{esc, LW}. We therefore consider our LW escape fractions to be lower limits.

Additionally, our UV source is placed in the centre of our 1D halo and all relative velocities of the contributing stars are neglected. For the slowly moving I-front cases, the LW near-field escape fractions are underestimated in this work. For outbreaking or returning I-fronts, the abundance of molecular hydrogen is either too low or too high to be influenced by relative motions at the centre of the halo.

We have shown in this study that the LW escape fraction from the first galaxies is high. Unless a protogalaxy harbours a star cluster with flat IMF or forms stars with a low ∼0.1 per cent SFE, a large fraction of its LW radiation can escape from the host halo. The LW escape fractions in the far-field exceeds 60 per cent in all cases. We thus conclude that the gas present in the host halo of a stellar cluster with 10⁴ M_⊙ to a few times 10⁶ M_⊙ in stars cannot prevent the build-up of the LW background.

The recently observed Lyα system CR7 at z = 6.6 (Sobral et al. 2015) is comprised of two star-forming components and one seemingly metal-poor component (however, see Bowler et al. 2016). Agarwal et al. (2016) showed that the star-forming components can easily produce the LW radiation field required for DCBH formation in the metal-poor clump. For an exponentially declining star formation history, they find that DCBH formation can occur in similar systems as early as z ∼ 20, if the star-forming component consists of a DM halo or haloes of M_DM ∼ 5 × 10⁸ M_⊙ and an SFE of 8 per cent or larger (Agarwal et al. 2017). This closely matches our case of Halo B-forming stars at its upper limit of SFE∼5 per cent with the same underlying IMF and stellar mass cut-offs. Given the separation between the star-forming component and the metal-poor component of CR7 over which the LW feedback allows for DCBH formation (Agarwal et al. 2016), we can apply our near-field limit. Thus, we can conclude that the star-forming component of CR7 type systems could have an escape fraction as high as 96 per cent at onset of DCBH formation.

In general, the near-field can be applied to neighbouring systems that move with no or only a small velocity relative to the halo-emitting LW radiation. The critical relative velocity separating the two regimes depends on the strength of the LW absorption lines in the emitting halo. When the lines are weak and dominated by thermal broadening, |$v_{\rm crit} \simeq v_{\rm th, H_{2}}$|⁠, the thermal velocity of the H₂, which in our simulated haloes is typically 7–15 km s⁻¹. On the other hand, if the main LW lines are strong and dominated by Lorentz broadening (i.e. if the near-field LW escape fraction is very small), then the linewidths and hence the critical relative velocity can both be much larger.

Acknowledgments

We would like to thank Mordecai-Marc Mac Low, Mattis Magg and Massimo Ricotti for useful discussions. Furthermore, we would like to thank the anonymous referee for his/her valuable comments. ATPS, BA, C-ER and DJW acknowledge support from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007 – 2013) via the ERC Advanced Grant ‘STARLIGHT: Formation of the First Stars’ (project number 339177). BA acknowledges support from a TCAN postdoctoral fellowship at Yale. SCOG and R S K. also acknowledge support from the Deutsche Forschungsgemeinschaft via SFB 881, ‘The Milky Way System’ (sub-projects B1, B2 and B8) and SPP 1573, ‘Physics of the Interstellar Medium’ (grant number GL 668/2-1). ML acknowledges funding from European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 656428. EZ acknowledges funding from the Swedish Research Council (project number 2011-5349). All simulations were run on the Sciama Supercomputer at the Institute of Cosmology and Gravitation at the University of Portsmouth.

REFERENCES