Understanding the escape of LyC and Ly$\alpha$ photons from turbulent clouds

TL;DR

This study uses radiation-hydrodynamic simulations to explore how turbulent molecular clouds' properties influence the escape of LyC and Lya photons, revealing factors that enhance photon escape and their observational signatures.

Contribution

It provides new insights into the dependence of photon escape fractions on cloud properties, feedback mechanisms, and stellar spectra, advancing understanding of galaxy reionization processes.

Findings

Escape fractions increase over time with cloud dispersal.

Lower metallicity and higher stellar mass limits enhance LyC escape.

Lya photons have higher escape fractions and distinctive spectral profiles.

Abstract

Understanding the escape of Lyman continuum (LyC) and Lyman alpha (Lya) photons from molecular clouds is one of the keys to constraining the reionization history of the Universe and the evolution of galaxies at high redshift. Using a set of radiation-hydrodynamic simulations with adaptive mesh refinement, we investigate how photons propagate and escape from turbulent clouds with different masses, star formation efficiencies (SFEs), and metallicities, as well as with different models of stellar spectra and supernova feedback. We find that the escape fractions in both LyC and Lya are generally increasing with time if the cloud is efficiently dispersed by radiation and supernova feedback. When the total SFE is low (1% of the cloud mass), 0.1-5% of LyC photons leave the metal-poor cloud, whereas the fractions increase to 20-70% in clouds with a 10% SFE. LyC photons escape more efficiently if gas metallicity is lower, if the upper mass limit in the stellar initial mass function is higher, if binary interactions are allowed in the evolution of stars, or if additional strong radiation pressure, such as Lya pressure, is present. As a result, the number of escaping LyC photons can easily vary by a factor of $\sim4$ on cloud scales. The escape fractions of Lya photons are systemically higher (60-80%) than those of LyC photons despite large optical depths at line centre ($\tau_0\sim10^6-10^9$). Scattering of Lya photons is already significant on cloud scales, leading to double-peaked profiles with peak separations of $v_{\rm sep}\sim400\,{\rm km\,s^{-1}}$ during the initial stage of the cloud evolution, while it becomes narrower than $v_{\rm sep} \le 150 \, {\rm km\,s^{-1}}$ in the LyC bright phase. Comparisons with observations of low-redshift galaxies suggest that Lya photons require further interactions with neutral hydrogen to reproduce their velocity offset for a given LyC escape fraction.

Tables

Table 1. Table 1: The eight photon groups used in our simulations. From left to right, each column describes the name of the photon group, minimum and maximum photon energy, dust opacity, and main purpose of each photon group.

Photon	${ϵ}_0$	${ϵ}_1$	$\kappa$	Main function
group	[eV]	[eV]	[${\mathrm{cm}}^2/g$]
EUVHeII	54.42	$\infty$	${10}^3$	HeII ionisation
EUVHeI	24.59	54.42	${10}^3$	HeI ionisation
EUVHI,2	15.2	24.59	${10}^3$	HI and $H_2$ ionisation
EUVHI,1	13.6	15.2	${10}^3$	HI ionisation
LW	11.2	13.6	${10}^3$	$H_2$ dissociation
FUV	5.6	11.2	${10}^3$	Photoelectric heating
Optical	1.0	5.6	${10}^3$	Radiation pressure
IR	0.1	1.0	5	Radiation pressure

Table 2. Table 2: List of simulations performed in this study. From left to right, each column shows the name of the simulation, the size of the simulated box, the maximum cell width, the mass and half-mass radius of the cloud, the total stellar mass, the initial gas metallicity, the type of SED, the upper mass limit to the IMF, the time of the last snapshot, the placement of stars in the initial conditions, and necessary remarks.

Name	$L_{\mathrm{box}}$	$\Delta x_{\min }$	$M_{\mathrm{cloud}}$	$r_{1/2}$	$M_{\mathrm{star}}$	$Z_{\mathrm{gas}}$	SED	$M_{\max }$	$t_{\mathrm{final}}$	SF	Remarks
	[pc]	[pc]	[$M_{\odot }$]	[pc]	[$M_{\odot }$]			[$M_{\odot }$]	[Myr]
M6_SFE10	512	0.25	${10}^6$	35	${10}^5$	0.002	BPASS	100	10	random
M6_SFE1	512	0.25	${10}^6$	35	${10}^4$	0.002	BPASS	100	20	random	Fiducial
M5_SFE10	256	0.25	${10}^5$	13	${10}^4$	0.002	BPASS	100	7	random
M5_SFE1	256	0.25	${10}^5$	13	${10}^3$	0.002	BPASS	100	20	random
M6_SFE10_sng	512	0.25	${10}^6$	35	${10}^5$	0.002	BC03	100	10	random	Single stellar SED
M6_SFE1_sng	512	0.25	${10}^6$	35	${10}^4$	0.002	BC03	100	20	random	Single stellar SED
M6_SFE10_300	512	0.25	${10}^6$	35	${10}^5$	0.002	BPASS	300	10	random
M6_SFE1_300	512	0.25	${10}^6$	35	${10}^4$	0.002	BPASS	300	20	random
M6_SFE10_Zsun	512	0.25	${10}^6$	35	${10}^5$	0.02	BPASS	100	10	random
M6_SFE1_Zsun	512	0.25	${10}^6$	35	${10}^4$	0.02	BPASS	100	20	random
M6_SFE1_noTurb	512	0.25	${10}^6$	35	${10}^4$	0.02	BPASS	100	20	random	No turbulence
M6_SFE10_noSN	512	0.25	${10}^6$	35	${10}^5$	0.002	BPASS	100	10	random	No SNe
M6_SFE1_noSN	512	0.25	${10}^6$	35	${10}^4$	0.002	BPASS	100	20	random	No SNe
M6_SFE1_dSF	512	0.25	${10}^6$	35	${10}^4$	0.02	BPASS	100	20	dense
M6_SFE1_dSF_PLya	512	0.25	${10}^6$	35	${10}^4$	0.02	BPASS	100	20	dense	\Lya pressure

Table 3. Table 3: Summary of the simulation results. From left to right, each column represents the name, the final epoch of the simulation, the luminosity-weighted, time-averaged LyC escape fractions, the total number of LyC photons escaped, the luminosity-weighted, time-averaged \Lya  escape fractions, and the total number of \Lya  photons that escaped from the cloud. Note that the last two columns do not include the contribution from the \Lya  photons created by escaping LyC photons.

Name	$t_{\mathrm{final}}$		$N_{\mathrm{esc}}^{\mathrm{LyC}}$		$N_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$
	[Myr]	[%]	[${10}^{50}$]	[%]	[${10}^{50}$]
M6_SFE10	10	45.4	132.2	76.7	41.7
M6_SFE1	20	4.8	1.6	47.0	8.8
M5_SFE10	7	71.6	18.7	61.1	3.9
M5_SFE1	20	23.1	0.8	51.0	1.6
M6_SFE10_300	10	53.1	240.9	82.6	51.5
M6_SFE10_sng	10	30.0	60.3	76.5	39.6
M6_SFE10_Zsun	10	5.2	15.1	19.9	11.4
M6_SFE10_noSN	10	47.8	139.1	75.2	38.1
M6_SFE1_300	20	7.1	3.5	56.8	13.7
M6_SFE1_sng	20	0.4	0.08	35.7	4.4
M6_SFE1_Zsun	20	1.9	0.6	20.5	2.3
M6_SFE1_noTurb	20	8.5	2.8	31.3	5.4
M6_SFE1_noSN	20	3.9	1.3	53.3	9.8
M6_SFE1_dSF	20	1.1	0.4	26.1	5.4
M6_SFE1_dSF_Lya	20	13.0	4.3	63.1	10.9

Equations

$$\displaystyle\mbox{$f_{\rm esc}^{\rm LyC}$}(t)=\frac{\int_{\nu_{0}}^{\infty}d\nu\int d\Omega\sum_{i}\left\{L_{i}(\nu,t)\exp\left[-\tau_{i}(\nu,t;\Omega)\right]\right\}/4\pi h\nu}{\int_{\nu_{0}}^{\infty}d\nu\left[\sum_{i}L_{i}(\nu,t)\right]/h\nu},$$

$$\frac{1}{c_{\mathrm{s},\mathrm{i}}}\frac{d{r}_{\mathrm{i}}(t)}{dt}={\left(\frac{r_{\mathrm{S}}}{r_i(t)}\right)}^{\frac{3}{4}}-{\left(\frac{c_{\mathrm{ext}}}{c_{\mathrm{s},\mathrm{i}}}\right)}^2{\left(\frac{r_{\mathrm{S}}}{r_i(t)}\right)}^{-\frac{3}{4}}+\frac{v_{\mathrm{ext}}}{c_{\mathrm{s},\mathrm{i}}},$$

$${\tau }_{\mathrm{ion}}={\int }_1^{r_{\mathrm{cloud}}/r_{\mathrm{S}}}\frac{r_{\mathrm{S}}/c_{\mathrm{s},\mathrm{i}}}{y^{-3/4}-{\left(c_{\mathrm{ext}}/c_{\mathrm{s},\mathrm{i}}\right)}^2{y}^{3/4}}dy$$

$${ϵ}_{\mathrm{Ly}\alpha }={ϵ}_{\mathrm{rec}}+{ϵ}_{\mathrm{coll}}$$

$${ϵ}_{\mathrm{rec}}=P_B(T){\alpha }_B(T)n_e{n}_{\mathrm{HII}}{e}_{\mathrm{Ly}\alpha }$$

$${ϵ}_{\mathrm{coll}}=C(T)n_e{n}_{\mathrm{HI}}{e}_{\mathrm{Ly}\alpha }$$

$$P_B(T)=0.686-0.106\log T_4-0.009T_4^{-0.44},$$

$${\alpha }_B=2.753\times 10^{-14}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}\frac{{\lambda }^{1.5}}{{\left[1+{\left(\lambda /2.74\right)}^{0.407}\right]}^{2.242}},$$

$$C(T)=\frac{2.41\times 10^{-6}\mathrm{c}{\mathrm{m}}^3{\mathrm{s}}^{-1}}{T^{0.5}}{T}_4^{0.22}\exp \left[\frac{-1.63\times 10^{-11}}{k_{B}T}\right].$$

$$\displaystyle v_{\rm sep}\approx 320\,\mbox{${\rm km\,s^{-1}}$}\left(\frac{N_{\rm HI}}{10^{20}\,{\rm cm^{-2}}}\right)^{1/3}\left(\frac{T}{10^{4}\,{\rm K}}\right)^{1/6},$$

$$n_{\mathrm{ext}}(r)=n_0{\left(r/r_0\right)}^{-w},$$

$$\frac{1}{c_{\mathrm{s},\mathrm{i}}}\frac{d{r}_{\mathrm{i}}(t)}{dt}=F(r,t)-{\left(\frac{c_{\mathrm{ext}}}{c_{\mathrm{s},\mathrm{i}}}\right)}^2\frac{1}{F(r,t)}+\frac{v_{\mathrm{ext}}(r,t)}{c_{\mathrm{s},\mathrm{i}}},$$

$$F(r,t)$$

$$={\left(\frac{r_{\mathrm{Strom}}}{r_i(t)}\right)}^{\frac{3}{4}-\frac{w}{2}}.$$

$$n_i={\left(\frac{3{\dot{N}}_{\mathrm{ph}}}{4\pi r_{\mathrm{Strom}}^3{\alpha }_B}\right)}^{1/2}.$$

$$r_{\mathrm{stall}}=r_{\mathrm{Strom}}{\left(\frac{c_{\mathrm{s},\mathrm{i}}}{c_{\mathrm{ext}}}\right)}^{\frac{4}{3-2w}}.$$

$$\frac{dy}{dt}=\frac{c_{\mathrm{s},\mathrm{i}}}{r_{\mathrm{Strom}}}\left(y^{-3/4}-{\left(\frac{c_{\mathrm{s},\mathrm{ext}}}{c_{\mathrm{s},\mathrm{i}}}\right)}^2{y}^{3/4}\right)$$

$${\tau }_{\mathrm{ion}}(r)={\int }_1^{r/r_{\mathrm{Strom}}}\frac{r_{\mathrm{Strom}}/c_{\mathrm{s},\mathrm{i}}}{y^{-3/4}-{\left(c_{\mathrm{ext}}/c_{\mathrm{s},\mathrm{i}}\right)}^2{y}^{3/4}}dy.$$

$$\displaystyle M_{\rm cloud}\approx 228\,\mbox{$\rm M_{\odot}$}\,(R_{\rm cloud}/{\rm pc})^{2.36}.$$

Full text

1 Introduction

In a $\Lambda$CDM paradigm, the initial density perturbations develop into the large-scale cosmic web through gravitational interactions. Dark matter haloes (DMHs) form at the intersection of the filaments inside of which gas collapses and forms stars via radiative cooling. During this process, massive OB stars in dwarf-sized galaxies produce a large number of Lyman Continuum (LyC) photons that ionize neutral hydrogen in the Universe (Madau et al., 1999). As more structures collapse, the ionized bubbles expand and percolate, and the Universe is believed to be fully re-ionized by $z\approx 6$. The detection of strong Lyman absorption troughs in the spectra of quasi-stellar objects around $z=6$ (Fan et al., 2001, 2006) confirmed that the Universe was indeed opaque to ionizing radiation at $z\ga 6$, supporting this picture.

The physics behind the propagation of ionizing radiation in an expanding Universe is straightforward, and can be modelled using contemporary cosmological radiation-hydrodynamics (RHD) techniques (Wise & Cen, 2009; Kimm & Cen, 2014; Wise et al., 2014; Gnedin & Kaurov, 2014; Pawlik et al., 2015; Ocvirk et al., 2016; Xu et al., 2016; Kimm et al., 2017; Rosdahl et al., 2018; Finlator et al., 2018). Studies show that the first ionized bubbles appear with the emergence of massive Pop III stars in haloes of mass \sim 10^{6-7}\,\mbox{\rm M_{\odot}}, and then expand as subsequent metal-enriched populations provide additional LyC photons. Because the star formation histories in low-mass galaxies are quite intermittent, the ionized hydrogen in the dense circumgalactic (CGM) and intergalactic medium (IGM) recombines quickly, sometimes reducing the volume filling fraction of the ionized regions. Once galaxies become massive enough to host a number of star-forming clouds, they provide LyC photons continuously, and the bubbles grow from the over-dense regions into the void, while relatively dense filamentary gas remains self-shielded from the background radiation fields (e.g. Faucher-Giguère et al., 2010; Rosdahl & Blaizot, 2012; Chardin et al., 2018). These excess ionizing photons photo-heat gas and prevent it from collapsing on to small DMHs, which have a virial velocity less than \sim 10\,\mbox{{\rm km,s^{-1}}}, possibly delaying the growth of stellar mass in dwarf galaxies that we observe in the local Universe (e.g. Efstathiou, 1992; Gnedin, 2000; Somerville, 2002; Okamoto et al., 2008; Geen et al., 2013).

An important conclusion from numerical experiments is that the majority of the ionizing photons arise from dwarf-sized haloes. Based on the stellar mass-to-halo mass relations and the escape fractions obtained from RHDs, Kimm et al. (2017) solved the simple equation for reionization (Madau et al., 1999) and found that LyC photons from the halos of mass 10^{8}\,\mbox{\rm M_{\odot}}\la\mbox{M_{\rm vir}}\la 10^{11}\mbox{\rm M_{\odot}} must be included to match the end of the reionization epoch as well as the Thompson optical depth measured from the polarization signals of the cosmic microwave background. Improving upon this, Katz et al. (2018, 2019) developed a photon tracer algorithm that can follow the sources of ionization directly inside RHD simulations, and found that metal-poor galaxies ($\sim 0.001-0.1\,Z_{\odot}$) embedded in halos of mass 10^{8}\,\mbox{\rm M_{\odot}}\la M_{\rm halo}\la 10^{10}\,\mbox{\rm M_{\odot}} are likely to be mainly responsible for the reionization of the Universe.

However, the detailed process of how LyC photons interact with the gas within the star-forming clouds remains elusive. Several studies point out that the optical depth on $\sim$ 100 pc scales is already quite significant (Dove & Shull, 1994; Kim et al., 2013a; Kimm & Cen, 2014; Paardekooper et al., 2015; Trebitsch et al., 2017), meaning that a large fraction of LyC photons are absorbed inside the star-forming clouds. Currently, resolving the turbulent structure of molecular clouds is still challenging in galactic-scale simulations, and only a few simulations that adopt very high resolution ($\la$ 1 pc) begin to reproduce the basic observed properties of the star forming clouds, such as the linewidth-size relation (Larson, 1981; Heyer & Brunt, 2004), in idealised disk simulations (Hopkins et al., 2012a; Grisdale et al., 2018). Considering that the typical resolution of a reionization simulations is even lower ($\sim 10-1000\,{\rm pc}$), one can expect that the complex turbulent structure and corresponding Stromgren sphere inside the star-forming clouds is under-resolved. Thus, the leakage of the ionizing radiation may have been crudely approximated, potentially affecting the conclusions on the escape of LyC photons from the DMHs.

Recently, Dale et al. (2012) performed smoothed particle hydrodynamics simulations with a photo-ionization code to model the evolution of metal-rich star-forming clouds with different masses. They concluded that the escape fractions of LyC photons can be as high as $\sim 90\%$ if clouds of mass 10^{4}\,\mbox{\rm M_{\odot}}<M_{\rm cloud}<10^{6}\,\mbox{\rm M_{\odot}} are efficiently dispersed. Dale et al. (2013) further showed that the escape fractions are usually very high (\mbox{f_{\rm esc}^{\rm LyC}}\sim 0.2-0.9) if partially unbound clouds convert a significant fraction ($10-30\%$) of mass into stars. Similarly, Howard et al. (2018) ran RHD simulations with the flash Eulerian code and argued that the leakage of LyC photons is very significant ($\sim 65\%$ in clouds of mass 5\times 10^{4}\,\mbox{\rm M_{\odot}}<M_{\rm cloud}<10^{5}\,\mbox{\rm M_{\odot}} with $\sim 20\%$ star formation efficiencies (SFEs)), although the escape fraction decreases to $\la 10\%$ in a $10^{4}$ or 10^{6}\,\mbox{\rm M_{\odot}} cloud. Given that the majority of the ionizing radiation is produced before SNe explode (Leitherer et al., 1999; Bruzual & Charlot, 2003), the overall high escape fractions indicate that photo-ionization heating plays a critical role in clearing channels for the LyC photons on cloud scales (e.g. Krumholz et al., 2007; Dale et al., 2012; Peters et al., 2017; Geen et al., 2016; Gavagnin et al., 2017; Kim et al., 2018; Kannan et al., 2018). Using high-resolution (0.7 pc) cosmological RHDs with and without photo-ionization heating, Kimm et al. (2017) also confirmed that star-forming gas in metal-poor ($Z\sim 0.003\,Z_{\odot}$), dwarf-sized haloes (M_{\rm halo}\sim 10^{8}\,\mbox{\rm M_{\odot}}) is disrupted rapidly due to photo-heating even before SNe explode, leading to \mbox{f_{\rm esc}^{\rm LyC}}\sim 0.4 (see also Wise et al., 2014; Xu et al., 2016). Note that such a high escape fraction is indeed observed in compact starburst galaxies where a copious amount of ionizing radiation is being emitted (de Barros et al., 2016; Shapley et al., 2016; Bian et al., 2017; Izotov et al., 2018a; Vanzella et al., 2018, c.f. Leitherer et al. 2016; Puschnig et al. 2017). These corroborate that the propagation of photons during the early stage of star formation needs to be better understood to make firm predictions on the escape of LyC photons from galaxies.

In principle, it would be best to directly observe the ionizing part of the stellar spectrum from galaxies at $z\ga 6$ to determine the contribution of dwarf galaxies to the reionization of the Universe. However, few LyC photons would survive from absorption by neutral IGM at $z\geq 6$. An alternative method has thus been put forward to select the LyC leaking candidates and study their properties based on the profile of Lyman $\alpha$ (\Lya) line (Verhamme et al., 2015; Dijkstra et al., 2016). The basic idea is that because \Lya photons resonantly scatter with neutral hydrogen, the emerging line width depends on how gas is distributed around the source. If there exists low-density channels in which few LyC photons would be absorbed, \Lya would escape from the medium with little shift in frequency, and thus LyC leaking galaxies are likely to show a velocity profile narrower than 300\,\mbox{{\rm km,s^{-1}}} in \Lya. Given that \Lya is one of the strongest lines observed in the spectra of galaxies at high redshifts (e.g. Shapley et al., 2003), the approach seems promising and may be used in large observational programmes (e.g. Marchi et al., 2017; Steidel et al., 2018).

The beauty of \Lya is that it can be used not only to pre-select the LyC leakers, but also to infer the kinematics of the interstellar medium (ISM) in galaxies. As is well established in the literature (e.g. Ahn et al., 2003; Verhamme et al., 2006; Dijkstra et al., 2006; Barnes et al., 2011), an expanding medium would preferentially absorb photons with frequencies shorter than the line centre, resulting in a spectrum with a pronounced red peak. This type of profile is often observed in star-forming, Lyman break galaxies (e.g. Steidel et al., 2010; Kornei et al., 2010), and may be used to estimate the amount of galactic outflows. However, because \Lya photons primarily arise from young star-forming regions through recombination in gas ionized by LyC radiation, it is necessary to understand how \Lya photons are created and propagate on cloud scales before they interact with the ISM. Unfortunately, current models make use of either idealized environments, such as uniform or clumpy distributions, or gas distributions from galactic scale simulations where ISM structures are under-resolved. Little work has been done thus far based on cloud simulations where internal turbulent structures are well resolved. Despite the success at reproducing the overall features of \Lya profiles (Verhamme et al., 2008; Gronke, 2017), the simple models may be improved to better interpret the line shape and to make more accurate predictions to find the LyC leaking candidates by studying the propagation of \Lya photons inside star-forming clouds.

To this end, we perform high-resolution RHD simulations of turbulent gas clouds with feedback from supernovae and stellar radiation. The aim of these experiments is three-fold. First, we attempt to understand the absorption and propagation of LyC photons in clouds with complex turbulent structures by varying the mass of the cloud and the total SFE. Second, there are several uncertainties regarding stellar evolution, such as the maximum mass of stars (Crowther et al., 2010) or the evolution of the spectral energy distributions (SEDs) due to binary interactions (e.g. Stanway et al., 2016). These can affect the predictions of the reionization of the Universe as well as the SFE inside the cloud (Geen et al., 2018), and thus it is necessary to understand what level of uncertainty is implicitly inherited from our limited understanding of stellar evolution. Finally, we aim to examine the emergent \Lya profiles from the clouds and compare them with previous results so that the LyC candidates are more efficiently pre-selected in observations. The \Lya profiles obtained from this work may be used to model the fraction of bright Lyman alpha emitters (LAEs) in semi-numerical approaches (e.g. Choudhury et al., 2015; Mesinger et al., 2015; Weinberger et al., 2018), and to infer the relevance of the demise in the fraction of strong LAEs during the epoch of reionization (e.g. Stark et al., 2010; Treu et al., 2013; Schenker et al., 2014).

This paper is organized as follows. In Section 2, we present the initial conditions and input physics of the simulations. Section 3 presents how the escape fractions evolve in a cloud with different SFE and stellar spectra. Section 4 discusses the effects of turbulent structures, the connection between the escape fractions of LyC and \Lya photons, the implications for reionization, and the impact of star formation and feedback schemes. The summary and conclusions are given in Section 5.

2 Simulations

We perform 15 RHD simulations with different cloud masses, SFEs, SEDs, metallicities, and feedback to examine the escape fraction of LyC and \Lya photons on cloud scales using ramses-rt (Teyssier, 2002; Rosdahl et al., 2013; Rosdahl & Teyssier, 2015). The Euler equations are solved using the HLLC method with the positivity-conserving slope limiter (Toro et al., 1994). We adopt a Courant number of 0.7. The Poisson equation is evolved using the multi-grid method (Guillet & Teyssier, 2011). A uniform UV background is turned on with the self-shielding approximation, such that gas denser than \mbox{n_{\rm H}}\simeq 0.01\,\mbox{{\rm cm^{-3}}} is not affected by heating (Rosdahl & Blaizot, 2012). In order to model the photoelectric heating as well as the transport of ionizing radiation, we use the GLF solver with eight photon groups, as detailed in Table 1. We adopt the frequency-dependent cross-sections (Katz et al., 2017) based on Rosdahl et al. (2013) and Baczynski et al. (2015). The evolution of seven chemical species (HI, HII, HeI, HeII, HeIII, H2, and $e^{-}$) is followed by solving photo-chemistry equations (see Katz et al., 2017; Kimm et al., 2017, for details). The speed of light is reduced to $10^{-3}\,c$ to reduce the computation, where $c$ is the full speed of light.

The initial conditions follow Gaussian density distributions with the maximum densities of \mbox{n_{\rm H}}=200\,\mbox{{\rm cm^{-3}}} and 72\,\mbox{{\rm cm^{-3}}} and a 1$\sigma$ radius of 10 pc and 30 pc for clouds of gas mass 10^{5}\,\mbox{\rm M_{\odot}} and 10^{6}\,\mbox{\rm M_{\odot}}, respectively. We then add Kolmogorov turbulence with a power spectrum of the form $\propto k^{-5/3}$, adopting the mixture of the solenoidal (60%) and compressive (40%) mode for one free-fall time (4 and 5 Myr, respectively). The resulting turbulent energy is about $\approx$80% of the gravitational binding energy inside the half-mass radius, and thus the simulated clouds are marginally gravitationally bound initially. At later epochs, turbulence is generated either by radiation feedback or SN explosions. The clouds are assumed to be metal-poor ($Z_{\rm gas}=0.1\,Z_{\odot}$), motivated by the fact that the metallicity of galaxies typically observed at high redshift is low ($\sim 0.2\,Z_{\odot}$, Pettini et al., 2000; Song et al., 2014; Bouwens et al., 2016; Tamura et al., 2018) and that reionization is likely to be driven by metal-poor dwarf galaxies (e.g. Katz et al., 2019). Note that this is one of the main differences of our work compared to previous studies where simulated clouds are preferentially metal-rich (e.g. Dale et al., 2012; Howard et al., 2018; Kim et al., 2018; Geen et al., 2018). However, we also test the solar metallicity case to determine the effects of metallicity on the escape of LyC and \Lya photons.

The simulations include five different forms of stellar feedback, i.e., photo-ionization heating by UV photons (Rosdahl et al., 2013), direct radiation pressure111We note that momentum transfer to the neighbouring cells may be under-estimated when the Stromgren sphere is not properly resolved. To alleviate this, we adopt a correction scheme that takes into account the isotropic flux (rt_isoPress) (Rosdahl & Teyssier, 2015). by UV and optical photons (Rosdahl & Teyssier, 2015), radiation pressure by multiple scatterings of IR photons (Rosdahl & Teyssier, 2015), photo-electric heating on dust (Katz et al., 2017; Kimm et al., 2017), and Type II supernova feedback. We model supernova feedback using the mechanical scheme (Kimm & Cen, 2014; Kimm et al., 2015, see also Hopkins et al. 2014, Smith et al. 2018, Lupi 2019 for a similar method), which is designed to ensure the correct momentum input to the surroundings (Thornton et al., 1998), with a realistic time delay between 4 and 40 Myr based on the lifetime of massive (8-100\,\mbox{\rm M_{\odot}}) main-sequence stars (Leitherer et al., 1999). We assume an SN frequency of 0.011\,\mbox{\rm M_{\odot}}^{-1}, appropriate for the Kroupa initial mass function (IMF, Kroupa, 2001). Each SN event ejects metals of mass 1.4 $\rm M_{\odot}$ (i.e., metallicity of the ejecta is set to 0.075), and we assume the initial ejecta energy to be $10^{51}\,{\rm erg}$. The effect of photo-electric heating on dust is included by explicitly following the propagation and absorption of the photons at the Habing band (see Table 1, Kimm et al., 2017). We assume that the amount of dust is proportional to the amount of metals with a dust-to-metal ratio of 0.4 (e.g. Dwek, 1998; Draine et al., 2007) at $T<10^{6}\,{\rm K}$. At higher temperatures, dust is assumed to be destroyed as a crude approximation for the thermal sputtering process.

The fiducial model uses the binary star SEDs from Stanway et al. (2016), but we also examine the SEDs with single stellar evolution (Bruzual & Charlot, 2003). These SEDs are generated with a lower (upper) limit of mass of 0.1\,\mbox{\rm M_{\odot}} (100\,\mbox{\rm M_{\odot}}) and are used to compute the instantaneous luminosity as a function of age, in 8 photon groups representing radiation energy (or frequency) intervals as shown in Table 1. We also test SEDs generated with the upper mass limit of 300\,\mbox{\rm M_{\odot}} with binary stellar evolution to examine the importance of the IMF on the predictions of the escape fractions. Note that the cumulative number of LyC photons from metal-poor populations can vary by a factor of $\sim 2$ depending on this assumption of the stellar population (Figure 1) (see also Topping & Shull, 2015, for different estimates of the total number of LyC photons arising from variations in stellar IMF, metallicity, and rotation rates). We summarize the initial conditions of the simulations in Table 2.

We consider two cases of SFE, 1% and 10% SFE, where we place star particles in the cloud comprising 1% and 10% of the cloud mass, respectively. Recent studies suggest that the total stellar mass formed in the simulations may depend on how radiation pressure is modelled (Hopkins & Grudic, 2018; Kim et al., 2018; Krumholz, 2018) and how turbulence is driven (Geen et al., 2018). It can also be affected by the inclusion of other feedback processes, such as radiation pressure by \Lya photons (Dijkstra & Loeb, 2008; Smith et al., 2017; Kimm et al., 2018) or re-processed infrared photons (Hopkins et al., 2012b; Krumholz & Thompson, 2012; Rosdahl et al., 2015; Skinner & Ostriker, 2015; Tsang & Milosavljević, 2015). In order to circumvent these uncertainties, we simply start with stars as an initial condition instead of modelling formation and accretion onto sink particles (e.g. Bate et al., 1995). Specifically, we randomly place star particles of mass 10^{3}\,\mbox{\rm M_{\odot}} in the inner region of the cloud following the Gaussian distribution of the dispersion of 5 pc. The resulting mean radius of the stellar distributions is $\approx$ 3 pc, and the average gas density of their host cells is initially \left<\mbox{n_{\rm H}}\right>\sim 10–20\,\mbox{{\rm cm^{-3}}}.

We stop the simulations when $\approx 98\%$ of the ionizing radiation is emitted from a simple stellar population with binaries ($t=20\,{\rm Myr}$) for the 1% SFE cases. This is chosen to minimize the computational cost, as covering 99% of the total LyC radiation would require the modelling of additional $\approx 20\,{\rm Myr}$. For the runs with 10% SFEs, we choose the final snapshot in which the instantaneous escape fraction reaches near 100%, which corresponds to 7 and 10 Myr for the M_{\rm cloud}=10^{5}\,\mbox{\rm M_{\odot}} and 10^{6}\,\mbox{\rm M_{\odot}} cloud, respectively. Note that a large number of snapshots ($\approx 70-100$) is generated with a time interval of $\Delta t=0.1\,{\rm Myr}$, $0.2\,{\rm Myr}$, and $0.25\,{\rm Myr}$ at $0\leq t\leq 1\,{\rm Myr}$, $1\leq t\leq 5\,{\rm Myr}$, $5\leq t\leq 20\,{\rm Myr}$, respectively, for an accurate determination of the luminosity-weighted time average of escape fractions.

3 Results

3.1 Evolution of the clouds

We begin our discussion by describing the general features of the evolution of the turbulent clouds. Figure 2 shows the temperature distributions of the clouds with different SFEs and cloud masses. Ionizing photons are initially completely absorbed by the host cells of young star particles, imparting momentum of $L_{\rm ion}/c$, where $L_{\rm ion}$ is the luminosity of ionizing radiation. The absorption of longer wavelength photons by dust also transfers momentum to the surroundings, but their contribution to the total direct radiation pressure is sub-dominant near young stellar populations (e.g. Leitherer et al., 1999). The Stromgren sphere develops shortly after several recombination time scales ($\sim 10^{-3}-10^{-2}\,{\rm Myr}$), which over-pressurises the gas surrounding the radiation source. Once photo-ionization heating creates low-density channels, the ionization bubble expands faster, as gas at larger radii is more tenuous. Between 4 Myr and 40 Myr, SNe explode intermittently, enhancing outflows that were originally accelerated by photo-ionization heating and direct radiation pressure.

We find that the cloud with a larger number of stars is disrupted more quickly than that with a low SFE. The turbulent cloud is destroyed as early as $\sim$ 1–3 Myr in the case of the model with a 10% SFE (M6_SFE10 and M5_SFE10), indicating that radiation feedback is strong enough to blow out the gas (Dale et al., 2012; Dale et al., 2014; Walch et al., 2012; Geen et al., 2016; Kim et al., 2017). Very dense, star-forming gas is all destroyed, and approximately 80% and 95% of the gas leaves the simulated domain by the end in the M6_SFE10 and M5_SFE10 runs, respectively (Figure 3). As discussed extensively in the literature, photo-ionization is mainly responsible for this process (e.g. Matzner, 2002), while direct radiation pressure is more significant in the dense regime (\mbox{n_{\rm H}}\ga 10^{5}\,\mbox{{\rm cm^{-3}}}) (e.g. Rosdahl et al., 2015; Kimm et al., 2017; Kim et al., 2018). Figure 3 indeed shows that the majority of the gas in the cloud is quickly ionized and accelerated to the velocity that significantly exceeds the escape velocity of the clouds (\sim 5-10\,\mbox{{\rm km,s^{-1}}}). Non-thermal radiation pressure from multi-scattered IR photons is not expected to play an important role in regulating the overall dynamics of the cloud, as they are not effectively trapped due to the low metallicity and low optical depths in the cloud. Skinner & Ostriker (2015) also show that even for the metal-rich environments, the opacity must be significantly higher ($\kappa_{\rm IR}>15\,{\rm cm^{2}/g}$) to unbind the cloud and reduce the amount of star formation as observed in some super star clusters.

Figure 2 also shows that some of the dense parts of the clouds survive the radiation feedback for a long period of time ($\sim 10-20\,{\rm Myr}$), creating comet-like structures in the runs with a 1% SFE (M6_SFE1 and M5_SFE1). This happens because the pressure of the dense gas is significantly higher ($P/k_{B}\sim 10^{5-7}\,{\rm cm^{3}\,K^{-1}}$) than the ram pressure from the warm/hot gas in the HII bubble ($P/k_{B}\sim 10^{3-4}\,{\rm cm^{3}\,K^{-1}}$). Once the comet-like structures move away from the centre where the bright sources are located, they become more difficult to photo-evaporate, as their solid angle on the sky becomes smaller and fewer LyC photons can interact with them. As a result, some of the clouds survive even though the majority of the gas is blown out from the system, and the velocity dispersion of the neutral gas (\sigma\la 10\,\mbox{{\rm km,s^{-1}}}) does not increase as significantly as that of the ionized gas (\sigma\ga 10-50\,\mbox{{\rm km,s^{-1}}}, Figure 3). In contrast, if the stellar population produces enough LyC photons ($\approx 10^{65}$) to keep the dense clumps ionized (i.e., 10% SFE), the dense clumps are quickly evaporated (the fourth column of the first and third rows in Figure 2).

3.2 Escape of LyC photons

Now that we understand the general features of the simulated turbulent clouds, we study where and when most LyC photons are absorbed in different environments. We also examine the effects of the shape of SEDs and the gas metallicity on the escape of LyC photons in this section.

3.2.1 Escape fractions of LyC radiation

In order to measure the escape of LyC photons, we post-process the simulation snapshots with a simple ray tracing method. This is done by casting 12288 ($=12\times 32^{2}$) rays per star particle using the healpix algorithm (Górski et al., 2005) and by measuring the remaining photons at the computational boundary after attenuation due to $\rm HI$, $\rm H_{2}$, $\rm HeI$, $\rm HeII$, and dust, as

[TABLE]

where $L_{i}(\nu,t)$ is the luminosity of the $i$-th star particle with age $t$, $\tau_{i}(\nu,t;\Omega)$ is the total optical depth due to dust and gas at a given frequency $\nu$ along the sight line $\Omega$, and $\nu_{0}=3.287\times 10^{15}\,{\rm Hz}$ is the frequency at the Lyman limit. Each ray carries an SED depending on the age and metallicity, and also on the stellar evolutionary model assumed (i.e., single vs binary). The optical depth to neutral hydrogen and singly ionized helium is computed using the photo-ionization cross-section of a hydrogenic ion with the nuclear charge $Z$ (Osterbrock & Ferland, 2006, Equation 2.4). For the photo-ionization cross-section of neutral helium and molecular hydrogen, we adopt the fitting formula from Yan et al. (1998, 2001). Also included is the absorption due to either small Magellanic cloud-type dust based on the fitting formula to the effective cross-section (Gnedin et al., 2008) or Milky Way-type dust (Weingartner & Draine, 2001) with a dust-to-metal ratio of 0.4 (Draine & Li, 2007) for the metal-poor ($Z_{\rm gas}=0.002$) or metal-rich ($0.02$) cloud, respectively. Note that the amount of dust is assumed to depend on the fraction of ionized hydrogen, as $f_{\rm dust}=(1-x_{\rm HII})+f_{\rm ion}x_{\rm HII}$, where $f_{\rm dust}$ is the relative mass fraction of dust with respect to the purely neutral case, $x_{\rm HII}$ is the mass fraction of ionized hydrogen, and $f_{\rm ion}=0.01$ is the free parameter that takes into account the observed abundance of dust in HII regions (see Laursen et al., 2009, for detailed discussion).222We note that this is not fully self-consistent with the modelling of dust in our RHD calculations. However, because our code does not follow the formation and destruction of dust explicitly, we post-process our simulation outputs using the simple dust model of Laursen et al. (2009), which better reflects the observed dust abundance in a variety of environments and reproduces \Lya properties of high-z galaxies.

In Figure 4, we show the escape of LyC photons in clouds of different SFEs. The escape fraction ($f_{\rm esc}^{\rm LyC}$) generally increases with time, reaching $f_{\rm esc}^{\rm LyC}\approx 100\%$ when the cloud is entirely disrupted. One can also see that the variation in $f_{\rm esc}^{\rm LyC}$ is more monotonic when a SFE is higher, as the disruption process is more efficient (c.f. Howard et al., 2018). In contrast, the run with a 1% SFE (producing only one star particle of mass 10^{3}\,\mbox{\rm M_{\odot}}, M5_SFE1) predicts noisy $f_{\rm esc}^{\rm LyC}$. This happens because, although it is initially placed in a relatively low-density environment, the only star particle encounters and is swallowed by adjacent dense clumps inside the giant molecular cloud (GMC). Ionizing radiation is temporarily blocked, but because photo-ionization heating and radiation pressure from the star are strong enough to destroy these local clumps, the high escape fraction is quickly recovered. The sign of the noisy feature can also be found in the M6_SFE1 run where ten star particles are emitting ionizing radiation.

As expected, we find that the time-averaged, luminosity-weighted escape fraction, $\left<f_{\rm esc}^{\rm LyC}\right>$, is a strong function of SFE (Dove et al., 2000). In the massive cloud with a large number of stars (M6_SFE10), $45.4\%$ of the LyC photons escape from the GMC, while only a small fraction ($4.8\%$) of LyC photons leave in the case of a 1% SFE (M6_SFE1)333We have examined the uncertainty due to randomizing the initial positions of star particles by running two additional simulations with a 1% SFE, and found that the luminosity-weighted escape fraction can change by $\sim 30\%$ ($\left<f_{\rm esc}^{\rm LyC}\right>=$3.2% and 5.2%).. The same trend is seen for the cloud with less mass (M_{\rm cloud}=10^{5}\,\mbox{\rm M_{\odot}}, $\left<f_{\rm esc}^{\rm LyC}\right>=71.6\%$ vs $22.9\%$). This can be attributed to the fact that the most important mechanism during the early phase, i.e., photo-ionization heating, can impart radial momentum in proportion to $N_{\rm ion}^{4/7}$ (Krumholz, 2015). This lends support to the claim that the escape fraction relies sensitively on the burstiness of star formation histories (Kimm et al., 2017; Trebitsch et al., 2017).

It is also interesting to note that less massive clouds show higher escape fractions for a given SFE. Because the initial average density is chosen to be higher by a factor of $\sim 2$ in the less massive cloud (\left<n_{\rm H}\right>\approx 160\,\mbox{{\rm cm^{-3}}}), the radial momentum from photo-ionization heating (p_{\rm rad}\propto\mbox{n_{\rm H}}^{-1/7}, Krumholz, 2015) as well as the momentum from SN explosion (p_{\rm rad}\propto\mbox{n_{\rm H}}^{-2/17}, Blondin et al., 1998; Thornton et al., 1998) would be smaller in the less massive cloud; hence, one may expect the escape fraction to be lower. However, this is opposite to our findings. This may simply be due to the fact that the less massive cloud happens to have low-density channels around young star particles, given that they would encounter a fewer number of neutral hydrogen atoms than in the run with the massive cloud. The high escape fractions can also be attributed to the fact that ionization fronts reach the edge of the cloud earlier in the less massive cloud. Geen et al. (2015a) show that, for a uniform medium, the propagation velocity of ionization fronts can be written as

[TABLE]

where $r$ is the radius, $r_{S}$ is the Stromgren sphere radius, $r_{i}$ is the radius of the ionization front, $t$ is the time, $c_{s,i}$ is the sound speed of the ionized medium, $v_{\rm ext}$ is the infall velocity of the ambient gas, and $c_{\rm ext}$ represents the velocity term due to the thermal and turbulent pressure. Then, the time for the front to reach the edge of the cloud ($r_{\rm cloud}$) may be written as

[TABLE]

where $y\equiv r_{\rm i}/r_{\rm S}$. As shown in Appendix (Figure 18), although the number of ionizing photons per atom is the same for a given SFE, the time required for the ionization front to propagate to the edge of the cloud is shorter in the cloud with M_{\rm cloud}=10^{5}\,\mbox{\rm M_{\odot}} compared to the one with M_{\rm cloud}=10^{6}\,\mbox{\rm M_{\odot}}. As a consequence, the bright phase of escaping ionizing radiation starts earlier in the less massive cloud, leading to a higher $\left<f_{\rm esc}^{\rm LyC}\right>$.

The escape fractions obtained from our simulations are useful to understand the reionization history of the Universe, but may not be practical to compare with observations. For example, because it is not possible to observe the whole wavelength range of LyC photons, $f_{\rm esc}^{\rm LyC}$ is often measured by comparing the average flux densities at two different wavelengths, such as $F_{900}$ (flux density at wavelength $880<\lambda<910$ Å ) and $F_{1500}$ (flux density at $1475<\lambda<1525$ Å) (e.g. Steidel et al., 2018). Motivated by this, we also present the escape fractions of the photons with wavelength $880<\lambda<910$ ($f_{\rm esc}^{900}$) in the middle panels of Figure 4. We find that $f_{\rm esc}^{900}$ is systemically smaller than $f_{\rm esc}^{\rm LyC}$ because the absorption cross-section to ionizing radiation due to neutral hydrogen is larger at longer wavelengths below the Lyman edge (Osterbrock & Ferland, 2006). The runs with the massive metal-poor cloud (M_{\rm cloud}=10^{6}\,\mbox{\rm M_{\odot}}) show $\left<f_{\rm esc}^{900}\right>$=3.1% and 37% for a 1% and 10% SFE, respectively, while 17% and 65% of the photons around 900Å escape from the less massive cloud. Therefore, one should keep in mind that the luminosity-weighted time average of $f_{\rm esc}^{900}$ can be $\sim 20-35\%$ lower than $\left<f_{\rm esc}^{\rm LyC}\right>$ when comparing the simulated escape fractions with observationally derived quantities. Corresponding luminosity-weighted flux density ratios ($F_{900}/F_{1500}$) are found to be 0.05, 0.56, 0.24, and 1.1 by the end of each simulation (M6_SFE1, M6_SFE10, M5_SFE1, and M5_SFE10, respectively). Note that the flux density ratios are somewhat higher than those observed in the compact starbursts (e.g. Shapley et al., 2016), likely because our clouds are very young and have no underlying population that preferentially produces the flux density around 1500 Å.

To compute the relative contributions to the absorption of LyC photons, we also measure the effective optical depth due to $\rm HI$, $\rm H_{2}$, $\rm HeI+HeII$, and dust, as $\tau_{\rm eff}\equiv-\ln f_{\rm esc,i}^{\rm LyC}$, where $f_{\rm esc,i}^{\rm LyC}$ is the instantaneous escape fraction after attenuation by each element $i$. Figure 5 (top panels) shows that the majority of the LyC photons are absorbed by neutral hydrogen in the metal-poor, massive clouds. Initially, the effective optical depth due to dust is also significant ($\tau_{\rm eff,dust}\sim 1$) but decreases steadily, as the fraction of ionized hydrogen increases and dust is assumed to be destroyed in the HII region. The absorption due to helium and molecular hydrogen is also minor compared to that due to neutral hydrogen, but we find that their total contribution is as equally important as dust even in the metal-rich environments (bottom panels).

3.2.2 Impact of SED models on the escape of LyC radiation

We find that the prediction of the escape fractions relies on the choice of the SEDs used in the simulations. Figure 6 (top two panels) shows that $\left<f_{\rm esc}^{\rm LyC}\right>$ in the run adopting a single stellar evolution model (Bruzual & Charlot, 2003) is lower than the run with binaries. This can be attributed to the fact that the former predicts an ionizing photon production rate that is a factor of $\sim 2$ smaller than the model with binaries at $t\ga 5\,{\rm Myr}$ (Figure 1). In particular, mergers and transfer of gas between binaries results in enhanced ionizing emission at $t\ga 3\,{\rm Myr}$ (Stanway et al., 2016, see also Figure 1). This allows LyC photons to escape more easily, as the clouds become significantly disrupted at this stage. In contrast, few ionizing photons are generated after $10\,{\rm Myr}$ in the single SED case (M6_SFE1_sng), leading to an order of magnitude smaller $\left<f_{\rm esc}^{\rm LyC}\right>$ of $\sim 0.4\%$ than $\left<f_{\rm esc}^{\rm LyC}\right>=4.8\%$ of the M6_SFE1 model. The difference becomes less notable if the SFE is high enough to disrupt the cloud in the early phase during which a large number of LyC photons are still produced. For example, the runs with a 10% SFE (M6_SFE10 and M6_SFE10_sng) yield $\left<f_{\rm esc}^{\rm LyC}\right>=45.4\%$ and $30.0\%$, respectively.

Uncertainties persist in the SED models as the maximum stellar mass remains unknown (e.g. Kroupa et al., 2013). Previous studies suggest that the initial mass functions do not extend to more than 120-200\,\mbox{\rm M_{\odot}} (Massey & Hunter, 1998; Oey & Clarke, 2005), but re-analysis of O stars in 30 Doradus suggests the possibility of 300\,\mbox{\rm M_{\odot}} stars (Crowther et al., 2010). To assess the impact of the variation in the upper mass limit on escape fractions, we perform two additional simulations (the $10^{6}$M*⊙* cloud with 1% and 10% SFE) adopting the binary SEDs with a maximum mass of 300\,\mbox{\rm M_{\odot}}. Note that this SED provides an additional 50% of LyC photons in the early phase ($t\la 5\,{\rm Myr}$) compared to our fiducial case (Figure 1). We find that the choice of the upper mass limit seems to have a secondary effect compared to the inclusion of binaries (Figure 6). With the higher upper mass, the escape fractions are slightly enhanced from $\left<f_{\rm esc}^{\rm LyC}\right>=4.8\%$ to $5.2\%$ in the run with a 1% SFE, and from $45.4\%$ to $53.1\%$ in the cloud with a 10% SFE. This indicates that the additional photons from very massive stars (100\leq M/\mbox{\rm M_{\odot}}\leq 300) are preferentially used to build up the initial HII bubbles and help to create the low-density channels instead of leaving the clouds.

3.2.3 Impact of gas metallicity on the escape of LyC radiation

The bottom two panels of Figure 6 show the effects of gas metallicity on the escape of LyC photons. Compared to the metal-poor case ($Z_{\rm gas}=0.002$), the escape from the massive cloud runs with solar metallicity (M6_SFEx_Zsol) is significantly reduced by a factor of 2.5 from $\left<f_{\rm esc}^{\rm LyC}\right>=4.8\%$ to $1.9\%$ for a cloud with 1% SFE or by an order of magnitude from $45.4\%$ to $5.2\%$ in the runs with a 10% SFE. Note that we use the same metallicity ($Z_{\rm star}=0.002$) for star particles in order to keep the number of ionizing photons the same for all simulations with the same cloud mass and SFE; thus, the decrease in $\left<f_{\rm esc}^{\rm LyC}\right>$ is due to the enhanced attenuation by dust and/or more efficient metal cooling. To identify the cause of the significant reduction, we further examine $\left<f_{\rm esc}^{\rm LyC}\right>$ from the metal-rich runs, assuming that the dust-to-metal ratio is smaller by an order of magnitude (i.e. $D/M=0.04$) in the post-processing step, and find that the LyC escape is still considerably smaller ($\left<f_{\rm esc}^{\rm LyC}\right>=1.9\%$ (SFE1) or $5.0\%$ (SFE10)) than in the low-metallicity clouds (see also Figure 5 for the relative contribution to the optical depth by dust). This demonstrates that the effects due to dust absorption are secondary on cloud scales and that the low escape fraction arises mainly because ionized hydrogen quickly recombines as a result of efficient metal cooling. Indeed, we find that the enhanced cooling limits the effects of photo-ionization heating and that the cloud is not dispersed as efficiently as in the metal-poor case (compare Figure 2 with Figure 20 in Appendix). These experiments imply that the observed escape fractions from the relatively metal-rich LBG galaxies at $z\sim 3$ (Steidel et al., 2001; Siana et al., 2010; Leitet et al., 2013; Mostardi et al., 2015; Grazian et al., 2016; Marchi et al., 2017; Steidel et al., 2018) are likely to be a lower limit of the escape fractions of high-redshift metal-poor dwarf galaxies that likely reionized the Universe.

3.2.4 Absorption scale of LyC radiation

In galactic scale simulations with finite resolution, the turbulent structure of the star-forming clouds is under-resolved, and the estimation of LyC escape from these simulations is often uncertain (e.g. Ma et al., 2016). In this regard, it is useful to compute where most of the LyC photons are absorbed in different environments. In Figure 7, we measure the absorption scale based on the ray tracing method described in the previous section and show that half of the LyC photons are absorbed on small scales, particularly when the SFE is low ($1\%$, M6_SFE1). More specifically, the photon number-weighted absorption scales for this cloud are 8, 33, and 83 pc for $f_{\rm abs}=50$, $90$, and $99\%$, respectively. When the star formation becomes more intense (10% SFE), 50%, 90%, and 99% of the LyC photons are absorbed at much larger distances (49, 140, and 191 pc). This suggests that in order to determine the instantaneous escape fractions with reasonable accuracy (within a factor of two), one should adopt a computational resolution better than $\Delta x_{\rm min}\la{\rm 10\,pc}$, provided that the formation and disruption of star-forming clouds are reasonably well captured in the simulations. Based on high-resolution (2 and 4 pc) cosmological RHD simulations adopting strong supernova feedback, Kimm & Cen (2014) also show that the optical depth to LyC photons in their atomic-cooling haloes is large on 100 pc scales ($\tau\sim 2-4$) (see also Paardekooper et al., 2015; Kim et al., 2013b) and that the luminosity-weighted escape fractions are converged at these resolutions.

3.2.5 Variability of the LyC escape fractions

Our simulations with SFEs greater than 1%, which is what observations appear to support on average (Lada et al., 2010; Evans et al., 2014; Vutisalchavakul et al., 2016; Lee et al., 2016, c.f. Leroy et al. 2017), suggest that $f_{\rm esc}^{\rm LyC}$ tends to increase monotonically once radiation blows the cloud gas away and develops low-density channels (see also Kim et al., 2018). This may sound contradictory to previous findings that $f_{\rm esc}^{\rm LyC}$ fluctuates rapidly over time in galaxies (Wise et al., 2014; Kimm & Cen, 2014; Paardekooper et al., 2015; Ma et al., 2016; Trebitsch et al., 2017). However, the variation observed in the galactic scale simulations is on the timescale of $\sim 10-30\,{\rm Myr}$, which is not inconsistent with the timescale of the variability in our simulations (c.f. the runs with a 1% SFE).

An exception is when stars form in deeply embedded environments where stellar feedback cannot disrupt the surrounding clumps early on (Howard et al., 2018). If this were the case, the escape fractions would be highly time-dependent even on Myr timescales, and the predictions from galactic scale simulations (e.g. Wise & Cen, 2009; Kimm & Cen, 2014; Xu et al., 2016; Trebitsch et al., 2017; Rosdahl et al., 2018), where small-scale structures are unresolved, might be quite uncertain. However, as we discuss later in Section 4.4, radiation feedback in dense regions needs to be properly addressed, especially when the Stromgren sphere is under-resolved. For example, none of the cloud simulations conducted thus far include \Lya feedback, which may disrupt the metal-poor dense clumps near young stars (Kimm et al., 2018). The inclusion of such strong early feedback leads to the efficient expansion of photoionized bubbles and thus results in a rather monotonic evolution of the escape fractions, as in our fiducial runs where star particles are randomly placed inside the cloud. For these models, we do not expect extremely variable escape fractions even if we resolve the detailed structure of the clouds (see Section 4.4), and the fluctuating escape fractions on $10$–$30\,{\rm Myr}$ are likely to persist as a ramification of the sporadic nature of star formation episodes distributed over the galaxy.

3.3 Properties of \Lya photons

We now turn to the scattering and absorption of \Lya photons in the simulated clouds. Note that the propagation of \Lya photons in a neutral medium is considerably different from LyC photons in the sense that the interaction with neutral hydrogen does not destroy \Lya photons but simply changes their frequency and direction.

To model the propagation of \Lya photons in star-forming clouds, we post-process our simulations using the Monte Carlo \Lya radiative transfer code, rascas (Michel-Dansac et al. in prep). We compute the \Lya emissivity by taking into account recombination and collisional radiation, as

[TABLE]

where $e_{\rm Ly\alpha}$ is the energy of \Lya photon (10.16 eV), $n_{e}$ and $n_{\rm HII}$ are the number density of electron and ionized hydrogen. Here, $P_{B}$ is the probability for an absorbed LyC photon to be re-emitted as a \Lya photon (Cantalupo et al., 2008),

[TABLE]

where $T_{4}=T/10^{4}\,{\rm K}$, $\alpha_{B}$ is the case B recombination coefficient (Hui & Gnedin, 1997),

[TABLE]

with $\lambda=315614\,{\rm K}/T$, and $C(T)$ is the coefficient for the cooling radiation (e.g. Callaway et al., 1987)

[TABLE]

Based on the emissivity of each cell, we randomly sample the initial position for $10^{4}$ \Lya photons. Once the position is determined, the initial frequency is drawn randomly from a Gaussian distribution with the thermal Doppler broadening set by the temperature of each cell. The Doppler parameter ($\left<b\right>=\sqrt{2k_{B}T/m_{\rm H}}$) at the source position is approximately $13$–15\,\mbox{{\rm km,s^{-1}}} for the runs that we examine in this work. We include scattering due to deuterium with a fixed abundance of $3\times 10^{-5}$, recoil effects, and the scattering and destruction due to dust based on Laursen et al. (2009). The loss of \Lya photons due to molecular hydrogen (Shull, 1978; Black & van Dishoeck, 1987) is neglected for simplicity. As such, molecular hydrogen is transparent to \Lya photons in this study, and the resulting escape fractions may be slightly over-estimated, especially in the warm ($10^{3}\la T\la 10^{4}\,{\rm K}$), molecular regions, although it is unlikely to be significant in the typical cold ($T\sim 100\,{\rm K}$) star-forming sites (see Figure 20 in Neufeld, 1990).

3.3.1 Production of \Lya radiation

Figure 8 shows an example of \Lya emissivity maps from the M6_SFE1 run at two different times (c.f. Figure 2). Initially, LyC photons are well confined within the cloud, and \Lya emissivity map closely follows the distributions of young stellar particles (i.e., LyC photon density). At this stage, most \Lya photons are produced at high densities, \left<n_{\rm H}\right>\sim 100\,\mbox{{\rm cm^{-3}}}, mostly via recombinative radiation. Once the pressure from photo-ionization pushes the dense gas away, more diffuse gas in the vicinity of young stars and the irradiated surfaces of the clumps with \left<n_{\rm H}\right>\sim 1-5\,\mbox{{\rm cm^{-3}}} become the main sites of \Lya production. The initial position of \Lya photons does not match the distribution of stellar particles precisely at the late stage of the cloud evolution, as LyC radiation that provides photoionized electrons is widespread over the cloud. To be more quantitative, we measure the optical depth to the line centre ($\tau_{0}$) for 3072 sight lines from each star particle and compare them with $\tau_{0}$ measured from the positions of the actual \Lya emitting gas to 30,000 random directions in the bottom panel of Figure 8. The plot shows that $\tau_{0}$ measured from the stars is systematically lower than $\tau_{0}$ measured from the gas, as the young stars photo-ionize the surrounding neutral hydrogen and create low density channels via stellar feedback. Note that the difference between the two measurements becomes more prominent in the latter stages of the cloud evolution, indicating that the actual line emitting gas distribution should be used to compute the profile of \Lya photons rather than the stellar distribution as a radiation source (e.g. Verhamme et al., 2012; Behrens & Braun, 2014) if the structure of the ISM is resolved.

Because collisional radiation (right panels) as well as recombination (left panels) requires electrons to produce \Lya, the emissivity maps from the two different mechanisms appear quite similar, although the former tends to better trace the dense structures and is thus less extended. Approximately 19% of the total \Lya radiation is produced via collisional radiation in the massive cloud with a 1% SFE. When the metallicity is increased by an order of magnitude (i.e., $Z_{\rm gas}=0.02$), the contribution from collisional radiation decreases to $\sim 5\%$ for a 1% SFE run. This is mainly because the enhanced metal cooling lowers the temperature of the \Lya-emitting, ionized gas from $\approx 12,700\,{\rm K}$ to $\approx 11,400\,{\rm K}$, and collisional ionization accordingly becomes less efficient.

Figure 9 shows the rate of \Lya photon production inside the cloud. When $f_{\rm esc}^{\rm LyC}$ is negligible (M6_SFE1 and M6_SFE1), the number of \Lya photons can be reasonably approximated by the simple calculation that 67% of all the LyC photons from the stellar population are converted into \Lya photons, as is often done in the literature (e.g. Verhamme et al., 2012; Dijkstra, 2014). This is possible i) because LyC photons are efficiently absorbed and ii) because the contribution from cooling radiation to the total \Lya photon budget is not very significant ($\sim 20\%$). The simple approximation ($N_{\rm Ly\alpha}\approx 0.67N_{\rm LyC}$) is no longer valid in the early phase ($t\la 1\,{\rm Myr}$) of bubble expansion during which a large number of stars shine simultaneously and the local surplus of LyC photons does not directly contribute to recombinative radiation (see M6_SFE10 for an example). In addition, if stellar feedback ejects a large amount of gas from the cloud, LyC photons would leave without interacting with ionizing neutral hydrogen, and a larger number of \Lya photons may be produced in the ISM/CGM of the galaxy compared with those produced inside the cloud.

3.3.2 Escape fractions of \Lya photons

Figure 10 shows the escape fractions of \Lya photons ($f_{\rm esc}^{\rm Lya}$) from the simulated clouds. We find that $f_{\rm esc}^{\rm Lya}$ generally increases with time as the covering fraction of dense dusty gas diminishes due to stellar feedback (Figure 2). Because this process occurs more rapidly in clouds with higher SFEs, more \Lya photons escape from the runs with 10% SFE (\left<\mbox{f_{\rm esc}^{\rm Lya}}\right>=76.7\% and $61.1\%$ for M_{\rm cloud}=10^{6}\,\mbox{\rm M_{\odot}} and 10^{5}\,\mbox{\rm M_{\odot}}) than the runs with 1% SFE (\left<\mbox{f_{\rm esc}^{\rm Lya}}\right>=47.0\% and $51.0\%$, respectively). Similar to the escape of LyC photons (Figure 4), $f_{\rm esc}^{\rm Lya}$ in the M5_SFE1 run becomes temporarily very low when the only star particle is enshrouded by a dense dusty clump (blue solid line in the bottom panel). This phase does not last long, as the clumps are destroyed by radiation feedback. High escape fractions ($\sim 40-80\%$) are commonly found in the clouds with $Z_{\rm gas}=0.002$, regardless of the choice of the SED (see Table 3), demonstrating that in metal-poor environments, most \Lya photons from young stellar populations are likely to escape from their birth clouds.

In contrast, \Lya photons are more efficiently destroyed in dusty environments. Neufeld (1990) showed that the escape of \Lya photons in a uniform medium may be described as a function of the optical depth to \Lya ($\tau_{0}$) and dust ($\tau_{d}$), as \mbox{f_{\rm esc}^{\rm Lya}}=1/\cosh\left[A\left(a\tau_{0}\right)^{1/6}\tau_{d}^{1/2}\right], where $A\approx 2$ is a fitting parameter and $a$ is the Voigt parameter. If we apply this formula to obtain \left<\mbox{f_{\rm esc}^{\rm Lya}}\right> in the metal-rich environment by replacing $\tau_{d}$ with $10\,\tau_{d}$, it should be $\approx 28\%$ and $9\%$ for the M6_SFE10 and M6_SFE1 run, respectively. However, post-processing of the metal-rich runs with rascas yields a somewhat different \left<\mbox{f_{\rm esc}^{\rm Lya}}\right> of $\approx 20\%$ both for 1% and 10% SFEs (Table 3). This is not unexpected, given that the metals are distributed inhomogeneously in the simulated domain and that the optical depth to individual photons can be different, as they are produced in relatively extended regions of the cloud.

Hayes et al. (2011) and Blanc et al. (2011) find that more dust reddened \Lya emitters exhibit lower escape fractions of \Lya photons by comparing the observed \Lya line with the intrinsic flux derived either from the UV or $H\alpha$. To examine whether the trend is established already on cloud scales, we compute the colour index of the simulated clouds by convolving the angle-averaged, dust-attenuated spectrum with the B and V band filter throughputs. Figure 11 demonstrates that the escape of \Lya photons is less efficient in more dust-reddened clouds, largely consistent with the observed trend. However, two interesting differences are found. First, when the radiation field from young stars does not permeate the cloud, dust reddening is very significant, despite that a few percent of \Lya photons still emerge from the cloud through low-density channels. As LyC photons ionize the neutral hydrogen around the young stars, dust is destroyed (by construction) while the majority of neutral gas remains intact (see Figure 20). Consequently, in this early phase ($t\la 3\,{\rm Myr}$), dust reddening decreases while $f_{\rm esc}^{\rm Lya}$ is kept nearly fixed, as is shown with smaller yellow symbols in Figure 11. Second, when dust reddening is negligible, i.e., $E(B-V)\la 0.1$, the majority of the \Lya photons escape, which is a factor of two higher than the measurement by Hayes et al. (2011). During the transparent phase, our simulated cloud is better represented by Calzetti et al. (2000), although attenuation due to the ISM/CGM can shift the sequence to the Hayes et al. (2011) line. We note that the two differences can potentially contribute to the scatter in the observed \mbox{f_{\rm esc}^{\rm Lya}}-E(B-V) relations, which is consistent with the results from the HETDEX pilot survey (Blanc et al., 2011).

3.3.3 Line profile of \Lya photons

We find that the velocity profile of \Lya photons is already complex on cloud scales. Figure 12 presents angle-averaged \Lya profiles as a function of time in the massive cloud (M_{\rm cloud}=10^{6}\,\mbox{\rm M_{\odot}}) with two different SFEs (1% and 10%). In the early phase, during which $f_{\rm esc}^{\rm LyC}$ is small and the optical depth is high, the velocity profile exhibits well-known symmetric double peaks (Neufeld, 1990; Ahn et al., 2001; Verhamme et al., 2006; Dijkstra, 2014). Once radiation feedback drives outflowing motions, the red peak becomes more pronounced (green line in the M6_SFE1 case) and eventually dominates the velocity profile (orange and pink lines in the M6_SFE1 run or green and orange lines in the M6_SFE10 run). Note that few photons with zero velocity shift escape from the cloud (c.f. Behrens et al., 2014).

We also note that the emergent spectrum is quite broad (\Delta v_{\rm peak}\sim 200-400\,\mbox{{\rm km,s^{-1}}}) in the early phase of the cloud evolution. Because the young stars are enshrouded by a large amount of neutral hydrogen, the initial average column density along the sight lines of star particles is large ($\left<N_{\rm H}\right>\sim 10^{22}\,{\rm cm^{-2}}$), and the line-centre optical depth is $\tau_{0}\sim 10^{8}-10^{9}$ (see the bottom right panel of Figure 8). However, since \Lya photons preferentially propagate along low-density channels due to their resonant nature, the velocity peaks are not separated as much as the uniform case with the given optical depth (v_{\rm sep}\equiv v_{\rm peak,red}-v_{\rm peak,blue}\sim 1000\,\mbox{{\rm km,s^{-1}}}), but this separation is certainly broader than the thermally broadened spectra by the warm ($T\sim 10^{4}\,{\rm K}$) ISM. At later stages ($t\ga 10\,{\rm Myr}$), the logarithmic mean of $\left<\tau_{0}\right>$ at the \Lya production sites becomes smaller ($\log\tau_{0}\sim 6-8$) in the case of M6_SFE1, resulting in the narrower offset of the velocity peak (\Delta v_{\rm peak}\la 100\,\mbox{{\rm km,s^{-1}}}). It is worth noting that \Lya photons scatter significantly at all times, possibly imparting a large amount of momentum to unbind the star-forming gas via resonant scattering (see Section 4.4, Kimm et al., 2018; Dijkstra & Loeb, 2008; Smith et al., 2017).

To quantify the asymmetry of the velocity profile (e.g. Erb et al., 2014), we present the ratio of the number of photons blueward of the line centre to the number of photons redward of the line centre ($=L_{\rm blue}/L_{\rm red}$) as a function of time in Figure 13 (top panel). When \Lya photons are efficiently trapped, the ratio is close to unity, although no simulated clouds exhibit a large velocity offset of v_{\rm peak,red}>300\,\mbox{{\rm km,s^{-1}}} with $L_{\rm blue}/L_{\rm red}\sim 1$ due to the presence of turbulent structures. Once the gas near the young stars is radially accelerated, the outflowing motion leads to a smaller ratio of $L_{\rm blue}/L_{\rm red}\approx 0.2-0.5$. Because the velocity of the outflows is not very large (Figure 3) and because the optical depth to the \Lya photon decreases (Figure 8), we find that the velocity peak of the spectrum redward of the line centre in the simulated clouds is no greater than 100$-$200 ${\rm km\,s^{-1}}$ at $t\ga 3\,{\rm Myr}$. This seems inconsistent with the observations that some galaxies show a large velocity offset (v_{\rm peak,red}\ga 400\,\mbox{{\rm km,s^{-1}}}) with a pronounced red peak ($L_{\rm blue}/L_{\rm red}\la 0.5$, Erb et al., 2014). However, the discrepancy may be reconciled if scattering with large-scale galactic outflows is included, as it can shift the frequency of the photons to even redward directions. It is also interesting to point out that the simulated spectra of the clouds cannot account for the small velocity offset with $L_{\rm blue}/L_{\rm red}\ga 1$, suggesting that a large fraction of \Lya photons may directly arise from an optically thin ISM with little outflowing motions in some galaxies. Not surprisingly, both cases demonstrate that the modelling of emergent \Lya spectra requires the scattering with birth clouds and the ISM/CGM (e.g. Smith et al., 2019).

Figure 12 also shows that the velocity profile of \Lya photons is highly time variable, which can potentially be an issue when modelling \Lya emission in under-resolved galactic-scale simulations. In this case, a possible option is to use the luminosity-weighted time average of the velocity profile as an input spectrum for star-forming regions and let them propagate from the cloud boundaries. This is certainly an approximation but would be a good alternative to using simple Gaussian profiles in simulations with under-resolved clouds especially if there is a large number of GMCs hosting stellar populations with different ages. To provide an idea of how the processed spectrum from the cloud would contribute to the total \Lya spectrum in a galaxy, we compute the luminosity-weighted, time average spectrum of \Lya photons in different clouds in Figure 14. Note that we assume 100% of escaping LyC photons are absorbed and re-processed into \Lya photons outside the cloud, i.e., in a volume-filling warm-ionized ISM with $T\sim 10^{4}\,{\rm K}$ (e.g. Kim et al., 2013b; Kimm et al., 2018). For simplicity, we do not take into account any existing macroscopic motions, such as outflows, which can make the spectrum broader (see the discussion section). The plot shows that when some fraction of the LyC photons leak out from the system (i.e., M6_SFE1 series), the average spectrum tends to have triple peaks. If we fit the velocity distributions with the single Gaussian profile, the FWHM would be $\approx 330$ and 410\,\mbox{{\rm km,s^{-1}}} for the binary and single stellar evolution SEDs, respectively, indicating that \Lya photons scatter more in the cloud with weaker radiation feedback (M6_SFE1_sng). When the absorption by dust increases (M6_SFE1_Zsun), photons that travel a longer distance are preferentially destroyed by dust, resulting in a narrower FWHM of \approx 220\,\mbox{{\rm km,s^{-1}}} (see also Laursen et al., 2009). Indeed, we confirm that if we artificially lower the gas metallicity by an order of magnitude ($Z_{\rm gas}=0.002$) before post-processing with rascas, we recover the large FWHM (\sim 400\,\mbox{{\rm km,s^{-1}}}) found in the M6_SFE1 run.

In contrast, if the SFE is high enough to disrupt the cloud very quickly (i.e. M6_SFE10 series), we find little dependence of the FWHM on the cloud properties. The velocity profiles of the runs with the binary star and single star SEDs look very similar essentially because both SEDs produce a similar number of LyC photons before the cloud is destroyed. An exception is the metal-rich case in which the cloud is not efficiently destroyed, as opposed to other high SFE runs. The different amplitude at the velocity centre between M6_SFE10 and M6_SFE10_sng is simply due to the fact that more ionizing photons are emitted in the binary SED than in the single SED. At the late stage of the evolution, the profiles become asymmetric and display signs of outflows (red peak), even though these weak signals are likely to be modified by the scattering due to neutral hydrogen in the ISM/CGM.

4 Discussion

In this section, we discuss how the presence of turbulent structures affects the escape of LyC and \Lya photons, and how these photons are related by comparison with observations. We also discuss the implication of the escape fractions measured from turbulent clouds to the reionization history of the Universe, and show possible effects derived from modelling of star formation and feedback.

4.1 Effects of turbulent structures

As mentioned previously, the typical resolution of galactic-scale simulations ($\Delta x\ga 10\,{\rm pc}$) is insufficient to capture complex turbulent structures in star-forming regions. In this case, low-density channels through which LyC and \Lya photons propagate may be under-resolved, leading to an underestimation in the true escape fractions. To examine the effects of turbulent structures, we run a control simulation without initial turbulence while keeping other parameters including positions, metallicity, age, and mass of the star particles fixed, as in M6_SFE1. This is done by fitting the radial gas distribution of the initial turbulent cloud to the analytic form of $n_{\rm H}=n_{\rm H,0}\left[1+\left(r/r_{s}\right)^{\alpha}\right]^{-\beta}$. We use the following set of parameters that best matches the gas distributions: n_{\rm H,0}=499.4\,\mbox{{\rm cm^{-3}}}, $r_{s}=40.1\,{\rm pc}$, $\alpha=1.343$, and $\beta=6.509$.

Figure 15 (top panel) shows that the propagation of the LyC photons is efficiently confined within the cloud in the early phase ($t\la 6\,{\rm Myr}$) in the absence of initial turbulence. Although SN explosions drive turbulence inside the cloud, the outer neutral gas shells that block LyC photons are relatively well maintained until the ionization front breaks out of the cloud. Later, the shell experiences hydrodynamic instabilities (e.g. Elmegreen, 1994), breaking into smaller clumps, similar to the final stage of the cloud evolution in the turbulent case (c.f. Figure 2). As such, the escape of LyC photons is efficiently suppressed until $\sim 6$ Myr and increases rapidly afterwards. The resulting luminosity-weighted escape fraction (\left<\mbox{f_{\rm esc}^{\rm LyC}}\right>=8.5\%) is found to be even higher than the turbulent case (\left<\mbox{f_{\rm esc}^{\rm LyC}}\right>=4.8\%), as the propagation of the ionization front is faster in the case with no turbulence where the mean density around young stellar populations drops steeply at large radii (c.f. Safarzadeh & Scannapieco, 2016). This indicates that the porous structure inside the cloud does not necessarily lead to the significantly enhanced escape fractions, as speculated by previous studies as a possible solution to reproduce an early reionization of the Universe (e.g. Kimm & Cen, 2014).

Interestingly, we find that the presence of a turbulent structure leads to more significant escape of \Lya photons, as suggested by the observations of 14 galaxies with strong \Lya emission (Herenz et al., 2016). The middle right panel of Figure 15 shows that the run without the initial turbulence begins with a high $f_{\rm esc}^{\rm Lya}$ value at $t\approx 0$, as the cloud is assumed to be initially half molecular and transparent to \Lya photons. However, molecular hydrogen is quickly dissociated by Lyman-Werner radiation from stars on a timescale of $\approx 1\,{\rm Myr}$, and the scattering of \Lya photons becomes significant again. Because there are no low-density channels due to the lack of a turbulent structure, \Lya photons are efficiently trapped inside the HI shells, resulting in larger velocity offsets compared to those in the turbulent clouds (see blue lines in the bottom panel). Furthermore, $f_{\rm esc}^{\rm Lya}$ is reduced as \Lya photons have a high probability of encountering dust and being destroyed. This phase stops and the high escape fractions observed in the turbulent case are recovered once the gas cloud fragments into smaller clumps ($t>5\,{\rm Myr}$). Because the majority of the \Lya photons are created in the early phase, the resulting net \Lya escape fractions become lower in the absence of turbulence (\left<\mbox{f_{\rm esc}^{\rm Lya}}\right>\approx 31\% vs. $47\%$).

4.2 Connection between LyC and \Lya photons

In Figure 16, we compare the escape fractions of ${\rm LyC}$ around 900Åand \Lya photons from different environments. We confirm the previous finding that \Lya photons leave the cloud more efficiently than ionizing photons (Yajima et al., 2014; Verhamme et al., 2015; Dijkstra et al., 2016; Verhamme et al., 2017). A large fraction ($\sim 20-80\%$) of \Lya photons leave the cloud regardless of the input stellar spectrum and size of the cloud, even when almost all ($f_{\rm esc}^{\rm 900}\la 10^{-3}$) of the ${\rm LyC}$ photons are absorbed by neutral hydrogen. Dijkstra et al. (2016) show that the escape fractions of \Lya photons tend to be low ($\la 10\%$) if there are more than $\sim 5$ clumps along the sight line from the centre to the edge of the cloud. In light of this, the high $f_{\rm esc}^{\rm Lya}$ observed in our turbulent cloud simulations suggest that \Lya photons do not interact with dense clumps and filaments very often but rather propagate through low-density channels. We also find that $f_{\rm esc}^{\rm Lya}$ decreases down to $\sim 10\%$ when the amount of dust increases (i.e., $Z_{\rm gas}=0.02$), although the difference is negligible once the cloud becomes optically thin to both \Lya and ${\rm LyC}$ photons. As a result, the positive correlation between $f_{\rm esc}^{\rm 900}$ and $f_{\rm esc}^{\rm Lya}$ becomes more pronounced in the metal-rich case.

Verhamme et al. (2015) argue that ${\rm LyC}$ leaking candidates may be pre-selected by one of the following two features. First, if the intervening medium is optically thin to ${\rm LyC}$ photons, the separation of double peaks in the \Lya profile should be smaller than v_{\rm sep}\la 300\,\mbox{{\rm km,s^{-1}}}. Second, if the ISM is clumpy, \Lya photons propagating through low-density channels will be seen at line centre, while some photons will escape on the blue side of the profile. The latter point is straightforward to understand, and indeed can be inferred from the top panel of Figure 14. Although the precise determination of the amplitude of the central peak would depend on the distribution of the ISM, the clouds with \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>\ga 2\% exhibit some photons near the line centre (e.g., M6_SFE1 and M6_SFE1_Zsun), while they nearly disappear in the cloud with a lower \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>=0.3\% (M6_SFE1_sng).

To explore the first possibility, we measure the peak separation ($v_{\rm sep}$) by fitting the \Lya line profile with two different components, i.e., two skewed Gaussian profiles on the left and right side of the line centre in the right panel of Figure 16. Our simulated velocity separation of the two peaks generally follows the simple analytic trends expected in a uniform medium (e.g. Neufeld, 1990; Dijkstra, 2017),

[TABLE]

but two differences are found. First, the predicted $v_{\rm sep}$ appears to be larger than the analytic calculation with $T\sim 10^{4}\,{\rm K}$, appropriate for our simulated clouds, which is due to the fact that LyC photons can escape through low-density channels, while the majority of the \Lya photons still scatter inside dense regions of the cloud. Second, $v_{\rm sep}$ does not approach zero even when the escape fraction becomes very high ($f_{\rm esc}^{\rm 900}\ga 50\%$), as most \Lya photons are generated near (a few surviving) relatively dense regions where the optical depth to \Lya is still high. More importantly, we find that when $f_{\rm esc}^{\rm 900}$ is intermediate or high ($\ga 1\%$), the separation of the peaks is small (v_{\rm sep}\la 200\,\mbox{{\rm km,s^{-1}}}), supporting the criterion proposed by Verhamme et al. (2015). However, the small separation does not guarantee pre-selection of the efficient LyC leakers on cloud scales (as pointed out by Verhamme et al., 2015) because the turbulent structure inside the cloud allows the \Lya photons to propagate more efficiently than the uniform or shell case. As a result, some clouds with low LyC escape fractions ($f_{\rm esc}^{\rm 900}\la 0.1\%$) show v_{\rm sep}\la 300\,\mbox{{\rm km,s^{-1}}}.

Figure 16 also shows that $f_{\rm esc}^{\rm Lya}$ from the simulated clouds has a weaker correlation with $f_{\rm esc}^{\rm 900}$ compared to that obtained from metal-poor ($\sim 0.1$–$0.2\,Z_{\odot}$) compact starburst galaxies (Vanzella et al., 2015; Izotov et al., 2016b, a; Verhamme et al., 2017; Vanzella et al., 2018; Izotov et al., 2018a). The two significant LyC emitters, Ion2 (Vanzella et al., 2015) and J1154+2443 (Izotov et al., 2018a), exhibit escape of the order of unity for both \Lya and LyC photons, similar to our findings. However, other observed galaxies with lower $f_{\rm esc}^{\rm 900}$ tend to have a smaller $f_{\rm esc}^{\rm Lya}$, but with a much more broadened \Lya spectrum (the right panel, Verhamme et al., 2017; Izotov et al., 2018a; Izotov et al., 2018b, see also the equation 2 of Izotov et al. 2018b). This suggests that the interaction of \Lya photons with neutral hydrogen and dust inside the star-forming clouds does not fully account for the observed features of the UV spectrum and that scattering with the ISM may significantly change the propagation of \Lya photons in galaxies. It may also be possible that multiple clouds in different stages of cloud evolution contribute in a complex way to reproduce the observations, but in this case a large degree of scatter may be present in the low $f_{\rm esc}^{\rm Lya}$ regime.

Recently, Vanzella et al. (2018); Rivera-Thorsen et al. (2017); Izotov et al. (2018a) observed LyC leakers with triple \Lya peaks. They claim that this is consistent with a simple model in which LyC and \Lya photons escape through the same cavity (Behrens et al., 2014; Verhamme et al., 2015). Our simulations also support this picture in the sense that the central peaks do not appear to arise from inside of the clouds, but only seem possible if we include \Lya photons that would be generated by LyC photons leaking into an ISM (Figure 14). An interesting difference from the results of Ion3 (Vanzella et al., 2018) is that the central peak in Ion3 extends to v\sim\pm 100\,\mbox{{\rm km,s^{-1}}}, whereas our thermally and turbulently broadened spectra do not reach more than v\sim\pm 50\,\mbox{{\rm km,s^{-1}}}. This indicates that the motions of \Lya emitting gas in the ISM from this compact starburst galaxy are likely to be dominated by strong turbulence corresponding to the Doppler parameter of b\ga 20\,\mbox{{\rm km,s^{-1}}} (see Figure 1 of Verhamme et al. 2015).

4.3 Implications for reionization of the Universe

In a dwarf galaxy-driven scenario for reionization, the key quantity that governs the expansion of ionized bubbles is the number of escaping photons ($N_{\rm esc}^{\rm LyC}$, e.g., Wise et al., 2014; Kimm et al., 2017; Koh & Wise, 2018; Rosdahl et al., 2018). However, making theoretical predictions for $N_{\rm esc}^{\rm LyC}$ is not trivial, not only because the turbulent structure for high-$z$ galaxies is unknown, but also because the input stellar spectra are not well constrained in the early Universe. We find that in the massive cloud with a 10% SFE, the runs adopting a single, binary, and binary SED with larger cut-off mass (i.e., M_{\rm upper}=300\,\mbox{\rm M_{\odot}}) yield total $0.6\times 10^{52}$, $1.3\times 10^{52}$, and $2.4\times 10^{52}$ number of escaping LyC photons during the first 10 Myr (Table 3). In the case of a lower SFE (1%), the total $N_{\rm esc}^{\rm LyC}$ within $20\,{\rm Myr}$ is found to be $0.08\times 10^{50}$, $1.6\times 10^{50}$, and $3.5\times 10^{50}$ for the single, binary, and binary SED with the larger cut-off mass, respectively. These results indicate that the total $N_{\rm esc}^{\rm LyC}$ can vary by a factor of $\sim 4-40$ depending on the choice of the stellar spectra, and that one must be careful with interpreting the results from previous simulations where the single stellar evolution is adopted for the photon production rate. Indeed, using high-resolution ($\sim$10 pc), cosmological radiation-hydrodynamic simulations, Rosdahl et al. (2018) show that the simulated volume of $(10\,{\rm Mpc})^{3}$ is fully ionized by $z_{\rm reion}\approx 7$ when the binary stellar evolution model is used, while single stars cannot ionize the volume in time ($z_{\rm reion}<6$) (see also Ma et al., 2016). Given that a large fraction of massive stars do live in binaries in the local environments (Sana et al., 2012), it is encouraging that the binary model manages to ionize the simulated universe early enough. But the question remains how efficiently star formation needs to be suppressed at high redshifts in order to precisely match the end of reionization (i.e. $z_{\rm reion}\approx 6$), as the maximum stellar mass of low-metallicity stars is not well constrained.

4.4 Impact of star formation and feedback schemes

In our simulations, we place star particles randomly inside the clouds instead of directly modelling star formation (e.g. Bate et al., 1995; Gong & Ostriker, 2013; Hubber et al., 2013; Bleuler & Teyssier, 2014) so that we can control the SFE per cloud while resolving the initial Stromgren sphere. Nevertheless, it is true that stars form in dense environments, and our experiments may under-estimate the interaction of LyC and \Lya photons with neutral hydrogen patches.

In order to understand the difference of adopting a more realistic star formation model, we run additional simulations by placing star particles preferentially in dense pockets of gas. This is done by computing the SFE per free-fall time based on the local thermo-turbulent conditions (Federrath & Klessen, 2012), as described in Kimm et al. (2017, Equation 2). We then assign star particles based on the probability of forming a star particle in each cell. The resulting average density of the host cell is \left<\mbox{n_{\rm H}}\right>\sim 2\times 10^{5}\,\mbox{{\rm cm^{-3}}} (M6_SFE1_dSF), which is $10^{4}$ times denser than that of the runs with a random sampling (i.e., M6_SFE1). Because young stars are deeply embedded in dense, high-pressure regions, LyC photons are efficiently absorbed and ionization fronts stall at smaller radii. Although some star particles destroy their local clump through a combination of radiation and supernova explosions, the ionization fronts tend to propagate slowly, as the HII bubble develops from over-dense regions. As a result, only $\approx 1\%$ of LyC photons escape from the cloud in the M6_SFE1_dSF run (Figure 17), which is a factor of $\sim 5$ lower than that of the fiducial model (M6_SFE1), demonstrating the importance of self-consistent modelling of star formation in simulations.

However, we find that the previous results rely sensitively on the presence of strong feedback processes, such as \Lya pressure. To examine the possible impact of \Lya pressure, we adopt the simple model developed in Kimm et al. (2018). Briefly, we compute the number of \Lya photons produced in each cell using Equations 4–8 and estimate the total neutral hydrogen column density by summing the neutral hydrogen from the host and the adjacent cell along the direction of the propagation. The direction of the momentum input is taken as the direction of the LyC flux, as it allows us to trace the position of the ionizing source. Note that this method neglects the long-range force due to \Lya that escapes from the cell of interest and propagates to the neighbouring cells (see Smith et al., 2019, for example). Also, we do not use the Sobolev approximation to conservatively estimate \Lya pressure. In this regard, the actual impact from \Lya could be even more significant. Despite these simplifications, we find that the average density of the host cell of the young star particles is reduced to \left<\mbox{n_{\rm H}}\right>\sim 0.1-1\,\mbox{{\rm cm^{-3}}} in about 0.1 Myr. Note that this is four orders of magnitude lower than that of the M6_SFE1_dSF run (\left<\mbox{n_{\rm H}}\right>\sim 4000\,\mbox{{\rm cm^{-3}}}). Consequently, a large fraction of LyC photons leave from the cloud even from $t\sim 2\,{\rm Myr}$, leading to \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>\approx 13\% (Figure 17, orange lines). This suggests that predictions from the simulations without strong radiation pressure should be taken with caution and that even our results based on the randomly distributed star particles without \Lya pressure are likely to be a lower limit of the true \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>. As a first step, we do not include \Lya pressure in our fiducial set of runs, but this issue needs to be addressed for a variety of conditions in the near future.

5 Conclusions

Motivated by the recent results that a large fraction of LyC photons are absorbed on small scales (Kim et al., 2013a; Kimm & Cen, 2014; Paardekooper et al., 2015; Trebitsch et al., 2017) and to give useful insights into how we should interpret the emergent \Lya spectra, we perform a suite of RHD simulations of turbulent star-forming clouds with stellar feedback processes, including direct radiation pressure, photo-ionization heating, photo-electric heating on dust, non-thermal pressure due to multiply scattering infrared photons, and Type II supernova explosions. By randomly placing star particles around the central region of the cloud, we follow the evolution of $f_{\rm esc}^{\rm LyC}$, $f_{\rm esc}^{\rm Lya}$, and the spectral shape of the \Lya photons from the runs adopting different cloud masses, SED shapes, gas metallicity, and turbulence. Our findings can be summarized as follows.

• (i)

The escape fractions of ${\rm LyC}$ tend to increase rapidly and rather monotonically over the cloud lifetime (Figure 4). Although the porous structures inside the turbulent cloud allow for LyC photons to propagate locally, the optically thin channels are not necessarily well aligned, and all of the LyC photons are absorbed by neighbouring gas in the early stage of the evolution (Figure 2). Once radiation feedback clears away the neighbouring neutral regions and blows out the dense clumps, the HII bubble expands and the covering fraction of optically thick regions becomes smaller, elevating the escape fractions. In the case of the runs with less efficient radiation feedback (i.e., M6_SFE1_sng), SN explosions help to disrupt the cloud, although the escape fractions are not as high as in the fiducial runs (Figure 6).

• (ii)

We find that the luminosity-weighted, time-averaged escape fractions of LyC photons (\left<\mbox{f_{\rm esc}^{\rm LyC}}\right>) are relatively low from a massive cloud (M_{\rm cloud}=10^{6}\,\mbox{\rm M_{\odot}}) with a 1% SFE, which is the typical value derived in the local GMCs (e.g. Heyer et al., 2009) or simulations (e.g. Grisdale et al., 2018). With binary star SEDs, the metal-poor clouds with $Z_{\rm gas}=0.002$ show \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>\sim 5\%. For the metal-rich run ($Z_{\rm gas}=0.02$), radiative cooling enhances the recombination and significantly reduces \left<\mbox{f_{\rm esc}^{\rm LyC}}\right> to $\sim 1\%$ (Figures 4 and 6). In contrast, when the SFE is higher or when the cloud mass is smaller, overpressure due to photo-ionization heating and SN explosions blow away the cloud more rapidly (Figure 2), leading to a very large \left<\mbox{f_{\rm esc}^{\rm LyC}}\right> of $30$–$70\%$.

• (iii)

The runs with binary star SEDs or a higher stellar mass upper limit in the IMF show a higher \left<\mbox{f_{\rm esc}^{\rm LyC}}\right>, as radiation feedback is enhanced due to the larger number of LyC photons produced via binary interactions or by very massive stars (100<M/\mbox{\rm M_{\odot}}<300) compared to the SED with single star evolution (Figure 6). As a result, the number of escaping LyC photons can easily be different by a factor of $\sim$4 depending on the choice of the SED (Table 3).

• (iv)

The majority of the absorption takes place on small scales especially when the SFE is low (1%). We find that 50%, 90%, and 99% of the LyC photons are absorbed within a distance of 8, 33, and 83 pc from each star particles, respectively. The scale becomes larger (49, 140, and 191 pc) if the SFE is higher (10%) as the cloud is dispersed due to stellar feedback.

• (v)

In the run with a $1\%$ SFE, most of the \Lya photons are generated via the recombinative process, while collisional recombination contributes to $\sim 20$–$30\%$ of the total \Lya (Figure 8). The latter fraction becomes lower ($\sim 10\%$) in the metal-rich cloud, as the \Lya-emitting gas becomes cooler. The resulting number of \Lya photons produced inside the cloud reasonably matches the simple estimate based on the assumption that 67% of the LyC photons available from young stars produce \Lya photons, provided that \left<\mbox{f_{\rm esc}^{\rm LyC}}\right> is low.

• (vi)

Our simple experiment without initial turbulence shows that fewer \Lya escape from the cloud until it becomes disrupted, while the majority ($\sim 40$–$80\%$) of the \Lya photons escape from the cloud with turbulent structures in the early phase of the evolution (Figure 15). Even clouds with a large amount of dust (i.e. $Z_{\rm gas}=0.02$) show a slightly lower $f_{\rm esc}^{\rm Lya}$ of $\approx 20\%$, which is systemically larger than $f_{\rm esc}^{\rm LyC}$ (Figure 16). This suggests that the low $f_{\rm esc}^{\rm Lya}$ observed in local and high redshift galaxies (Deharveng et al., 2008; Cowie et al., 2010; Ono et al., 2010; Hayes et al., 2010) may be largely due to the substantial absorption occurring in the ISM.

• (vii)

We find that emergent \Lya spectra can be broad even on cloud scales. When ionizing radiation is effectively confined (i.e., \mbox{f_{\rm esc}^{\rm LyC}}\approx 0), the velocity spectrum from the 10^{6}\,\mbox{\rm M_{\odot}} cloud shows symmetric double peaks separated by \sim 400\,\mbox{{\rm km,s^{-1}}} (Figure 12). The peak separation becomes smaller ($v_{\rm sep}\sim 100$–300\,\mbox{{\rm km,s^{-1}}}) if \mbox{f_{\rm esc}^{\rm LyC}}\approx 0.1-1\%, but does not become closer to \Delta v\la 50\,\mbox{{\rm km,s^{-1}}} even when $f_{\rm esc}^{\rm LyC}$ becomes very high ($\ga 30\%$), as there exists residual neutral hydrogen in the relatively dense regions around which \Lya photons are produced (Figure 16). Consequently, the luminosity-weighted \Lya profiles over the cloud lifetime are found to be more complex than the simple Gaussian profile that is often used as an input spectrum for \Lya photons in the galactic scales simulation (Figure 14).

• (viii)

Finally, LyC leaking clouds ( \mbox{f_{\rm esc}^{\rm LyC}}\ga 1\%) show the separation of peaks less than v_{\rm sep}\approx 150\,\mbox{{\rm km,s^{-1}}}, consistent with Verhamme et al. (2015), as the turbulent structure allows \Lya photons to escape more efficiently (Figure 16). However, we find that the predicted $v_{\rm sep}$ for a given $f_{\rm esc}^{\rm LyC}$ is a factor of two smaller than the observed in compact metal-poor systems (e.g. Vanzella et al., 2018; Izotov et al., 2018a; Verhamme et al., 2017), again suggesting that the interaction of \Lya photons with the ISM is likely to be crucial to determine the emergent spectrum.

We note that the number of simulations performed in this study is limited, hence it may be difficult to generalize our results to the GMCs under a wide variety of conditions (e.g. Heyer et al., 2009). However, these experiments clearly demonstrate that the escape fractions of LyC photons are driven by radiation and SN feedback to steadily establish the low-density channels, although this process is highly dependent on the small-scale physics, such as star formation and input SEDs. Our work also emphasizes the complexity in predicting the \Lya line profiles, necessitating a more comprehensive understanding of the dynamics in the star-forming clouds and the ISM. Future simulations that can resolve these processes in galactic scales will be the natural step forward to make firm predictions on the escape of LyC and \Lya photons in the high-redshift Universe.

Appendix A Analytic calculations of the escape of ionizing radiation

In this section, we present a simple model for the propagation of an ionization front in a spherically symmetric cloud and compute the timescale for which the ionization front reaches the edge of the GMC. One may then calculate the photon number-weighted average escape fractions, assuming that all of the ionizing radiation leaves the cloud once the ionization front reaches the edge. Note that this experiment is highly idealized but provides useful insights into understanding the dependence of the escape fraction on the basic properties of simulated clouds.

Geen et al. (2015b, Appendix A) derive the propagation of the ionization front (see also Raga et al., 2012) in a cloud with a power-law density profile,

[TABLE]

as

[TABLE]

where $r$ is the radius, $t$ is the time, $c_{s,i}$ is the sound speed of the ionized medium, $v_{\rm ext}$ is the infall velocity of the ambient gas, $c_{\rm ext}$ represents the velocity term due to the thermal and turbulent pressure, and

[TABLE]

Here we assume that $r_{i}(t=0)=r_{\rm Strom}$ and that the density profile is quasi-static. This neglects the timescale necessary to develop the Stromgren sphere, i.e., a roughly recombination timescale ($\tau<{\rm 0.01\,Myr}$), but this is usually negligible compared to the lifetime of a GMC ($\sim\rm 1-10\,Myr$). The density of the ionized region is related to the source as

[TABLE]

The ionization front stalls at the radius where $\dot{r}_{\rm i}=0$, and for a pressure-dominated region without infalling motion, one finds

[TABLE]

This indicates that in a dynamically cold medium ($c_{\rm ext}<c_{\rm s,i}$), the over-pressure due to photo-ionization heating is significant only if the cloud is not too compact ($w<3/2$). If this condition is met, the ionized bubble would expand with time and can reach the edge of the cloud even though it may take a long time.

For the uniform profile ($w=0$), Equation 11 can be written as

[TABLE]

where $y\equiv r_{\rm i}(t)/r_{\rm Strom}$. The time required for the ionization front to propagate to some radius $r$ in a pressure-dominated region (i.e., $v_{\rm ext}\sim 0$) can then be computed as

[TABLE]

Note that $f_{\rm esc}^{\rm LyC}$ is zero if $r_{\rm stall}<r_{\rm cloud}$.

The radius of a cloud is taken from the empirical relation between the mass and radius (Roman-Duval et al., 2010), as

[TABLE]

Figure 18 shows that the ionization front increases steadily, as the LyC photon production rates are nearly constant until $5\,{\rm Myr}$. The ionization front can in principle shrink back if the external pressure is strong enough to counterbalance the pressure due to photo-ionization heating. This occurs mostly after $\sim$ 10 Myr when the emissivity drops rapidly, which is well after the ionization propagates to the outer region in this simple setup. Although massive clouds are more compact and that recombination is more efficient, the resulting escape fractions are higher because they are smaller in size.

Appendix B Resolution test and temperature distributions of the turbulent clouds

Figure 19 shows the escape fractions of LyC photons in our fiducial run with different resolutions. Although the Stromgren sphere around individual star particle is initially well resolved, the propagation of the LyC photons are affected by the maximum AMR resolution, as star particles encounter dense clumps and become enshrouded by neutral hydrogen before stellar feedback entirely destroys the cloud. As a result, 4.8, 5.1, and 7.3% of the total LyC photons escape from the cloud with 0.25, 0.5, and 1.0 pc resolution. However, the general trends in the evolution of the clouds are very similar, indicating that our conclusions are little affected by the resolution.

In Figure 20, we show the projected temperature distributions of turbulent clouds with different SED, metallicity, and input physics.

Acknowledgements

We thank Julien Devriendt, Kwang-Il Seon, and an anonymous referee for helpful comments. This work was supported by the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2017-C3-0073) and in part by the ERC Advanced Grant 320596 “The Emergence of Structure during the Epoch of Reionization". This study was performed under the umbrella of the joint collaboration between Yonsei University Observatory and the Korean Astronomy and Space Science Institute. TK was supported in part by the Yonsei University Future-leading Research Initiative (RMS2-2018-22-0183) and in part by the National Research Foundation of Korea (No. 2017R1A5A1070354 and No. 2018036146). HK thanks Brasenose College and the support of the Nicholas Kurti Junior Fellowship as well as the Beecroft Fellowship. JR and JB acknowledge support from the ORAGE project from the Agence Nationale de la Recherche under grand ANR-14-CE33-0016-03. TG is grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon for its financial support within the programme ‘Investissements d’Avenir’ (ANR-11-IDEX-0007) of the French government operated by the National Research Agency (ANR). TG acknowledges support from the European Research Council under grant agreement ERC-stg-757258 (TRIPLE). This work made use of the following DiRAC facilities (www.dirac.ac.uk): the Data Analytic system at the University of Cambridge (funded by BIS National E-infrastructure capital grant ST/K001590/1, STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1) and the Complexity system at the University of Leicester (funded by BIS National E-infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1). DiRAC is part of the National E-Infrastructure.