\mathrm{Ly}\alpha emitters in a cosmological volume I

 emitters in a cosmological volume I: the impact of radiative transfer

Siddhartha Gurung-López, Álvaro A. Orsi, Silvia Bonoli, Carlton M. Baugh, and Cedric G. Lacey.
Centro de Estudios de Física del Cosmos de Aragón, Plaza San Juan 1, piso 2, Teruel, 44001, Spain.
Institute for Computational Cosmology, Department of Physics, Durham University, Science Laboratories, South Road, Durham, DH1 3LE.
E-mail: sidgurung@cefca.es
Accepted XXX. Received YYY; in original form ZZZ

Lyman- emitters (LAEs) are a promising target to probe the large scale structure of the Universe at high redshifts, . However, their detection is sensitive to radiative transfer effects that depend on local astrophysical conditions. Thus, modeling the bulk properties of this galaxy population remains challenging for theoretical models. Here we develop a physically-motivated scheme to predict LAEs in cosmological simulations. The escape of  photons is computed using a Monte Carlo radiative transfer code which outputs a  escape fraction. To speed-up the process of assigning escape fractions to individual galaxies, we employ fitting formulae that approximate the full Monte Carlo results within an accuracy of 10% for a broad range of column densities, gas metallicities and gas bulk velocities. We apply our methodology to the semi-analytical model GALFORM on a large -body simulation. The  photons escape through an outflowing neutral gas medium, implemented assuming different geometries. This results in different predictions for the typical column density and outflow velocities of the LAE population. To understand the impact of radiative transfer on our predictions, we contrast our models against a simple abundance matching assignment. Our full models populate LAEs in less massive haloes than what is obtained with abundance matching. Overall, radiative transfer effects result in better agreement when confronting the properties of LAEs against observational measurements. This suggest that incorporating the effects of  radiative transfer in the analysis of this galaxy population, including their clustering, can be important for obtaining an unbiased interpretation of future datasets.

keyword1 – keyword2 – keyword3
pubyear: 2018pagerange:  emitters in a cosmological volume I: the impact of radiative transferB

1 Introduction

During the past two decades, surveys targeting the  emission in star-forming galaxies, the so-called  emitters (LAEs), have detected objects out to redshift (e.g. Steidel et al., 1996; Hu et al., 1998; Rhoads et al., 2000; Malhotra & Rhoads, 2002; Taniguchi et al., 2005; Kashikawa et al., 2006; Guaita et al., 2010; Konno et al., 2016; Sobral et al., 2017).The study of this galaxy population has allowed us to explore the kinematics of the interstellar medium (ISM) in high redshift galaxies (Shapley et al., 2003; Steidel et al., 2010, 2011; Kulas et al., 2011; Guaita et al., 2017; Chisholm et al., 2017), the large scale structure (Gawiser et al., 2007; Orsi et al., 2008; Ouchi et al., 2010; Bielby et al., 2016; Kusakabe et al., 2018; Ouchi et al., 2018), the epoch of reionization (Santos et al., 2004; Kashikawa et al., 2006; Dayal et al., 2011; Inoue et al., 2018) and to test galaxy formation models (Le Delliou et al., 2006; Kobayashi et al., 2007; Nagamine et al., 2010; Orsi et al., 2012).

Despite the success in detecting progressively larger samples of LAEs, their physical interpretation has proven to be a difficult challenge (see Dijkstra, 2017, for a review).  photons are easily scattered by neutral hydrogen, causing a large increase in the path that the photon needs to travel through neutral hydrogen clouds (e.g. Harrington, 1973; Neufeld, 1990). This results in an increased probability of interaction with dust grains, and thus, absorption. Hence, the  radiative transfer through a neutral medium reduces the  flux that escapes the galaxy and also modifies the line profile, since each scattering event changes the frequency of the photons. These physical processes also take place in the surrounding intergalactic medium (IGM) of galaxies and can also modify the observed  flux and line profile (Santos et al., 2004; Dijkstra et al., 2011).

Analytical approximations for  radiative transfer have been derived for over-simplistic neutral gas configurations (e.g. Harrington, 1973; Neufeld, 1990; Dijkstra et al., 2006). More realistic configurations can be explored with a Monte Carlo algorithm. Individual  photons are generated inside a neutral hydrogen cloud with a given geometry, kinematics and temperature. The path of  photons is tracked including their interactions, which produce scattering events, until the photons escape or are absorbed by dust. This approach has been studied in several scenarios (Ahn et al., 2000; Zheng & Miralda-Escudé, 2002; Ahn, 2003; Verhamme et al., 2006; Gronke et al., 2016). Most notably, Monte Carlo radiative transfer has shown to reproduce the diversity of observed  line profiles by allowing photons to escape through an outflowing medium (e.g. Schaerer & Verhamme, 2008; Orsi et al., 2012).

Theoretical models of galaxy formation have introduced the effect of radiative transfer in different approximate ways to predict the properties of the LAE population. The first model of LAEs in a hierarchical galaxy formation framework implemented a constant escape fraction of  photons to reproduce their observed abundance and clustering (Le Delliou et al., 2005; Le Delliou et al., 2006; Orsi et al., 2008). Further attempts introduced radiative transfer effects over simple geometries in semi-analytical models (Orsi et al., 2012; Garel et al., 2012). Cosmological hydrodynamical simulations also incorporated  radiative transfer in post-processing. One approach has been to track  rays to simulate different lines of sight (e.g. Laursen & Sommer-Larsen, 2007; Laursen et al., 2009, 2011) over small volumes. With a Monte Carlo radiative transfer code, Zheng et al. (2010) showed that the proper treatment of  photons radiative transfer has dramatic effects on the clustering of LAEs. However, recently, Behrens et al. (2017) found no significant change in the clustering of LAEs after implementing  radiative transfer in the Illustris simulation (Nelson et al., 2015), and attribute the claims of Zheng et al. (2010) about the clustering of LAEs to resolution effects.

In the next years many ground-based large surveys such as HETDEX (Hill et al., 2008), J-PAS (Benitez et al., 2014) and space missions like ATLAS-Probe (Wang et al., 2018), will aim to detect LAEs over large areas to trace the large scale structure (LSS) at high redshifts. Such measurements could potentially deliver cosmological constraints in redshift ranges well above those currently targeted by Multi-Object Spectroscopic surveys. With progressively larger and more accurate datasets, it becomes crucial to improve our theoretical understanding of galaxies as tracers of the underlying matter distribution (Orsi & Angulo, 2018). One of our aims in this work is to understand the impact of radiative transfer effects on clustering measurements.

The model for the  luminosity of star-forming galaxies presented here is based on a fast implementation of a Monte Carlo radiative transfer. To avoid the prohibitively long time that it would take to run a Monte Carlo code over millions of galaxies, we develop fitting formulae that reproduce the full Monte Carlo results accurately. To illustrate the potential of our model, we apply this methodology to the semi-analytic model GALFORM run over an -body simulation. This is a first paper in a series that explores the properties of galaxies selected by their  luminosity. Here we focus on the impact of the  RT in defining the properties of the LAE galaxy population. In a forthcoming paper we implement the impact of the intergalactic medium (IGM) and the effects of reionization on the LAE population.

The outline of this paper is as follows: in §2 we develop fitting formulae to predict the escape fraction of Lyman alpha photons through outflows. In §3, we describe our model for LAEs that combines galaxy formation physics and  radiative transfer in addition to the implementation of the  RT in a galaxy formation model is presented. We analyze the LAE population predicted by our model in §4. We discuss our results in §5. Finally, conclusions and future work are summarized in §6.

2 Model ingredients

Figure 1: A schematic illustration of the different outflow geometries implemented in this work: Thin Shell (left), Wind (middle) and Biconical Wind (right). The gas density is represented by the gray colour scale. Different possible trajectories of photons are labeled from to . The red cross over photon illustrates the point where this photon is absorbed by the medium.

In this section we describe our model ingredients and the methodology we follow to predict the properties of LAEs in a cosmological simulation.

2.1 Ly radiative transfer

We track the scattering, absorption and escape of  photons making use of the Monte Carlo radiative transfer code described in Orsi et al. (2012), which has been made publicly available 111https://github.com/aaorsi/LyaRT. This code is similar to others in the literature (e.g. Zheng & Miralda-Escudé, 2002; Ahn, 2003, 2004; Dijkstra et al., 2006; Verhamme et al., 2006; Laursen & Sommer-Larsen, 2007; Barnes & Haehnelt, 2010, and references therein). A detailed review of  radiative transfer can be found in Dijkstra (2017). Below we summaries the main features of the Orsi et al. (2012) code that are most relevant to this work.

The code receives as input a configuration of a 3D neutral gas geometry, temperature, expansion velocity , neutral hydrogen column density and optical depth of dust . For a given gas distribution, the code generates a  photon with a random direction and follows its interactions with hydrogen and dust until it is either absorbed by dust or escapes from the neutral gas medium. Every interaction with a hydrogen atom results in a scattering event that changes the direction and frequency of the photon. Interactions with dust, on the other hand, can change the direction of the photon or result in absorption depending on the assumed albedo of the dust grains. The process is repeated for photons, recording in the end the frequency of every photon that escaped and those that were absorbed by dust grains. This allows us to compute the escape fraction  and wavelength distribution (i.e. the  line profile) for every outflow geometry over which both the neutral gas and the dust are distributed. In this work we implemented three different outflow geometries, which are illustrated in Fig. 1.

  1. Thin Shell. This geometry consists of an expanding isothermal homogeneous spherical shell. This spherical shell is thin and can be described by an inner and an outer radius, and respectively, which satisfy . The shell is expanding outwards, thus it has a radial macroscopic velocity . The neutral hydrogen column density is given by:


    where is the total neutral hydrogen mass and is the mass of a hydrogen atom.

    The empty cavity in the center of the shell produces photon backscatterings, i.e. photons can bounce back into the empty cavity multiple times, as illustrated by photons and in Fig. 1.

  2. Galactic Wind. This geometry consists of an expanding spherical gas distribution with a central empty cavity of radius . The gas is isothermal and is expanding radially at a constant velocity . Unlike the Thin Shell, the gas is distributed with a radial density profile given by:


    where is the ejection neutral hydrogen mass rate. Thus, the column density in the Wind geometry is


    This geometry is illustrated in the middle panel of Fig. 1.

    We define a large outer radius where the computation is forced to end and any photon that have reached this radius is considered to have escaped. We have checked that for greater values of the code provides the same line profile and escape fraction. Thus, we conclude that our results converge for our choice of .

    Figure 2: (Left) versus the dust optical depth for different geometries in outflows with the same physical properties ( and ), as indicated in the figures. The output of the radiative transfer code is represented by green circles, blue diamonds and red squares for the Thin Shell, galactic wind and biconical geometries respectively. Additionally, our analytical fit is represented by solid lines with the same color code as the code’s output. (Right)  line profile for different geometries with the same physical properties. In colored lines the radiative transfer code output is plotted for the Thin Shell geometry (green), the galactic wind (blue) and the biconical galactic wind (red).
  3. Biconical Wind. This geometry shares the same properties of theWind, but additionally it features an aperture angle, , which defines the volume of gas and dust. In particular we arbitrarily set . The resulting polar asymmetry is thus the main difference between the two previous geometries and this one. This is shown in the right panel of Fig. 1.

    Furthermore, in this geometry we force photons to be emitted from the center of the geometry (as in the other geometries) and within the aperture of the bicone, i.e. no photons are emitted outside the bicone. Additionally, due to the empty regions in this geometry, photons that scatter off the internal cavity and escape off the the bicone are considered absorbed by the external medium (e.g. photon in Fig. 1). This is equivalent to assuming that there is a dusty optically thick medium surrounding the bicone.

Fig. 2 illustrates the difference between the  escape fraction (left panel) and line profile (right panel) predicted by each geometry, for a particular choice of column density and expansion velocity. As expected, the escape fraction decreases towards higher values of in all geometries, as greater amounts of dust absorb more photons. However, the impact on the geometry of the medium is evident: even if the three configurations have the same and , photons have the highest escape fractions from the Wind geometry, and the lowest from the Bicone. This is due to the complicated  RT. For example, as in the Bicone configuration photons that leak through the empty cavity are considered absorbed, the escape fraction does not reach 1 even if there is no dust in the outflow, making a great difference with respect to the other two geometries. Additionally, even if the Wind and Thin Shell configurations share spherical symmetry (unlike the Bicone) the dependence of  on , and is different due to the distinctive hydrogen density radial profiles of the two configurations. This dependence on the geometry does not only affect the  but also the line profile of the  emission. The predicted line shape changes dramatically from a geometry to another: in the case of the Wind it is a broad line, for the Bicone it is a narrow line and for the Thin Shell it assumes a double-peak profile. We use these three different outflow geometries to estimate the variance in the LAEs population depending on the geometry.

2.2 Fitting formulae for  radiative transfer

Thin Shell
Galactic Wind
Biconical Wind
Table 1: Constant parameter values used to derive the escape fraction of the different geometries

As discussed in §2.1, the Monte Carlo radiative transfer code can take a long time to run for a given configuration of parameters. For a single photon, the average number of scatterings, and thus, calculations, grows as a power-law function of the column density of the medium (Harrington, 1973). In the parameter space explored here, the completion time of the code can vary from a few seconds up to a few hours in the most extreme cases. Applying this directly in a catalog of millions of objects would result in prohibitively long execution times.

To overcome this, we develop empirical (measured from the radiative transfer Monte Carlo code) expressions that approximate the results of the Monte Carlo runs. We start by constructing a grid to scan the parameter space with configurations spanning the ranges , and . We run the Monte Carlo code with photons and obtain the  escape fraction,   as a function of , and, .

To construct an analytic expression for  we start from a generalized form of the expression for the  in a homogeneous, static slab derived in Neufeld (1990):


where and are functions of and for all geometries. Additionally, is set to in the Thin Shell and Wind geometries, but is a function, , in the Bicone, since, in this geometry, the escape fraction is always less than 1 (see section 2.1). We perform a Monte Carlo Markov Chain (MCMC) with the emcee222http://dfm.io/emcee/current/ code (Foreman-Mackey et al., 2013a) to determine the functional form of , and , by minimizing the function


where corresponds to the escape fraction of photons obtained with the MC code over each configuration in the grid, and is the error in the calculation of the escape fraction, given by the dispersion in a binomial distribution with probability of success :


where is the -th percentile of the standard normal distribution. In particular we use the quantile 95, i.e. . Additionally, is the number of generated photons in each configuration.

The functional form and parameter values of the fits for , and for each geometry are shown in Table 1.

Fig. 2 compares the  computed analytically with Eq. (4) and with the free parameters obtained with the MCMC (lines), and that obtained with the full MC RT code (symbols) for a given values of and and the three different geometries. The analytical expression reproduce remarkably well the results of the full MC RT code.

The accuracy of our analytic expressions varies with , , and the geometry. In particular, there is a strong dependence on : for every geometry we find that the accuracy decreases with increasing . We find that, in general, the discrepancy with the full MC RT code in configurations with becomes greater than 10%. Galaxies with such a large dust absorption, in general, will not be observed as a LAE so we are not concerned about the low accuracy at high . Additionally, we checked that, after calibration of our LAEs model (see §LABEL:sssec:Calibrating), less than 2% of the galaxies in every geometry have , making the contribution of these galaxies negligible. For galaxies with , the discrepancy is just a few percents for between and and between and . Moreover, for the discrepancy is typically below the 1% in the same parameter range.

A detailed assessment of the accuracy of the analytical expressions for  is presented and discussed in Appendix A

redshift Geometry
Thin Shell 4.440 4.911 -12.367 -11.839
Wind 4.857 4.914 -7.065 -5.338
Bicone 4.982 4.258 -8.140 -7.249
Thin Shell 4.337 4.549 -12.465 -11.915
Wind 4.691 4.769 -7.440 -5.166
Bicone 4.896 4.338 -8.436 -6.404
Thin Shell 4.737 4.428 -13.906 -11.808
Wind 4.660 3.782 -8.292 -6.180
Bicone 4.612 3.590 -8.078 -7.614
Thin Shell 4.659 4.279 -13.81 -11.934
Wind 4.589 3.871 -8.073 -5.910
Bicone 4.455 3.561 -7.848 -7.647
Table 2: Free parameters as defined in equations 7 and 8 after the calibration with the observed luminosity function for different geometries and redshifts.

2.3 Simulation and semi-analytical model.

We combine the radiative transfer code described above with the semi-analytical model of galaxy formation GALFORM (Lacey et al., 2016) run on the P-Millennium -body simulation (Baugh et al., in prep.).

The P-Millennium is a state-of-the-art dark matter only N-body simulation using the Plank cosmology: , , , (Planck Collaboration et al., 2016). The box size is and the particle mass ( dark matter particles). Between the initial redshift, , and the present, , there are 272 snapshots. In this work we use snapshots 77, 84, 120 and 136 corresponding to redshifts 6.7, 5.7, 3.0, 2.2, respectively.

A full review on semi-analytical models of galaxy formation can be found in Baugh (2006). The variant of GALFORM used in this work is based on earlier versions described in Cole et al. (2000); Baugh et al. (2005) and Bower et al. (2006). In brief, GALFORM computes the properties of the galaxy population following the hierarchical growth of dark matter halos. Halo merger trees are extracted from an -body simulation (the P-Millennium  in our case), so the model can also predict the spatial distribution and peculiar velocities of galaxies.

In GALFORM, galaxies are formed and evolve as a result of the following processes: i) the radiative cooling and the shock-heating of gas inside halos; ii) the subsequent cooling of gas forming a disk at the bottom of the potential well; iii) quiescent star formation in the disk and starbursts in bulges resulting from disk instabilities and galaxy mergers; iv) feedback processes (supernovae, AGN and photoionization) regulating the star formation, and v) the chemical enrichment of stars and gas that results from star-formation and feedback episodes. Additionally, the variant of GALFORM used in this work assumes different initial mass functions (IMFs) for quiescent and starburst modes of star-formation (see Lacey et al., 2016, for more details).

GALFORM generates a composite spectral energy distribution (SED) for each individual galaxy based on its star-forming history and computes the rate of emission of hydrogen ionizing photons, , by integrating the galaxy SED over wavelengths bluer than the Lyman break at Å. All ionizing photons are assumed to be absorbed by the neutral medium. Then case B recombination (Osterbrock, 1989) is used to compute the intrinsic line luminosity of , where a fraction of of ionizing photons contribute to generating  photons.

2.3.1 Radiative transfer parameters

To combine the Monte Carlo radiative transfer code with GALFORM, we need to derive the parameters that define the neutral gas configuration from the galaxy output properties. In particular, the column density , expansion velocity and the optical depth of dust are key to determine the escape fraction. The expansion velocity is computed for the three geometries as:


where the index denotes the galaxy component (disk or bulge), and are the SFR and half mass radius of each galaxy component, is the total stellar mass of the galaxy and are two (one per galaxy component) free parameters.

The neutral hydrogen column density is computed in different ways depending on the geometry (see section 2.1) :


where and are, respectively, the cold gas mass and a free parameter of the galaxy component .

All the free parameters linking GALFORM properties to and are calibrated by fitting the observed LAE luminosity function at different redshifts. For further details see §3.1.

Finally, the is computed for every geometry as:


where is the albedo at the Ly wavelength, is the ratio for solar metallicity, (Granato et al., 2000) and is the cold gas metallicity of the galaxy component .

The intrinsic  LF predicted by GALFORM (see Figure 3) results from two populations: normal star forming galaxies (populating the low luminosity range) and galaxies with an ongoing star formation burst (populating the high luminosity range). Consequently, the values of and control the shape of the faint-end LF, whereas and control the bright end of the LF. In both regimes, increasing (decreasing) leads to an increase (decrease) of the distribution. This leads to a decreasing (increasing) in the resulting  distribution and thus lowers (increases) the number of galaxies with higher luminosities. Also, increasing (decreasing) leads to a increase (decrease) of the distribution, increasing (decreasing)  and the number of galaxies with high luminosities.

Figure 3: LAE LF at redshift 2.2 (top left), 3.0, (top right), 5.7 (bottom left) and 6.7 (bottom right). The LF computed for different geometries is plotted as colored continuum lines, in blue for the Wind geometry, in red for the Bicone geometry and in green the Thin Shell geometry. In continuum black we show the intrinsic  LF. The black dashed lines show the combined LF that is fitted that, at the same time, is the AM-noRT LF (detailed in §4) LF. At redshift 2.2 we also show the LF observed by Kono et al 2016 (blue dots), Sobral et al. 2016 (purple diamonds) and Cassata et al 2011 (green squares). At redshift 3.0 we show the LF observed by Cassata et al 2015 (green squares) and Ouchi et al. 2008 (blue dots). At redshift 5.7 and 6.7 we show the LF observed by Ouchi et al. 2008 (blue dots) and Konno et al. 2018 (purple diamonds).

3 Implementing  radiative transfer in a semi-analytical model.

In this section we describe how we incorporate the  radiative transfer processes inside the semi-analytical galaxies from GALFORM. We make use of the fitting formula described above to predict the  escape fraction and line profiles. The strategy to fit the value of the free parameters of Eqs. 7 and 8 is described below.

3.1 Calibrating the model.

In order to calibrate the model and compute the values of the free parameters for each geometry, we fit our model to the observed LAE luminosity function at redshifts and . We run emcee (Foreman-Mackey et al., 2013b) to perform an MCMC to find the values of . The dynamical range of each free parameter is determined by limiting the expansion velocity and column densities of each component to lie within and for at least 90% of the resulting galaxy population with  rest frame equivalent width Å and  luminosity . These limits are imposed by the range of validity of the fitting formulae to derive the escape fraction (see §2.2).

This calibration is done independently for each outflow geometry and individual redshift bin. To combine multiple observed LFs at redshift 2.2 and 3.0 we compute a 5th-order polynomial fit (in logarithm of  luminosity - logarithm LF space) taking into account the uncertainties of each survey to obtain a single curve that represents the observational measurements. We choose to use a 5th-order polynomial at these redshifts as some recent works suggest that the typical Schechter function is not able to reproduce the observe LF (Konno et al., 2016; Sobral et al., 2017). Additionally, at redshift 5.7 and 6.7 we use the best fitting Schechter function to the observed LAE LF computed by Konno et al. (2018). The LF used to calibrate our model are shown in Fig.3 in black dashed lines.

The model  luminosity of galaxies, for each geometry and choice of , , , is computed as follows: i) we compute the intrinsic  luminosity of each component, , of each galaxy, which is directly proportional to the ionizing photon production predicted by GALFORM; ii) we compute , and using Eqs. (8) and (9); iii) we obtain for each galaxy component using Eq. (4); iv) the observed  luminosity of each component is obtained by multiplying the intrinsic luminosities by their respective ; and v) the total  luminosity for each galaxy is the sum of the observed luminosity of each component (disk bulge).

Fig. 3 shows the observed LAE LF (points), the full GALFORM intrinsic  LF (thin black line), the predictions for each geometry (thick colored lines) using the free parameters that result from the MCMC (listed in table 2) at the different redshifts implemented in this work.

The intrinsic  LF in divided into two populations: normal SFR galaxies in the low luminosity range and starburst galaxies in the high luminosity range. In general, in GALFORM the galaxy disk component in dominated by a quiescent SFR while in bulges the main mode of star formation is starburst, although quiescent star formation is also included. Additionally, in GALFORM the quiescent SFR and the starburst have different IMFs, which produces the bumps in the LF. On one hand, at lower redshifts, the predicted intrinsic LF is above the observations at all luminosities, thus galaxies at these redshifts require a significant in order to reduce the amplitude of the LF. On the other hand, at redshifts 5.7 and 6.7, the intrinsic LF at low (disk-dominated region) matches observations, implying that galaxies in this range must have . Additionally, the intrinsic high redshift LF at high luminosities (bulge-dominated regime) requires .

In general, the MCMC approach finds good matching solutions for the models including the  radiative transfer. First, we find that the Thin Shell is consistent with the measured LF at at all redshifts. Secondly, the Wind geometry performs quite well at and 5.7 while at it underpredicts the number density of LAE. However, we have checked that by allowing to be slightly higher, the observed LF is matched at redshift 6.7 as well. In the third place, the Bicone geometry matches the observed LF at and 3.0 while at and 6.7 it fails. The low abundance of LAEs predicted with the Bicone geometry arises due to the low escape fractions predicted by this geometry. In fact, at high redshifts, faint  emitters require escape fractions close to 1 to match the observed LFs, and this is not possible in the Bicone geometry by construction, as shown in Fig. 2.

3.2 A simplified model with no  radiative transfer

In order to highlight how radiative transfer changes the properties of LAEs, we compare the properties of our model with an abundance matching approach. We perform a simple SFR- mapping where no  radiative transfer is taken into account. We refer to this model variant as ’AM-noRT’.

To construct the AM-noRT model, we rank galaxies by their SFR. We assign a  luminosity to each galaxy based on their total SFR in a monotonic way. Objects with the highest SFR are assigned the brightest  luminosity. We compute the  equivalent width using the assigned  luminosity and continuum luminosity around the  frequency provided by GALFORM.  luminosities are assigned recursively towards lower luminosities such that the Ly observed luminosity function (using the cut of each survey) is recovered at each redshift. The resulting  luminosity distribution is shown in Fig. 3 as dashed black line. We compute a , which corresponds to the ratio between the assigned luminosity and the intrinsic one.

In contrast with our RT models, the  in the SFR-only model does not depend on properties such as the cold gas mass or the galaxy metallicity. Due to the way that  luminosities are computed, the resulting  can be higher than 1 in some cases.

4 Results.

In this section we describe the main predictions of our radiative transfer model when applied to GALFORM with the different outflow geometries.

4.1 The and distributions.

Figure 4: Ouflow expansion velocity and neutral hydrogen column density distributions for each redshift ( and from left to right) and for each geometry color coded as stated in the legend. The dark and light shaded contours enclose the 40 and 80 percentiles of the galaxy population, respectively.

Since the parameters in our model are calibrated to match the observed LFs for each geometry independently, the resulting distributions of and are different for each configuration. Though this work unless it is different stated, we define LAE as a galaxy with a  restframe equivalent width Å as typically in the literature (e.g. Ouchi et al., 2018). In this section, we use a subsample of the full LAE population obtained from each model by imposing a number density cut in  luminosity of .

Fig. 4 shows the distribution of and for each geometry. Since each quantity is computed for the disk and bulge component of each galaxy separately, we weight each component by their observed  luminosity to build the distributions shown in Fig. 4.

Overall, the distribution is relatively compact at redshifts (z=2.2,3.0) and more extended at higher redshift (z=5.7,6.7). The Thin Shell tends to have lower and than the Wind geometry. Additionally, there is a strong difference between low and high redshift for these two distributions, while, in the case of the Bicone, remains generally unchanged across cosmic time. Additionally, most of the galaxies lie within the  analytic expression optimal accuracy region defined in §3. Moreover, we have checked that the fraction of galaxies outside the this region is lower than a 7% for every geometry and redshift.

Typical values the are found to be around and for the Thin Shell and Wind geometries respectively at . Meanwhile, is found at for the Thin Shell and for the Wind. Notably, at higher redshifts, and 6.7, the distributions acquire a ’V’ shape (especially visible for the Thin Shell) due to the division of each GALFORM galaxy into a disk and bulge and the significant difference in  for starburst and normal SFR galaxies at these redshifts. Lower column densities are favored by disk-dominated galaxies, requiring a higher  in order to fit the LF. the distribution of these galaxies peak around and . Bulge-dominates starbursts require a lower  to fit the LF, thus they favor high and low distributions centered around and respectively.

The Bicone geometry displays noticeable differences with respect to the other two geometries. The Bicone distributions are very similar across the different redshifts used in this work and present the available highest and lowest distributions (peaking around and respectively), maximizing as much as possible the escape of  photons. This is due to the fact that the typical  is always lower in the Bicone compared to the other geometries, and it never reaches . Thus, the  LF with this geometry is not able to fit the observed LF, as shown above.

4.2 Breaking down the  Lf

Figure 5: Left panels : Break down of the  LF or galaxies with  Å in bins of halo mass (upper left), stellar mass (upper right), star formation rate (lower left) and metallicity (lower right). The bins are indicated in the legends. In each quantity the bins are represented in lighter colors for low values and darker as they increase. The total LF is plotted in thick gray line. Right panels : The probability distribution function of the different properties. In black we show the bin cuts.

To illustrate the properties of LAEs, Fig.5 shows the  LF obtained with the Thin Shell geometry at , split by the contribution of different ranges of halo and stellar mass, star formation rate and gas metallicity. We note that other redshifts and geometries show a similar behavior to what is shown in Fig. 5. Here we are analyzing a subsample composed of every LAE (Å) with  luminosity .

When splitting the LF based on the halo mass of LAEs (upper-left panel), we find that the majority of LAEs are hosted by haloes of moderate mass, which dominates the bright and moderate luminosities. LAEs with host halo masses below dominate the very faint end of the LF, with . Finally, the most massive haloes host galaxies do not contribute significantly to the LF shape. Furthermore, we have checked that there is no clear correlation between halo mass and  luminosity.

In the upper right panel in Fig. 5 the LF is split according to the stellar mass of the emitting galaxy. The whole body of the LF is dominated by LAEs with stellar mass about . Moreover, galaxies with a very low () or a very high () stellar mass do not contribute to bright or the faint ends. As in the case, we do not find any clear correlation between stellar mass and  luminosity.

The star formation rate, as expected, contributes in a roughly monotonic way to the  LF. The faint-end of the  LF is dominated by galaxies with low . Additionally, the intermediate luminosities are dominated by moderate while the bright end is populated by galaxies with the highest (although with a significant scatter). Note that this trend only means that the of LAEs scales with SFR, but not that every galaxy with high SFR would result in a LAE. Finally, we note that typically, galaxies with do not contribute to the LF.

The break down of the  LF in terms of gas metallicity is less intuitive. Naively one would expect to find an anti-correlation between metallicity and  luminosity, since  decreases with increasing dust, and thus, metallicity. However, we find the opposite: for LAEs with , the low metallicity bins contribute to the lower luminosities and vice versa. This trend is broken for due to the low  at this metallicity range. The galaxies with highest do not contribute anymore to the bright end but to low and average luminosities . This leads to the bulk of the  emitter population being dominated by galaxies with average metallicities, spanning the range . We dig deeper in this relation in §4.3

4.3 The bulk properties of LAEs.

Figure 6: a) The stellar mass - halo mass distribution at The gray shaded region shows the distribution for the full GALFORM sample. The solid yellow and blue lines and correspond to the median of GALFORM central galaxies disk and bulges properties respectively. The shade regions show the 10-90 percentiles. The red dots show the Thin Shell LAE sample median, 10-90 percentiles (vertical) and the bin size (horizontal). b) Same as a) but for the stellar mass - metallicity distribution. c) Same as a) but for the stellar mass - star formation distribution. The top panels show the distributions of the halo mass, star formation and metallicity, respectively, for the full GALFORM (yellow and blue for disk and bulge dominated respectively) and the Thin Shell model (red). The stellar mass distribution is shown in the right vertical panel.
Figure 7: The   as a function of SFR (left panels) and stellar mass (right panels) at . Gray points are from Oyarzún et al. (2017). Each panel displays our model predictions with a different outflow geometry, as shown in the legend. The bottom-right corner displays the predictions of the model with no radiative transfer. The solid line in each panel is the median of  predicted by our models. The dark and light coloured shaded regions display the and percentiles of the models predictions, respectively.

In this section we analyze the galaxy properties of our simulated LAE, focusing on the results at redshift and for the Thin Shell geometry (we checked that different geometries and redshifts give similar results). We restrict our analysis only to central LAEs with a number density cut in  luminosity (we check that different number density cuts produce similar results), and we compare it with the properties of the underlying population of central galaxies, i.e., the full population of galaxies predicted by GALFORM with .

Figure 6 shows some physical properties of the LAEs (red dots) and for the general population of galaxies from GALFORM selected using the same number density cut as the LAEs (yellow for disk properties and blue for bulge properties). Each panel includes the distribution of halo mass , star formation rate , metallicity and stellar mass and the correlation between , and .

The distribution in the LAE sample peaks at intermediate and spans between . LAEs halos trace the massive end of the disk-dominated distribution while avoiding the most massive dark matter halos, even if they host the strongest starburst episodes. This is caused by the predicted by GALFORM that associates high metallicites (low ) to high .

The metallicity and the of the LAE sample behave in a similar way due to the tight relation. The bulk of the LAE sample peaks at intermediate values of and , avoiding the extremes of the full GALFORM distribution. In particular, the galaxies with the highest SFR are not selected as LAE as the metallicity is also too high, causing a lower . Additionally, the galaxies with extreme low are not selected either as their in too low in these galaxies.

The relations (Fig. 6) for disk and bulge-dominated galaxies behaves in the same way. On the other hand, in the LAE sample this relation is the same as in the underlying galaxy population up to the peak of the and distributions, where the relation flattens for higher halo masses. In the high halo mass regime, LAEs typically have lower stellar masses than the overall average. This behavior is given by the tight relation causing to be lower for galaxies with higher as they become more dust rich.

In the LAE sample, the relation is consistent with the bulk of the disk-dominated galaxies for . After a transition around , is consistent with starburst galaxies. At metallicities below that transition the LAE relation is slightly above the overall relation.

In the LAE sample the relation is below the full GALFORM relation. This implies that for a fixed stellar mass, galaxies with higher SFR are selected, as the intrinsic correlates directly with the SFR.

4.4 The predicted   against observational estimates

Figure 8: Top: the halo occupation distribution (HOD) at redshift 2.2 , 3.0 , 5.7 and 6.7 from left to right. Model with radiative transfer show as blue, red and green solid lines for the Wind, Bicone and Thin Shell geometry respectively. The LAE sample AM-noRT is plotted as dashed black line. Bottom: fraction of galaxies that are considered LAE times the bias of the hosting dark matter halo. This quantifies the contribution of the different to the overall bias of the population.

In this section we compare our model predictions for the  against observational estimates from Oyarzún et al. (2017) at . In order to mimic their sample selection function we select galaxies with , and .

Fig. 7 shows the relation between the   and the SFR and stellar mass. The  computed in Oyarzún et al. (2017) displays a noticeable anti-correlation between SFR and  . In the models including RT galaxies with higher SFR have lower values of , in remarkable agreement with the observational estimates. The scatter in the observational data of Oyarzún et al. (2017) is consistent with the spread predicted by our models. This anti-correlation is caused by intrinsic link between SFR and . Even if the is higher for greater SFR (equation 7), dust plays the mayor role in the escape of  photons and reduces .

The stellar mass is also anti-correlated with the , as shown in the right panel of Fig. 7. This is due to the known correlation between and . Although our models reproduce the observationally inferred trend, the stellar masses predicted by GALFORM are systematically larger by .

Interestingly, the abundance matching model AM-noRT does not display the same trends found in Oyarzún et al. (2017), highlighting the importance of considering radiative transfer effects to predict LAE galaxy properties consistent with observational datasets.

4.5 The dark matter haloes hosting LAEs

Figure 9: Top panels : Monopole (3D auto-correlation function) for the AM-noRT sample (black), for the Thin Shell (green), galactic wind (blue) and biconical galactic wind (red) for redshift 2.2 , 3.0 , 5.7 and 6.7 from left to right. Middle panels : The ratio between the different LAE sample and the dark matter correlation function. Bottom panels : relative difference between the different samples and the AM-noRT monopole correlation function.

In the following we study the properties of dark matter halos hosting LAEs. To compare different model predictions, we select the brightest LAEs with a number density cut of .

Fig. 8 shows the halo occupation distribution (HOD) at and . This is constructed by computing the mean number of galaxies within different halo mass bins. All models including radiative transfer display a similar HOD at and . Central galaxies have a peak abundance in haloes of mass . Satellite galaxies start dominating the abundance of haloes of mass . None of the HODs at these redshifts reach . Even at the peak of occupation, less than of haloes host a LAE, regardless of radiative transfer effects.

At the HOD of the Bicone model falls significantly below that from the Thin Shell and Wind models. This reflects the differences in the LFs at these high redshifts. As the Bicone model is not able to reproduce the observed LF the resulting LAE population have quite different properties to the other RT samples.

The model with no radiative transfer systematically places LAEs in higher mass haloes compared to the radiative transfer models at low redshift. The occupation peak for centrals in the AM-noRT model is shifted to slightly more massive halos at and 3.0. Additionally, at these redshifts, the occupation of dark matter halos with is much greater in the AM-noRT model than in the models including RT. At redshifts and 6.7, the trend is inverted as LAE (Thin Shell and Wind geometry) populate halos slightly more massive than the AM-noRT model. Also, the occupation of halos with is greater in the RT models.

The bottom panels of Fig. 8 show the quantity


where is the number of sources in our LAEs samples in a halo mass bin, is the number of galaxies in the same bin and the galaxy bias is defined as


where and are the two point correlation functions for the galaxies and dark matter. This exhibits the contribution of different mass bins to the overall clustering bias of the LAE population. There is an evolution in the that contributes to the bias, being greater at lower redshifts and lower at higher redshift. In particular, the peak values varies from at to at

At low redshift ( and 3.0) the greater contribution to the bias come from lower mass halos in the RT models than in the AM-noRT model. However, this trend is inverted at z=5.7. Additionally, at the main contribution to the bias comes from the same halo mass for all the models.

4.6 The clustering of LAEs.

In this section we study how  radiative transfer impacts the clustering of LAEs for each of the outflow geometries implemented. The sample used in this section is the same as the one used in §4.5.

In Fig. 9 the top panel shows the spherically-averaged 2-point auto-correlation function (2PCF) in real space at and . The middle panel shows the bias, defined as in Eq.10. Moreover, in order to highlight the differences in the RT samples and the AM-noRT we show in the bottom panel of Fig.9 the relative difference of the 2PCF of the LAE samples and the AM-noRT, i.e., , where is the AM-noRT 2PCF.

Overall, the clustering amplitude increases towards higher redshifts regardless of the LAE model variant. In detail, each model predicts a slightly different clustering bias. There is a strong scale-dependence of the clustering bias in all models and at all redshifts for separations below .

At and the clustering amplitude of the AM-noRT sample is about 10% above the one predicted by the RT models. This is a consequence of LAEs being hosted by higher mass dark matter halos for this model, as shown in previous sections. At and , the clustering amplitude of the Thin Shell and Wind LAE samples are above that of the AM-noRT and Bicone models. Interestingly, as shown in the bottom panels of Fig. 9, towards redshifts the AM-noRT sample features a slightly different slope with respect to the RT models.

In summary, the predicted clustering of LAEs at is overall slightly lower when radiative transfer is included, and slightly higher towards . The relative differences in the amplitude of clustering, with respect to the AM-noRT model, are of the order of . These differences result from the non-trivial relation between the  luminosity of galaxies and the dark matter halo population hosting these objects.

4.7 The clustering in mock catalogs of LAE surveys

Survey Thin shell Wind Bicone AM-noRT
Kusakabe et al. (2018) 2.2 0.0773 0.93 104.9 93.6 448 1248
Bielby et al. (2016) 3.0 0.0633 1.07 60.0 119.1 468 643
Ouchi et al. (2018) 5.7 0.0954 7.67 43.5 401.5 18 734
Ouchi et al. (2018) 6.7 0.1078 21.2 41.0 696.5 19 873
Table 3: Mock catalog characteristics including the redshift , the redshift width , sky coverage (Area), the size along the line of size , the distance perpendicular to the line of sight , the number of mocks sliced from the simulation box and the median number of LAEs the mocks with the 32 and 68 percentiles.

In this section we compare our clustering prediction against several measurements of the clustering of LAEs at different redshifts from Kusakabe et al. (2018) at , Bielby et al. (2016) at and Ouchi et al. (2010, 2018) at and , respectively. We build LAE mock catalogs mimicking the properties of the different surveys to allow a close comparison with the observational datasets. These surveys use narrow band photometry to detect LAEs over a restricted redshift range. The main difference in the mock catalogs comes from the specific area, flux depth and equivalent width limit () of the individual survey.

Figure 10: Comparison between different model mocks (Thin Shell, Wind, Bicone and AM-noRT in rows from top to bottom) and the observed 2-point projected correlation function (Kusakabe et al., 2018; Bielby et al., 2016; Ouchi et al., 2010, 2018) at redshifts 2.2, 3.0, 5.7 and 6.7 in each column from left to right. The observational data is shown by dots and the best fitting power law extracted from their original work are plotted as dashed black lines. The solid lines correspond to the median for the mocks and the darker and lighter shades to the 32-68 and 5-95 percentiles respectively.

To build the mock catalogues, we choose a direction as line of sight (LoS). Assuming a distant observer, a galaxy coordinate is transformed in redshift space using


where is the galaxy coordinate along the LoS, is the galaxy peculiar velocity along the LoS and and are the scale factor and the Hubble parameter, respectively, at the  pivot redshift, , of the NB filter. Additionally, we conserve the periodicity of the box along the LoS direction.

Although some surveys have complicated footprints due to multiple pointings, our mocks are constructed as squares comprising an area equal to that of the target survey. Thus, the simulation box is simply split in slices along the LoS. The size perpendicular to the LoS is computed as


where is the survey sky coverage. The thickness (along the LoS) of the slice is computed as


where is the comoving distance at the geometric redshift . Additionally,


where and are the pivot wavelength and the full width half maximum of the narrow band filter and is the Lyman wavelength.

We calculate the limiting luminosity and the minimum rest frame equivalent width for each survey by matching the LAE number density, of the surveys to the one in the whole simulation box (see Appendix B). Then, our mock catalogs consist of galaxies with luminosity above and above . Table 3 lists the properties of the mocks, including the parallel and transverse sizes along the LoS, the redshift window , the number of mocks, , sliced from the simulation box and the number of LAE in each survey, and the median with 32-68 percentiles of the number of LAEs in the mocks.

The value of for narrow-band surveys is typically very small compared to the box length of the simulation. This allows for a big fragmentation of the simulation box along the LoS. On the other hand, can vary significantly between surveys. While, at low redshift () is relatively small and allows a large number of mock surveys, at only one cut is possible due to the large size required for the mock surveys. As a result of this, the number of mocks at (448 and 468 respectively) is much larger than that at (18 and 19 respectively).

Since in the simulation box is set to match the observed of each survey (see Appendix B), the observed number of LAE and the median number of LAE in our mocks, , are compatible within 1 sigma. Additionally, the dispersion of is higher (lower) at and 3.0 (5.7 and 6.7), since the comoving volume is smaller (larger). Hence, the impact of cosmic variance on clustering measurements is stronger (weaker).

We construct mock catalogs of LAE surveys from Kusakabe et al. (2018) at , Bielby et al. (2016) at , Ouchi et al. (2010) at and Ouchi et al. (2018) at . Figure 10 shows the comparison between the observed angular 2-point correlation function of these surveys, , and that computed from the mock catalogues, .

Overall, is very similar among our different model variations, including the AM-noRT model. The differences in the clustering due to the different bias of the samples are small in comparison with the scatter due to cosmic variance, making all models indistinguishable from each other.

At redshift 2.2 there is a good agreement between the the mocks and the clustering measurements in Kusakabe et al. (2018). At is significantly below the . However, the slope of the different samples are very similar to observations. At higher redshifts the LAE clustering predicted by the mocks is overestimated in our models. In particular, at , for angular distances , overestimates the clustering, while at larger the mocks match very well . Additionally, at redshift 6.7 the bias is significantly (about 2-sigma) overestimated in comparison with . This discrepancy could be caused by multiple reasons. The moderate contamination of interlopers () in the Ouchi et al. (2010) sample could decrease the measured clustering amplitude. Also, the observed LAE population at this redshift might contain a significant contribution of objects at the mass resolution limit of our simulation (), thus making our predictions biased towards higher masses and clustering amplitudes.

5 Discussion.

Here we discuss some of the results found in previous sections. In particular, in subsection 5.1 we discuss how the different outflow geometries impact the predicted properties of the LAE populations. Then, in subsections 5.2 and 5.3 we discuss the limitations of our methodology.

Model Thin Shell Wind Bicone AM-noRT
2.2 Thin Shell 1.000 0.814 0.555 0.229
Wind 0.814 1.000 0.592 0.189
Bicone 0.555 0.592 1.000 0.197
AM-noRT 0.229 0.189 0.197 1.000
3.0 Thin Shell 1.000 0.805 0.401 0.188
Wind 0.805 1.000 0.427 0.151
Bicone 0.401 0.427 1.000 0.227
AM-noRT 0.188 0.151 0.227 1.000
5.7 Thin Shell 1.000 0.413 0.798 0.108
Wind 0.413 1.000 0.322 0.063
Bicone 0.798 0.322 1.000 0.160
AM-noRT 0.108 0.063 0.160 1.000
6.7 Thin Shell 1.000 0.354 0.663 0.104
Wind 0.354 1.000 0.229 0.076
Bicone 0.663 0.229 1.000 0.259
AM-noRT 0.104 0.076 0.259 1.000
Table 4: Fraction of shared galaxies between pairs of models at the same redshift.

5.1 Differences between the RT models.

In this work we have used three different gas outflow geometries (Thin Shell, spherical galactic wind and biconical galactic wind) to model the  radiative transfer inside galaxies. The galaxy properties predicted for LAEs are very similar. The only significant difference between the predictions of different geometries is on the required distributions of column density and expansion velocity.

In Table 4 we list the fraction of galaxies shared by pairs of LAE models imposing Å and a number density cut of in . We find that the Wind and Thin Shell geometries share a high fraction of galaxies () at redshifts 2.2 and 3.0. However, at high redshift these geometries select different galaxies as the shared fraction in relatively low ( overlap). This might be due to the fact that there is a necessity of and the recipes to compute and are different. However, quite the opposite relation is seen between the Thin Shell and Bicone, as at low redshift they share a relatively low percentage of galaxies () and this increase at higher redshifts ().

Finally, when comparing the galaxies in the Wind and Bicone geometry we surprisingly find a low overlap between them. In particular, the maximum overlap happens at () and it drops down to only at . This shows the impact of the gas geometry on how the RT shapes the LAE selection function; even though the intrinsic galaxy population and the recipes to derive and are the same, the two geometries predicts different populations (although with similar characteristics).

We conclude that the RT LAE samples, in general, share a big fraction of galaxies () although the implemented gas geometries are very different. This is due to the fact that  behaves similarly for all of them. In particular, even if the exact dependence is different for each geometry, decreasing , increasing and decreasing increase   thus the visibility of the object for all of them. This makes the RT LAE samples very similar, as galaxies with properties that maximize and  are selected.

5.2 Limitations of the simple AM-noRT model.

We have also used a very simplistic LAE model where were radiative transfer effects are not taken into account and depends monotonically on the SFR. In Table 4 we also list the overlap between the radiative transfer and AM-noRT LAE sample. We find that the fraction of galaxies shared between the AM-noRT and RT catalogs is low, reaching its maximum value at z=2.2 () and then decreasing to at redshift 6.7.

As shown in Fig. 7 the AM-noRT sample fails to match not only the observed and relations but also the overall trend where anti-correlates with these two properties due to the RT (as described above). Additionally, the dark matter halo population, and thus the clustering, is different in comparison with the RT samples.

This work highlights the importance of taking into account the  RT inside galaxies when modeling LAEs. In particular, unlike in RT LAE samples, the galaxy properties model AM-noRT differ from observations, making them less attractive to study galaxy formation and evolution.

5.3 Limitations of the RT models.

The IGM also plays a mayor role in the detectability of galaxies based on  flux (Dijkstra et al., 2006; Zheng et al., 2011; Behrens et al., 2017). The IGM opacity becomes more important at higher redshifts () where the universe is denser and colder. However, the IGM might already also have an impact on the LAE selection function at as, even if the universe is highly ionized, the cross-section of neutral hydrogen atoms for scattering  photons is very high. The IGM impact might alleviate some of the tension that we find when we compare LAE models with observations. We will implement the effect of the IGM opacity in future work.

In Fig. 7 we found that, although the observed relation is perfectly reproduced by our RT models, the relation is not. Even if the overall trend is similar, we find a significant difference (about 0.5 dex) in the stellar mass. This is probably not caused by our implementation of RT in a semi-analytic model, but by GALFORM itself, as we note that full GALFORM  relation at redshift 3.0 is overestimated (also about 0.5 dex) in comparison with the observed one (Behroozi et al., 2010). Another possible source for this discrepancy is the different stellar population synthesis models used by Oyarzún et al. (2017) and GALFORM.

Another limitation of the RT models is that they predict very similar galaxy properties for the three different geometries. This degeneracy makes it difficult to determine from observations which geometry is the one driving the  photons escape. Nonetheless, the three gas geometries used in this work have very different  line profiles (as shown in figure 2) which might break the degeneracies and lead to a better understanding of the escape channels of  radiation. We will implement line profiles in a upcoming work.

6 Conclusions and future work.

Lyman- emitters are a promising galaxy population to trace the large scale structure of the Universe at high redshifts, . One of the main advantages of LAEs is their high luminosity at the  rest frame wavelength, making them easy to detect. Additionally, due to the Hubble expansion (Hubble, 1929), the  line is observable in the optical from to , allowing ground-based measurement of these galaxies. However, their selection function is quite complex as it depends upon  radiative transfer, which is sensitive to local astrophysical conditions.

We have designed a theoretical model of LAE based on a Monte Carlo Radiative Transfer code that can be applied to huge cosmological volumes. In particular, we have applied our model the N-body only-dark-matter simulation P-Millennium and the semi-analytical model of galaxy formation and evolution GALFORM (Lacey et al., 2016).

Monte Carlo Radiative Transfer codes have demonstrated to be a powerful tool to understand how  photons escape from galaxies. Unfortunately, the high computational cost prohibits the capability of being directly run over cosmological volumes. In order to avoid this problem we have developed analytical expressions for the  escape fraction  that are quite accurate for a wide range of outflow expansion velocities , neutral hydrogen column densities and metallicities .

Our methodology computes  for each galaxy as a function of , and , which characterise the gas outflows from which  photons escape. We compute these quantities using galaxy properties such as the size, SFR or halo mass. Free parameters to compute these quantities are chosen to fit the observed luminosity function over a wide range of redshifts. After calibration we find that every geometry reproduces well the observed LAEs LF at low redshift while only the Thin Shell and Wind manage to match them at high redshift. We conclude that our Bicone geometry (as described in this work), at high redshift, is less favoured with respect to the others.

We have analysed the relative abundance of  emitters by breaking down their LF in terms of several properties. Halo or stellar masses are not significantly correlated with  luminosities. The LF is actually mostly dominated by relatively low mass galaxies. However, when the LF is split in SFR bins we find a clear positive correlation with  luminosity. Finally, when the LF is divided into metallicity bins we find a scattered correlation for . Moreover, the contribution of high metallicities () to the bright end of the LF is small.

We also compared the properties of a  selected sample to the bulk of the galaxy population at high redshifts. We find that LAEs lie in relatively low mass halos. Additionally, the galaxies with the strongest starburst episodes are not selected as LAE since these galaxies typically have higher metallicities, and thus their  is low.

To validate our predicted , We have compared our LAE samples to the observational data from Oyarzún et al. (2017). We find a remarkable good agreement between our predictions and the observationally measured  - SFR relation. The LAE samples including RT reproduce successfully this anti-correlation and the scatter found between these quantities. However, the predicted  - plane is offset by in with respect to the data from Oyarzún et al. (2017). This difference can be due to the different assumptions about the stellar population synthesis models used by Oyarzún et al. (2017) and GALFORM, the impact of a different IMF in GALFORM, or simply that GALFORM predicts significantly more massive star-forming galaxies at these higher redshifts with respect to observational estimates. Finally, we find that our LAE AM-noRT sample based on assuming a monotonic relation between SFR and is not able to reproduce any of the observed trends. This highlights the crucial role of RT in shaping the LAE selection function.

We have also studied the dark matter halo population hosting LAEs in our models. We find differences between the samples including RT and the sample without RT. At low redshift, in comparison with the AM-noRT, the RT models predicts lower mass dark matter halos host LAE. This trend reverses at high redshift, as LAEs lie in more massive halos in the RT samples. We also find that the satellite fraction is low at all redshifts () and similar for all of the model variants.

The difference in the DM halo populations is directly translated into clustering discrepancies between the AM-noRT and RT samples. At low redshift, as a consequence of LAEs modeled with RT lying in lower mass DM halos, we find that they have a lower galaxy bias than the AM-noRT sample. This trend is reversed at high redshifts, when RT LAEs lie in more massive dark matter halos. Thus, we find that the RT models have a steeper galaxy bias evolution than the model excluding RT.

Finally, we have compared our model clustering predictions with observations finding some tension. While at redshifts 2.2 and 5.7 the observed clustering is well reproduced, at redshifts 3.0 and 6.7 the galaxy bias is poorly constrained. As studied in previous works (Zheng et al., 2011) the IGM transmission could have an impact on  selected samples that might alleviate this tension.

We have demonstrated the importance of RT in shaping the selection function of LAEs for galaxy properties as metallicity, SFR or DM halo properties. On one hand, the peculiar observational trends found can not be reproduce with a simple monotonic relation between SFR and . On the other hand, the inclusion of RT changes in a very particular way the clustering of  selected samples. All this make extremely important to construct models with  RT in order to understand the galaxy properties, formation and evolution of LAEs. Moreover, future surveys tracing the large scale structure of the Universe through LAEs will require a deep understanding of the channels through which  photons escape in order to obtain unbiased cosmological constrains.

In future work we plan to implement the transmission of  photons through the IGM, which is especially important at high redshifts. In order to do so we will develop analytic expression for the  line profile and a model to compute the IGM transmission in large cosmological volumes. These tools will enable us to explore how the IGM shapes the LAE galaxy properties and clustering.


The authors acknowledge the useful discussions with Zheng Zheng, Mark Dijkstra, Anne Verhamme in addition to the whole CEFCA team. We also thank Grecco Oyarzun for making their observational data available to us. The authors also acknowledge the support of the Spanish Ministerio de Economiaa 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 Data Centric system at Durham University, operated by the Institute for Com- putational Cosmology on behalf of the STFC DiRAC HPC Fa- cility (www.dirac.ac.uk). This equipment was funded by BIS Na- tional E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Opera- tions grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.


  • Ahn (2003) Ahn S., 2003, Journal of Korean Astronomical Society, 36, 145
  • Ahn (2004) Ahn S., 2004, ApJ, 601, L25
  • Ahn et al. (2000) Ahn S.-H., Lee H.-W., Lee H. M., 2000, Journal of Korean Astronomical Society, 33, 29
  • Barnes & Haehnelt (2010) Barnes L. A., Haehnelt M. G., 2010, MNRAS, 403, 870
  • Baugh (2006) Baugh C. M., 2006, Reports on Progress in Physics, 69, 3101
  • Baugh et al. (2005) Baugh C. M., Lacey C. G., Frenk C. S., Granato G. L., Silva L., Bressan A., Benson A. J., Cole S., 2005, MNRAS, 356, 1191
  • Behrens et al. (2017) Behrens C., Byrohl C., Saito S., Niemeyer J. C., 2017, preprint, (arXiv:1710.06171)
  • Behroozi et al. (2010) Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  • Benitez et al. (2014) Benitez N., et al., 2014, preprint, (arXiv:1403.5237)
  • Bielby et al. (2016) Bielby R. M., et al., 2016, MNRAS, 456, 4061
  • Bower et al. (2006) Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645
  • Cassata et al. (2011) Cassata P., et al., 2011, A&A, 525, A143
  • Chisholm et al. (2017) Chisholm J., Orlitová I., Schaerer D., Verhamme A., Worseck G., Izotov Y. I., Thuan T. X., Guseva N. G., 2017, A&A, 605, A67
  • Cole et al. (2000) Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
  • Dayal et al. (2011) Dayal P., Maselli A., Ferrara A., 2011, MNRAS, 410, 830
  • Dijkstra (2017) Dijkstra M., 2017, preprint, (arXiv:1704.03416)
  • Dijkstra et al. (2006) Dijkstra M., Haiman Z., Spaans M., 2006, ApJ, 649, 14
  • Dijkstra et al. (2011) Dijkstra M., Mesinger A., Wyithe J. S. B., 2011, MNRAS, 414, 2139
  • Foreman-Mackey et al. (2013a) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013a, PASP, 125, 306
  • Foreman-Mackey et al. (2013b) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013b, PASP, 125, 306
  • Garel et al. (2012) Garel T., Blaizot J., Guiderdoni B., Schaerer D., Verhamme A., Hayes M., 2012, MNRAS, 422, 310
  • Gawiser et al. (2007) Gawiser E., et al., 2007, ApJ, 671, 278
  • Granato et al. (2000) Granato G. L., Lacey C. G., Silva L., Bressan A., Baugh C. M., Cole S., Frenk C. S., 2000, ApJ, 542, 710
  • Gronke et al. (2016) Gronke M., Dijkstra M., McCourt M., Oh S. P., 2016, ApJ, 833, L26
  • Guaita et al. (2010) Guaita L., et al., 2010, ApJ, 714, 255
  • Guaita et al. (2017) Guaita L., et al., 2017, A&A, 606, A19
  • Harrington (1973) Harrington J. P., 1973, MNRAS, 162, 43
  • Hill et al. (2008) Hill G. J., et al., 2008, in Kodama T., Yamada T., Aoki K., eds, Astronomical Society of the Pacific Conference Series Vol. 399, Panoramic Views of Galaxy Formation and Evolution. p. 115 (arXiv:0806.0183)
  • Hu et al. (1998) Hu E. M., Cowie L. L., McMahon R. G., 1998, ApJ, 502, L99+
  • Hubble (1929) Hubble E., 1929, Proceedings of the National Academy of Science, 15, 168
  • Inoue et al. (2018) Inoue A. K., et al., 2018, preprint, (arXiv:1801.00067)
  • Kashikawa et al. (2006) Kashikawa N., et al., 2006, ApJ, 648, 7
  • Kobayashi et al. (2007) Kobayashi M. A. R., Totani T., Nagashima M., 2007, ApJ, 670, 919
  • Konno et al. (2016) Konno A., Ouchi M., Nakajima K., Duval F., Kusakabe H., Ono Y., Shimasaku K., 2016, ApJ, 823, 20
  • Konno et al. (2018) Konno A., et al., 2018, PASJ, 70, S16
  • Kulas et al. (2011) Kulas K. R., Shapley A. E., Kollmeier J. A., Zheng Z., Steidel C. C., Hainline K. N., 2011, preprint, 11074367 (arXiv:1107.4367)
  • Kusakabe et al. (2018) Kusakabe H., et al., 2018, PASJ,
  • Lacey et al. (2016) Lacey C. G., et al., 2016, MNRAS, 462, 3854
  • Laursen & Sommer-Larsen (2007) Laursen P., Sommer-Larsen J., 2007, ApJ, 657, L69
  • Laursen et al. (2009) Laursen P., Razoumov A. O., Sommer-Larsen J., 2009, ApJ, 696, 853
  • Laursen et al. (2011) Laursen P., Sommer-Larsen J., Razoumov A. O., 2011, ApJ, 728, 52
  • Le Delliou et al. (2005) Le Delliou M., Lacey C., Baugh C. M., Guiderdoni B., Bacon R., Courtois H., Sousbie T., Morris S. L., 2005, MNRAS, 357, L11
  • Le Delliou et al. (2006) Le Delliou M., Lacey C. G., Baugh C. M., Morris S. L., 2006, MNRAS, 365, 712
  • Malhotra & Rhoads (2002) Malhotra S., Rhoads J. E., 2002, ApJ, 565, L71
  • Nagamine et al. (2010) Nagamine K., Ouchi M., Springel V., Hernquist L., 2010, PASJ, 62, 1455
  • Nelson et al. (2015) Nelson D., et al., 2015, Astronomy and Computing, 13, 12
  • Neufeld (1990) Neufeld D. A., 1990, ApJ, 350, 216
  • Orsi & Angulo (2018) Orsi Á. A., Angulo R. E., 2018, MNRAS, 475, 2530
  • Orsi et al. (2008) Orsi A., Lacey C. G., Baugh C. M., Infante L., 2008, MNRAS, 391, 1589
  • Orsi et al. (2012) Orsi A., Lacey C. G., Baugh C. M., 2012, MNRAS, 425, 87
  • Osterbrock (1989) Osterbrock D. E., 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  • Ouchi et al. (2008) Ouchi M., et al., 2008, ApJS, 176, 301
  • Ouchi et al. (2010) Ouchi M., et al., 2010, ApJ, 723, 869
  • Ouchi et al. (2018) Ouchi M., et al., 2018, PASJ, 70, S13
  • Oyarzún et al. (2017) Oyarzún G. A., Blanc G. A., González V., Mateo M., Bailey III J. I., 2017, ApJ, 843, 133
  • Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  • Rhoads et al. (2000) Rhoads J. E., Malhotra S., Dey A., Stern D., Spinrad H., Jannuzi B. T., 2000, ApJ, 545, L85
  • Santos et al. (2004) Santos M. R., Ellis R. S., Kneib J., Richard J., Kuijken K., 2004, ApJ, 606, 683
  • Schaerer & Verhamme (2008) Schaerer D., Verhamme A., 2008, A&A, 480, 369
  • Shapley et al. (2003)</