Ly{\alpha} emission from galaxies in the Epoch of Reionization

TL;DR

This study uses advanced simulations of a high-redshift galaxy to explore how dust, gas, and galaxy structure influence Lyα emission, revealing complex attenuation effects and their implications for galaxy evolution during reionization.

Contribution

It combines hydrodynamical and radiative transfer simulations to analyze Lyα escape in a prototypical high-redshift galaxy, providing new insights into dust effects and emission variability.

Findings

Lyα escape fractions are low (<1%) with large scatter.

Dust in clumps significantly attenuates Lyα emission.

Anti-correlation between Lyα-[CII] velocity shift and Lyα luminosity.

Abstract

The intrinsic strength of the Ly$\alpha$ line in young, star-forming systems makes it a special tool for studying high-redshift galaxies. However, interpreting observations remains challenging due to the complex radiative transfer involved. Here, we combine state-of-the-art hydrodynamical simulations of 'Althaea', a prototypical Lyman Break Galaxy (LBG, stellar mass $M_{\star}$$\simeq$ $10^{10}{\rm M}_{\odot})$ at $z=7.2$, with detailed radiative transfer computations of dust/continuum, [CII] 158 $\mu$m, and Ly$\alpha$ to clarify the relation between the galaxy properties and its Ly$\alpha$ emission. Althaea exhibits low ($f_\alpha< 1\%$) Ly$\alpha$ escape fractions and Equivalent Widths, EW $\lesssim 6$ Angstrom for the simulated lines of sight, with a large scatter. The correlation between escape fraction and inclination is weak, as a result of the rather chaotic structure of high-redshift galaxies. Low $f_\alpha$ values persist even if we artificially remove neutral gas around star forming regions to mimick the presence of HII regions. The high attenuation is primarily caused by dust clumps co-located with young stellar clusters. We can turn Althaea into a Lyman Alpha Emitter (LAE) only if we artificially remove dust from the clumps, yielding EWs up to $22$ Angstrom. Our study suggests that the LBG-LAE duty-cycle required by recent clustering measurements poses the challenging problem of a dynamically changing dust attenuation. Finally, we find an anti-correlation between the magnitude of Ly$\alpha$-[CII] line velocity shift and Ly$\alpha$ luminosity.

Equations

$$L_{\alpha }=0.68L_{\mathrm{ion}}(1-f_{esc}).$$

$$R_S^3=\frac{3Q}{4\pi n^2{\alpha }_B}$$

$$P=p(\theta )\exp (-{\tau }_o)$$

Full text

††pubyear: 2019††pagerange: Ly$\alpha$ emission from galaxies in the Epoch of Reionization–A.4

1 Introduction

Understanding the formation and evolution of young galaxies via the Ly$\alpha$ line produced by massive stars has become a standard tool of extragalactic observations. However, such a strategy is still affected by a number of uncertainties that can be ultimately associated with the resonant nature of this line, and that result in a complex relation between gas distribution, kinematics, and the observables.

This complexity fostered a large number of attempts to simulate radiative transfer (RT) in galaxies, using both simplified - but yet successful - models (Verhamme et al., 2006; Dijkstra et al., 2006; Behrens et al., 2014) and more detailed models that follow the RT in an isolated disk galaxy (Verhamme et al., 2012; Behrens & Braun, 2014), or even galaxies in their cosmological environments (Laursen, 2010; Gronke & Bird, 2017).

While considerable progress has been made, many puzzles remain. An outstanding one concerns the decreasing abundance of Ly$\alpha$ emitting galaxies towards high-redshift and in the Epoch of Reionization (EoR). The luminosity function of Lyman-$\alpha$ emitters (LAEs) shows only mild evolution between redshift 2 and 6, but rapidly declines above $z=6$ (Pentericci et al., 2011; Stark et al., 2010; Hayes et al., 2011; Schenker et al., 2012; Ono et al., 2012). In particular, the LAEs abundance drops more strongly than the UV luminosity function.

Several proposals have been made to explain this discrepancy. One possibility is related to the evolution of the neutral hydrogen fraction in the intergalactic medium (IGM), which might rapidly increase at $z>6$ as a result of the cosmic reionization process. The problem with this explanation is that it requires a dramatic change of the IGM neutral fraction $x_{HI}$ of $\Delta x_{HI}\sim 0.5$ in a short redshift range, $\Delta z=1$ (Pentericci et al., 2011; Ono et al., 2012; Schenker et al., 2012; Dijkstra et al., 2011), requiring both a late and very rapid reionization. While the latest Planck results (Planck Collaboration et al., 2018) do allow for a late reionization scenario, the rapidness of the transition is hard to model. However, e.g. Mesinger et al. (2015) find that a joint evolution of the ionization field and the neutral fraction is moderately consistent with the data. Other possible explanations of the vanishing of LAEs focus on changes in the circumgalactic medium (CGM). Bolton & Haehnelt (2013); Weinberger et al. (2018) argue that in the final stages of reionization, a decrease in the mean free path length of ionizing photons lowers the transmissivity of the CGM to Ly$\alpha$ photons; Sadoun et al. (2017) invoke an explanation related to the ionization state of the infalling gas around LAEs .

While this possibility relates the LAEs dismissal at high redshifts to the cosmological evolution of the Universe, a different idea focuses on a parallel evolution of the interstellar medium (ISM) properties. Some encouraging hints come from observations of the line shift of the Ly$\alpha$ line deduced from other non-resonant lines, like e.g. [C $\scriptstyle\rm II$]. Typically, the observed shift is $\sim 400$ km/s for Lyman break galaxies (LBGs) at intermediate redshift (e.g. at $z=3$ Erb et al., 2014), while Pentericci et al. (2016) find a typical shift of the Ly$\alpha$ line of only 100-200 km/s at redshift $z\approx 7$ (see also Carniani et al., 2017, for an analysis of currently available line shifts at high redshift).

Thus, a change in the ISM properties causing a reduction of the Ly$\alpha$ line velocity shift might represent a viable solution to this puzzle, as the reduced line shift would render the Ly$\alpha$ radiation more susceptible to attenuation by the intergalactic medium (IGM). If this is the correct explanation, one expects a positive correlation between the velocity shift of the Ly$\alpha$ line and its escape fraction.

Moreover, the connection between LAEs and LBG has also remained elusive (Dayal & Ferrara, 2012). Clustering analysis at intermediate redshifts suggests an overlap of the populations, taking the form of an effective duty-cycle, with LBG for some fraction of time attaining a much larger escape fraction and turning into a LAE (Kovač et al., 2007). Such a duty-cycle scenario is typically quantified in terms of the fraction of dark matter halos hosting a LAE. For example, from the SILVERRUSH survey data, Ouchi et al. (2018) infer duty cycles of $<1\%$ at $z\sim 6$ using halo occupation models, while Sobacchi & Mesinger (2015) estimate it to be less than few percent from combining observational data with modelling of the EoR. However, the mechanism that governs this transition of a galaxy from Ly$\alpha$-dark to Ly$\alpha$-bright (or vice versa) remains unclear.

Theoretically, Pallottini et al. (2017b) performed high-resolution zoom simulations of a prototypical high-redshift LBG called “Althæa”. By $z=6$, Althæa has a stellar mass of $M_{\star}\simeq 10^{10}{\rm M}_{\odot}$ and a $\rm SFR\simeq 100\,{\rm M}_{\odot}/{\rm yr}$. It shows an exponentially rising SFR and a SFR-stellar mass relation compatible with that derived from high-$z$ observations (e.g. Jiang et al., 2016). It also closely follows the Schmidt-Kennicutt relation (Kennicutt & Evans, 2012; Krumholz et al., 2012). In Behrens et al. (2018a) we studied the dust Far Infra-Red (FIR) emission from Althæa. The resulting high dust temperatures we found might solve the puzzling low infrared excess recently deduced for the dusty galaxy A2744_YD4 observed by Laporte et al. (2017) (also see Katz et al., 2019).

In this paper, we study the Ly$\alpha$ line transfer in Althæa, focusing on observables like the Ly$\alpha$ escape fraction, line shift, and spectra along different lines of sight. We augment the original simulation with detailed RT calculations of the Ly$\alpha$ radiation emitted by young stars, complementing the analysis of the UV and FIR continuum (Behrens et al., 2018a), and of the FIR emission lines (Pallottini et al., 2017b; Vallini et al., 2018). By combining these observables, we aim at understanding the relative importance of geometry, dust content, velocity field and viewing angle in determining the fraction of escaping radiation.

2 Methods

First we summarize the main characteristics of the adopted hydrodynamical simulations (Sec. 2.1), then we detail the assumption for the stellar emission model, dust absorption, and treatment of hydrogen ionization (Sec. 2.2), and finally we describe the RT simulations (Sec. 2.3).

2.1 Hydrodynamical simulations

We use the simulations described in Pallottini et al. (2017b). The simulations are based on a modified version of the publicly available ramses code (Teyssier, 2002). The simulations evolve a comoving cosmological volume of $(20\,{\rm Mpc/h})^{3}$, focusing on a halo of mass $10^{11}{\rm M}_{\odot}$ hosting the galaxy Althæa in a zoom-in fashion. The gas mass resolution in the zoomed region is $2\times 10^{4}{\rm M}_{\odot}$ and the adaptive mesh refinement (AMR) allows us to reach a resolution of $\simeq 25\,\rm pc$ by $z=7$. Outside of the zoom-in regions, our grid has a spatial resolution of 14 kpc.

In the simulation stars form from molecular hydrogen, whose abundance is computed from non-equilibrium chemistry using the krome package111https://bitbucket.org/tgrassi/krome (Grassi et al., 2014; Bovino et al., 2016). Stellar feedback is modeled as described in Pallottini et al. (2017a) and includes supernova explosions, winds from massive stars, radiation pressure and accounts for the sub-grid evolution of the blastwave within molecular clouds. In the simulation stellar clusters are assumed to have a Kroupa (2001) Initial Mass Function (IMF).

Most of the results presented in this work refer to the snapshot at $z=7.2$, i.e. the same one used in the fiducial model in Behrens et al. (2018a)222Note, that in Behrens et al. (2018a) the galaxy was shifted to a higher redshift for better comparison with observations from Laporte et al. (2017).. At this redshift Althæa has $\mathrm{SFR}\simeq 78\,{\rm M}_{\odot}/{\rm yr}$, $M_{\star}\simeq 10^{10}{\rm M}_{\odot}$, and an age of about $513\,{\rm Myr}$.

2.2 Emission and absorption properties

2.2.1 Stellar emissivity

In Althæa we keep track of age ($t_{\star}$) and metallicity ($Z_{\star}$) of all its stellar clusters. Therefore, we can use population synthesis codes to derive the Spectral Energy Distribution (SED) of each stellar cluster. In particular, we use Bruzual & Charlot (2003) prescribing the same Kroupa (2001) IMF assumed for the feedback in the hydrodynamical simulations.

These assumptions yield the intrinsic continuum flux directly. Integrating the ionizing part of each SED, we obtain the total ionizing power $L_{\mathrm{ion}}$. Assuming Case-B recombination and a gas temperature of 104 K, the latter can be related to the intrinsic Ly$\alpha$ radiation via (Dijkstra, 2014)

[TABLE]

We will assume in the following that the escape fraction of ionizing photons, $f_{esc}=0$. This maximizes the Ly$\alpha$ output. We note that the factor of 0.68 is only strictly valid at a temperature of $10^{4}$ K. However, as the temperature dependence is very weak (at 103 K it is 0.77), we neglect it.

Practically, we use a table and interpolate age and metallicity to obtain the rate of ionizing photons. In Fig. 1, we show this rate for different $Z_{\star}$ as a function of $t_{\star}$. We assume that the stars in a cluster form in a single burst, and treat the stellar clusters as point emitters of Ly$\alpha$ radiation.

We do not consider any additional source of Ly$\alpha$ radiation. In particular, we do not include cooling radiation from hydrogen around temperatures of $\sim 10^{4}$ K. The uncertainties in estimating the intrinsic luminosity of Ly$\alpha$ in cooling radiation are large (e.g., see Goerdt et al., 2010; Rosdahl & Blaizot, 2012; Faucher-Giguère et al., 2010). As an upper limit, we can calculate the intrinsic Ly$\alpha$ luminosity assuming that all the gas with temperature $10^{4}<T/{\rm K}<5\times 10^{4}$ cools only due to Ly$\alpha$ radiation. In this case, we find that the contribution from cooling radiation is 1% of the intrinsic luminosity coming from stellar ionizing flux. Hence we have decided to neglect this contribution.

2.2.2 Dust model

The evolution of the gas and stellar metallicity is traced in Althæa. However, we stress that the simulation does not keep track of dust formation and/or destruction processes, and therefore does not follow the evolution of dust content/properties directly. Therefore, we need to make an assumption for the conversion between metallicity and dust content. We take the same route as in Behrens et al. (2018a) by assuming a linear relation between metal and dust mass, with a constant conversion factor $f_{\mathrm{d}}=0.08$, and use the dust model developed by Weingartner & Draine (2001); Li & Draine (2001) for the Milky Way (MW). Our value for $f_{\mathrm{d}}$ is lower than the one determined for the MW ($\sim$ 0.3-0.5). The choice is motivated by the constraints obtained by Behrens et al. (2018a) in order to fit the ALMA observations from Laporte et al. (2017). With this value, the dust mass in Althæa is about $1.6\times 10^{6}{\rm M}_{\odot}$. As we will see, the exact value does not affect our results substantially.

We note that Laursen et al. (2009) have developed a model explicitly taking into account that dust may have a different abundance in ionized regions (like H$\scriptstyle\rm II$ regions). This translates into reducing $f_{\mathrm{d}}$ in such regions in our notation. We do not follow this route, but will analyze an extreme case (setting $f_{\mathrm{d}}$ to zero in H$\scriptstyle\rm II$ regions), see the next section.

2.2.3 Hydrogen ionization

In order to do the Ly$\alpha$ RT, we need to estimate the neutral hydrogen content, i.e. the fraction of hydrogen in ionized and molecular state must be calculated. The hydrodynamical simulation provides the molecular hydrogen fraction of the gas. To obtain the ionized fraction of gas, we assume the prescription of Chardin et al. (2017), which extends and refines the model of Rahmati et al. (2013) for $z>6$. Such prescription accounts for collisional excitation and a Haardt & Madau (2012) UV background radiation, that is suppressed at high gas densities. At typical ISM densities, the UV background plays essentially no role (e.g. Gnedin, 2010), as on these scales local sources are the dominant producers of ionizing radiation.

To recover the effect of local sources, we post-process the hydrodynamical simulation by considering H$\scriptstyle\rm II$ regions around young stars. The size of an H$\scriptstyle\rm II$ region can be estimated by the radius of a Strömgren sphere:

[TABLE]

with $Q$ being the ionizing photon rate, $n$ the gas number density, and the Case-B recombination rate $\alpha_{B}=2.6\times 10^{-13}$ cm3/s. Consider that stellar clusters in the simulation have masses of about $2\times 10^{4}$ ${\rm M}_{\odot}$; with the chosen IMF and stellar SED (see Sec. 2.2.1), we find $Q\sim 4\times 10^{50}$ s*-1*, in turn yielding $R_{S}\sim 10$ pc for $n=10^{2}{\rm cm}^{-3}$. Such density is selected since it is the typical of star forming regions in the simulation and it is consistent with the expected one in H$\scriptstyle\rm II$ regions (e.g. Dopita & Evans, 1986; Kewley & Dopita, 2002; Kewley et al., 2013). Note that the Strömgren solution is reached on timescales of $\sim 10^{4}$ yr while pressure equilibrium is marginally achieved on time scales of 106 yr, when photoevaporation effects start to play a role (Decataldo et al., 2017), whereas the stars responsible of the ionizing flux have lifetime of around 107 yr (e.g. see Lejeune & Schaerer, 2001).

This simple estimate neglects the complex evolution of H$\scriptstyle\rm II$ regions, the history of each star-forming region, the interaction of overlapping region, and the effect on gas dynamics. In our simulations, the linear size of the most-refined cells is 25 pc, and we therefore chose to treat each of the cells hosting stellar cluster as ionized bubble, setting the gas neutral fraction to zero333If a cell is not refined to the finest level $L_{max}$, we set its partial ionization by ionizing a physical volume of (25 pc)3 in the cell, that is, we set the ionization level to 1/$(L_{max}-L)^{3}$, where $L$ is the refinement level of the cell in question..

The adopted method gives a crude approximation of the radial profile within H$\scriptstyle\rm II$ regions, however it is justified since we do not have the spatial resolution to resolve the ionized bubbles. Setting the neutral fraction to zero in regions around young stars maximize their effect on the Ly$\alpha$ RT, since the gas becomes transparent to the Ly$\alpha$ radiation. Because of this, we do not need to modify the temperature or velocity fields in the H$\scriptstyle\rm II$ regions. In total, our prescription turns about $7\times 10^{8}{\rm M}_{\odot}$ of neutral gas into ionized gas, i.e. about 20% of the gas mass in the ISM.

As stated above, this prescription is part of our fiducial setup; below, we will sometimes compare to the results we obtain when we instead completely ignore local sources of ionizing photons. We will refer to this model as “no H$\scriptstyle\rm II$ bubbles” model.

2.3 Radiative transfer

For the Ly$\alpha$ RT, we use our own code, iltis (Sec. 2.3.1). For the continuum UV RT (Sec. 2.3.2), we make use of the public skirt444version 8, http://www.skirt.ugent.be code (Baes et al., 2003; Baes & Camps, 2015; Camps & Baes, 2015; Camps et al., 2016). While these calculations rely on Monte-Carlo methods, the [C $\scriptstyle\rm II$] emission can be modeled using an semi-analytic approach (Sec. 2.3.3). We describe the three approaches in more detail below.

2.3.1 Ly$\alpha$ line

We use an updated version of the code iltis already used in Behrens & Niemeyer (2013); Behrens & Braun (2014); Behrens et al. (2014, 2018b). We refer the reader to these papers for more details on the implementation, and to the public code repository555https://github.com/cbehren/Iltis. In the following, we briefly summarize the RT algorithm.

iltis makes use of the usual Monte-Carlo approach for Ly$\alpha$ RT; it follows the escape of a significant number of tracer photons from their emission sites through the simulation volume666Note that iltis can directly handle Ramses datasets, using a modified version of the ramsesread++ library originally written by Oliver Hahn (see https://bitbucket.org/ohahn/ramsesread).. Photons can be scattered by gas, or scattered/absorbed by dust. In case they are absorbed by dust, they are considered “lost”, since reemission will take place in the IR regime. The optical depth a photon experiences is determined by the density of the neutral gas along its path, the gas temperature, the bulk velocity of the gas, the density of dust, and the photon frequency. Initial frequencies, points of interaction, thermal velocities of scattering atoms, and scattering angles are drawn from appropriate probability distributions. On top, we also consider the redshifting of photons due to the Hubble flow, whose effect is to shift photons at wavelengths shorter than the line center (“blue” photons) back into the line center, rendering them subject to a large attenuation even in the diffuse IGM, far away from the ISM.

In order to produce meaningful surface brightness maps, we use the peeling-off technique (e.g. Zheng & Miralda-Escude, 2002; Dupree & Fraley, 2004). At each scattering event, we evaluate the probability that the photon directly escapes into the direction of the observer, which can be understood as weighting this contribution with the total flux escaping in the specified direction. This probability can be written as

[TABLE]

where $p(\theta)$ is the probability to be scattered in the direction of the observer, and $\tau_{o}$ is the optical depth the photon would penetrate before reaching the observer. Technically, this means that we need to integrate the optical depth from the point of interaction throughout the simulation box.

We perform radiative transfer at the spatial resolution of the hydro simulation (25 pc). For comparison, Smith et al. (2015) performed Ly$\alpha$ RT at considerably higher resolution (up to 830 au or 0.004 pc), but in a smaller box, as they considered halos of mass $\sim$ 10${}^{7}{\rm M}_{\odot}$. Trebitsch et al. (2016) reached 434 pc; (Verhamme et al., 2012, also see their table 1 for more references) and Behrens & Braun (2014) reached 18/30 pc, but with an idealized setup for an isolated disk galaxy and no cosmological initial conditions. Recently, Smith et al. (2019) presented high-resolution simulations ($<10pc$) for a LAE with a stellar mass of about $5\times 10^{8}{\rm M}_{\odot}$.

More technical details about the numerical setup can be found in App. A. In brief, we launch $10^{3-4}$ photon packages per source from the center of the hydro cell. The intrinsic spectrum is a Gaussian of width 10 km/s, and we use a standard acceleration scheme to avoid core scatterings (Dijkstra et al., 2006).

When presenting the results, note that our definition of the Ly$\alpha$ escape fraction $f_{\alpha}$ is related to the flux actually observable along an individual line of sight (los). It therefore combines two very different mechanisms: (a) direct loss of photons due to absorption by dust in ISM on the one hand, and (b) indirect losses due to diffuse scatterings out of the los in the CGM/IGM. In the former case, we expect the photons to be re-emitted in the IR; in the latter, we expect them to contribute to a diffuse background or an extended halo. On top, we do consider here the escape along one specific los; even in a dust-free galaxy, and without considering the IGM, we expect anisotropic escape to take place due to the relative distribution of gas/stars and peculiar velocity fields.

To recover observables like the line flux from the galaxy, the escape fraction, or spectra, we integrate over a circular aperture of radius 2 kpc around the center of the galaxy. We chose this size of the aperture to include only the escaping radiation from Althæa, which has a size of about 1 kpc, and not the diffuse background coming from Ly$\alpha$ scattered in the CGM/IGM.

2.3.2 UV and IR continuum

We make use of skirt, a publicly available code for detailed simulations of UV continuum, absorption by dust, and re-emission of the absorbed energy in the IR. The code is very flexible thanks to its ability to load e.g. data from AMR or Smoothed-particle hydrodynamics (SPH) simulations, a variety of implemented, SED (Sec. 2.2.1), dust models (Sec. 2.2.2), and physical mechanisms.

Since we are mostly interested here in the escaping UV continuum rather than the mid/FIR emission, we do not take into account self-absorption of dust or stochastic heating. Note that the modeling of the UV/IR emission used here is the same as in Behrens et al. (2018a).

2.3.3 [C $\scriptstyle\rm II$] line emission

We compute [C $\scriptstyle\rm II$] with a semi-analytical method detailed in Pallottini et al. (2017a, see appendix C). The basic features are as follows. First we compute the column density of C $\scriptstyle\rm II$ in each cell by using a grid of photoionization models obtained with cloudy777version 13.03, see https://www.nublado.org/ (Ferland et al., 2013). For the cloudy calculations, we assume a uniform UV interstellar radiation field with a MW-like spectral shape and an intensity scaled to the Althaea SFR. Then, the [C $\scriptstyle\rm II$] emissivity is computed as a function of gas temperature and density following Dalgarno & McCray (1972), and Vallini et al. (2013, 2015). Finally, we account for CMB flux suppression as in Vallini et al. (2015); Pallottini et al. (2015). As the [CII] line is observed against the CMB, the spin temperature of the transition must differ from the CMB one to be detectable. As this can be obtained essentially only via collision with other species, the emission for low density gas is largely suppressed. This is an important effect that needs to be included.

In the present paper the [C $\scriptstyle\rm II$] emission is used only to calculate its velocity shift with respect to the Ly$\alpha$ line. A full analysis of the kinematic/morphology of the [C $\scriptstyle\rm II$] line from Althæa is presented in Kohandel & et al. in prep. (2018). In order to acquire meaningful line shifts, in this work the systemic velocity of the galaxy is taken to be the one of the [C $\scriptstyle\rm II$]  which is determined by fitting the [C $\scriptstyle\rm II$] spectrum of Althæa for each los with a Gaussian.

We note that [C $\scriptstyle\rm II$] is optically thin, and therefore solely affected by the local conditions of the emitting gas.

3 Results

3.1 Morphology and Ly$\alpha$ escape fraction

We start by analyzing the galaxy morphology resulting from our fiducial model ($z=7.2$). This is shown in Fig. 2 for two lines of sight, chosen to be face-on and edge-on, respectively.

We first focus on the face-on case. As Althæa is very compact, the mean neutral hydrogen column density is $N_{HI}=6\times 10^{21}\mathrm{cm}^{-2}$; in the proximity of the star-forming clumps, the typical neutral column density is larger ($\sim 3\times 10^{22}\mathrm{cm}^{-2}$); however, since the plot takes into account the ionized bubble model outlined above, some of the clumps display a core of reduced column density, owing to the ionization from local sources.

The los component of the velocity of the neutral gas, $v_{\rm los,HI}$, shows a very complex pattern, owing to the dynamics of accretion through filaments, tidal stripping of gas streams from nearby satellites, and SN-driven outflows. The distribution of stellar mass has a half-mass radius of 0.5 kpc and shows a diffuse stellar component extending to about 1 kpc, also featuring multiple clumps in the disk. Typical velocities within the disk are few $\sim 100$ km/s. The diffuse component mainly consists of relatively old stars, as star formation in Althæa takes place mostly in the central region and in the (relatively small) clumps (see Behrens et al., 2018a). Such star forming clumps are characterized by high stellar mass surface densities ($\Sigma_{\star}\sim 10^{11}{\rm M}_{\odot}/{\rm kpc}^{-2}$). Young stars producing most of the intrinsic Ly$\alpha$ radiation are predominantly found in these regions, as clearly shown by the corresponding map of the intrinsic Ly$\alpha$ emission $S_{\alpha,\rm int}$. Clumps can reach surface brightness of $S_{\alpha,\rm int}\sim 10^{12}{\rm L}_{\odot}/{\rm kpc}^{-2}$.

By comparing $S_{\alpha,\rm int}$ with the surface brightness of the escaping Ly$\alpha$ emission, $S_{\alpha}$, one sees that Ly$\alpha$ radiation escapes virtually only from the south-west side of the disk, with a peak located in an inter-arm region at about 1 kpc from the center. Such region is characterized by a low gas column density ($N_{H}\sim 10^{19}{\rm cm}^{-2}$), and shows the signature of an outflow expanding at a velocity $v_{\rm los}\sim 100$ km/s along the los.

In general the UV emission is more diffuse than the Ly$\alpha$ one. While the central region is Ly$\alpha$-dark, it is relatively bright in UV; the same is true for most of the star-forming clumps. However, the UV continuum map also shows a surface brightness peak co-located with the Ly$\alpha$ peak. Both observations can be explained as a consequence of the larger attenuation of the resonantly scattered Ly$\alpha$ radiation: In regions that allow Ly$\alpha$ to escape, UV will escape as well.

The galaxy-averaged Ly$\alpha$ escape fraction is $f_{\alpha}=3\times 10^{-4}$; while the intrinsic Ly$\alpha$ luminosity of Althæa  is $2\times 10^{44}$ erg/s, only $5\times 10^{40}$ erg/s reaches the observer888We have checked that the global morphology and escape fraction do not vary appreciably if we drop the fiducial model assumption that the gas around stars is ionized.. Recall that in the computation of $f_{\alpha}$ we include the effects of both the ISM and CGM/IGM (see Sec. 2.3.1). For comparison, the UV escape fraction (evaluated at 1500 Å) in this case is $f_{\mathrm{UV}}=0.05$.

For the fiducial case, given the FWHM of $\sim$ 2 Å, we obtain a total Ly$\alpha$ flux of $\sim 10^{-19}$ erg s*-1cm-2*. The detection of these Ly$\alpha$ fluxes is challenging for current telescopes (but see Wisotzki et al., 2018), but still achievable with future instruments. For example, the Multi-Object Spectrograph planned for E-ELT, is expected to achieve a sensitivity of $\sim 10^{-19}$erg s*-1cm-2* in 20 hr of observing time (Evans et al., 2015).

We show the same set of maps also for the edge-on direction999The orientation of the plots has been chosen such that the face-on direction discussed before is pointing upwards (along the y axis) in these panels. in Fig. 2. The column densities in this case are obviously higher. We can see that while the disk is relatively thin, it is far from being flat, owing to the complex gas motions resulting in a bent and warped disk structure. This is not only visible in the gas density and kinematical structure, but also in the distribution of stars.

As expected, the escape of Ly$\alpha$ radiation is much more difficult in this direction: in fact, we find $f_{\alpha}=5\times 10^{-6}$, i.e. about 60 times lower than in the face-on case. The Ly$\alpha$ transmission in the edge-on case exclusively comes from photons that are scattered into this los at several 100 pc above the disk. We also note a strong above/below asymmetry, that is, photons preferentially come from above the disk. This asymmetry is confirmed by the more systematical variation of the inclination angle discussed in Sec. 3.2. For comparison, $f_{\mathrm{UV}}=0.02$ in this case.

Finally, we derive the Ly$\alpha$ line EW. While the intrinsic value obtained by considering the attenuated SED of the stars in Althæa is 103 Å, the equivalent width (EW) after RT decreases to 0.7 (0.04) Å for the face-on (edge-on) case.

3.2 Inclination effects

Several authors (Laursen et al., 2009; Verhamme et al., 2012; Behrens & Braun, 2014; Behrens et al., 2014) have pointed out that due to its resonant nature, the Ly$\alpha$ line is particularly susceptible to geometrical effects or viewing angles, owing to the large optical depth Ly$\alpha$ photons experienced from neutral gas compared to UV continuum photons. However, this conclusion is geometry-dependent. If photons travel in a diffuse ISM in which dense, optically-thick clouds are embedded, it can be shown that under certain conditions Ly$\alpha$ the escape of Ly$\alpha$ photons is enhanced as they simply bounce off the surface of the dense clouds instead of penetrating them like continuum UV photons. However, this so-called Neufeld effect (Neufeld, 1991) is only valid in a very narrow regime of parameters for e.g. the velocity dispersion of the clouds (Laursen et al., 2013; Duval et al., 2014). Hence, in general, it is expected that the inclination effect leaves a stronger imprint on the Ly$\alpha$ properties. We investigate the possible inclination effect on our results in the following.

To this aim, we have generated 3072 los by using an equal area and iso-latitude tessellation (Gorski et al., 1999), and ran Ly$\alpha$ and UV/IR simulations for all of them. For illustration, the Ly$\alpha$ spectra of some los are shown in Fig. 3. All the spectra show a red peak as expected, but the total line luminosity, and position of the peaks vary from los to los.

In the left panel of Fig. 4, we show the escape fraction from each los in a Mollweide plot. $f_{\alpha}$ varies by about six orders of magnitude and the brightest los reaches a luminosity of about $5\times 10^{41}$ erg/s. The line width (FWHM) varies largely from 60 km/s to 2000 km/s.Observationally, line widths of few hundred km/s are typical at both lower redshift (e.g. see Verhamme et al., 2018, figure 2) and higher redshift (e.g. Vanzella et al., 2011; Matthee et al., 2017; Pentericci et al., 2018). In Fig. 5, we show the Ly$\alpha$ luminosity, inclination, and line full-width half-maximum (FWHM) for each of the 3072 lines of sight considered. A value $\cos\theta=0$ ($\cos\theta=1$) corresponds to a edge-on (face-on) orientation. In contrast to e.g. Verhamme et al. (2012); Behrens et al. (2014), the preferential Ly$\alpha$ escape along face-on los is less pronounced; the behavior is instead much more random. This can be understood by considering that Behrens et al. (2014); Verhamme et al. (2012) used an idealized, non-cosmological setup for a disk galaxy, whereas Althæa is a complex object in its cosmological environment, with gas flowing in along filaments, a warped disk with many small-scale features, and satellites swirling around the main disk. Comparing with the results of Laursen et al. (2013), we do not find a significant increase of the EW in face-on directions, similar to what Yajima et al. (2012) found in a simulated galaxy at intermediate redshift ($z=3.1$). In the latter case, inclination effects become prominent only at later times ($z\approx 0$), when the galaxy has settled into an ordered disk.

Consistent with the results in the above Section, the dependency of the UV continuum on the viewing angle (Fig. 4, center) is shallower compared to the Ly$\alpha$; the emerging UV luminosity varies with inclination by a factor of 4 at most (see also Fig. B1 in Behrens et al., 2018a) and it is essentially invariant with respect to the azimuthal angle. The Ly$\alpha$ line EW for every individual los is shown in the right panel of Fig. 4. Its angular distribution broadly correlates with the $f_{\alpha}$ map, and shows a maximum value of 5.9 Å.

3.3 What quenches the Ly$\alpha$ line?

To get a deeper insight on the physics determining the low $f_{\alpha}$ values found, it is useful to inspect the line spectra along different los. We build spectra for the two los shown in Fig. 2.

In Fig. 6, we show the resulting Ly$\alpha$ spectra for the face-on (top) and edge-on (bottom) cases. Both the fiducial run and the “no H$\scriptstyle\rm II$ bubble” case are shown, together with some tests discussed below. In the following, we try to disentangle different effects determining the resulting escape fraction and the shape of the spectrum.

IGM attenuation

The blue wing of the line is almost completely blanketed by the intervening IGM. This is a well known effect (Dijkstra et al., 2006; Laursen, 2010) caused by the Hubble flow: photons leaving the galaxy blue-wards of the Ly$\alpha$ line center will redshift into the line center where they are scattered out of the los. Indeed of the Ly$\alpha$ radiation escaping the ISM is blue-shifted, indicating that infall (rather than outflow) dominates the overall dynamics. This is not surprising, as Althæa is by far the most massive object in its environment.

IGM quenching can be distilled by switching off Hubble expansion in the simulation. We stress that we do not render the surrounding gas outside of the galaxy, but only switch off the Hubble expansion. In this case we find a “classical” double-peaked spectrum with a flux of $2\times 10^{41}$ erg/s in the blue peak (dotted line in Fig. 6). This blue peak is completely erased by the IGM in the presence of the Hubble flow, while the red peak is not strongly affected (given that at $z\approx 7.2$ the mean IGM ionization is $x_{\rm HI}\approx 10^{-3}$ in the simulation). Note that switching off the Hubble flow changes the velocity field that photons experience everywhere; this is the reason why the red part of the spectrum changes when the Hubble flow is turned off. In the case of this particular los, the effect on the red wing is quite strong, probably owing to the specific density structure above the disk along this los; for other los, it is at the 10% level, so this has to be seen as a peculiarity of the los.

We can factorize the total escape fraction as $f_{\alpha}=f_{\mathrm{ISM}}f_{\mathrm{IGM}}$, i.e. the fraction of photons escaping the galaxy times the fraction of photons reaching us through the IGM. For the face-on case, the fiducial model predicts $f_{\mathrm{ISM}}=0.001$, $f_{\mathrm{IGM}}=0.23$. If we instead of switching off the Hubble flow remove all the gas outside a radius of 30 pkpc to assess the influence of the IGM, $f_{\mathrm{IGM}}$ goes up by a factor of 1.8, and $f_{\mathrm{ISM}}$ goes down accordingly.

We conclude that while the IGM does reduce the escaping flux by a factor $\leq 5$ it is not the major cause of line quenching.

Neutral hydrogen in the ISM

The presence of neutral gas increases the scattering rate and ultimately the probability for the Ly$\alpha$ photon to be absorbed by dust. To check its importance in determining the escape fraction, in Fig. 6 we also plot the spectra from the model lacking H$\scriptstyle\rm II$ regions around stellar clusters (dashed lines). These spectra deviate moderately from the fiducial ones. They exhibit a larger line shift ($\sim 1.8$ Å vs. $\sim 1.0$ Å), and the flux is reduced by a factor of $\sim$ 2.101010Line shifts are obtained by fitting a skewed Gaussian to the Ly$\alpha$ line profile.. However, in terms of the escape fraction, both the fiducial and the no bubble case differ only by factor of $\sim 2$ . In other words, increasing the HI fraction close to stars does not affect the escape fraction significantly. Stated differently, the gas close to the star forming sites alone does not appreciably increase the path length of Ly$\alpha$ photons, which is mainly set by the residual gas in the galaxy.

Infall/outflows

The ratio of the fluxes bluewards/redwards of the line center depends on the peculiar velocity field. Relatively mild outflows, with mean velocities around $150\,{\rm km}/{\rm s}$, are present in the galaxy; infalling streams also share the same velocity scale (Gallerani et al., 2018). Infall motions are dominating the radiative transfer, as can be seen by the large amplitude of the blue peak in the spectrum when we ignore the Hubble flow (Fig. 6, dotted line). In order to check the effect of the peculiar velocities, we re-run the Ly$\alpha$ RT by down-scaling their norm by a factor of 4.

The resulting spectrum is shown in Fig. 6 (dashed-dotted lines). As now less photons are transferred to the blue side of the line center (and more are transferred to the red side), less flux is scattered out of the los by the CGM/IGM Although $f_{\alpha}$ and EW are enhanced by a factor of $\simeq 2$ with respect to the fiducial model, such increase is not sufficient to promote Althæa into a LAE.

In Fig. 7, we show $f_{\alpha,\mathrm{emitter}}$ for each individual stellar cluster in the fiducial simulation along the brightest los as a function of the dust density at the location of the cluster, and the velocity differential, $\Delta v_{los}$, along the los measured 25 pc above the source. Positive values of $\Delta v_{los}$ correspond to outflows. In addition to the strong negative correlation between dust density and escape fraction, there is a clear, increasing trend of $f_{\alpha}$ with $\Delta v_{los}$. The effect imprinted by peculiar velocities is nevertheless sub-dominant with respect to dust absorption, as we discuss further below.

Absorption by dust

The vast majority of Ly$\alpha$ photons in the simulation are absorbed on the spot, that is, in the gas element they are spawned in. About 60% of all absorptions occur in the production cell in the fiducial model, i.e. within $\sim 15$ pc from their emission spot. In the “no bubble model”, this fraction increases to 95% due to the increased absorption probability driven by an increased pathlength due to more scatterings on hydrogen. It is therefore an interesting question to ask whether the attenuation of the Ly$\alpha$ is mainly affected by (a) the total amount of dust, or (b) by its distribution, given the fact that dust is found to be highly concentrated in the star-forming knots discussed above.

To answer this question, we first down-scaled the dust mass in the whole simulation volume, that is, we reduced $f_{\mathrm{d}}$ by a factor of 10. As a result, we find that $f_{\alpha}$ increases from $3\times 10^{-4}$ to $4\times 10^{-3}$ in the face-on case. However, the UV continuum raises as well, with the net result that the EW is lower ($EW=0.3\;\AA{}$ for the face-on case).

Next, we ran a different set of simulations in which we selectively removed the dust from each star-hosting cell, that is, we rendered H$\scriptstyle\rm II$ regions dust-free. We refer to this model as the “dust-free bubbles” model. Since the dust is highly concentrated in the star forming knots, this procedure removes about 60% of the total amount of dust from Althæa. Note that the total amount of dust is still 4$\times$ higher than in the previous case in which we decreased the dust content everywhere. In this case, the escape fraction increases by a factor of 100, yielding $f_{\alpha}=0.03$ (face-on); this corresponds to a Ly$\alpha$ luminosity of $6\times 10^{42}$ erg/s. However, as the UV increases as well, the boost of the EW is more modest (4.2 Å).

This finding is also backed by a simple estimate of the optical depth due to dust along the line of sight that Ly$\alpha$ photons suffer from. Averaging over the (face-on) lines of sight from young stellar cluster which we consider a source of Ly$\alpha$ radiation, we find a mean optical depth of $\tau_{D}=17$, albeit with a wide distribution skewed to small values (for example, 40% of the emitters have $\tau_{D}<5$). The optical depth from dust in the source cells alone is about 8 on average, again with a large scatter.

These two experiments show that to transform Althæa into an object with significant transmission of Ly$\alpha$a selective depletion of the dust contained in H$\scriptstyle\rm II$ regions is required. Quite naturally, one may wonder whether the dust-free bubbles model is based on some physical grounds. While it is plausible that H$\scriptstyle\rm II$ regions near young stellar clusters might be dust-depleted to some extent (e.g. through dust destruction in the ionized regions, e.g. Pavlyuchenkov et al., 2013), we cannot uniquely individuate a physical mechanism that completely destroys the dust in the vicinity of the star-forming regions. However, it is worth mentioning the results by Draine (2011) who first studied radiation pressure effects on the internal density structure of H$\scriptstyle\rm II$ regions. He finds that if the H$\scriptstyle\rm II$ region is powered by stars with a strong ionizing flux surrounded by an initially dense medium, the gas and dust in the cavity get evacuated and pile up into a shell marking the ionized boundary.

This effect is not fully caught by our simulations. In spite of their high resolution ($<30$ pc), our simulations might be insufficient to properly describe both the Ly$\alpha$ RT in the inherently turbulent medium characterizing molecular cloud interiors, and the above radiative feedback effects: In order to do so, sub-pc resolution would be required.

With the above caveats, we compare in more detail the fiducial and the dust-free bubbles model. For the dust-free bubbles, we have run simulations with 48 lines of sights. In these simulations, we find qualitatively similar dependencies of the Ly$\alpha$ escape on inclination as for the fiducial model. However, the Ly$\alpha$ and UV luminosities are higher in the dust-free bubbles model (see Fig. 8), with EW reaching 22 Å. In this case, Althæa would be classified as a LAE by the EW cuts typically adopted by observational practice ($\sim 20$ Å). The maximum escape fraction is $f_{\alpha}=0.15$. For one los with a high EW, we show in Fig. 8 both the resulting spectrum (left panel) and the morphology (right panel) of the escaping Ly$\alpha$ radiation; we indicate the continuum flux level with a dotted line. The line widths in the dust-free case are similar to the fiducial case; for the lines of sight with EW > 20 Å, we find a width of $\approx 600$ km/s. In the dust-free bubbles model, the total Ly$\alpha$ flux predicted by our model ($\sim 2\times 10^{-17}$ erg s*-1* cm${}^{-}2$, for a FWHM$\sim$3 Å) is consistent with the ones typically observed in $z\sim 7$ LAEs (e.g. Vanzella et al., 2011; Caruana et al., 2014; De Barros et al., 2017), thus being detectable in $\sim$ 15 hr with the FORS2 spectrograph on the ESO Very Large Telescope.

In Fig. 9 we report the relation between $f_{\alpha}$ and $E\rm(B-V)$ for both the fiducial and dust-free bubbles models at $z=7.2$111111To calculate $E\rm(B-V)$, we used the optical depth in the V/B-band at 0.552/0.442 $\mu\mathrm{m}$ as calculated in the UV continuum RT.. The main effect of the dust-free bubbles model is to shift horizontally rightwards the points in the plot.

For comparison, in Fig. 9 we also plot available observational data from Kornei et al. (2010); Hayes et al. (2011) for low-redshift galaxies ($z=2-3$). The lower (upper) limits taken from Hayes et al. (2011) correspond to objects detected in Ly$\alpha$ (Hα), but not in Hα (Ly$\alpha$).

While in this paper, we focused on Althæa at redshift $z=7.2$, we show the escape fraction $f_{\alpha}$ for three different evolutionary stages and 48 lines of sight per stage in Fig. 10 together with the median (crosses) and the lower/upper quartiles (bars). At the additional redshifts of 6.5 and 6.1, Althæa has a stellar mass of 1.4 and 2.0 $\times$ $10^{10}{\rm M}_{\odot}$, respectively. We find only a mild variations from snapshot to snapshot in terms of the escape fraction; interestingly, the extinction of the UV is the highest at $z=7.2$, whereas the total dust mass is the lowest at this evolutionary stage. This can be related to the fact that in this stage, recent burst of star formation have deposited metals in their immediate surroundings, leading to large attenuation.

Summary

We conclude that the clumpy distribution of dust is the key factor inducing Ly$\alpha$ line quenching, with the peculiar velocity field playing a sub-dominant role. In turn, the dust is clumpy because grains tend to cluster around their production factories, at these high redshifts predominantly massive stars (Todini & Ferrara, 2001).

3.4 Ly$\alpha$–[CII] line shift

In Fig. 11, we show the relative shift between the Ly$\alpha$ and [CII] lines as a function of the Ly$\alpha$ luminosity for different los and the fiducial/dust-free bubbles model. Within each of the two samples, we find a weak anti-correlation, i.e. the shift increases towards lower Ly$\alpha$ luminosities (and EW). In our case, this is readily explained by the low escape fractions. As already discussed in the literature (see e.g. Laursen et al., 2009), Ly$\alpha$ photons far in the wings are more subject to absorption, since they typically have traveled a longer path through the dusty medium; they originate from regions of high gas densities, and undergo a large number of scatterings, which in turn induce changes of direction and enhance their path lengths. Enhanced pathlengths boost the attenuation of these photons, thus dust damps preferentially the outer parts of the spectrum, even though the dust absorption cross section itself is only weakly depending on frequency. Accordingly, we find the average number of scatterings prior to escape to be anti-correlated with the observed luminosity along a given los.

Verhamme et al. (2018) reported a positive correlation between line shift and Ly$\alpha$ line width. We do not find such a relation; this holds also in the dust-free bubbles case, i.e. when Althæa is a LAE. However, we note that even in this favorable case, the escape fraction is at the percent level due to severe dust attenuation. This suggests the relationship might only valid in low-attenuation systems like the ones Verhamme et al. (2018) consider, where the escape is largely driven by the outflow of gas. In contrast, in our case, it is the dust content that dominates the escape, and outflows only play a minor role.

4 Summary and Conclusions

Performing radiative transfer, we have post-processed Althæa, a simulated prototypical LBG at $z=7.2$ (Pallottini et al., 2017b), to obtain its Ly$\alpha$, continuum, and [C $\scriptstyle\rm II$] line properties. We used the publicly available dust continuum code skirt (Baes et al., 2003; Camps & Baes, 2015), a semi-analytic model for the [C $\scriptstyle\rm II$] line (Vallini et al., 2015; Pallottini et al., 2017a), and iltis, an updated version of the Ly$\alpha$ code used by Behrens et al. (2018b); Behrens & Braun (2014); Behrens et al. (2014). We included both the interaction of radiation with gas and with dust, and employed a simple model for H$\scriptstyle\rm II$ bubbles forming around star-forming regions, as we do not directly simulate ionizing radiation. We performed our simulations using up to 3072 lines of sight.

In our fiducial simulation, Ly$\alpha$ radiation escapes Althæa solely from one side and one quadrant on the outskirt of the disk, characterized by a low column density gap in between two filamentary structures. Both the star-forming core and the $\sim$ 100 pc-scale star-forming clumps remain dark in Ly$\alpha$. Althæa is a resilient LBG with low Ly$\alpha$ EWs (<6 Å) compared to an intrinsic Lya EW of 103 Å; the implied escape fractions are less than 0.1%.

In contrast to studies investigating highly idealized, isolated disk galaxies (Verhamme et al., 2012; Behrens et al., 2014), we only find a weak correlation between inclination of the los and the escape fraction/Ly$\alpha$ luminosity because of the complex, anisotropic structure of gas/dust in and around Althæa. However, the escape fraction $f_{\alpha}$ strongly depends on the specific los, with variations of up to six orders of magnitude.

The emerging spectra typically show a well-defined asymmetric, single-peaked shape located redwards of the Ly$\alpha$ line center. However, the spectrum prior to entering the IGM is double-peaked. The dominant blue peak is washed out by the IGM, as bluer photons are progressively redshifted into resonance and scatter on residual hydrogen out of the line of sight. We have investigated the cause for the very small escape fraction found in our fiducial simulation, and find that transfer is dominated by the clumpiness of dust close to star-forming regions, absorbing Ly$\alpha$ photons very efficiently. Outflows do not generally play a major role for the escape mechanism, apart for some specific los.

To turn our LBG into a LAE, we have artificially removed dust from the galaxy. While reducing the overall dust content does not make Althæa a LAE, by selectively removing the dust content from star-forming regions we reach EWs up to 22 Å and escape fractions of up to $f_{\alpha}$=15%. This further justifies our conclusion that it is the dust within the star-forming clumps that dominates the RT. As a by-product, neither the precise value for the dust-to-metal ratio, nor the inclusion/exclusion of ionized regions around stellar clusters change our results dramatically. While in this paper we have mostly focused on one evolutionary stage of Althæa, we have also run similar simulations on two snapshots at different redshifts that exhibit lower star formation rates, without dramatic changes in the escape fraction.

Taken at face-value, our results raise the question of what mechanism could turn LBGs into LAEs and vice versa in the popular duty-cycle scenario. Such a mechanism would require to lower the attenuation by orders of magnitude within a relatively short time. Stellar clusters that are less attenuated in our simulations have typically moved out of their birth clouds, and are therefore too old to contribute to the ionizing flux generating Ly$\alpha$. A stronger feedback mechanism dispersing the molecular clouds more quickly, however, might be a viable solution if it can also get rid of the dust.

Finally, we investigated the relation between the Ly$\alpha$ luminosity and the line shift of the Ly$\alpha$ line, measured with respect to the [CII] line. For the face-on fiducial model, the Ly$\alpha$ shift with respect to the [CII] line is 1 Å ($\sim 250$ km/s). We find a negative correlation between the two, i.e. the shift increases towards lower Ly$\alpha$ luminosities (and EW). This does not support the suggestion by Pentericci et al. (2016) that smaller line shifts might correlate with a more strongly damped Ly$\alpha$ line due to an enhanced attenuation in the IGM. This difference might be related to the overall small value of $f_{\alpha}$ in our case.

Since our simple model of ionized bubbles around stellar clusters can only be considered a rough estimate, direct RT of the ionizing photons emerging from stars will be required to get the full picture. Another caveat of our work is the lack of substructure of molecular clouds in our simulation due to our resolution limit of 25 pc. Sub-grid models for the escape of photons from the multiphase molecular clouds will be necessary (also see the work by Hansen & Oh, 2006; Gronke et al., 2016; Kimm et al., 2019), including e.g. the evaporation of the cloud by young stars. We plan to address these issues in future work.

Acknowledgements

AF acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. LV acknowledges funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 746119. We acknowledge use of the Python programming language (Van Rossum & de Boer, 1991), and use of the packages IPython (Perez & Granger, 2007), Matplotlib (Hunter, 2007), Numpy (van der Walt et al., 2011), Pymses (Labadens et al., 2012), and Jupyter (Kluyver et al., 2016). This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013).

Appendix A Details on the Ly$\alpha$RT

A.1 Peeling-off algorithm

As presented in Sec. 2.3, the peeling-off algorithm requires to integrate the optical depth along the los to each scattering event. In order to speed up this calculation, we discard contributions whenever the line of sight optical depth exceeds $\tau_{o}=20$; since the weight of each contribution scales exponentially, this means we discard contributions that add less than $O(10^{-9})$, which translates to $10^{33}$ erg/s in our units.

A.2 Hubble flow

Photons redshift in-between scatterings. In order to make sure we resolve the passage of a photon on the blue side of the spectrum through the line center, we invoke a limit on the maximum step size of $\sim 1$ kpc, which corresponds to a velocity shift of 0.9 km/s at $z=7.2$.

A.3 Dust properties

For the Ly$\alpha$, we used the same dust model as in skirt, that is, the Weingartner & Draine (2001) model. However, as we are interested in the line only, we assumed the cross sections, the albedo, and the asymmetry parameter to be the constant with respect to frequency. Their values changes only modestly (e.g., the cross section varies by $\sim 3\%$ within $\pm 10\;\AA{}$ of the Ly$\alpha$ line center), and thus we chose to set them to their values at the Ly$\alpha$ line center.

We explicitly checked that UV continuum RT in skirt and Ly$\alpha$ RT in iltis are consistent with each other, that is, we checked that the Ly$\alpha$ RT and the UV continuum RT yield the same results at the Ly$\alpha$ line center if we only consider dust (by removing all gas artifically).

A.4 Parameters of the RT

We launch Ly$\alpha$ photons from each stellar cluster that has an intrinsic Ly$\alpha$luminosity of above 1038 erg/s in order to reduce the number of sources we need to consider. The intrinsic spectrum is a Gaussian with a width of 10 km/s; As we expect our Ly$\alpha$ photons to originate from recombinations in the vicinity of young stars, the width is motivated from the typical velocity dispersion in such regions. We arrive at a total number of 2590 sources. We launch a minimum of $10^{4}$ ($5\times 10^{3}$ for the runs using 48 lines of sight) tracer photons (equivalently; photon packages) per source. While this number is small given the fact that the escape fraction is very low, we acquire of the order $O(10^{6})$ contributions from the peeling-off scheme and find good convergence. For the runs with 3072 lines of sight, we only launch $10^{3}$ photon packages in the Ly$\alpha$  transfer. We have verified that these runs have sufficiently converged using the simulations with a lower number of lines of sight. We make use of the usual acceleration scheme, skipping core scatterings by cutting the thermal distributions of scattering atoms if the dimensionless frequency of the photon is $x<3$(e.g. Dijkstra et al., 2006). For the continuum photons, we use $10^{6}$ photon packages, as this value is sufficient to achieve convergence.