Lyman-Werner Escape Fractions from the First Galaxies

TL;DR

This study uses radiation-hydrodynamical simulations to determine how the escape fraction of Lyman-Werner UV photons from the first galaxies varies with different ionisation front propagation modes, impacting black hole seed formation.

Contribution

It provides detailed analysis of LW escape fractions based on ionisation front types, revealing dependencies on microphysics and spectral properties not previously quantified.

Findings

LW escape fraction exceeds 95% for R-type I-fronts

Near-field escape fraction can be as low as 3% for certain I-fronts

Far-field escape fraction can reach 61% or 100% depending on conditions

Abstract

Direct collapse black holes forming in pristine, atomically-cooling haloes at $z \approx 10-20$ may act as the seeds of supermassive black holes (BH) at high redshifts. In order to create a massive BH seed, the host halo needs to be prevented from forming stars. H$_2$ therefore needs to be irradiated by a large flux of Lyman-Werner (LW) UV photons in order to suppress H$_2$ cooling. A key uncertainty in this scenario is the escape fraction of LW radiation from first galaxies, the dominant source of UV photons at this epoch. To better constrain this escape fraction, we have performed radiation-hydrodynamical simulations of the growth of HII regions and their associated photodissociation regions in the first galaxies using the ZEUS-MP code. We find that the LW escape fraction crucially depends on the propagation of the ionisation front (I-front). For an R-type I-front overrunning the halo, the LW escape fraction is always larger than 95%. If the halo recombines later from the outside--in, due to a softened and weakened spectrum, the LW escape fraction in the rest-frame of the halo (the near-field) drops to zero. A detailed and careful analysis is required to analyse slowly moving, D-type I-fronts, where the escape fraction depends on the microphysics and can be as small as 3% in the near-field and 61% in the far-field or as large as 100% in both the near-field and the far-field.

Tables

Table 1. Table 1: Summary of all spectral energy distributions used in this study.

$M_{\min }$ [M⊙]	$M_{\max }$ [M⊙]	runtime [Myr]	shape
50	500	3.6	Salpeter
9	500	20.2	flat
1	500	20.2	log-normal
0.1	100	20.2	Kroupa
0.1	100	20.2	Kroupa + nebula

Table 2. Table 2: Star formation efficiencies and the corresponding stellar masses for both haloes. The fiducial SFEs are calculated using a formula relating the mass in stars to the total mass from O’Shea et al . ( 2015 ) ; see Equation 1.

	$M_{\mathrm{vir}}$ [M⊙]	$M_b$ [M⊙]	$M_{\star }$ [M⊙]	$f_{\star }$ [%]	SFE [%]
Halo A, fiducial	$5.6\times {10}^7$	$9.4\times {10}^6$	$4.6\times {10}^4$	0.082	0.5
Halo A, lower limit			$9.4\times {10}^3$	0.015	0.1
Halo A, upper limit			$9.4\times {10}^4$	0.15	1.0
Halo B, fiducial	$4.0\times {10}^8$	$6.7\times {10}^7$	$1.4\times {10}^6$	0.35	2.1
Halo B, lower limit			$6.7\times {10}^5$	0.16	1.0
Halo B, upper limit			$3.3\times {10}^6$	0.79	5.0

Table 3. Table 3: Time-averaged LW escape fractions in the near-field, averaged over the runtime of the simulation. All values are given in percent. The upper six rows refer to Halo A, the lower six rows to Halo B. The upper three rows of haloes A and B, respectively, are calculated with H 2 shielding only, the lower three rows with H and H 2 shielding. The letter in parentheses refers to the case of I-front behaviour: (o) stands for outbreaking I-front, (r) for returning I-front and (s) for slowly moving I-front.

			Salpeter	flat	log-normal	Kroupa	Kroupa + nebula	shielding
	SFE		50–500 M⊙	9–500 M⊙	1–500 M⊙	0.1–100 M⊙	0.1–100 M⊙
Halo A	0.1		63 (s)	3 (s)	38 (s)	58 (s)	63 (s)	H2 only
	0.5		93 (r)	11 (r)	26 (r)	88 (s)	90 (s)
	1.0		99 (o)	13 (r)	44 (r)	90 (s)	92 (s)
	0.1		59 (s)	3 (s)	37 (s)	47 (s)	51 (s)	H2 and H
	0.5		93 (r)	11 (r)	26 (r)	84 (s)	86 (s)
	1.0		99 (o)	13 (r)	44 (r)	88 (s)	89 (s)
Halo B	1.0		100 (o)	16 (r)	61 (r)	66 (r)	70 (r)	H2 only
	2.1		100 (o)	18 (r)	78 (r)	86 (r)	88 (r)
	5.0		100 (o)	20 (r)	94 (r)	96 (r)	97 (r)
	1.0		99 (o)	16 (r)	61 (r)	64 (r)	68 (r)	H2 and H
	2.1		100 (o)	18 (r)	78 (r)	84 (r)	86 (r)
	5.0		100 (o)	20 (r)	94 (r)	95 (r)	96 (r)

Table 4. Table 4: Time-averaged LW escape fractions in the near-field, averaged over the first two Myr. All values are given in percent. The structure of the rows is as in Table 3 .

			Salpeter	flat	log-normal	Kroupa	Kroupa + nebula	shielding
	SFE		50–500 M⊙	9–500 M⊙	1–500 M⊙	0.1–100 M⊙	0.1–100 M⊙
Halo A	0.1		54 (s)	27 (s)	29 (s)	41 (s)	48 (s)	H2 only
	0.5		87 (r)	94 (r)	45 (r)	57 (s)	64 (s)
	1.0		97 (o)	100 (r)	85 (r)	66 (s)	72 (s)
	0.1		54 (s)	27 (s)	28 (s)	40 (s)	47 (s)	H2 and H
	0.5		87 (r)	93 (r)	44 (r)	56 (s)	63 (s)
	1.0		97 (o)	99 (r)	84 (r)	65 (s)	71 (s)
Halo B	1.0		99 (o)	100 (r)	97 (r)	46 (r)	58 (r)	H2 only
	2.1		100 (o)	100 (r)	99 (r)	70 (r)	76 (r)
	5.0		100 (o)	100 (r)	100 (r)	92 (r)	93 (r)
	1.0		99 (o)	99 (r)	96 (r)	45 (r)	58 (r)	H2 and H
	2.1		99 (o)	99 (r)	98 (r)	69 (r)	75 (r)
	5.0		99 (o)	99 (r)	99 (r)	91 (r)	92 (r)

Table 5. Table 5: LW flux averaged LW escape fractions in the near-field, averaged over the total LW flux of the simulation. All values are given in percent. The upper three rows refer to Halo A, the lower three rows refer to Halo B. Shielding by H 2 only and by the combination of H and H 2 yields the same results.

			Salpeter	flat	log-normal	Kroupa	Kroupa + nebula	shielding
	SFE		50–500 M⊙	9–500 M⊙	1–500 M⊙	0.1–100 M⊙	0.1–100 M⊙
Halo A	0.1		66	27	50	55	60	H2 and H
	0.5		92	92	61	82	85
	1.0		99	99	82	86	88
Halo B	1.0		100	100	91	74	78	H2 and H
	2.1		100	100	96	89	90
	5.0		100	100	99	97	97

Table 6. Table 6: Time-averaged LW escape fractions in the far-field. All values are given in percent. The upper six rows refer to Halo A, the lower six rows to Halo B. The upper three rows of haloes A and B, respectively, are calculated with H 2 shielding only, the lower three rows with H and H 2 shielding.

			Salpeter	flat	log-normal	Kroupa	Kroupa + nebula
	SFE		50–500 M⊙	9–500 M⊙	1–500 M⊙	0.1–100 M⊙	0.1–100 M⊙
Halo A	0.1		98	73	96	99	100	H2 only
	0.5		99	85	96	100	100
	1.0		100	84	98	100	100
	0.1		79	61	83	68	68	H2 and H
	0.5		91	78	88	84	84
	1.0		94	79	90	86	86
Halo B	1.0		100	85	99	99	100	H2 only
	2.1		100	86	100	100	100
	5.0		100	88	100	100	100
	1.0		97	76	91	87	87	H2 and H
	2.1		98	78	93	90	90
	5.0		99	80	95	93	93

Table 7. Table 7: Time-averaged ionising escape fractions in the near-field. All values are given in percent. The upper three rows refer to Halo A, the lower three rows to Halo B. All values are derived from the position of the I-front and provide lower limits.

			Salpeter	flat	log-normal	Kroupa	Kroupa + nebula
	SFE		50–500 M⊙	9–500 M⊙	1–500 M⊙	0.1–100 M⊙	0.1–100 M⊙
Halo A	0.1		0	0	0	0	0
	0.5		85	10	18	0	0
	1.0		96	13	36	0	0
Halo B	1.0		98	17	54	26	26
	2.1		99	19	66	42	42
	5.0		100	22	84	58	58

Equations

$$f_{\star }=1.26\times 10^{-3}{\left(M_{\mathrm{vir}}/10^8{M}_{\odot }\right)}^{0.74},$$

$$\mathrm{SFE}=\frac{M_{\star }}{M_{\mathrm{b}}},$$

$$f_{\mathrm{esc},\mathrm{LW}}^{(\mathrm{o})}\ge 0.95$$

$$f_{\mathrm{esc},\mathrm{LW}}^{(\mathrm{r})}=\frac{t_{\mathrm{return}}}{t_{\mathrm{total}}},$$

Full text

††pagerange: Lyman-Werner Escape Fractions from the First Galaxies–References††pubyear: 2002

1 Introduction

In recent years, several highly luminous quasars have been observed, for example a $1.2\times 10^{10}$ M*⊙* BH at redshift $z\approx 6.3$ (Wu et al., 2015) and a $9\times 10^{9}$ M*⊙* BH at redshift $z\approx 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$\alpha$ emitter CR7 (Sobral et al., 2015), recently discovered at $z=6.6$, is found to have a unique structure. A deep HST image shows it consists of three clumps, separated by $\approx$ 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 (Hartwig et al. 2016; Agarwal et al. 2016; Pallottini et al. 2015; Dijkstra, Gronke & Sobral 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 the one hand, SMBH seeds have been suggested to be remnants from Pop III stars (Madau & Rees, 2001; Fryer, Woosley & Heger, 2001; Inayoshi, Haiman & Ostriker, 2016) with relatively low masses ($10^{2-3}$ M*⊙) that accrete close to the Eddington limit for large fractions of their lifetime. On the other hand, very massive seeds ($10^{4-6}$ 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^{4}$ K) halo in which gas contracts quasi-isothermally without cooling and fragmenting. The central clump, to which the halo compresses, contains $\geq 10\%$ 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 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 prep.) 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 place of such a direct collapse BH.

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: 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 built up of the LW background?

We have recently performed a study on LW escape fractions from minihaloes (Schauer et al. 2015, from here on S15; see also Kitayama et al. 2004), finding the escape fraction varying with halo- and stellar mass from 0 to 95%. 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^{7}$ - $10^{8}$ 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\approx 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 of the LW escape fractions for all halo, SED and SFE combinations over time in both 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 ionising UV transport self consistently to non-equilibrium primordial chemistry and hydrodynamics to model the growth of cosmological ionisation fronts (I-fronts; Whalen & Norman, 2006). We include a network of nine chemical species (H, H*+, He, He+, He2+, H-, H${}^{+}_{2}$, H2* and e*-*). Heating and cooling processes like photoionisation, collisional ionisation, excitation and recombination, inverse Compton scattering from the CMB, bremsstrahlung emission and H2 cooling mostly following 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 eV 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 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 20000 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% 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, H2 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 H2 self-shielding in the low density and LTE limits, respectively. They provide simple fitting functions describing the shielding as a function of H2 column density and temperature. LW photons that do not dissociate H2 molecules or are not absorbed and re-emitted as lower-energy photons, can escape the halo. For high column densities, H2 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_{\mathrm{H2}}>10^{13}$ cm*-2*), we use the shielding functions from WH11, which depend on the column density $N_{\mathrm{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 H2 contributes, and another where we additionally account for shielding from neutral hydrogen. In the H2-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_{\mathrm{H2}}<10^{13}$ cm*-2*), 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 restframe 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 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 H2 or H in the halo. If we denote this fraction as $f_{\mathrm{abs}}$, then the far-field escape fraction is simply $f_{\mathrm{esc}}=1$ - $f_{\mathrm{abs}}$.

To calculate $f_{\mathrm{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_{\mathrm{abs}}$, from which $f_{\mathrm{esc}}$ follows trivially.

3 Parameter Space

In this section, we present our parameter space. We consider two haloes with masses of $5.6\times 10^{7}$ (Halo A) and $4\times 10^{8}$ 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).

3.1 Haloes

We perform our simulations with data from cosmological simulations from Latif & Volonteri (2015). All three dimensional quantities are mapped onto one dimensional radial profiles. In Halo A, we smooth out a region at radius $\approx$ 20 pc, where a substructure is falling into the halo. In our one dimensional picture, this provides a more physical result. In Figure 1, we show the temperature, density and velocity profiles as well as the initial H2 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.

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 IMF of first stars is still very uncertain (Clark, Glover & Klessen, 2008; Stacy, Greif & Bromm, 2010; Clark et al., 2011; Stacy & Bromm, 2014; Hirano et al., 2014; Stacy, Bromm & Lee, 2016), we consider different descriptions that cover the range currently discussed in the literature. A total of 5 initial mass functions (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 Starburst 99; 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 Figure 2 for the first and last moment of the simulation for all five SEDs in terms of $\mathrm{d}Q/\mathrm{d}E$, where $\mathrm{d}Q/\mathrm{d}E$ is the number of photons emitted per second and per eV. Depending on the SED, the simulation runtimes differ. The SED with a Salpeter slope ranging from 50 – 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, 20.2 Myr (Schaerer, 2002). A summary of the SEDs can be seen in Figure 3, where we show the number of ionising photons and LW photons produced as a function of time.

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 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, tend = t${}_{9\mathrm{M}_{\odot}}$ = 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 – 500 M*⊙* with a typical mass of 100 M*⊙* (Schaerer, 2002). The second is a log-normal distribution IMF with a typical mass of 10 M*⊙* and the standard deviation of the underlying normal distribution $\sigma=1.0$, ranging from 1 M*⊙* 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 log-normal 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, Schaerer & Fosbury (2010) for absolute metallicities $Z=10^{-7}$ and $Z=10^{-5}$ 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 photoionisation 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 initial mass function.

We point out here that there are some physical limitation of the Kroupa plus nebular emission case. In principle, the radiation in the halo would be a mixture of starlight and nebular emission. No datasets 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 maximised. Excluding the Kroupa plus nebular emission case from our work does not qualitatively change our results.

3.3 Star Formation Efficiencies

The SEDs from the previous section are normalised to a total number of stars, calculated for each halo. For the fiducial case, we adopt the stellar to virial mass ratio from O’Shea et al. (2015)

[TABLE]

which is valid in the mass range $10^{7}$ M*⊙$\leq M_{\mathrm{vir}}\leq 10^{8.5}$ M⊙. Our most massive halo, Halo B with $M_{\mathrm{vir}}=10^{8.6}$ M⊙*, is slightly outside this range and we overestimate the stellar mass by a small amount.

The stellar mass, $M_{\star}=f_{\star}\times M_{\mathrm{vir}}$, relates to a star formation efficiency:

[TABLE]

where $M_{\mathrm{b}}$ is the baryonic mass in a halo approximated by $M_{\mathrm{b}}=\Omega_{\mathrm{b}}/\Omega_{0}\times M_{\mathrm{vir}}$. Using the values from Jarosik et al. (2011), $\Omega_{\mathrm{b}}=0.0456$ and $\Omega_{0}=0.2726$, we derive SFEs of 0.5% (Halo A) and 2.1% (Halo B). In addition, we adopt an upper and a lower SFE limit on each halo, ranging from 0.1% to 5.0%. All properties are listed in Table 2.

For our lower limit SFE of 0.1% in Halo A we get a total of $\approx 10^{4}$ M*⊙* in stars. Mas-Ribas, Dijkstra & Forero-Romero (2016) show that for a total stellar mass of $10^{3}$ M*⊙, there is a factor of few difference in the LW flux for different stochastically sampled IMFs, but this is much reduced for $\approx 10^{4}$ M⊙* in stars.

4 Results

Consistent with S15, we find that the LW escape fraction depends on the evolution of the ionisation front (I-front). In a completely ionised halo, the escape fraction is 100%, 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 subsections we will tabulate the time-averaged escape fractions in both near-field (Section 4.2) and 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 ionising escape fractions based on the position of the I-front in Section 4.5. The escape fractions we provide are lower limits: our one dimensional 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 ionisation front behaviour

4.1.1 Outbreaking ionisation front

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

The abundances of the main chemical species can be seen in Figure 5 for four different output times. All species are colour coded. The black solid line shows the number density in cm*-3*. Initially (upper left panel), the halo is mostly neutral with only a small amount of ionised hydrogen outside 200 pc. Already at 0.02 Myr (upper right panel), a very steep I-front has advanced out to 9 pc, indicated by the turquoise H ii region that has a sharp transition to the light blue neutral H outer region. The I-front continues to proceed as an R-type I-front through the halo and crosses the virial radius at 0.08 Myr (lower left panel). Helium is doubly ionised out to 1000 pc. The density still shows its initial profile. At the end of the simulation, at 3.6 Myr (lower right panel), the halo is still completely ionised. The expansion of the ionised gas has started to drive a shock front outwards, but this has only reached 200 pc. It is therefore well within the H ii region and hence is not relevant for our LW escape fraction calculation. Throughout the simulation, the H2 abundance remains much smaller than $10^{-4}$ and the halo is optically thin to LW radiation. The time-averaged LW escape fraction in the near-field therefore is 99% in this case. In general,

[TABLE]

for both the near-field and far-field limit in all models with I-front breakout (indicated by the superscript (o)).

4.1.2 Outbreaking and returning ionisation front

Similar to the case just mentioned, here, the ionisation front quickly advances through the halo, leading to a sudden increase in the LW near-field escape fraction from 0 to 100% at the beginning of the simulation. However, after some time, the radiation source weakens and produces fewer ionising 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% when recombination sets in.

This can be clearly observed in the lower limit SFE case of a flat IMF in Halo B. In Figure 6, a fast outbreak of the I-front yields a LW escape fraction of 100% almost immediately after the simulation starts. The abundances and number density are displayed in Figure 7. At 2.6 Myr (upper left panel), the halo is still completely ionised, but as the spectrum weakens, at 3.0 Myr (upper right panel), the outer part of the halo partly recombines. One can observe that at 3.4 Myr (middle left panel), the recombination proceeds outside-in. At 4.0 Myr (middle right panel), the I-front has moved back to the shock-front position. Here, the high density and the large fraction of free electrons produce the ideal environment for the formation of molecular hydrogen. The H2 abundance peaks close to 0.01 and hinders all LW photons from escaping in the near-field. Molecular hydrogen continues to form outside of the shock-front of the halo at later times (8.1 Myr, lower left panel) until the end of the simulation (at 20.2 Myr, lower right panel). The time-averaged near-field escape fraction is only 16%, as after the outside-in recombination in the halo, LW photons are blocked completely. The LW near-field escape fraction can be approximated in simulations of returning I-fronts (indicated by the superscript (r))as

[TABLE]

where $t_{\mathrm{return}}$ is the point in time when the I-front starts moving backwards through the halo and $t_{\mathrm{total}}$ is the total runtime of the simulation.

4.1.3 Slowly moving ionisation 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 H2 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 H2 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% of the simulations are run with our standard resolution, introducing a maximum error of 5% in the time-averaged LW escape fraction.

As an example, in Figure 8 we show the I-front position and the escape fraction as a function of time for a 0.1% SFE and Salpeter IMF in Halo A. The initial 5% of the simulation had to be run without the refinement around the I-front and therefore shows a LW near-field escape fraction that is probably too low. In the upper left panel of Figure 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^{-6}$ forms at the position of the I-front that is able to shield some of the LW radiation. The thin H2 shell continues moving with the I-front and shock front (upper right panel and middle left panel). From 2.3 Myr on (middle right 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 panel), the He2+ 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 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%. For a slowly moving I-front, we cannot give a simple approximation, in contrast to the outbreaking or returning I-front cases.

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, 20.2 Myr (or 3.6 Myr for the Salpeter IMF). We list the corresponding escape fractions in Table 3.

Shielding by neutral hydrogen only makes a slight difference: the LW escape fractions are at most 16% higher if we only consider H2 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 $\geq$ 95%. For slowly moving or returning I-fronts, the values can be as small as 3%. 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 ionising 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 other simulations, the change in spectrum is not as extreme and ionising 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. Second, 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%. 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). Two 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%. This change occurs because the initial few timesteps until outbreak with escape fractions smaller than 100% now carry a slightly larger weight. For 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 Kroupa IMF quite late) 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 Figure 3) by factors of up to $4\times 10^{4}$ 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 Figure 10. Except for a SFE of 0.1%, the escape fraction is always larger than 75%. 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%. Therefore we only provide the escape fraction obtained by considering both H2 and H shielding together.

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% to 100%. Our values are listed in Table 6. Including shielding by neutral hydrogen can reduce the LW escape fraction by up to 29% (comparing 0.1% 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% 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 ionising escape fractions

Since we track the position of the I-front through the halo, we are able to calculate the escape fractions of ionising photons as well. The values need to be taken as a first order approximation only, since our one dimensional setup prevents us from seeing the outbreak of ionised cones like authors recently do with high-resolution three dimensional studies (Ritter et al., 2012; Paardekooper, Khochfar & Dalla Vecchia, 2015; Stacy, Bromm & Lee, 2016).

The escape fraction of ionising 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_{\mathrm{esc,ion}}(t)=$ 1; otherwise, $f_{\mathrm{esc,ion}}(t)=$ 0. Averaging these values over the runtime of the simulation yields the ionising escape fractions $f_{\mathrm{esc,ion}}$, listed in Table 7.

The values of $f_{\mathrm{esc,ion}}$ are similar to $f_{\mathrm{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 ionising escape fractions of zero while the LW escape fractions vary from 3% to 85% in the near-field. In these cases, the outer halo is optically thick to ionising 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\times 10^{7}$ and $4.0\times 10^{8}$ M*⊙, that 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% and 5.0% 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)}}\geq$ 95%. This is obtained for SFEs greater than 0.5% with a Salpeter slope IMF ranging from 50 to 500 M*⊙*: the I-front overruns the halo and all gas is ionised 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_{\mathrm{return}}$ is the time when the recombination sets in and $t_{\mathrm{total}}$ is the runtime of the simulation.

• •

For a slowly moving I-front, the LW escape fraction can be as small as 3% and as large as 88%, 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, Gnedin & Shull 2001, discussing H2-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$\beta$ and higher, that 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%. In general, all LW escape fractions in the far-field are higher than 75% in the fiducial case and the halo only slightly hinders the build-up of the LW background.

We also tabulated values of the ionising 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 ionising escape fractions are larger than 95%. If the I-front first breaks out and later returns due to outside-in recombination of the halo, the escape fraction ranges from 10% to 85%.

The biggest drawback of this study is our 1D approach. In 3D simulations, ionising 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_{\mathrm{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 $\sim 0.1$% SFE, a large fraction of its LW radiation can escape from the host halo. The LW escape fractions in the far-field exceeds 60% in all cases. We thus conclude that the gas present in the host halo of a stellar cluster with $10^{4}$ M*⊙* to a few times $10^{6}$ M*⊙* in stars cannot prevent the build-up of the LW background.

The recently observed Ly$\alpha$ 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\sim 20$, if the star-forming component consists of a DM halo or haloes of $M_{\mathrm{DM}}\sim 5\times 10^{8}$ M*⊙* and a SFE of $8\%$ or larger (Agarwal et al. in preparation). This closely matches our case of Halo B forming stars at its upper limit of SFE$\sim 5\%$ 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% 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 H2, which in our simulated haloes is typically 7–15 km s*-1*. 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. A. T. P. S., B. A., C.-E. R. and D. J. W. 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). B. A. acknowledges support of a TCAN postdoctoral fellowship at Yale. S. C. O. G. 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). M. L. acknowledges funding from European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 656428. E. Z. 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.