Lyman-alpha emitters in a cosmological volume II: the impact of the intergalactic medium

TL;DR

This paper models how the intergalactic medium affects the observability and clustering of Lyman-alpha emitters, revealing redshift-dependent couplings that impact cosmological measurements like baryon acoustic oscillations.

Contribution

It introduces a new model incorporating Lyman-alpha radiative transfer in galaxy formation simulations, showing its effects on LAE properties and large-scale clustering.

Findings

LAEs tend to have low metallicities and intermediate star formation rates.

At high redshift, LAEs are linked to IGM velocity and density gradients.

The model predicts modifications to the BAO peak shape and position.

Abstract

In the near future galaxy surveys will target Lyman alpha emitting galaxies (LAEs) to unveil the nature of the dark energy. It has been suggested that the observability of LAEs is coupled to the large scale properties of the intergalactic medium. Such coupling could introduce distortions into the observed clustering of LAEs, adding a new potential difficulty to the interpretation of upcoming surveys. We present a model of LAEs that incorporates Lyman-alpha radiative transfer processes in the interstellar and intergalactic medium. The model is implemented in the GALFORM semi-analytic model of galaxy of formation and evolution. We find that the radiative transfer inside galaxies produces selection effects over galaxy properties. In particular, observed LAEs tend to have low metallicities and intermediate star formation rates. At low redshift we find no evidence of a correlation between the spatial distribution of LAEs and the intergalactic medium properties. However, at high redshift the LAEs are linked to the line of sight velocity and density gradient of the intergalactic medium. The strength of the coupling depends on the outflow properties of the galaxies and redshift. This effect modifies the clustering of LAEs on large scales, adding non linear features. In particular, our model predicts modifications to the shape and position of the baryon acoustic oscillation peak. This work highlights the importance of including radiative transfer physics in the cosmological analysis of LAEs.

Tables

Table 1. Table 1: Values of the best fitting K a subscript 𝐾 𝑎 K_{a} used to compute the ionization field in our model.

redshift	$\log K_a\left[\mathrm{cMpc}{h}^{-1}\right]$
2.2	2.79
3.0	3.65
5.7	4.48

Table 2. Table 2: Free parameters as defined in equations 1 and 2 after the calibration with the observed luminosity function for different geometries and redshifts.

redshift	Geometry	$\log {\kappa }_{V,disk}$	$\log {\kappa }_{V,bulge}$	$\log {\kappa }_{N,disk}$	$\log {\kappa }_{N,bulge}$
$z=2.2$	Thin Shell	4.43	4.27	-12.33	-12.11
	Galactic Wind	4.14	4.80	-8.01	-5.58
$z=3.0$	Thin Shell	3.78	4.23	-12.73	-12.06
	Galactic Wind	4.04	4.47	-8.10	-5.69
$z=5.7$	Thin Shell	4.55	3.48	-14.03	-12.24
	Galactic Wind	4.35	3.28	-9.62	-6.85

Table 3. Table 3: Parameters of the the ZZ11 analytic expression (Eq. 6.2 ) to the clustering measurements of our model at redshift 5.7.

Redshift	Geometry	$b$	$\log \widetilde{{\alpha }_4}/r_H$
$z=5.7$	Thin Shell	4.82	-0.44
	Galactic Wind	4.87	-0.2

Table 4. Table 4: Properties of our mock catalogs. In particular we list the redshift z 𝑧 z , the redshift bin width δ ​ z 𝛿 𝑧 \delta z , the size along and parallel to the line of sight ( L ∥ subscript L parallel-to \rm L_{\parallel} and L ⟂ subscript L perpendicular-to \rm L_{\perp} respectively), the Number of mock catalogs within our simulation volume and the median number of LAE in the catalogs for different models as well as the ± 1 ​ σ plus-or-minus 1 𝜎 \pm 1\sigma dispersions.

Authors	$z$	$\Delta z$	$L_{\parallel }$	$L_{⟂}$	$N_{\mathrm{mocks}}$	$⟨N_{\mathrm{LAE}}⟩$
			$\left(\mathrm{cMpc}\right)$	$\left(\mathrm{cMpc}\right)$		Survey	Thin shell	Galactic Wind
Kusakabe et al. (2018)	2.2	0.0773	104.9	93.6	448	1248	${1197}_{-86}^{+96}$	${1193}_{-76}^{+96}$
Bielby et al. (2016)	3.0	0.0633	60.0	119.1	468	643	${633}_{-52}^{+49}$	${633}_{-47}^{+55}$
Ouchi et al. (2018)	5.7	0.0954	43.5	401.5	18	734	${714}_{-28}^{+53}$	${723}_{-51}^{+53}$

Table 5. Table 5: Rest frame equivalent cut EW 0 subscript EW 0 \rm EW_{0} and Ly ​ α Ly 𝛼 \mathrm{Ly}\alpha  luminosity cut L Ly ​ α , cut subscript L Ly 𝛼 cut \rm L_{Ly\alpha,cut} in the different surveys and FLAE mocks.

Authors		${\mathrm{EW}}_{0,\mathrm{cut}}[$Å$]$				$L_{\mathrm{Ly}\alpha ,\mathrm{cut}}\left[\mathrm{erg}{s}^{-1}\right]$
		Survey	Thin Shell	Galactic Wind		Survey	Thin Shell	Galactic Wind
Kusakabe et al. (2018)		$20.0$	$19.52$	$20.3$		${\mathrm{1.62\; 10}}^{42}$	${\mathrm{1.47\; 10}}^{42}$	${\mathrm{1.77\; 10}}^{42}$
Bielby et al. (2016)		$65.0$	$38.45$	$46.27$		${\mathrm{1.62\; 10}}^{42}$	${\mathrm{1.32\; 10}}^{42}$	${\mathrm{1.52\; 10}}^{42}$
Ouchi et al. (2018)		$20.0$	$20.06$	$20.06$		${\mathrm{6.3\; 10}}^{42}$	${\mathrm{6.98\; 10}}^{42}$	${\mathrm{6.64\; 10}}^{42}$

Equations

$${\mathrm{V}}_{\exp ,\mathrm{c}}$$

$${\rm N_{H,c}}=\left\{\begin{array}[]{c l}\kappa_{N,c}\frac{\rm M_{\rm cold,c}}{\rm r^{2}_{c}}&{\rm Thin\;Shell}\\ \kappa_{N,c}\frac{\rm M_{{\rm cold},c}}{\rm r_{c}V_{{\rm exp},c}}&{\rm Galactic\;Wind}\end{array}\right.,$$

$${\tau }_{\mathrm{a},\mathrm{c}}$$

$${\tau }_{IGM}(\lambda )={\sigma }_0{\int }_{R_{vir}}^{\infty }{n}_{\mathrm{H}}(s){\chi }_{{\mathrm{H}}_{\mathrm{I}}}(s)\varphi (\lambda ,V_{\mathrm{shift}},T_{\mathrm{gas}})ds,$$

$$n_{\mathrm{H}}(\overrightarrow{x})=x_{\mathrm{H}}\frac{{\Omega }_{\mathrm{b}}}{{\Omega }_{\mathrm{DM}}}{\rho }_{\mathrm{DM}}(\overrightarrow{x})/m_{\mathrm{H}},$$

$${\delta }_a(\overrightarrow{x})=\frac{{\rho }_a(\overrightarrow{x})-⟨{\rho }_a⟩}{⟨{\rho }_a⟩},$$

$${\delta }_{\mathrm{halos}}(\overrightarrow{x})=b(M){\delta }_{\mathrm{DM}}(\overrightarrow{x}),$$

$${\rho }_{\mathrm{DM}}(\overrightarrow{x})=\left(\frac{{\delta }_{\mathrm{halos}}(\overrightarrow{x})}{b_{\mathrm{eff}}(\overrightarrow{x})}+1\right)⟨{\rho }_{\mathrm{DM}}⟩,$$

$$b_{\mathrm{eff}}(\overrightarrow{x})=\frac{\int b(M)\frac{dN}{dM}(\overrightarrow{x})dM}{\int \frac{dN}{dM}(\overrightarrow{x})dM},$$

$$\Gamma (\overrightarrow{x})=\sum _i\frac{{\sigma }_0{\dot{Q}}_{\mathrm{H},i}{f}_{\mathrm{esc}}^{\mathrm{ion}}}{4\pi |\overrightarrow{x}-{\overrightarrow{x}}_i{|}^2}\frac{\beta }{\beta -3}G(|\overrightarrow{x}-{\overrightarrow{x}}_i|),$$

$$G(|\overrightarrow{x}-\overrightarrow{x}_{i}|)=\left\{\begin{array}[]{ll}0&\mathrm{if\ }|\overrightarrow{x}-\overrightarrow{x}_{i}|>R_{ion,i}\\ 1&\mathrm{if\ }|\overrightarrow{x}-\overrightarrow{x}_{i}|<R_{ion,i}\end{array}\right.,$$

$$R_{ion,i}=K_a({\dot{Q}}_{H,i}/10^{55}{s}^{-1}{)}^{1/3},$$

$$X_{\rm H_{I}}(\overrightarrow{x})=1+a(\overrightarrow{x})/2-\big{[}a+(a(\overrightarrow{x})/2)^{2}\big{]}^{1/2},$$

$$\frac{T_{\mathrm{gas}}(\overrightarrow{x})}{T_0}={\left(\frac{{\rho }_{\mathrm{DM}}(\overrightarrow{x})}{⟨{\rho }_{\mathrm{DM}}⟩}\right)}^{\gamma -1}$$

$$\begin{array}{rl}{\tau }_{\mathrm{IGM}}(\lambda )={\sigma }_0{L}_{\mathrm{cells}}\sum _{l=l_{\mathrm{gal}}}& [n_{{\mathrm{H}}_{\mathrm{I}}}(i,j,l)\times X_{{\mathrm{H}}_{\mathrm{I}}}(i,j,l)\times \\ & \times \varphi (\lambda ,V_{\mathrm{shift}}(i,j,l),T_{\mathrm{gas}}(i,j,l))],\end{array}$$

$$V_{\mathrm{shift}}(i,j,l)=V_{\mathrm{H}}(i,j,l)+H_{\mathrm{Hubble}}(z)\times (l_{\mathrm{gal}}-l)\times L_{\mathrm{cells}},$$

$${\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{c}}={\mathrm{L}}_{\mathrm{Ly}\alpha }^{0,\mathrm{c}}{f}_{\mathrm{esc}}^{\mathrm{c}}.$$

$$f_{\mathrm{esc}}^{\mathrm{IGM},\mathrm{c}}=\frac{\int {\mathrm{\Phi }}^{\mathrm{c}}(\lambda )\times {\mathrm{T}}^{\mathrm{gal}}(\lambda )d\lambda }{\int {\mathrm{\Phi }}^{\mathrm{c}}(\lambda )d\lambda },$$

$${\mathrm{L}}_{\mathrm{Ly}\alpha }={\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}\times f_{\mathrm{esc}}^{\mathrm{IGM},\mathrm{Disk}}+{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}}\times f_{\mathrm{esc}}^{\mathrm{IGM},\mathrm{Bulge}}.$$

$$\mathrm{EW}={\mathrm{L}}_{\mathrm{Ly}\alpha }/{\mathrm{L}}_{\mathrm{continuum}},$$

$${\mathrm{L}}_{\mathrm{Ly}\alpha }={\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}+{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}},$$

$$\mathrm{Z}=\frac{{\mathrm{Z}}^{\mathrm{Disk}}{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}+{\mathrm{Z}}^{\mathrm{Bulge}}{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}}}{{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}+{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}}}.$$

$$f_{\mathrm{esc}}^{\mathrm{IGM}}=\frac{{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}\times f_{\mathrm{esc}}^{\mathrm{IGM},\mathrm{Disk}}+{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}}\times f_{\mathrm{esc}}^{\mathrm{IGM},\mathrm{Bulge}}}{{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Disk}}+{\mathrm{L}}_{\mathrm{Ly}\alpha }^{\mathrm{Bulge}}},$$

$$b(r)={\left[\frac{\xi (r)}{{\xi }_m(r)}\right]}^{1/2},$$

$${\delta }_g=b{\delta }_m\times [1+{\tilde{\alpha }}_1{\delta }_m+{\tilde{\alpha }}_2\frac{1}{aH}{\partial }_z{\mathrm{V}}_z+$$

$$+{\tilde{\alpha }}_3\frac{1}{aH}\left({\partial }_x{\mathrm{V}}_x+{\partial }_y{\mathrm{V}}_y\right)+$$

$$+{\tilde{\alpha }}_4\frac{1}{aH}\frac{{\mathrm{V}}_z}{r_{\mathrm{H}}}+{\tilde{\alpha }}_5{r}_{\mathrm{H}}{\partial }_z{\delta }_m],$$

$$P_0(k)=\left[{\gamma }_1^2+\frac{2}{3}{\gamma }_1{\gamma }_2+\frac{1}{5}{\gamma }_2^2+\frac{1}{3}{\gamma }_3^2\right]b^2{P}_m(k),$$

$${\gamma }_1=\left(1+\frac{{\tilde{\alpha }}_1-{\tilde{\alpha }}_{3}f}{b}\right),$$

$${\gamma }_2=\left({\tilde{\alpha }}_3-{\tilde{\alpha }}_2\right)\beta ,$$

$${\gamma }_3=\left({\tilde{\alpha }}_4\beta \frac{1}{k{r}_{\mathrm{H}}}+{\tilde{\alpha }}_5\frac{k{r}_{\mathrm{H}}}{b}\right),$$

Full text

††pubyear: 2018††pagerange: $\mathrm{Ly}\alpha$ emitters in a cosmological volume II: the impact of the intergalactic medium–A

1 Introduction

Galaxies exhibiting strong Lyman-$\alpha$ emission, the so-called LAEs, are one of the most important tracers of the high redshift Universe in modern astrophysics. After their first detection about 20 years ago (e.g. Steidel et al., 1996; Hu et al., 1998; Rhoads et al., 2000; Malhotra & Rhoads, 2002) studies have focused on finding LAEs at low (Orlitová et al., 2018; Henry et al., 2018) and high redshift (Ouchi et al., 2008; Oyarzún et al., 2017; Matthee et al., 2017; Caruana et al., 2018) including at the epoch of reionization (Sobral et al., 2015; Ouchi et al., 2018; Shibuya et al., 2018).

In the coming years, cosmological surveys such as HETDEX (Hill et al., 2008a) and J-PAS (Benitez et al., 2014) aim to unveil the mystery of the dark energy. These surveys will scan the sky chasing LAEs to trace the underlying dark matter density fluctuations and make clustering measurements. Additionally, these surveys will contribute notably to our knowledge of galaxy formation and evolution. Thus, it is becoming timely to understand properly the selection function of LAEs and how it might affect the apparent spatial distribution and galaxy properties.

The large cross-section of neutral hydrogen atoms around $\mathrm{Ly}\alpha$ wavelengths makes these photons suffer multiple scattering events that modify their frequency and direction. These frequency changes produce the characteristic $\mathrm{Ly}\alpha$ line profiles widely studied (Verhamme et al., 2006; Gronke et al., 2016). The typical distance covered by $\mathrm{Ly}\alpha$ photons inside atomic hydrogen is increased drastically due to multiple scattering events, making these photons very sensitive to dust absorption.

In order to understand the complex escape of $\mathrm{Ly}\alpha$ photons consider the illustration shown in Fig.1. The journey of $\mathrm{Ly}\alpha$ photons begins when a star formation episode takes place inside a galaxy. During these events, hot massive stars (mostly O-type and B-type) emit high energy photons capable of ionizing neutral hydrogen. A fraction of these photons dissociates the H i in the surroundings of the star forming region. Another part is absorbed by dust grains, while the rest escapes from the galaxy and ionizes H i in the intergalactic medium.

Then, the free electrons within the ionized region surrounding the young stars recombine with H ii ions in excited energy levels, causing a cascade until electrons reach the ground energy level. $\mathrm{Ly}\alpha$ photons are emitted when an electron decays from the first excited level to the ground level, an event that occurs with probability $\sim 2/3$ (Spitzer, 1978)) per ionizing photon.

$\mathrm{Ly}\alpha$ photons have to get through the intricate interstellar medium (ISM) before escaping the galaxy (Fig. 1 top left panel)(Neufeld, 1991). The ISM morphology includes dusty gas rich regions such as bars (Spinoso et al., 2017), arms (Kormendy, 2013), H ii bubbles, outflows (Cazzoli et al., 2016) and other structures that complicate the radiative transfer of $\mathrm{Ly}\alpha$ photons. Resonant scattering inside the ISM enhances dust absorption, thus only a fraction of the emitted $\mathrm{Ly}\alpha$ photons manages to escape from the galaxy. Additionally, in this process the $\mathrm{Ly}\alpha$ line profile is modified due to consecutive H i interactions (Harrington, 1973). The final $\mathrm{Ly}\alpha$ escape fraction and line profiles depend strongly on the ISM topology and kinematics (Verhamme et al., 2006; Gurung-Lopez et al., 2018a).

After emerging from galaxies (Fig.1 zoom panel), $\mathrm{Ly}\alpha$ photons enter the IGM and interact with the H i atoms within, producing further $\mathrm{Ly}\alpha$ scattering events (Fig.1 main panel). While inside galaxies the flux reduction is due to dust absorption, in the IGM the $\mathrm{Ly}\alpha$ photons are scattered out of the line of sight, as illustrated in the bottom panel of Fig.1. The Hubble flow redshifts the emitted photons, causing further IGM absorptions at wavelengths bluewards of $\mathrm{Ly}\alpha$  while redder photons travel freely.

The RT processes, occurring in both the ISM and IGM, complicate the selection function of galaxy surveys using $\mathrm{Ly}\alpha$ line detection as a tracer. Galaxies selected with this technique tend to have low metallicity and high specific star formation rates (Sobral et al., 2018a). Observational evidence also suggests that this galaxy population lies in low density environments (Shimakawa et al., 2017). Further observations are needed to cast light on this matter.

The RT in the ISM has been modeled and explored using Monte Carlo Radiative Transfer codes (Zheng & Miralda-Escudé, 2002; Ahn, 2003; Verhamme et al., 2006; Orsi et al., 2012; Gronke et al., 2016; Gurung-Lopez et al., 2018a). These tools generate photons in H i structures and track the subsequent interactions, changes in direction, frequency and possible absorptions. Monte Carlo Radiative Transfer codes have been implemented in cosmological simulations in the $\Lambda$CDM scenario to understand the effect of the ISM on the selection function of LAEs (Orsi et al., 2012; Garel et al., 2012; Gurung López et al., 2018b).

Radiative transfer inside the IGM has been implemented in several ways. For example, Dijkstra et al. (2007) made use of analytic expressions to determine the velocity field, density field and ionization state of the IGM around galaxies as a function of some galaxy properties like the host halo mass or the circular velocity. Meanwhile, Laursen et al. (2011) studied the IGM transmission around the $\mathrm{Ly}\alpha$ wavelength in a hydrodynamic simulation tracking $\mathrm{Ly}\alpha$ rays along different lines of sight. A different approach was taken by Zheng et al. (2010), where they run a Monte Carlo radiative transfer code (similar to the algorithms used to model the ISM) to study the observability of LAEs as a function of the IGM large scale properties. In particular, Zheng et al. (2010) found that the resonant nature of $\mathrm{Ly}\alpha$ might introduce new clustering features in galaxy samples selected by $\mathrm{Ly}\alpha$ detection. However, using a higher resolution simulation, Behrens et al. (2017) claimed to find only a marginal coupling between LAEs and the IGM.

If the IGM impacts significantly on the observed spatial distribution of LAEs, this could introduce dramatic biases (e.g. Wyithe & Dijkstra, 2011) into the cosmological interpretation of clustering data from surveys like HETDEX (Hill et al., 2008b), and future space mission concepts, such as the ATLAS Probe (Wang et al., 2018) and the Cosmic Dawn Intensity Mapper (Cooray et al., 2016). Hence, understanding the role of the IGM in shaping the observed properties of high redshift LAEs is of crucial importance.

This is the second in a series of papers that tackle the selection effects on LAE caused by the RT of $\mathrm{Ly}\alpha$ resonant scattering nature. In our first paper (Gurung López et al., 2018b) we focused in the RT inside the ISM and how it determines the properties of galaxies observed as LAEs. Here, we expand our model and include the RT in the nearby IGM. We study the coupling between $\mathrm{Ly}\alpha$ observability and different IGM large scale properties and how this modifies the clustering of LAEs.

This work is structured as follows : In §2 we present our model and its calibration (§3). Then, in §4 we briefly study the selection function of galactic properties while in §5 we focus on the selection effects of the IGM and its impact on LAE clustering (§6). Finally, we compare our work with the literature (§7) and present our conclusions (§8).

2 A comprehensive model for LAEs.

The model presented here is a follow up of the work presented in Gurung López et al. (2018b). We combine a wide range of physical scenarios in order to assemble a realistic LAE model. Our model is built upon four main pillars:

The P-Millennium N-Body simulation (Baugh et al., 2019), a state-of-the-art dark matter $N$-body simulation with box size $\rm(542.16cMpc/h)^{3}$ with $5040^{3}$ dark matter particles of mass ${\rm M_{p}=1.061\times 10^{8}M_{\odot}}h^{-1}$ and Planck cosmology : $\rm H_{0}=67.77\;km\;s^{-1}Mpc^{-1}$, $\Omega_{\Lambda}=0.693$, $\Omega_{M}=0.307$ , $\sigma_{8}=0.8288$ (Planck Collaboration et al., 2016). The P-Millennium models the hierarchical growth of structures in the $\rm\Lambda CDM$ scenario in our model. 2. 2.

GALFORM, a semi-analytic model of galaxy formation and evolution. The GALFORM version used in this work is detailed in Lacey et al. (2016) and Baugh et al. (2019) and it is based on Cole et al. (2000).

In short, initially GALFORM populates the dark matter halos extracted from a high redshift output of the P-Millennium with gas. Then, the gas is evolved tracking the merger histories of halos and several physical mechanisms are included, such as i) shock-heating and radiative cooling of gas inside halos; ii) formation of galactic disk with quiescent star formation ; iii) the triggering of starburst episodes due to disk instabilities and mergers in bulges; iv) AGN, supernovae and photoionization feedback to regulate the star formation rate; v) the chemical evolution of gas and stars. Additionally, the GALFORM version used in this work implements different stellar initial mass functions (IMFs) for quiescent star formation and starburst episodes (for further details in Lacey et al. (2016)). To ensure a proper resolution, in our model we consider only galaxies with stellar masses higher than $10^{7}{\rm M_{\odot}}h^{-1}$, which roughly corresponds with dark halo masses around $10^{10}{\rm M_{\odot}}h^{-1}$ ($\sim 100$ dark matter particles). 3. 3.

FLaREON111 https://github.com/sidgurun/FLaREON (Gurung-Lopez et al., 2018a), an open source python code to predict Lyman $\alpha$ escape fractions and line profiles in minimal computational time. FLaREON is based on the radiative transfer Monte Carlo code LyaRT222 https://github.com/aaorsi/LyaRT (Orsi et al., 2012) that fully tracks the trajectory of $\mathrm{Ly}\alpha$ photons in outflows with different gas geometries, hydrogen column densities ($\rm N_{H}$), macroscopic expansion velocities ($\rm V_{exp}$) and dust optical depths ($\rm\tau_{a}$). Briefly, FLaREON combines precomputed grids of $f_{\mathrm{e}sc}^{\rm Ly\alpha}$ and $\mathrm{Ly}\alpha$ line profiles in the $\rm N_{H}$-$\rm V_{exp}$-$\rm\tau_{a}$ space for several outflow geometries and different algorithms, such as multi-dimensional interpolation and machine learning, achieving high accuracy at a very low computational cost. We make use of FLaREON to include the RT inside galaxies in our model (see §2.1). 4. 4.

Radiative transfer of $\mathrm{Ly}\alpha$ photons in the intergalactic medium. We estimate the IGM transmission for every galaxy depending on the local environment properties, such as the density, velocity and ionization state of the IGM (see §2.2).

In reality there is an important interface between the ISM and the IGM: the circumgalactic medium (CGM). In principle, the CGM is also populated by neutral hydrogen, implying that $\mathrm{Ly}\alpha$ photons are prone to continue experiencing the resonant scatter. Some observational studies have shown that LAEs present extended $\mathrm{Ly}\alpha$ halos around them (e.g. Leclercq et al., 2017). Also, theoretical works have studied the impact of the CGM in the RT of $\mathrm{Ly}\alpha$ photons (e.g. Zheng et al., 2011; Wyithe & Dijkstra, 2011; Behrens et al., 2017). In this work, our simulation lacks the enough spatial resolution to include properly a CGM contribution to our analysis. However, we expect that the changes introduced by the CGM in our model would be small. In fact, the CGM could modify slightly the line frequency distribution of $\mathrm{Ly}\alpha$ photons entering into the IGM, which in principle could increase the IGM-LAE coupling found in this work (see §5.2) if the CGM is infalling into the galaxy or it could decrease it if the CGM was being ejected from the galaxy.

In the following subsections we describe in detail the design of our LAE models.

2.1 Modeling the radiative transfer of ${\bf Ly\alpha}$ photons inside galaxies.

The physics of $\mathrm{Ly}\alpha$ photons escaping galaxies through galactic outflows are implemented with FLaREON. We focus on two outflow geometries: i) expanding homogeneous Thin Shell (e.g. Verhamme et al., 2006; Orsi et al., 2012; Gurung López et al., 2018b); ii) expanding Galactic Wind (Orsi et al., 2012; Gurung López et al., 2018b) with a density gradient. Both geometries exhibit an empty cavity in the center of the geometry, where monochromatic $\mathrm{Ly}\alpha$ photons are generated. We assume a constant temperature $T=10^{4}K$. These outflow geometries are detailed in Gurung-Lopez et al. (2018a).

GALFORM galaxies are divided into two components, disk and bulge, with distinct galactic properties, such as metallicity, cold gas mass or star formation rate. Therefore, each component is assigned a unique $f_{\mathrm{e}sc}^{\rm Ly\alpha}$ and $\mathrm{Ly}\alpha$ line profile. This assumes that photons generated in a certain galaxy component only interact with that galaxy component.

We calculate outflow properties as in Gurung López et al. (2018b). In particular, motivated by observational studies (e.g. Cazzoli et al., 2016), the outflow expansion velocity is computed as

[TABLE]

where the subscript $c$ denotes the galaxy component (disk or bulge), $\rm SFR_{c}$ is the star formation rate in ${\rm M_{\odot}Gyr^{-1}}h^{-1}$ units, $\rm r_{c}$ is the half stellar mass radius in ${\rm pMpc\;}h^{-1}$, $\rm M_{*}$ is the total stellar mas of the galaxy in ${\rm M_{\odot}}h^{-1}$ units. Additionally, $\rm\kappa_{V,c}$ are free dimensionless parameters regulating the efficiency of gas ejection.

In FLaREON, the Thin Shell and the Galactic Wind geometries present different density profiles, thus different column densities. The outflow neutral hydrogen column density of each component is computed as

[TABLE]

where $\rm M_{cold,c}$ is the cold gas mass in ${\rm M_{\odot}}h^{-1}$ and $\rm\kappa_{N,c}$ are dimensionless free parameters. These are calibrated later in §3. Note that in the Galactic Wind we have used $M_{\rm cold}$ as proxy for the ejected cold gas mass. Therefore, $\kappa_{N,c}$ absorbs the relation between these two galaxy properties and has dimension of $\rm[Time]^{-1}$ .

The free parameters presented above outline the outflow properties, regulating the $\mathrm{Ly}\alpha$ luminosity distribution. In general, the parameters related to the quiescent star formation (galactic disks) give form to the faint end of the $\mathrm{Ly}\alpha$ luminosity function. Meanwhile, the bright end is shaped by the free parameters regulating the starburst episodes (galactic bulges) (Lacey et al., 2016).

Finally, the dust absorption optical depth is simply computed as

[TABLE]

where $\rm E_{\odot}=1.77\times 10^{-21}{\rm cm^{-2}}$ is the ratio $\rm\tau_{a}/N_{H}$ for solar metallicity, $A_{\rm Ly\alpha}=0.39$ is the albedo at the Ly$\alpha$ wavelength, the solar metallicity is $\rm Z_{\odot}=0.02$ (Granato et al., 2000) and $\rm Z_{c}$ is the metallicity of the cold gas in $\rm Z_{\odot}$ units.

2.2 Modeling the radiative transfer of $\mathrm{Ly}\alpha$ photons inside the IGM.

While inside galaxies the losses of $\mathrm{Ly}\alpha$ flux are due to dust absorption, in the IGM photons are scattered out of the line of sight by the neutral hydrogen. The total opacity of the IGM is given by (Dijkstra et al., 2007)

[TABLE]

where $\lambda$ is the wavelength, $s$ is the proper distance to the galaxy where the photon is emitted, $n_{\rm H}$ is the IGM hydrogen number density, $\chi_{\rm H_{I}}$ is the fraction of neutral hydrogen in the IGM and $\phi(\lambda,V_{\rm shift},T_{\rm gas})$ is the Voigt profile Doppler-shifted by the velocity between the emitting galaxy and the IGM at the temperature $T_{\rm gas}$.

We compute the IGM transmission of each galaxy, which depends on the local environment. To do so we compute $n_{\rm H}(\overrightarrow{x})$ , $\chi_{\rm H_{I}}(\overrightarrow{x})$ and the hydrogen velocity field, $V_{\rm H}(\overrightarrow{x})$ from our simulation. Due to disk storage limitations, the dark matter particles of P-Millennium at the snapshots used in this work were not saved. Hence, these quantities are computed from the halo catalogs as shown in the next subsections.

2.2.1 Hydrogen number density field.

We assume that density of hydrogen, $\rho_{\rm H}(\overrightarrow{x})$, is coupled to the dark matter density, $\rho_{\rm DM}(\overrightarrow{x})$. Thus,

[TABLE]

where $x_{\rm H}=0.74$ is the hydrogen fraction of baryonic matter in the Universe, $\Omega_{\rm b}$ and $\Omega_{\rm DM}$ are the densities of baryons and dark matter respectively and $m_{\rm H}$ is the hydrogen mass.

In order to compute $\rho_{\rm DM}(\overrightarrow{x})$, it is useful to define the overdensity field of a given quantity, $\delta_{a}(\overrightarrow{x})$, as

[TABLE]

where $\rho_{a}(\overrightarrow{x})$ is the density field of some quantity and $\langle\rho_{a}\rangle$ is its average. Additionally, the definition of bias for a given dark matter halo mass, $b(M)$, is

[TABLE]

where $\delta_{\rm halos}(\overrightarrow{x})$ is the overdensity field of the dark matter halos and $\delta_{\rm DM}(\overrightarrow{x})$ is the dark matter overdensity field.

By combining eq. 6 and eq. 7, $\rho_{\rm DM}(\overrightarrow{x})$ can be expressed as

[TABLE]

where we have defined the effective bias $b_{\rm eff}(\overrightarrow{x})$ in each cell as

[TABLE]

where ${{dN}\over{dM}}(\overrightarrow{x})$ is the halo mass distribution in each cell.

2.2.2 Velocity field.

We assume that the motion of the dark matter and gas are the same, i.e., $V_{\rm H}(\overrightarrow{x})=V_{\rm DM}(\overrightarrow{x})$. In practice, we divide our simulation box into smaller volumes and compute $V_{\rm DM}(\overrightarrow{x})$ from the halo catalog simply as the median of the velocity distribution of halos within the subvolume.

In Fig.2 we show the matter velocity along the line of sight (chosen arbitrarily as the $Z$ coordinate of the simulation box) at redshift 3.0 in a slice of width 2${\rm cMpc}h^{-1}$ of our simulation box. The volume is divided in big chunks with coherent positive motion along the line of sight (red) and negative motion (blue). The typical scales of these areas are hundreds of comoving megaparsecs. Meanwhile, the transition between these regions is relatively small. This causes great contrasts of velocity along the line of sight on scales below $\sim 10{\rm cMpc}h^{-1}$. We have checked that this behaviour is also found at the other redshifts studied in this work.

2.2.3 Fraction of neutral hydrogen field.

Star forming galaxies and active galactic nuclei (AGN) are the main sources of ionizing photons (photons wavelength shorter than 912Å) in the Universe (Kimm & Cen, 2014). Thus the photoionization field, $\rm\Gamma(\overrightarrow{x})$, is coupled to the number density field of these objects. In principle, regions with a high rate of ionizing photons, i.e. close to these sources, will be more ionized, thus lowering the IGM opacity to Ly$\alpha$ photons.

In order to compute $\Gamma(\overrightarrow{x})$ we assume that the ionizing radiation is produced only by galaxies. Whether or not AGNs and QSO contribute significantly to $\Gamma(\overrightarrow{x})$ at high redshift is still an open debate. In fact, some works in the literature (e.g. Parsa et al., 2018) have shown that AGNs or QSOs, are not sufficiently abundant to significantly contribute to $\Gamma(\overrightarrow{x})$, meanwhile, other works in the literature suggest the opposite (e.g. Giallongo et al., 2019) . However, we do not expect this to significantly affect our result since, as described later, $\Gamma(\overrightarrow{x})$ is calibrated to reproduce the observed mean $\Gamma$ at the different epochs studied here.

We compute $\Gamma(\overrightarrow{x})$ as the superposition of every ionizing field generated by each galaxy, i.e,

[TABLE]

where the sum is over all the galaxies in the box, $\overrightarrow{x}_{i}$ is the location of each galaxy, $\overrightarrow{x}$ is the position where ionizing field in evaluated. Additionally, we assume a global escape fraction of ionizing photons $f^{\rm ion}_{\rm esc}=0.1$ (Kimm & Cen, 2014) and $\dot{Q}_{{\rm H},i}$ is the total luminosity of ionizing photons given by GALFORM for each galaxy. We assume that the SED of the galaxy takes the form $J(\nu)=\nu^{\beta}$ in bluer parts of the hydrogen ionization frequency threshold, $\nu_{\rm H}$. We also assume that the photoionization cross section as $\sigma_{\rm H}(\nu)=\sigma_{0}(\nu/\nu_{\rm H})^{-3}$ with $\sigma_{0}=6.3\times 10^{-18}cm^{-2}$ . Finally, the function $G(|\overrightarrow{x}-\overrightarrow{x}_{i}|)$ takes into account the fact that photons emitted by a single galaxy do not reach every point in space and is given by

[TABLE]

where $R_{ion,i}$ is the radius of the sphere centered in the location of the galaxy $i$. We use a similar expression to the Strömgren radius to compute $R_{ion,i}:$

[TABLE]

where $K_{a}$ is a free parameter. Increasing (decreasing) $K_{a}$ leads to greater (lower) $R_{ion,i}$. Hence, at each point, more (less) galaxies contribute to the the ionizing radiation field, which augments (lowers) the $\braket{\Gamma}$ of the simulation box.

We determine the $K_{a}$ values at each redshift by fitting the $\braket{\Gamma}$ of our simulations to the $\braket{\Gamma}$ given by Haardt & Madau (2012). The best fitting values are listed in Table 1. We find that $K_{a}$ decreases with redshift. Since the number density of ionizing sources decreases from $z=2.2$ to 5.7, the volume ionized by a single galaxy must be increased with redshift in order to fit observations.

Fig.3 displays a comparison between the observed $\braket{\Gamma}$ and the $\Gamma(\overrightarrow{x})$ distribution (percentiles 2.5, 16, 50, 84 and 97.5) of our model. The right panel shows the spatial variations of $\Gamma(\overrightarrow{x})$ over cosmic time in the same slice of width $\sim 2{\rm cMpc}h^{-1}$. By construction, $\Gamma(\overrightarrow{x})$ reproduces Haardt & Madau (2012) quite well.

Additionally, in Fig.4 it is shown the PDF (left) and spatial distribution of the photoionization rate of the IGM. At $z=2.2$ $\braket{\Gamma}$ peaks and decreases towards higher redshifts. The dispersion of $\Gamma(\overrightarrow{x})$ evolves with redshift too. While at high redshift the photoionization field exhibits a complex structure with high contrasts (broader PDF), as the Universe evolves, it becomes smoother (tighter PDF).

Finally, we compute the neutral hydrogen fraction field as in Dijkstra & Loeb (2009) :

[TABLE]

with $a(\overrightarrow{x})=\Gamma(\overrightarrow{x})/n_{\rm H}(\overrightarrow{x})\alpha_{\rm rec}(\overrightarrow{x})$ where the case-A333Case-A recombination assumes that the medium is optically thin to ionizing photons. recombination coefficient is given by $\alpha_{\rm rec}(\overrightarrow{x})=4.2\times 10^{-13}(T_{\rm gas}(\overrightarrow{x})/10^{4}K)^{0.7}cm^{3}s^{-1}$. We assume that the IGM temperature behave as

[TABLE]

where $T_{0}=2\times 10^{4}K$ and $\gamma-1=0.45$ (Rudie et al., 2012). We have tested that our results do not depend on our choice of $T_{0}$ by repeating our analysis for several values of $T_{0}$ spawning between $10^{2}K$ and $10^{5}K$. Additionally, we have implemented different models for $T_{\rm gas}(\overrightarrow{x})$ and we find that our results depend very weakly (or not at all) on how $T_{\rm gas}(\overrightarrow{x})$ is computed. In fact, we have also made the whole analysis presented in this work assuming $T_{\rm gas}(\overrightarrow{x})=T_{0}$, finding very similar results, even for the different values of $T_{0}$ discussed above.

2.2.4 IGM transmission.

We compute $n_{H}(\overrightarrow{x})$, $V_{H}(\overrightarrow{x})$ and $\chi_{H_{I}}(\overrightarrow{x})$ on a grid of $250^{3}$ cells with side ${\rm\sim 2cMpc}\;h^{-1}$. Then, by linear interpolation, we reevaluate these fields on a grid with higher spatial resolution along an arbitrary direction chosen as the line of sight. In particular, the new grid is composed of 250$\times$250$\times$3000 cells, where the size of the cells along the line of sight is $\rm L_{cells}\sim 0.2cMpc/h$. We choose this grid size to ensure good signal to noise in the computed field and the IGM transmission curves with high enough resolution in frequency space.

Following eq. 4, the IGM optical depth in each cell in the simulation box rest frame is given by

[TABLE]

where the cell indices {$i$, $j$} and {$l$} are perpendicular and parallel, respectively, to the line of sight. Moreover, $l$ starts at the cell where the galaxy lies, $l_{\rm gal}$, and iterates towards the observer along the line of sight direction. Additionally, $V_{\rm shift}$ is the relative velocity along the line of sight between the cell containing the galaxy ($i$,$j$,$l=l_{\rm gal}$) and the IGM in the iterated cell ($i,j,l$), given by

[TABLE]

where $H_{\rm Hubble}(z)$ is the Hubble constant at redshift $z$.

Finally, the IGM transmission at each position is computed as $T(\lambda)=e^{-\tau_{\rm IGM}(\lambda)}$. In Fig. 5 we show the median and dispersion of the IGM transmission at different redshifts for a sample of the 150000 GALFORM galaxies with the highest specific star formation rate (sSFR) as representative of an emission line galaxy population. In general, we find a good qualitative agreement between the shape of our mean transmission curves and the mean transmission curves obtained by Laursen et al. (2011) computed from hydrodynamical simulation with RT physics implemented. The IGM absorbs photons bluer than $\lesssim 1216$Å. Moreover, as galaxies lie in overdense regions, the IGM opacity is higher close to the galaxy, causing the drop in the transmission close to $\mathrm{Ly}\alpha$ wavelength. Then the IGM transmission flattens to the IGM cosmic transmission.

Additionally, the optical depth of the IGM evolves with redshift, producing higher transmissions at lower redshifts. This becomes dramatic at $z=5.7$, where the IGM transmission goes below 1% at bluer frequencies than $\mathrm{Ly}\alpha$.

The large dispersion of the IGM transmission reflects the complex variety of environments surrounding galaxies. In order to test this idea, we rank our LAE samples in IGM density $\rho_{\rm H}$ and split them into 3 subsamples: underdense (below the percentile 33rd of density), intermediate (between the 33 and 66 percentiles) and overdense (above the 66 percentile). In Fig. 6 we show the median IGM transmission at $z=3.0$ for these subsamples. In general the behaviour of the three population is the same. At redder wavelengths than $\mathrm{Ly}\alpha$  the transmission is 1. Meanwhile, at bluer wavelengths far from $\mathrm{Ly}\alpha$ ($\rm\lambda\sim 1210$Å) the transmission converges to the mean IGM transmission. However, at blue wavelengths around $\mathrm{Ly}\alpha$ ($\rm\lambda\sim 1214$Å) the density has a great impact on the IGM transmission. We find that LAEs hosted in denser environments exhibit lower IGM transmission than their counterparts in low density regions. The transmission is $\sim 0.4$, $\sim 0.2$ and $\sim 0.1$ for the underdense, intermediate and overdense environments. For completeness, we also did this analysis at the other snapshots of our simulations, finding the same trend. The typical transmission at $z=2.2$ around $\mathrm{Ly}\alpha$ is $\sim 0.9$, $\sim 0.85$ and $\sim 0.8$ for the underdense, intermediate and overdense environments. Meanwhile, at $z=5.7$ the transmission remains below 1% even in the underdense regions.

2.3 The observed $\mathrm{Ly}\alpha$ luminosity.

The observed $\mathrm{Ly}\alpha$ luminosity is a convolution of galactic and IGM properties. We compute the observed $\mathrm{Ly}\alpha$ luminosity as follows.

Physical processes taking place inside galaxies are implemented as in (Gurung López et al., 2018b). For each galactic component: i) the intrinsic $\mathrm{Ly}\alpha$ luminosity, ${\rm L_{Ly\alpha}^{0,c}}$ (where the superindex "$c$" denotes the disk or bulge component), is predicted by GALFORM from the instantaneous star formation rate; ii) the outflow properties are computed with equations 1, 2 and 3; iii) the outflow $\mathrm{Ly}\alpha$ escape fraction $f_{\rm esc}^{c}$ and line profile $\rm\Phi^{c}(\lambda)$ are predicted with FLaREON; iv) the $\mathrm{Ly}\alpha$ luminosity escaping the ISM is computed as

[TABLE]

The IGM transmission is calculated in the simulation rest frame, as discussed in the previous section. However, the galaxy rest frame mismatches, in general, the simulation rest frame due to galaxy peculiar velocities. Hence, the IGM transmission for a given galaxy is Doppler shifted by the peculiar velocity of the galaxy along the line of sight. Then, for each component, the fraction of photons that travel unscattered through the IGM, $f_{esc}^{\rm IGM,c}$, is computed as

[TABLE]

where ${\rm T^{gal}(\lambda)}$ is the IGM transmission at the galaxy position and rest frame and $\rm\Phi^{c}$ is the $\mathrm{Ly}\alpha$ line profile of each galaxy component.

In Fig. 7 we illustrate how $f_{\rm esc}^{\rm IGM,c}$ is computed. In this particular case the IGM and the galaxy are approaching along the line of sight. As shown previously, in the rest frame, the IGM absorbs photons bluer than $\mathrm{Ly}\alpha$. However, due to the relative motion between the galaxy and the IGM along the line of sight, the IGM transmission (blue curve) is redshifted in the galaxy rest frame. Hence, the $\mathrm{Ly}\alpha$ photons emerging from the ISM (black curve) are partially absorbed. As a result the observed $\mathrm{Ly}\alpha$ line profile (red curve) is modified. $f_{\rm esc}^{\rm IGM,c}$ is calculated as the ratio between the integrals of the observed and the emerging line profiles.

Then, for each galaxy, the observed $\mathrm{Ly}\alpha$ luminosity is computed as

[TABLE]

Moreover, we estimate the rest frame $\mathrm{Ly}\alpha$ equivalent width, $\rm EW$, as

[TABLE]

where $\rm L_{continuum}$ is the continuum luminosity per unit of wavelength around $\mathrm{Ly}\alpha$ computed by GALFORM. This quantity is based on the evolution of the composite stellar population of each modeled galaxy. From now on, unless stated otherwise, we define LAEs in our model as galaxies with $\rm EW>20$Å.

3 Lyman-$\alpha$ emitters

3.1 Model calibration.

In these section we briefly discuss how the free parameters described in Eq.1 and 2 are calibrated. For further details we refer the reader to Gurung López et al. (2018b).

In short, for each redshift and outflow geometry, we perform an MCMC analysis using the open source Python library emcee (Foreman-Mackey et al., 2013) to estimate the values of $\rm\kappa_{N,c}$ and $\kappa_{V,c}$ that best reproduce the observed LAE luminosity function (LF) at different redshifts. In particular, at $z=2.2$ we fit the observed LF from Cassata et al. (2011), Konno et al. (2016) and Sobral et al. (2017), while at $z=3.0$ we fit the LFs from Ouchi et al. (2008) and Cassata et al. (2011). Finally, at $z=5.7$ we model the LFs of Ouchi et al. (2008) and Konno et al. (2018). In order to fit the different observed LFs at the same same redshift, we combine them at each redshift by fitting a 5th-order polynomial (in \log$$\rm L_{{\rm Ly}\alpha}-$\log$LF space) taking into account the uncertainties of each data set. The functional form of a 5th-order polynomial was chosen due to the fact that some recent works suggest that observed LFs of LAEs are more complex than the typical Schechter function (Konno et al., 2016; Sobral et al., 2018b). The free parameter values are listed in Table 2. Additionally, the resulting $\rm N_{H}$, $\rm V_{\rm exp}$ and $\tau_{a}$ distributions are discussed in Appendix A.

In Fig. 8 we compare the observed LF at several redshifts with our models. For illustration we show the GALFORM intrinsic LAE LF (black thin line), which over predicts the number of LAEs over the full luminosity range at $z=2.2$ and $z=3.0$. Hence, the total $\mathrm{Ly}\alpha$ escape fraction $f_{\rm esc}^{\rm Ly\alpha}$ (galaxy+IGM) must be $<1$. However at redshift $z=5.7$, while bright intrinsic LAE surpass observations, the intrinsic number counts of faint LAEs resembles observations, implying $f_{\rm esc}^{\rm Ly\alpha}<1$ and $f_{\rm esc}^{\rm Ly\alpha}\sim 1$ respectively.

In general, after calibration, our models (colored solid lines) match the observed LAE LF at all redshifts by construction. The good agreement is remarkable at $z=5.7$. The Thin Shell geometry matches slightly better observations than the Galactic Wind. Additionally, we show the LAE LF of our calibrated model excluding the IGM absorption (colored dashed lines). Although the LAE number counts of the model without IGM exceeds the abundance of LAE in the complete models at every $\mathrm{Ly}\alpha$ luminosity, their LAE LF are very similar. This points to the fact that the RT inside galaxies is the main driver shaping the observed LAE LF.

3.2 LAE samples

Throughout this work we analyze and compare the properties of different LAE samples to highlight the RT selection effects. Here we describe how the samples are built.

3.2.1 Full Ly$\alpha$ emitters samples.

Full $\mathrm{Ly}\alpha$ emitters (FLAE). These samples represent the observed LAE population. They include all the radiative transfer processes explained above (ISM + IGM). The FLAE samples are derived from the full GALFORM population. The outflow properties of the GALFORM galaxies are computed with the calibrated free parameters (listed in Table 2). Then we assign to each galaxy a $\mathrm{Ly}\alpha$ luminosity (as described in §2.3). Note that these samples, by construction, reproduce the observed LF. We rank these populations by $\mathrm{Ly}\alpha$ luminosity and perform a number density cut of $4\times 10^{-3}({\rm cMpc}\;h^{-1})^{-3}$. The chosen number density cut is arbitrary. We obtain similar results for higher and lower number density cuts.

3.2.2 Partial Ly$\alpha$ emitters samples.

Partial $\mathrm{Ly}\alpha$ emitters (PLAE). These samples include galactic RT physics but lack the IGM absorption. They are also subsamples of the full GALFORM galaxy population. PLAEs can be seen as the LAE population that would be observed if the IGM was completely transparent. In particular, $\rm N_{H}$ and $\rm V_{exp}$ are computed with the same calibration as the models with full RT (FLAE samples). However, in contrast to FLAE populations where the $\mathrm{Ly}\alpha$ luminosity is computed through Eq. 19, in the PLAE samples, the observed $\rm L_{{\rm Ly}\alpha}$ is computed simply as

[TABLE]

where the ISM RT is included in $\rm L_{Ly\alpha}^{Disk}$ and $\rm L_{Ly\alpha}^{Bulge}$.

Finally, we rank galaxies by $\rm L_{{\rm Ly}\alpha}$ and make the same number density cut as in the FLAE samples. Note that, by construction, the intrinsic galaxy population and properties are identical for FLAE and PLAE samples. However, the FLAE samples include the IGM selection effects. Therefore, the comparison between these samples sheds light on the impact of the IGM.

3.2.3 Shuffled Ly$\alpha$ emitters samples.

Shuffled $\mathrm{Ly}\alpha$ emitters (SLAE) samples are built to exhibit the same clustering as the FLAE samples if no LAE-IGM coupling is found. Throughout this work we analyze the clustering and the galactic properties of LAEs and how the IGM affects them. Since the FLAE and PLAE are computed from the same galaxy population they are useful to understand how the IGM shapes the galactic properties. As seen later on, the FLAE and PLAE exhibit different mass functions. This causes that each population displays a different clustering, and in particular, different bias. Additionally, if the large scale IGM properties are coupled to the $\mathrm{Ly}\alpha$ observability, further clustering distortions are expected (Zheng et al., 2011). Therefore, in order to study the LAE-IGM coupling it would be desirable to use samples with the same bias.

SLAE samples are derived from the P-Millennium halo catalog mimicking the FLAE halo mass functions. In detail:

We separate our full RT LAE population into centrals and satellites. We find that at all redshifts and gas geometries, central galaxies constitute $\sim 98\%$ of the full RT LAE samples. 2. 2.

Comparing with the dark matter halo catalogs of the simulation we compute the fraction of halos occupied by LAEs as a function of mass (HOD) individually for the central and satellite population. 3. 3.

From the halo catalog we select central halos reproducing the HOD of the central LAE sample using a uniform random distribution. 4. 4.

We determine the number of LAE satellites hosted in each dark matter halo with a random Poisson distribution with mean equal to the satellite’s HOD evaluated at the mass of the halo. Following Jímenez et al. (in prep.), we assign to the satellite the same location as the central halo and restrict our clustering analysis to the 2-halo-term scales. 5. 5.

We combine the new satellite and central population in a single SLAE sample.

By construction, SLAE  are free of the IGM selection effects and exhibit the same halo mass distribution as their full RT progenitors. Therefore, if the IGM is not shaping the LAE spatial distribution, the clustering in the SLAEs and full RT LAE samples should be the same444 Note that, by construction, the SLAE samples do not isolate the IGM-LAE coupling, but it also includes the effects of assembly bias (Contreras et al., 2019). This might cause differences between the clustering of the FLAE and SLAE populations. To check the assembly bias impact on the clustering of our LAE populations we built SLAE samples from the PLAE samples. No significant difference was found between the clustering of these populations, i.e., there was no assembly bias evidence. Therefore, we assume that the difference between the FLAE and SLAE (computed from the FLAE) samples are due to the IGM-LAE coupling. on scales greater than the 1-halo term.

4 LAE galaxy properties.

In this section we briefly study the selection function of LAEs. First, we analyze which galaxies would be observed as LAEs if the IGM was completely transparent to photons around $\mathrm{Ly}\alpha$. For this goal, we contrast the full GALFORM galaxy population with the PLAE samples, which includes RT only in the ISM. In the second case, we characterize the IGM impact by directly comparing the PLAE and the full RT LAE sample. The variations among these samples are caused by the IGM, since the only difference between the PLAE and FLAE samples is that FLAE also include RT in the IGM.

In Fig.9 we compare the galaxy property distributions between our $\mathrm{Ly}\alpha$ flux selected samples and the full galaxy population predicted by GALFORM at different redshifts. We define a $\mathrm{Ly}\alpha$-weighted average gas metallicity for a galaxy as

[TABLE]

We find very little difference between the FLAE and PLAE populations. This indicates that the IGM does not induce significant selection effects on the galaxy properties of LAEs. Moreover, independently of the outflow geometry and cosmic time, galaxies with strong $\mathrm{Ly}\alpha$ emission present moderate stellar mass and SFR and low metallicity. We find that models implementing different outflow geometries behave in a similar fashion. Still there are tiny differences between them.

The typical stellar mass of galaxies observed as LAEs also evolve over cosmic history. In general, we find the same trend in the LAEs samples as in the full galaxy population: LAEs at lower redshift exhibit higher stellar content than their homologous at higher redshift. In detail, the $\rm M_{*}$ distributions peak around moderate masses; $10^{8.5}{\rm M_{\odot}}h^{-1}$, $10^{9}{\rm M_{\odot}}h^{-1}$ and $10^{9.2}{\rm M_{\odot}}h^{-1}$ at redshift 5.7, 3.0 and 2.2 respectively. The shapes of the distributions are very similar between $z=2.2$ and $z=3.0$. The dynamical range of $\rm M_{*}$ at these redshifts is very similar to the observed LAEs (Oyarzún et al., 2017). Additionally, very small differences can be found between the two outflow geometries implemented in this work. At $z=2.2$ and 3.0 the Thin Shell predicts a greater abundance of massive galaxies in comparison with the Galactic Wind. This trend disappears at $z=5.7$, when both stellar mass distribution are almost identical.

We note that at $z=5.7$ the galaxy stellar mass distributions for our LAEs samples truncates abruptly at $\rm M_{*}=10^{7}{\rm M_{\odot}}h^{-1}$. This is caused by the cut in GALFORM at this stellar mass imposed to ensure a good resolution of galaxies. This suggests that a fraction of the LAE population at $z=5.7$ would inhabit galaxies with $\rm M_{*}<10^{7}{\rm M_{\odot}}h^{-1}$. However, as our galaxy population lacks these low mass galaxies, other more massive galaxies are selected as LAEs. Computing the precise number of LAEs that should be hosted by galaxies with $\rm M_{*}<10^{7}{\rm M_{\odot}}h^{-1}$ is challenging. As a simple calculation, we rescale the stellar mass distribution found in our model at $z=3.0$ towards smaller $\rm M_{*}$ such as the peak of the distribution matches the one at $z=5.7$. Then, the fraction of galaxies below the resolution limit is $\sim 0.04$. This hints that only a small part of the LAE population a would lie in galaxies not resolved in our model.

Meanwhile, the SFR of the full galaxy population also evolves with redshift. GALFORM predicts a progressive increase in the SFR from redshift 5.7 to 2.2 (middle column of Fig.9). However, LAE populations exhibit a moderate SFR distribution that remains almost frozen through cosmic time. For all redshifts the SFR peaks at $10^{0.25}{\rm M_{\odot}}h^{-1}yr^{-1}$, well below the galaxies with the highest SFR ($\sim 10^{2.5}{\rm M_{\odot}}h^{-1}yr^{-1}$). We find no significant evolution from $z=2.2$ to 3.0, while at $z=5.7$ the distribution becomes a little bit broader.

The outflow geometry has a significant impact on the SFR distribution of LAEs. In particular, at all redshifts, the LAE samples characterized by the Thin Shell exhibit lower SFR than those using the Galactic Wind. This is caused by several reasons. First, the recipes to link galaxy proprieties to outflow properties are different between both geometries (see Eq.1 and Eq.2). Second, the escape fraction of $\mathrm{Ly}\alpha$ photons depends strongly on the geometry (Gurung López et al., 2018b; Gurung-Lopez et al., 2018a).

GALFORM also predicts an evolution through cosmic time of the metal abundance in galaxies (right column of Fig.9). In fact, galaxies at lower redshifts exhibit higher metallicity. This is a consequence of the consecutive events of SFR that pollute the initially pristine interstellar medium. The typical maximum metallicities are $\rm Z\sim 10^{-0.75}Z_{\odot}$, $\rm\sim 10^{-0.6}Z_{\odot}$ and $\rm\sim 10^{-0.5}Z_{\odot}$ at redshift 5.7, 3. and 2.2 respectively.

The galaxies observed as LAEs exhibit a low metallicity independently of redshift or cosmic time. Additionally, LAEs also show the $Z$ evolution present in the full galaxy population. Their metallicity distributions peak around $\rm\sim 10^{-2.5}Z_{\odot}$ at $z=2.2$ and 3.0, while at $z=5.7$ it peaks at $\rm\sim 10^{-2.75}Z_{\odot}$. There also differences between the Thin Shell and Galactic Wind models at low redshift. At $z=2.2$ and 3.0 the Thin Shell predicts a higher number of LAEs with $\rm Z\sim 10^{-1.25}Z_{\odot}$. Finally, very little differences are found between the LAE samples including RT in the IGM and excluding it.

5 LAE IGM properties.

In this section we study how large scale properties of the IGM affect our full RT LAE samples (FLAE). Later, in §7, we compare our results with previous works.

We analyze the $\mathrm{Ly}\alpha$ IGM escape fraction computed as the ratio of the observed $\mathrm{Ly}\alpha$ emission and the total $\mathrm{Ly}\alpha$ flux escaping galaxies, i.e,

[TABLE]

as a function of the density $\rho$ (Eq. 8), the IGM line of sight (LoS) velocity $\rm V_{z}$, the IGM gradient along the LoS of the velocity $\partial_{z}V_{z}$555 Throughout this paper we use the Einstein notation for partial derivatives, i.e., $\partial_{z}\coloneqq{{\partial}\over{\partial z}}$. and density $\partial_{z}\rho$. The gradients are computed from their respectively fields by a 3-point derivative method.

5.1 The Transmission of the IGM

In Fig.10 and Fig.11 we show the $f_{\rm esc}^{\rm IGM}$ behaviour of the Thin Shell and Galactic Wind models against the IGM properties at different redshifts. Overall, we find that the Thin Shell and Galactic Wind exhibit the same global trends. Since most of the $\mathrm{Ly}\alpha$ photons are redshifted by the ISM, the IGM absorbs only a small fraction of the $\mathrm{Ly}\alpha$ flux that escaped from galaxies. For both geometries, $f_{\rm esc}^{\rm IGM}$ is close to unity. The median $f_{\rm esc}^{\rm IGM}$ increases with the age of the Universe, as the IGM becomes more transparent to $\mathrm{Ly}\alpha$ photons. We find that at $z=2.2$ and 3.0 the IGM absorbs more photons in the Thin Shell model than in the Galactic Wind, while at $z=5.7$ the absorption is comparable. In particular, for the Thin Shell, $f_{\rm esc}^{\rm IGM}\gtrsim 0.99$ , $\sim 0.96$ , $\sim 0.92$ at redshift $2.2$, $3.0$ and $5.7$ respectively. Meanwhile, the median $f_{\rm esc}^{\rm IGM}$ for the Galactic Wind rounds $\sim 0.9999$ , $\sim 0.999$ and $\sim 0.92$ at redshift 2.2, 3.0 and 5.7 respectively.

The differences between the Thin Shell and Galactic Wind are mainly due to their different $\mathrm{Ly}\alpha$ line profiles. Our model is calibrated by fitting the luminosity function of our model with observations. There are two actors converting the intrinsic $\mathrm{Ly}\alpha$ LF into the observed one: the RT in the ISM and in the IGM. As seen in Fig.8 the IGM effect in the LF is small at any redshift, which puts most (but not all) of the weight of the fitting in the ISM component. Roughly speaking, the shape of the observed LF determines the $\rm V_{exp}$ and $\rm N_{H}$ distributions (see Appendix) in our model through equations 1 and 2. Moreover, these distributions determine the properties of the $\mathrm{Ly}\alpha$ line profiles and, therefore, the IGM absorption. Galaxies emitting more flux bluewards of $\mathrm{Ly}\alpha$ would be more obscured than if they were emitting only redwards of $\mathrm{Ly}\alpha$. In practice, we find that the $\rm V_{exp}$ and $\rm N_{H}$ distributions make the Thin Shell to be more coupled to IGM at $z=2.2$ and $3.0$ while the opposite happens at $z=5.7$, where the Galactic Wind is more affected.

The median $f_{\rm esc}^{\rm IGM}$ varies over the dynamical range of the large scale IGM properties studied in this work. We note that the $f_{\rm esc}^{\rm IGM}$ variation is smaller than the dispersion around the median. However, the scatter is not caused by uncertainties, but by the great diversity of combinations of $\rho$, $\rm\partial_{z}\rho$, $\rm V_{z}$ and $\rm\partial_{z}V_{z}$. This points in the same direction as Fig. 5, where the great variety of the IGM transmission curves is shown. These trends are statistically significant, i.e., they are not caused by noise or sample variance in our data set. In general, we find the same trends in the Thin Shell and Galactic Wind models. The strength of the correlations evolves with redshift. As the IGM becomes more transparent the trends become weaker. In fact, the lower the redshift, the weaker are the dependencies on $f_{\rm esc}^{\rm IGM}$. In particular, for the Thin Shell, the typical changes in $f_{\rm esc}^{\rm IGM}$ are $\lesssim 1\%$, $\sim 2\%$ and $\sim 5\%$ at $z=2.2$, 3.0 and 5.7 respectively. Meanwhile, in the Galactic Wind modes, the $f_{\rm esc}^{\rm IGM}$ variations are $\lesssim 0.01\%$, $\lesssim 0.1\%$ and $\sim 10\%$.

We find an anti-correlation between the IGM transmission and the local density. This is more apparent at $z=2.2$ and dilutes at higher redshift. In particular, at $z=5.7$ the median $f_{\rm esc}^{\rm IGM}$ is quite flat and the dispersion becomes greater at higher densities.

In general, $f_{\rm esc}^{\rm IGM}$ correlates with $\partial_{z}\rho$ at all redshifts. In addition to this trend, at redshift $2.2$ and $3.0$, $f_{\rm esc}^{\rm IGM}$ peaks at $\partial_{z}\rho=0$, where the IGM transmission is slightly higher. However this peak is not present at $z=5.7$.

We find a clear correlation between the IGM transmission and the IGM velocity along the line of sight at all redshifts. Galaxies in IGM regions moving towards the observer ($\rm V_{z}<0$) suffer, statistically, greater absorption than galaxies moving away from the observer ($\rm V_{z}>0$). The amplitude of the correlation augments towards higher redshifts, when the IGM becomes more optically thick. In particular, the IGM transmission variation in the $\rm V_{z}$ dynamical range rounds $<1\%$, $\sim 1\%$ and $\sim 5\%$ in the Thin Shell at redshift 2.2, 3.0 and 5.7 respectively.

Finally, at $z=2.2$, $f_{\rm esc}^{\rm IGM}$ correlates with $\rm\partial_{z}V_{z}$. However, at $z=3.0$ no correlation is found. At $z=5.7$ we find a small anti-correlation that produces variations of $\sim 2\%$ in the IGM transmission in the $\rm\partial_{z}V_{z}$ dynamical range.

5.2 The IGM-LAE coupling.

In Zheng et al. (2011) (from now on ZZ11) the authors present a model that features the LAEs coupling with the large scale properties of the IGM. The LAE models presented in this work exhibit similar trends to the relations found by ZZ11. In this section we discuss the origin of the IGM-LAE coupling, that matches ZZ11 interpretation for some IGM properties ($\rho$, $\partial_{z}\rho$, $\rm V_{z}$) while differs in others ($\rm\partial_{z}V_{z}$).

5.2.1 IGM Density.

At low redshift, $f_{\rm esc}^{\rm IGM}$ anticorrelates with the density. This is caused by the column density along the line of sight between the galaxy and the distance where the IGM becomes transparent due to the Hubble flow. In regions of high density the amount of H i that $\mathrm{Ly}\alpha$ photons have to go through is greater than in environments with lower density. This causes a greater number of scattering events and absorptions, in the high density regions. Actually, this can be seen in Fig.6, where underdense regions exhibit a higher transmission than denser regions.

However, at $z=5.7$ the trend reverses, and $f_{\rm esc}^{\rm IGM}$ is higher the greater the IGM density is. It is still true that the H i column density is systematically higher in over dense regions. However, independently of the environment the IGM transmission at wavelengths bluer than $\mathrm{Ly}\alpha$ is below 1%. This causes that the IGM selection effect on the density less important. At the same time, the impact of the IGM increases with redshift, causing strong selection effects on other IGM properties. Therefore, the correlation found between $f_{\rm esc}^{\rm IGM}$ and density is caused by the other IGM properties.

5.2.2 IGM density gradient.

Our models predict a correlation between $f_{\rm esc}^{\rm IGM}$ and $\partial_{z}\rho$. This trend is caused by the difference in H i column density between the galaxy and the distance where the IGM becomes transparent. This scenario is illustrated in the right panel of Fig.12. On one hand, for a source of $\mathrm{Ly}\alpha$ photons at a fixed $\rho$, if the density gradient is positive (galaxy ’c’), then the IGM density decreases as the photons travel through it, allowing more photons to escape. On the other hand, in the case that $\partial_{z}\rho<0$ (galaxy ’a’), the IGM density increases as photons go through the IGM, causing a higher number of lost photons.

Additionally, there are two main different regimes with $\partial_{z}\rho=0$: i) the peak of overdensities (galaxy ’d’ of Fig.12) with low IGM transmission due to the high column density; and ii) the bottom of underdensities (galaxy ’b’) with high $f_{\rm esc}^{\rm IGM}$. The relative occupancy of these regions changes due to the RT in the IGM and ISM. LAE are preferentially observed in underdense environments. This is caused by two main reasons. First, the low IGM transmission in overdense regions. Second, the RT in the ISM prevents very massive galaxies (which would lie in overdense regions) to be observed as LAEs (see Fig.9). Therefore, most of the galaxies observed in environments with $\partial_{z}\rho=0$ are in underdense regions. Meanwhile, the IGM transmission is higher the lower the density is. This causes the transmission peak found at $\partial_{z}\rho=0$ for $z=2.2$, $3.0$.

5.2.3 IGM kinematics

When the relative velocity between the galaxy and the IGM is zero only photons bluewards of $\mathrm{Ly}\alpha$ are absorbed, as seen in Fig. 5. However, if the IGM and the galaxy are moving towards each other, the $\mathrm{Ly}\alpha$ line is blueshifted in the IGM rest frame, causing greater absorption. In the opposite scenario, where the galaxy and the IGM are moving away from each other, the $\mathrm{Ly}\alpha$ photons are redshifted and escape the nearby IGM more easily.

The relative motion between the IGM and galaxies causes selection effects. On average, galaxies moving away from the observer are more likely to be observed as LAE. In Fig. 13 we list the different combinations of velocities between IGM and galaxies along the LoS.

• •

Case 1) Both the galaxy and the IGM are moving away form the observer. However, the galaxy velocity is greater. The IGM only absorbs bluer photon than $\mathrm{Ly}\alpha$.

• •

Case 2) Both the galaxy and the IGM are moving away form the observer. Additionally, the IGM moves faster than the galaxy. In this case, redder wavelengths than $\mathrm{Ly}\alpha$ are absorbed.

• •

Case 3) The galaxy goes away while the IGM approaches the observer. $\mathrm{Ly}\alpha$ photons are redshifted when they reach the IGM. Only bluer frequencies than $\mathrm{Ly}\alpha$ suffer absorption.

• •

Case 4) Both the IGM and the galaxy approaches the observer. In this case the galaxy moves faster and the IGM absorbs redder wavelengths than $\mathrm{Ly}\alpha$.

• •

Case 5) Both the galaxy and IGM come towards the observer, but the IGM travels faster. The absorption only happens in bluer frequencies than $\mathrm{Ly}\alpha$.

• •

Case 6) The IGM moves away from the observer while the galaxy approaches it. The $\mathrm{Ly}\alpha$ photons are blueshifted when they reach the IGM. Wavelengths redder than $\mathrm{Ly}\alpha$ get absorbed.

Additionally, we have computed the relative abundance of these six scenarios in the galaxy population of GALFORM. We find that the cases 1) and 4) are equally probable, as well as 2) with 5) and 3) with 6). In particular, cases 1) and 4) constitute the 25% of the galaxy population each. Moreover, 2) and 4) represent 20% each, while 3) and 6) only 5% each.

On one hand, among the three scenarios where galaxies are moving away from the observer (cases 1, 2 and 3), in two of them (cases 1 and 3) the relative velocity between the IGM and galaxy is positive. In these cases the $\mathrm{Ly}\alpha$ line profile is received redshifted in the IGM frame, causing low absorption. On the other hand, when the galaxy is approaching the observer (cases 4, 5 and 6) in two scenarios (4 and 6) the IGM sees $\mathrm{Ly}\alpha$ blueshifted and absorbs $\mathrm{Ly}\alpha$ photons. In addition, 60% of galaxies moving away from the observer are redshifted in the frame of the IGM. Meanwhile the 60% of galaxies approaching the observer are seen blueshifted by the IGM. This asymmetry causes that galaxies with $\rm V_{z}>0$ (getting away from the observer) are more likely to be observed as LAEs.

5.2.4 IGM velocity gradient.

A further level of complexity is given by the cosmological velocity structure of the Universe (see Fig.2). The IGM is divided into regions of coherent motion of hundreds of $\rm cMpc$ that can be collapsing or moving away from each other. Additionally, the areas with $V_{\rm LoS}\sim 0$ between these regions are small (${\rm\sim 5cMpc}h^{-1}$), causing a great contrast of velocity on scales critical666 Typically, the Hubble flow redshifts the $\mathrm{Ly}\alpha$ line 1Å per ${\rm\sim 3cMpc}$ at $z=3.0$. to the IGM absorption.

There are two opposite effects controlling the IGM transmission dependence on $\partial_{z}{\rm V}_{z}$. On the one hand, the transition between coherent motion regions facilitates the escape of $\mathrm{Ly}\alpha$ photons emitted in clouds with $\partial_{z}V_{z}>0$, as illustrated in the left panel of Fig.12. For example, if a given galaxy lies in the border of the IGM cloud with $\rm V_{z}>0$ (bottom panel), $\mathrm{Ly}\alpha$ photons are observed strongly redshifted in static ($\rm V_{z}\sim 0$) regions, and even more in clouds with $\rm V_{z}<0$. In this case the $\mathrm{Ly}\alpha$  escapes more easily. Meanwhile, in the opposite scenario, where the galaxy lies in a region with $\partial_{z}V_{z}<0$, the neutral hydrogen increases the velocity towards the galaxy as the photon travels. This would result in a greater absorption. On the other hand, in regions with $\partial_{z}V_{z}<0$ the Hubble flow effect is enhanced and photons escape more easily.

The Hubble parameter determines which of these effects prevail. If $H(z)$ is low, the typical distance that $\mathrm{Ly}\alpha$ photons have to travel before the IGM becomes transparent is greater. In this case, this distance is compatible with the transition region between $\rm V_{z}<0$ and $\rm V_{z}>0$. Therefore, the $\mathrm{Ly}\alpha$ transmission correlates with $\partial_{z}{\rm V}_{z}<0$. We find that in our model this happens at z=2.2 and 3.0. Meanwhile, if $H(z)$ is high, then the distance at which the IGM becomes transparent is smaller than the transition region. In this scenario regions with $\partial_{z}{\rm V}_{z}>0$ exhibit greater transmission. We find that this last scenario dominates $z=5.7$ in our models.

6 The clustering of LAEs

In this section we analyze how the coupling of $\mathrm{Ly}\alpha$ detectability and the IGM large scale properties modify the clustering of LAEs. To do so we compare the full RT LAE sample (FLAE) and the SLAE. The SLAE populations inherit the mass function from the FLAE population. In contrast with the FLAE galaxies, by construction, the SLAE positions do not depend on the IGM transmission.

6.1 2-point 3D Correlation function.

In this section we analyze the clustering of our samples in real space. Note that similar result are found in redshift space.

In Fig. 14 we show the real-space 3D 2-point correlation function $\xi(r)$ and the bias computed as

[TABLE]

where $\xi_{m}(r)$ is the matter correlation function.

By comparing LAE samples that include the full RT processes (FLAE) and their shuffled samples (SLAE) we can understand the IGM impact. We note that, by construction, our SLAE samples do not isolate the IGM effects but also includes assembly bias processes. However, we have checked that if we create shuffled samples from the PLAE populations (only RT in the ISM) the clustering measurements are identical, i.e., no evidence of assembly bias is found. Therefore, we attribute the differences in the clustering between FLAE and SLAE samples to only the IGM impact.

Overall, we find that the SLAE populations (both, Thin Shell and Galactic Wind geometry) behave in the same way in every redshift bin. Below $\sim 5{\rm cMpc}h^{-1}$ the bias of SLAE is not constant and decays with distance. From $\sim 5{\rm cMpc}h^{-1}$ on, their bias becomes scale-independent at a very similar value for both geometries. The bias increases with redshift. In detail, the bias of the SLAE on scales larger than $10{\rm cMpc}h^{-1}$ is $\sim 1.8$, $\sim 2.4$ and $\sim 5.2$ at redshift 2.2, 3.0 and 5.7 respectively.

Meanwhile, the FLAE samples exhibit different behaviours at different epochs. On one hand, at redshifts 2.2 and 3.0, the FLAE and SLAE clustering exhibit the same trends and they are almost indistinguishable. Therefore, we find that the IGM does not shape the LAE clustering at these redshifts.

On the other hand, at larges distances, we find that the IGM increases the clustering of FLAEs at $z=5.7$. At scales smaller than $\sim 5{\rm cMpc}h^{-1}$ the FLAE and SLAE clustering are identical. However, the FLAE samples including RT in the IGM exhibit a scale dependent clustering excess on scales ${\rm>20cMpc}h^{-1}$. The boost is present in both outflow geometries. However, its amplitude changes with the geometry. We find that the Galactic Wind geometry exhibit a more powerfull boost than the Thin Shell. We attribute this to the greater coupling with the IGM in the Galactic Wind than in the Thin shell. In particular, the clustering at $50{\rm cMpc}h^{-1}$ is boosted a factor of $\sim 1.1$ and $\sim 1.3$ in the Thin Shell and Galactic Wind geometry respectively. We study the origin of the clustering boost in next section.

6.2 The LAE clustering at large scales.

In this section we interpret our clustering measurements with the physical model presented in ZZ11. They described the overdensity field of galaxy population correlated with the large scale properties of the IGM as

[TABLE]

where $a$ is the scale factor, $\delta_{m}$ is the overdensity of matter in the universe, $r_{\rm H}$ is a length scale set the coefficients dimensionless and the parameters $\tilde{\alpha_{i}}$ are free parameters that quantify the coupling with the IGM properties.

From this expression, the monopole of the galaxy power spectrum can be expressed in real space as a function of the matter power spectrum $P_{m}$ :

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

where $f$ the growth factor and $\beta=f/b$. Meanwhile in redshift space, due to the redshift-space distortions (Kaiser, 1987) $\gamma_{2}$ is rewritten as

[TABLE]

In this description the monopole ($P_{0}$) exhibits other additional dependencies on scale rather than $P_{m}$. There are two terms modifying the shape of the LAE power spectrum. First, a term proportional to $\tilde{\alpha}_{4}/k$ that enhances the clustering at large scales if the velocity field and LAE distribution are coupled. Second, a term proportional to $\tilde{\alpha}_{5}k$ that amplifies the clustering in small scales if the LAEs are coupled to the gradient of density along the line of sight. Finally, the terms proportional to $\tilde{\alpha}_{1}$, $\tilde{\alpha}_{2}$ and $\tilde{\alpha}_{3}$ are scale independent and they only modify the clustering amplitude, leaving the monopole shape unaffected.

We perform a MCMC fit to Eq.26 using our model predictions at z=5.7 in real space. We assume that the bias $b$ of the halo population hosting the FLAE samples is the same as their derived SLAE sample. In practice we compute $b$ for each geometry outflow comparing the SLAE clustering to the matter clustering through Eq.24 and averaging between ${\rm 10cMpc}h^{-1}$ and ${\rm 30cMpc}h^{-1}$. Additionally, we set $\tilde{\alpha}_{1}$, $\tilde{\alpha}_{2}$, $\tilde{\alpha}_{3}$ and $\tilde{\alpha}_{5}$ to zero, while $\rm\tilde{\alpha}_{4}/r_{H}$ is the only free parameter. In fact, we checked that the quality of the fit does not improve when the other parameters are included in the analysis. Finally, we restrict our analysis to the linear regime, i.e., between 5 and $60{\rm cMpc}h^{-1}$.

The bias measurements and fitting results are listed in Table 3. In general, we find that the Thin Shell and Galactic Wind exhibit similar bias, as seen in Fig.14. Additionally, $\tilde{\alpha}_{4}$ is greater in the Galactic Wind than in the Thin Shell, indicating that the LAE model with the first one is more coupled to the IGM.

In Fig. 15 we compare our FLAE models (dots) with the MCMC output (solid lines). The ZZ11 description matches our simulations remarkably well for both outflow geometries. The extension of our models to larger scales highlights the complicated shape of the clustering at cosmological scales. In general, the ratio between the $\xi_{\rm LAE}(r)$ and $\xi_{\rm DM}(r)$ shows a hill at $\sim 80{\rm cMpc}h^{-1}$ and a valley at $\sim~{}100{\rm cMpc}h^{-1}$, while it increases at even larger scales.

Also, in Fig. 15 we compare the matter 2-point correlation function (2PCF) with the analytic clustering description (calibrated with the MCMC analysis) of our Thin Shell and Galactic Wind FLAE samples. For the sake of a better comparison, we have renormalized the clustering of the different samples so that the 2PCF maximum between 90 and 110${\rm cMpc}h^{-1}$ matches in all cases. In the matter 2PCF we find the baryon acustic oscilation (BAO) peak, which is produced by the balance between gravity and pressure in the early Universe (e.g. Chaves-Montero et al., 2018). However, in the analytical description calibrated using our FLAE model, the shape of the BAO peak is distorted. In particular, it becomes broader in the Thin Shell geometry and even wider in the Galactic Wind geometry, as it is more affected by the IGM. Not only that, but also, the position of the maximum in our fits shifts by up to ${\rm~{}1cMpc}h^{-1}$ and ${\rm~{}3cMpc}h^{-1}$ respectively.

6.3 Clustering parallel and perpendicular to the line of sight.

In this section we study the clustering of LAEs in the perpendicular and parallel directions to the line of sight. In order to understand the role of the IGM we compare directly our FLAE and SLAE samples.

First, we make our analysis in real space to quantify the impact of the IGM on the LAE apparent spatial distribution. We perform our analysis at $z=5.7$, where the coupling with the IGM is strong. Note that we find similar results at the other redshifts too. In Fig.16 we show the real space clustering of our LAE populations split in parallel and perpendicular to the line of sight of the IGM. We find that the LAE parallel and perpendicular clustering are symmetric. No angular dependence is apparently found in any of our models.

We also analyze the clustering of our LAE samples in redshift space. In order to convert from real space to redshift space we modify the galaxy positions with their peculiar motion along the LoS of the IGM ($\rm Z$ axis of our simulation) following

[TABLE]

where $x_{\rm Z}$ is the galaxy position along the LoS and $\rm V_{\rm Z}$ is the galaxy peculiar velocity along the same direction.

In the bottom panel of Fig.16 we show the 2D clustering of our LAE samples at $z=5.7$. In general, we find that in redshift space the clustering of all our models is suppressed along the LoS and conserved in the direction perpendicular to the LoS in every redshift bin. Additionally, the suppression along the LoS increases at lower redshift. In general, the peculiar velocity of galaxies points towards overdensities. Hence, the fluctuations along the line of sight are enhanced, which translates into the distortion of the clustering along the LoS.

Furthermore, we also analyze the clustering quadrupole of our FLAE and SLAE samples. In particular, the quadrupole is computed as

[TABLE]

where $\ell=2$, $\mathbcal{L}_{\ell}$ is the Legendre’s polynomial of degree $\ell$ and $\mu=\cos\theta$, where $\theta$ is the angle between the line of sight and the sources.

In Fig.17 we show the quadrupole of the FLAE (continuum lines) and SLAE (dashed lines) samples at different redshifts in real space (top panels) and in redshift space (bottom panels). Overall, we find that at all redshifts the FLAE and SLAE quadrupoles are compatible in both, real and redshift space. Additionally, in real space (top panels) the quadrupole is null at all scales. This implies that the LAE-IGM coupling affects the parallel and perpendicular clustering in the same fashion, and no asymmetry is created. Then, in redshift space, we find a non-null quadrupole that evolves with redshift. The quadrupole becomes steeper on scales lower than $30{\rm cMpc}h^{-1}$ at higher redshifts.

7 Discussion.

In this section we compare our results to other theoretical works and current LAE observations.

7.1 Comparison with previous theoretical works.

In ZZ11the authors studied the coupling of the LAE observability and the IGM large scale properties. Later, Behrens et al. (2017) (from now on BC17) made the same analysis with similar techniques implemented in a simulation with higher spatial resolution. In this context, we compare our approach and results to these works.

7.1.1 Simulations

ZZ11 used a hybrid approach, where first a high resolution dark matter only N-body simulation was evolved while, on the fly, the hydrodynamic physics were run at lower resolution. The size of their simulation was $({\rm 100cMpc}h^{-1})^{3}$ with a constant spatial resolution of the neutral hydrogen density field close to ${\rm 0.13cMpc}h^{-1}$. For further detail about the simulation we refer the reader to Zheng et al. (2010).

Meanwhile, BC17 implemented the IGM RT in the Illustris simulation (Nelson et al., 2015), a high resolution full hydrodynamic simulation using an adaptive mesh refinement approach. BC17 rebinned Illustris to a uniform grid of resolution ${\rm 0.02cMpc}h^{-1}$. The size of this simulation was $({\rm 107cMpc}h^{-1})^{3}$. For further information see their original work (Behrens et al., 2017).

Neither ZZ11 or BC17 were able to resolve the complex ISM structure due to its sub-$\rm kpc$ nature. Additionally, the volume of their simulation was not large enough to trace the large scale variation of the velocity field. In fact, their simulation could be enclosed within one of the many coherent motion regions populating our simulation (see Fig.2).

In contrast, our work uses a dark matter N-body simulation and implements the baryons in a post-processing flavour. Our resolution for the neutral hydrogen density is ${\rm 0.2cMpc}h^{-1}$. However, our simulation size is ${\rm 542.16cMpc}h^{-1}$, which translates into more than 125 times more volume than previous studies. This allow us to make accurate clustering predictions up to ${\rm\sim 60cMpc}h^{-1}$ and resolve the large scale structure of the velocity field of the Universe.

7.1.2 IGM radiative transfer methodology.

Both, ZZ11 and BC17  implement the RT in the IGM using post-processing Monte Carlo approaches similar to other works in the literature studying the $\mathrm{Ly}\alpha$ RT in the ISM (e.g. Verhamme et al., 2006; Orsi et al., 2012; Gronke et al., 2016; Gurung-Lopez et al., 2018a). In particular, in galaxy locations they generate photons in random directions and frequencies assuming a Gaussian line profile centered on $\mathrm{Ly}\alpha$. Then they track the photon’s trajectory and changes in frequency. Finally, photons are collected and an $f_{\rm esc}^{\rm IGM}$ is computed for each galaxy comparing the number of emitted photons and the number of photons received within a given aperture centered on the galaxy position. Additionally, the $\mathrm{Ly}\alpha$ luminosity emitted to the IGM by each galaxy is assumed to be directly proportional to the SFR (although in a different way in each work). For further details, we refer the reader to ZZ11 and BC17 original works.

Meanwhile, here we use a different approach also explored by Weinberger et al. (2018). We divide our simulation into thin cells along the line of sight and analytically compute the IGM transmission of each cells. Then we compute an absorption profile for each galaxy by summing all the small contribution of each cell. Finally, we convolve the IGM transmission with the $\mathrm{Ly}\alpha$ line profile that depends on galaxy ISM properties such as the cold gas mass or the metallicity, as well as the $\mathrm{Ly}\alpha$ luminosity emitted to the IGM (Gurung López et al., 2018b).

7.1.3 IGM transmission

In comparison, ZZ11 and BC17  predict an IGM escape fraction well below that in our models. ZZ11 and BC17 RT approaches greatly overestimate the IGM absorption in comparison to our work for two main reasons:

The assumption of a Gaussian line profile centered on $\mathrm{Ly}\alpha$ in contrast to much observational evidence supporting that the $\mathrm{Ly}\alpha$ line profile is modified (and normally redshifted) by the ISM (e.g. Verhamme et al., 2008; Gronke, 2017; Sobral et al., 2018b). In this way, too much flux is put at bluer frequencies than $\mathrm{Ly}\alpha$, where the IGM is more efficient at absorbing photons (see Fig. 5). 2. 2.

The assumption that the $\rm L_{{\rm Ly}\alpha}$ emitted by a galaxy is directly proportional to the SFR. In fact, in our previous work (Gurung López et al., 2018b) we found that the $\mathrm{Ly}\alpha$ RT in the ISM breaks this relation. Their assumption places LAEs in more massive DM halos, hosted in denser environments, thus with lower median transmission (see Fig.6).

However, our model also has limitations. For example, we assume that every photon that interacts with the IGM is lost. In this way, we underestimate $f_{\rm esc}^{\rm IGM}$ since we do not take into account photons that are scattered out of the LoS and thanks to consecutive scatters are sent back in the LoS again. The contribution of these photons to the received $\rm L_{{\rm Ly}\alpha}$ is predicted to be small (Zheng et al., 2010).

We note that at $z=5.7$ the galaxy stellar mass distributions for our LAE samples truncates abruptly at $\rm M_{*}=10^{7}{\rm M_{\odot}}h^{-1}$. This is caused by the cut in GALFORM at this stellar mass imposed to ensure a good resolution of galaxies. This suggests that a fraction of the LAE population at $z=5.7$ would inhabit galaxies with $\rm M_{*}<10^{7}{\rm M_{\odot}}h^{-1}$. However, as our galaxy population lacks these low mass galaxies, other more massive galaxies are selected as LAEs. Computing the precise number of LAEs that should host galaxies with $\rm M_{*}<10^{7}{\rm M_{\odot}}h^{-1}$ is challenging. As a simple calculation, we rescale the stellar mass distribution at $z=3.0$ towards smaller $\rm M_{*}$ so that the peak of the distribution matches the one at $z=5.7$. Then, the fraction of galaxies below the resolution limit is $\sim 0.04$. This hints that only a small part of the LAE population would lie in galaxies not resolved in our model. However, this limitation of our model might cause an overestimation of the clustering bias at this epoch (see below).

On one hand, ZZ11 found a strong correlation between the IGM and the LAE population. Their model predicted a relatively small dependence on $\rm V_{z}$, $\rho$ and $\rm\partial_{z}\rho$ and a tight relation between $\rm\partial_{z}V_{z}$ and $f_{\rm esc}^{\rm IGM}$, causing differences in $f_{\rm esc}^{\rm IGM}$ greater than 1 order of magnitude across the $\rm\partial_{z}V_{z}$ dynamical range.

In their model, the strong coupling between LAEs and $\partial_{z}V_{z}$ has dramatic consequences in the clustering of LAEs at $z=5.7$. Their model predicts a scale independent enhancement in the clustering along the line of sight, which created an asymmetry between the clustering parallel and perpendicular to the LoS in real space. Moreover, in their model, the amplitude of this effect outpowers the Kaiser boost. Hence, even in redshift space the clustering along the LoS is more powerful than in the perpendicular direction.

On the other hand, BC17 claimed to find only a marginal coupling between the LAE and the IGM large scale properties at $z=2.0$, 3.0, 4.0 and 5.85. Additionally, they studied the clustering of LAE samples and did not find any asymmetry between the directions parallel and perpendicular to the LoS clustering, nor any other strange clustering feature. Additionally, BC17 claimed ZZ11 results to be a consequence of poor spatial resolution. In particular, they recovered the strong selection on $\partial_{z}V_{z}$ after lowering their resolution to match ZZ11 simulation.

In this work we find a coupling between the $\mathrm{Ly}\alpha$ observability and the large scale IGM properties. The amplitude of the coupling depends on redshift and on the outflow structure assumed in the ISM. Indeed, typical variations in $f_{\rm esc}^{\rm IGM}$ are $<1\%$ , $\sim 2\%$ and $\sim 5\%$ at $z=2.2$, 3.0 and 5.7 respectively. In other words, we detect a LAE-IGM coupling but it is small, in particular at $z=2.2$ and 3.0 where it has negligible effect on the clustering.

We do not find the dramatic asymmetry of clustering parallel and perpendicular to the line of sight (see Fig. 16) that ZZ11 measured in their original work. Our models do not predict any other clustering modification on scales where ZZ11 or BC17 simulations allowed them to measure the clustering accurately (${\rm\sim 10cMpc}h^{-1}$). Hence our clustering predictions are in agreement with BC17 and differ from ZZ11.

Our model predicts a strong boost of the LAE clustering at scales ${\rm>20cMpc}h^{-1}$. This feature comes from the coupling with the IGM velocity field (see from eq.6.2 to eq.30). We attribute the detection of this effect in our work and the non-detection in BC17 to the difference in volume probed by the different simulations. Meanwhile, our simulation resolves the velocity large scale structure of the Universe, ZZ11 and BC17 simulations are small enough to be enclosed in one of the huge regions of coherent motion (see Fig. 2).

7.2 Comparison with an analytical model

In Wyithe & Dijkstra (2011) (from now on WD11), authors studied the impact of the LAE-IGM coupling in the accuracy of LAE cosmological surveys. Here we make a brief comparison between this work and WD11 and we encourage the readers to visit their original article for full details.

The approach taken by WD11 is very different to the one developed through this work. On one side, we have created an LAE model over a cosmological dark matter simulation (P-Millennium) using a model of galaxy formation and evolution (GALFORM) and then we have implemented the RT in the ISM (through FLaREON) and in the IGM by computing the IGM transmission a long the line of sight for each galaxy. On the other side, WD11 designed an analytical physically motivated model of the LAE power spectrum that took into account the coupling between the observability of $\mathrm{Ly}\alpha$ and the IGM conditions (similar to the one presented in ZZ11). In particular, their model included the dependence with the density of the IGM, its velocity gradient along the line of sight and its ionization rate. In this work we have found that the coupling between the $\mathrm{Ly}\alpha$ observability and the IGM velocity along the line of sight is the most relevant on large scales (see §6.2). Meanwhile we report that our model does not present a significant coupling with the ionization state of the IGM in the redshift range studied here.

There are two main conclusions on WD11 regarding the accuracy of cosmological constrains using LAEs :

• •

On one hand, the BAO wiggles are not heavily affected by the coupling of IGM and $\mathrm{Ly}\alpha$ observability. In fact, they reported that precision on cosmology given by BAO studies using LAEs should be similar to the one achieved using non-IGM-coupled galaxy populations.

• •

On the other hand, they argue that in such galaxy surveis most of the constrain power in the cosmology comes from the full shape of the power spectrum (e.g. Shoji et al., 2009) . As a consequence, they found that the constrains in cosmology given by an LAE surveys are weaker than in a normal galaxy survey with the same volume and number density. Moreover, the impact on the accuracy would depend on the level of coupling with the IGM and it could become highly significant if the IGM coupling was strong.

Overall, we find a good agreement with the finding of WD11. First, we find that the BAO peak is shifted only $\sim$2Å (depending on the outflow geometry) from its original position at redshift 5.7, while smaller displacement are found at lower redshifts. This would indicate that the BAO analysis would be affected, but not heavily by the LAE-IGM coupling predicted by our model. Additionally, our model includes different models for the escape of $\mathrm{Ly}\alpha$ photons from galaxies, which also might increase the uncertainty in BAO studies.

Second, in contrast to our work, their model do not include the coupling with the velocity field, which, in our case is the one injecting more power to clustering of LAEs. Therefore, the impact of the LAE-IGM coupling in cosmological contains could be stronger in our models, since we find that our clustering is much more disrupted than the one explored in their work. However, this analysis needs a high level of detail and we will leave it to a future study.

7.3 Comparison with current observations.

In this section we compare our model’s clustering with current LAE clustering observations across cosmic time. In order to make a fair comparison we mimic the observations performed by Ouchi et al. (2010, 2018), Bielby et al. (2016) and Kusakabe et al. (2018). These surveys select LAE candidates by combining narrow and broad band photometric filters. The main differences between the different surveys are the flux depth, the sky coverage and the redshift of LAE candidates. The different properties of the observations used in this work are listed in Table 5 and Table 4.

Here we summarize the mock construction and refer the readers to Gurung López et al. (2018b) for a deeper explanation. In short:

We convert the galaxy coordinates to redshift space along the same line of sight (LoS) used to compute the IGM absorption (Eq.31). 2. 2.

We divide our box into rectangular subvolumes. The faces perpendicular to the LoS are squares with the same area as the survey sky coverage. The depth along the LoS is determined by the narrow band full width half maximum of each survey. 3. 3.

We make a rest frame equivalent width ($\rm EW_{0,cut}$) and luminosity cut ($\rm L_{Ly\alpha,cut}$). We choose the closest $\rm EW_{0,cut}$ and $\rm L_{Ly\alpha,cut}$ to the survey cuts that best match the observed LAE number density. These cuts are listed in Tab. 5.

The mock properties are listed in Table 4. At $z=2.2$ and $3.0$ the survey volumes are small in comparison with our simulation and there are about 450 mocks in each redshift bin. Hence, we get a good estimate of the cosmic variance of these observations. However, at $z=5.7$ the number of mock catalogs decreases to 18 and the constrain on cosmic variance is weaker. Additionally, the observed number density is well reproduced in our mocks. In particular, the median number of LAEs in each mock is within 1$\sigma$ of the observed value.

In Fig. 18 we compare our full RT model mocks clustering with measurements at several redshifts. Overall, we find a good agreement with observations. In particular, at redshift $2.2$, our RT models match very well observation at all scales. Meanwhile, at $z=3.0$ the clustering is perfectly reproduced within the cosmic variance at angular separations higher than $\theta\gtrsim 10\;arcsec$, while at smaller scales is underpredicted. Finally, at $z=5.7$ we find that the full RT LAE samples exhibit a higher bias at all scales. This small disagreement is not caused by LAE-IGM coupling predicted by our model. If this was the case, only the large scales would disagree. In fact, the clustering excess here is partially caused by the overpredicted stellar masses at $z=5.7$ in our FLAE samples, as described above (see Fig.9).

The characteristic excess of power in the clustering of FLAEs populations has been tentatively detected in Ouchi et al. (2018). For example, in Fig. 18 we show the LAE clustering measurements of SILVERRUSH (Ouchi et al., 2018) (in blue dots). The clustering at angular separations $>10^{2}arcsec$ is boosted. However, this trend is not found in Ouchi et al. (2010) (green dots). We attribute this to the different volume traced by surveys. On the one hand, Ouchi et al. (2010) only covered $\rm\sim 1\;deg^{2}$, an area small enough to fall inside one region of coherent motion (see Fig.2). On the other hand (Ouchi et al., 2018) explored $\rm\sim 13\;deg^{2}$ with contiguous patches of up to $\rm\sim 5\;deg^{2}$, big enough to begin to resolve these regions.

8 Conclusions

We have created a cosmological model of Lyman-$\alpha$ emitter galaxies that includes Lyman $\alpha$ radiative transfer physics in both the interstellar and intergalactic medium. The cosmological background is imprinted by the P-Millennium N-body simulation. Meanwhile, the semi-analytic model of galaxy of formation and evolution GALFORM populates the P-Millennium DM halos with galaxies. For the ISM transmission we used FLaREON (Gurung-Lopez et al., 2018a), an open Python package based on a Monte Carlo RT code (Orsi et al., 2012) that predicts the $\mathrm{Ly}\alpha$ line profiles and escape fractions of photons in outflows of different characteristics. Meanwhile, the RT in the IGM is implemented by computing the $\mathrm{Ly}\alpha$ transmission at each position of our simulation. Our main conclusions are:

The RT in the ISM produces a strong selection effect over galaxy properties such as metallicity or SFR. Meanwhile, the RT in the IGM leave the galaxy property distributions nearly unchanged. In fact, LAEs tend to have a low-intermediate metallicity, moderate SFR and intermediate stellar mass. For further analysis on the galaxy properties of LAEs we refer the reader to Gurung López et al. (2018b). 2. 2.

Our models predict that the $\mathrm{Ly}\alpha$ IGM escape fraction depends on the large scale properties of the IGM such as the IGM density, motion, the density and velocity gradient along the line of sight, as first studied by Zheng et al. (2011). While at low redshift ($z=2.2$ and 3.0) the correlations are weak and lead to variations in $f_{\rm esc}^{\rm IGM}$ of $\sim 1\%$, it intensifies at higher redshifts ($z=5.7$), reaching variation on $f_{\rm esc}^{\rm IGM}$ around $\sim 5\%$. 3. 3.

The level of coupling between the LAE distribution and the large scale IGM properties depends on the RT inside the ISM. Our model predicts that if the outflows driving the $\mathrm{Ly}\alpha$ photons escape from galaxies have a Thin Shell geometry the coupling is greater than if it is driven by a Galactic Wind at low redshift. However, the opposite is found at $z=5.7$. 4. 4.

The IGM-LAE coupling can have an impact on their clustering. At redshift 5.7 the shape of the 2-point correlation function is modified, introducing extra power at scales larger than $\sim{\rm 20cMpc}h^{-1}$ in a scale-dependent fashion. Meanwhile, at $z=2.2$ and z=3.0 we do not find any modification in the 2PCF, as the IGM-LAE coupling is too weak. 5. 5.

In order to study the LAE clustering at large scales we fit our model 2-point correlation function (2PCF) to an analytic clustering model including the IGM-LAE coupling (Eq.26) introduced by Zheng et al. (2011). The IGM-LAE coupling disrupts the LAE clustering and modify the shape of the 2PCF at the baryon acustic oscilation (BAO) peak scales, which becomes broader and modify the position of the maximum by $\sim 1{\rm cMpc}h^{-1}$. 6. 6.

We have made a comparison between current LAE clustering measurements and our LAE models. Overall we find good agreement between our model and observations, even at $z=5.7$. This suggests that a greater sky coverage is necessary to detect the clustering excess presented in this work. However, hints of a scale dependent bias at $z=5.7$ can be found in the literature (Ouchi et al., 2018).

In future works we plan to implement the model presented in this work to a larger simulation to determine the IGM-LAE coupling impact on the clustering at low redshifts. Additionally, we will also implement the physics of reionization to understand the clustering of LAEs during this interesting epoch.

Acknowledgements

The authors acknowledge the useful discussions with Zheng Zheng, Mark Dijkstra and Anne Verhamme in addition to the whole CEFCA team. The authors also acknowledge the support of the Spanish Ministerio de Economia y Competividad project No. AYA2015-66211-C2-P-2. We acknowledge also STFC Consolidated Grants ST/L00075X/1 and ST/P000451/1 at Durham University. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Appendix A The $\rm N_{H}$, $\rm V_{exp}$ and $\rm\tau_{a}$ distributions

In our model, the $\mathrm{Ly}\alpha$ radiative transfer physics inside galaxies are integrated through FLaREON. FLaREON is open source code that predicts escape fraction and the emerging $\mathrm{Ly}\alpha$ line profile from different outflow configurations among several gas geometries. In particular, FLaREON is based on pre-computed {$\rm N_{H}$, $\rm V_{exp}$, $\rm\tau_{a}$ } grids of the full radiative transfer Monte Carlo code LyaRT (Orsi et al., 2012). Then, different algorithms such as multidimensional interpolation are used to obtain the line profile and $f_{\mathrm{e}sc}^{\rm Ly\alpha}$. Hence, the high performance of FLaREON is limited to the space covered by the grid. In the following we study the fraction of LAEs in our model that fall within the FLaREON range.

In Fig. 19 we show the $\rm N_{H}-V_{exp}$ distributions at different redshifts (solid colored lines) and the FLaREON accuracy boundaries (black lines). Over all, the Thin Shell and the Galactic Wind models behave likely and there is always some overlap between them. Additionally, as described above, there is little variation between the properties of the samples including the full RT (FLAE) and galaxies with only RT in in the ISM (PLAE). At low redshifts the $\rm N_{H}-V_{exp}$ distributions are quite compact. In particular, at redshift 2.2 both geometries have a similar $\rm N_{H}$ distribution that peaks at $\sim 10^{20.5}cm^{-2}$, while the velocities of the Thin Shell ($\sim 10^{2}km\;s^{-1}$) are slightly above the Galactic Wind distribution ($\sim 10^{1.8}km\;s^{-1}$). Meanwhile, $z=3.0$ the velocity rounds $\sim 10^{1.5}km\;s^{-1}$ and $\sim 10^{1.8}km\;s^{-1}$ in the Thin Shell and Galactic Wind respectively, while $\rm N_{H}$ are close to $\sim 10^{20.2}cm^{-2}$ and $\sim 10^{20.8}cm^{-2}$.

Finally, at $z=5.7$ the distributions become broader. The velocities are higher and similar in both geometries, rounding $\sim 10^{2.4}km\;s^{-1}$, while the Thin Shell have higher $\rm N_{H}$ ($\sim 10^{19.8}cm^{-2}$) than the Galactic Wind ($\sim 10^{18.8}cm^{-2}$). We find that at that the fraction of LAEs outside the FLaREON range is less than a 3% in any of our samples.

In Fig. FLaREON we show the resulting $\tau_{a}$ distribution after calibration for our LAE samples. The $\tau_{a}$ distribution both geometries at $z=2.2$ are narrow and centered at $\sim 10^{-1.5}$. At redshift 3.0 the distribution becomes slightly broader. Additionally, the Thin Shell exhibits lower dust optical depth ($\sim 10^{-1.8}$) than the Galactic Wind ($\sim 10^{-1.5}$). Meanwhile, at $z=5.7$ the distribution become wider and the Galactic Wind have lower $\tau_{a}$ ($\sim 10^{-3.5}$) than the Thin Shell ($\sim 10^{-2.5}$). We also show the accuracy border of FLaREON at $\log\tau_{a}=0.0$. We find that less than 1% of the LAE in our sample lie outside FLaREON $\tau_{a}$ dynamical range. All together we find that the fraction of LAEs outside FLaREON accuracy range is negligible.