Effect of Intergalactic Medium on the Observability of Lyman Alpha Emitters during Cosmic Reionization

We perform a systematic study of how the inhomogeneities in the Inter-Galactic Medium (IGM) affect the observability of Lyman Alpha Emitters (LAEs) around the Epoch of Reionization. We focus on the IGM close to the galaxies as the detailed ionization distribution and velocity fields of this region could significantly influence the scattering of Ly-alpha photons off neutral H atoms as they traverse the IGM after escaping from the galaxy. We simulate the surface brightness (SB) maps and spectra of more than 100 LAEs at z=7.7 as seen by an observer at z=0. To achieve this, we extract the source properties of galaxies and their surrounding IGM from cosmological simulations of box sizes 5-30 Mpc/h and follow the coupled radiative transfer of ionizing and Ly-alpha radiation through the IGM using CRASH-alpha. We find that the simulated SB profiles are extended and their detailed structure is affected by inhomogeneities in the IGM, especially at high neutral fractions. The detectability of LAEs and the fraction of the flux observed depend heavily on the shape of the SB profile and the SB threshold (SB_th) of the observational campaign. Only ultradeep observations (e.g. SB_th ~ 10^-23 ergs/s/cm^2/arcsec^2) would be able to obtain the true underlying mass-luminosity relation and luminosity functions of LAEs. The details of our results depend on whether Ly-alpha photons are significantly shifted in the galaxy to longer wavelengths, the mean ionization fraction in the IGM and the clustering of ionizing sources. These effects can lead to an easier escape of Ly-alpha photons with less scattering in the IGM and a concentrated SB profile similar to the one of a point source. Finally, we show that the SB profiles are steeper at high ionization fraction for the same LAE sample which can potentially be observed from the stacked profile of a large number of LAEs.


INTRODUCTION
predicted that high redshift galaxies would emit copious amounts of Lyα photons (∼ 10 per cent of the galaxy luminosity) due to the reprocessing of ionizing photons from the young massive stars by the interstellar medium (ISM). This stimulated a large number of studies aimed at observing the high-redshift galaxies which mostly lead to non-detections (e.g. Davis & Wilkinson 1974;Partridge 1974;Djorgovski & Thompson 1992, and the references there in). Models were made to reason the null de-E-mail: akila@mpa-garching.mpg.de tections such as absorption of Lyα radiation by dust (e.g. Spitzer 1978;Meier & Terlevich 1981;Hartmann et al. 1988;Charlot & Fall 1993) and lower number of massive stars due to ageing of stellar populations (Valls-Gabaud 1993).
The first detections of high redshift Lyα emitting galaxies (which we also refer to as 'Lyα Emitters', or LAEs) were reported in the mid-1990s (e.g. Hu & McMahon 1996). Since then a large number of LAEs 1 have been observed using both narrow band photometry and spectroscopy up to z 7 (e.g. Rhoads & Malhotra 2001;Dawson et al. 2004;Rhoads et al. 2004;Taniguchi et al. 2005;Iye et al. 2006;Nilsson et al. 2007Nilsson et al. , 2011Ouchi et al. 2009Ouchi et al. , 2010Wang et al. 2009;Guaita et al. 2010;Lehnert et al. 2010;Wolf et al. 2010;Cassata et al. 2011;Gronwall et al. 2011;Kashikawa et al. 2011;Pentericci et al. 2011;Ono et al. 2012;Schenker et al. 2012). In addition to being one of the main methods of detecting high-redshift galaxies to study the evolution of galaxy properties with redshift, LAEs can also be used to probe the epoch of reionization (EoR) (e.g. Miralda-Escude & Rees 1998; McQuinn et al. 2007) and to constrain dark energy properties using baryonic acoustic oscillations in the LAE clustering signal (e.g. Eisenstein et al. 2005;Hill et al. 2008). Since this paper is motivated by the use of LAEs for the detection of EoR, we will discuss this in more detail.
Due to its resonant nature, the Lyα line is scattered by even very small amounts of neutral hydrogen along its path. The optical depth τ due to scattering away from the line-ofsight to the observer (at z = 0) for a LAE at redshift zs is given by τ 6.02 × 10 5 xHI 1 + zs 10 for photons entering the intergalactic medium (IGM) with wavelength λ λα = 1215.67Å (Lyα rest wavelength; e.g. Gunn & Peterson 1965;Barkana & Loeb 2001) and by τ = 2.9 xHI ∆v 600 km s −1 −1 1 + zs 10 for photons entering the IGM with λ = λα (1 + ∆v/c) for an equivalent velocity shift of ∆v (e.g. Miralda-Escude & Rees 1998; Dijkstra & Wyithe 2010). Here, xHI is the neutral fraction of the IGM, c is the speed of light in vacuum and the cosmology used to derive the above equations and for the rest of the work in this paper is the standard model -Ω0,m=0.3, Ω0,Λ=0.7, Ω 0,b =0.04, h=0.7, n=1 and σ8=0.9. The effect of scattering on the shape of the intrinsic spectrum can be used to estimate the ionization fraction of the IGM around the source and thus study the EoR. LAEs luminosity functions (LFs), number density and clustering are the main methods used to study the EoR history. It is expected for example that an increment in the neutral fraction of the universe leads to a reduction of the observed number of sources, thus supressing the LFs towards higher redshifts (e.g. Haiman & Spaans 1999). Following this idea, Ouchi et al. (2010) used the drop in LAEs LF at z = 5.7 − 6.6 to constrain xHI 0.2 ± 0.2 at z = 6.6. A less stringent limit of xHI 0.5 has been found at the same redshift by Ouchi et al. (2010) using measurements of LAE clustering. These are marginally consistent with the estimate by Taniguchi et al. (2005) that at least 20-50 per cent of the volume of the universe needs to be ionized to account for the observed number of LAEs at z = 6.5. The use of number density (e.g. Malhotra & Rhoads 2006) and clustering properties (e.g. Furlanetto, Zaldarriaga, & Hernquist 2006;McQuinn et al. 2007) of LAEs have also been includes ultraviolet (UV) selected galaxies for which follow-up spectroscopy revealed the presence of a strong Lyα emission line.
extensively discussed in the literature as tools to probe the EoR.
A large number of authors have modelled LAEs at different redshifts using analytic (e.g. Haiman 2002; Kobayashi, Kamaya, & Yonehara 2006;Dijkstra, Wyithe, & Haiman 2007;Kobayashi, Totani, & Nagashima 2007;Samui, Srianand, & Subramanian 2009;Tilvi et al. 2009) and seminumeric (e.g. Le Delliou et al. 2005Delliou et al. , 2006Furlanetto, Zaldarriaga, & Hernquist 2006;Tasitsiomi 2006;McQuinn et al. 2007;Iliev et al. 2008;Dayal et al. 2009;Zheng et al. 2010;Dijkstra, Mesinger, & Wyithe 2011;) methods. It has been shown that properties like dust content (e.g. Kobayashi, Totani, & Nagashima 2007;Dayal et al. 2009) and gas inflows (e.g. Santos 2004;Dijkstra, Lidz, & Wyithe 2007) reduce the observability of LAEs, whereas top-heavy initial mass function (IMF; e.g. Malhotra & Rhoads 2002;, clustering of ionizing sources, 'patchy' reionization (e.g. Wyithe & Loeb 2005;Furlanetto, Zaldarriaga, & Hernquist 2006;McQuinn et al. 2007;Iliev et al. 2008;Mesinger & Furlanetto 2008a,b;Dijkstra, Mesinger, & Wyithe 2011) and gas outflows (e.g. Santos 2004;Dijkstra & Wyithe 2010) improve their observability at each redshift. Most works simulate the scattering of Lyα photons away from the line of sight by suppressing the intrinsic spectrum of an LAE by the optical depth calculated along the line of sight (see equations 1 and 2). Only a few studies perform a full 3D RT (RT) of Lyα photons (e.g. Cantalupo et al. 2005;Tasitsiomi 2006;Laursen & Sommer-Larsen 2007;Laursen, Sommer-Larsen, & Andersen 2009;Kollmeier et al. 2010;Zheng et al. 2010; Barnes et al. 2011;. In this work, we look at the structure in the IGM at z = 7.7, where the current constraints point towards nonnegligible neutral fraction in the volume averaged IGM, specifically focusing on the effects of the gas close to the object. This region is important because the inhomogeneities in the gas density play a crucial role in determining the size of the ionized region, especially at very low ionization fractions. The residual neutral hydrogen close to the source, together with the velocity field in this same region strongly affects the scattering pattern of the Lyα photons and their escape towards an observer. Here we investigate these issues for the first time in detail for a large (i.e. > 100) sample of objects. LAEs have been detected up to z ∼ 7 (e.g. Iye et al. 2006;Ono et al. 2012). Because observations in the near future are aiming for z ∼ 7.7 (e.g. Tilvi et al. 2010;Clément et al. 2012), in this study we have chosen to concentrate on objects at z = 7.7. Zheng et al. (2010) previously simulated a large sample of LAEs with Lyα RT at moderate spatial resolution in a highly ionized IGM at z = 5.7. Similarly, Laursen, Sommer-Larsen, & Razoumov (2011) studied scattering in the postreionization IGM using hydrodynamical simulations that had much higher spatial resolution for RT. Their work is complementary to our study.
In this work we simulate LAEs using galaxies from hydrodynamical cosmological simulations and RT using CRASHα (Pierleoni, Maselli, & Ciardi 2009). We make surface brightness (SB) maps of the objects and study the effect of IGM structure on the SB profiles. In the following Section, we describe the cosmological simulations and Section 3 briefly introduces the RT code CRASHα. Section 4 explains the pipeline used to simulate a statistically significant sample of LAEs. The results are described in Section 5 and the parameter study to understand the uncertainties in our work is shown in Section 6. In Section 7, we discuss the effect of IGM structure on the methods used in estimating the mean ionization fraction of the universe. Finally the conclusion are drawn in Section 8.

SIMULATIONS OF GALAXY FORMATION
To study high redshift galaxies and their surrounding IGM, we need simulations with a large box size to provide a statistical sample of halos in a wide mass range in different environments. At the same time, high resolution to resolve the halo and the structure in the surrounding IGM is required. A compromise was achieved by using medium range box sizes of 5-30 h −1 Mpc comoving. The simulations used in this paper are described in Maio et al. (2010), although additional ones have been run for different box sizes. The simulations were run using the TREE-PM SPH code Gadget-2 (Springel 2005) modified to include primordial Hydrogen, Helium and Deuterium based chemistry [e − , H, H + , He, He + , He ++ , H2, H + 2 , H − , HeH + , D, D + , HD] (Yoshida et al. 2003;Maio et al. 2007), stellar evolution and metal pollution ) and fine structure metal transition cooling (O, C + , Si + , Fe + ) at T < 10 4 K (Maio et al. 2007(Maio et al. , 2009. Metals produced by AGB stars and supernovae (SNII, SNIa) are spread by supernovae/wind feedback. The IMF is chosen according to the metallicity of the stellar particles, Z. We assume that a transition from metal-free/very metal-poor Population III stars to metalenriched, more standard Population II/I stars takes place whenever the gas reaches a critical metallicity Zcrit, which is determined by the cooling and fragmentation properties of the gas. Zcrit quoted by different authors range between 10 −6 to 10 −3.5 Z (e.g. Schneider et al. 2003;Bromm & Loeb 2003;Schneider et al. 2006;Santoro & Shull 2006). In these simulations we choose Zcrit = 10 −4 Z . Below Zcrit we assume a Salpeter IMF in the mass range [100,500] M , while for metallicities above Zcrit a Salpeter IMF in the mass range [0.1,100] M is used.
The simulations analyzed in this paper have box sizes of L = 5, 10, 20 and 30 h −1 Mpc comoving with 320 3 particles each in dark matter and gas. The mass of dark matter particles is 3.93 × 10 5 M (L/5 h −1 Mpc) 3 and the gas particle mass is 6.04 × 10 4 M (L/5 h −1 Mpc) 3 . The comoving softening length is 0.78 h −1 kpc (L/5 h −1 Mpc) 3 which is ∼ 1/20 of the mean inter-particle distance. The details of the simulations are given in Table 1.

THE RADIATIVE TRANSFER CODE CRASHα
CRASHα (Pierleoni, Maselli, & Ciardi 2009) is the first RT code for cosmological application where the propagation of Lyα and ionizing photons are coupled (refer to  for a similar code). The ionizing part is based on CRASH (Ciardi et al. 2001;Maselli, Ferrara, & Ciardi 2003;Maselli, Ciardi, & Kanekar 2009) which is a 3D ray-tracing grid-based RT code using Monte-Carlo (MC) techniques to sample the probability distribution of several parameters like spectrum of sources, emission direction and optical depth. The ionizing radiation is propagated through an arbitrary static H/He gas density field. Both radiation from point sources and diffuse radiation from ultraviolet background or recombination of H/He gas can contribute to the ionizing flux. The total ionizing photon count, Es [phot], from each source 2 s with photon rateṄion,s [phot s −1 ] over the whole time tsim of the RT simulation, is distributed among Np photon packets containing ionizing photons of different frequencies depending on the source/background spectrum. The emission of photon packets from each source happens at equally spaced time intervals dt = tsim/Np. The direction of emission is assigned by the MC sampling of the angular characteristic of the source. As the packets are propagated through each cell i, photons are absorbed according to the cell optical depth τ i c calculated combining the contributions from HI, HeI and HeII. The probability that a photon in a packet is absorbed is calculated using P (τ i c ) = 1 -e τ i c . This changes the temperature and ionization conditions of the cell and the photon distribution in a packet. The trajectory of a packet is followed till all the photons are absorbed or in case of non-periodic boundary conditions, till it reaches the edge of the box.
CRASH is modified to follow the time evolution of Lyα photons through an evolving ionization configuration of gas. A statistical approach to the Lyα RT is adopted using pre-compiled tables from MCLyα (Verhamme, Schaerer, & Maselli 2006). Lyα radiation is emitted by both point sources and recombination of ionized gas. Just like in the ionizing case, the Lyα photons produced by each point source s with photon rateṄLyα,s over the entire simulation time tsim is divided between Lyα photon packets. Np,Lyα Lyα photon packets are emitted by each source s at regular time intervals dtLyα = tsim/N em,l in random directions. N em,l is the number of times the Lyα photons are propagated during the RT simulation. When a packet enters a cell, depending on the conditions in the cell, the photons are either absorbed or scattered. In case of scattering, a new wavelength and the time of escape from this cell is obtained from the precalculated tables. The typical time of escape from a cell of size 1 kpc, temperature 10 K and number density 0.01 cm −3 giving an optical depth of ∼ 6 × 10 7 is about 3 × 10 5 years.
(e.g. Adams 1975;Dijkstra & Loeb 2008) Photons from the recombination in the gas are added to the packets. Dust is an important factor which determines Lyα radiation transport (e.g. Neufeld 1991). Dust optical depth is calculated as a fraction fτ of the Lyα optical depth. fτ is determined as fτ = m p/dust ×f H/dust ×σ dust ×ρ cell ×d cell where proton-to-dust mass ratio m p/dust = 5.0 × 10 −8 , gas-to-dust ratio f H/dust = 5 × 10 −3 (Verhamme, Schaerer, & Maselli 2006), dust absorption cross section σ dust = 2 × π × r 2 dust for a dust grain of radius r dust = 2.0 × 10 −6 cm. ρ cell is the gas density in the cell and d cell is the distance traveled within the cell by the photon. In high density regions, where Lyα photons have a high probability of scattering multiple times within the cell volume, the probability of dust absorption gets high, with the consequence that the effects due to dust scattering on Lyα RT are negligible. For example, assuming that the probability of scattering of Lyα photons by dust is equal to that of absorption, in 20 photon-dust interactions, only 1 in 10 6 photons would continue to scatter without being absorbed. Thus in high density regions, we neglect dust scattering.
For more details, we refer the reader to the original papers.

SIMULATING LYα EMITTERS
Our plan is to study how the observability of LAEs is affected by transmission through the IGM. For this, we calculate Lyα SB profiles of a large number of simulated LAEs at z = 7.7 covering a wide range of dark matter halo masses, thus sampling an equally wide range in IGM environments. The redshift was chosen as observational efforts are underway to investigate LAEs at z = 7.7 (e.g. Tilvi et al. 2010;Clément et al. 2012). As at this relatively high z the IGM is expected to be still substantially neutral, a full 3D RT approach is necessary to investigate the observability of LAEs.
The method adopted to simulate LAEs is to extract a cube around each dark matter halo from the snapshot of the simulation of galaxy formation, grid the density and velocity fields, get the spectrum produced by its stellar population from STARBURST99 (Leitherer et al. 1999) using the halo properties, input all the details to CRASHα along with the default values (refer to Section 4.3) for temperature and ionization fields and run the RT simulation to obtain Lyα photons escaping the simulation cube. We can use the details of the photons exiting each of the six faces of the cube to make SB maps. Each step is explained in detail in the following Sections.

Extracting Halos and Gridding
A snapshot of the hydrodynamical simulation is taken and the properties of the dark matter haloes are obtained using the Friends-of-Friends (FoF) routine (Davis et al. 1985). The centre-of-mass (CoM) of the dark matter (DM) halo, the DM, gas and stellar mass, and mean gas metallicity are used in the rest of the analysis. The mean gas metallicity is defined as the fraction of gas mass in metals. A summary of the DM halo properties for the four different simulations used in this paper are given in Table 2. In our study, we choose a subset of the resolved haloes (i.e. with mass MDM > 100 mDM; e.g. Trenti et al. 2010) from L05-30 to get a fair sampling of the DM mass range 10 9−11 M . This results in ∼ 130 haloes.
In our simulations, we resolve the structure in the IGM, but we do not properly resolve the ISM of the individual galaxies. To mimic the absorption/scattering of radiation within the ISM, we use the commonly adopted prescription based on the escape fraction of ionizing photons, fesc,ion,ISM, and the escape fraction of Lyα photons, fesc,Lyα,ISM (see Sec. 6.1). Thus we need to exclude from our RT calculations those cells which represent the ISM, to avoid accounting for their effect twice. We obtain this by removing all cells with a density larger than 0.03 cm −3 . Even values up to 0.05 cm −3 are acceptable, but might show slight dependence on the gridded density resolution especially for high resolution gridding in the case of small haloes (see discussion in the next paragraph). For the case of 0.03 cm −3 , the effect of density resolution is negligible in the range of simulations we deal with in this work. Values < 0.03 cm −3 give a similar ionization structure but lead to loss of some of the IGM gas. As a precaution to avoid removal of gas from high density regions not associated with the specific dark matter halo, we restrict the removal of cells to a distance of ∼ 0.7 × r200 from the source. The precise choice of this radial distance does not affect our main results, as long as it is less than the nearest massive dark matter halo in the cube.
A cube around the CoM of the halo is extracted with a side of 35 × r200, where r200 is the radius in comoving units at which the mean DM density inside the sphere is 200 times the critical density of the universe at that redshift. This is the smallest box containing the HII region produced by the galaxy at its equilibrium, so that we get the highest possible spatial resolution for the gridding of density and velocity Figure 1. Example of the extraction of a cube around a DM halo and the gridding of the gas density and velocity fields for CRASHα. Top: The gas distribution within a thin slice of the simulation box L05 centred on the CoM (shown as a green dot) of the most massive halo at z = 7.7 in comoving units of length. The comoving thickness of the slice is 10 h −1 kpc. The cube to be extracted around the object for the RT simulation is shown with a solid red square. Bottom: The gridded gas number density field (in logarithmic scale) corresponding to the cube in the top panel.
fields for our RT simulations. The exact value of the cube dimension though does not affect our main results. A cube size in the range 25 − 50 × r200 gives similar results, but the lower bound is too close to the edges of the ionized region, while larger boxes would have lower spatial resolution at fixed grid size, which is undesirable.
Once the cube is extracted, the gas density and velocity fields are gridded as input for CRASHα. For gridding, we use GadgettoGrid (Pakmor 2010). The code distributes the particle mass and peculiar velocities using an SPH kernel to 64 neighboring particles and then grids the fields. The default grid size used is 256 3 , which is set by memory and runtime constraints of CRASHα, although different grid sizes have been used for testing purposes. An important caveat is that, because the grid resolution is a factor of the virial radius of the halo, this results in a lower physical resolution for higher mass objects. Figure 1 shows a slice through the simulation box. The top panel refers to the gas number density in the box with the to-be-extracted cube around a halo marked by a solid red line and the CoM of the halo marked by a green dot. The extracted density field is gridded and shown in the bottom panel. Note that the gas density in the IGM is not uniform and the gridded density field has a range in values spanning many orders of magnitude.

Luminosity of Stellar Sources
As mentioned in the previous Section, we obtain the halo properties using the FoF and use the associated IMF, stellar mass and metallicity to calculate the corresponding ionizing and Lyα photon rate. To this purpose, we use STAR-BURST99 (Leitherer et al. 1999), which is a population synthesis code that provides the spectrum of a stellar population given an IMF, stellar mass, metallicity and age.
The lowest metallicity available in STARBURST99 is Z = 0.0004 ∼ 0.02 Z for Padova tracks. We choose the Padova original stellar evolution tracks instead of Padova AGB tracks as the Hubble age at z = 7.7 is less than 1 Gyr, allowing us to ignore the contribution of AGB stars to the spectrum. The average ionizing flux for the same stellar mass is mildly higher for lower metallicities, i.e. ∼ 17% higher for Z = 3 × 10 −3 Z (refer to Schaerer 2003 for low metallicity values). Since our metallicities are above Zcrit, we use a Salpeter IMF in the mass range [0.1,100] M in the STARBURST99 models, consistently with our cosmological simulations. The spectrum is scalable with respect to the stellar mass of the halo. Therefore, we choose M * = 10 6 M as the default stellar mass for the STARBURST99 model and normalize the photon count to this stellar mass. We also choose instantaneous star formation mode for computing the spectrum. The mean ionizing photon rateṄion was obtained by averaging the number of photons emitted over t = 2 × 10 8 years in the wavelength range [91, 912]Å. The time t was chosen as the cumulative number of ionizing photons reaches convergency. Thus, we obtain a normalized ionizing photon rateṄion = 7.77× 10 50 (M * /10 6 M ) phot s −1 . Only a fraction of these ionizing photons reaches the IGM, the rest is converted to Lyα photons in the ISM. Therefore, the ionizing photon rate reaching the IGM,Ṅ esc ion , is: where fesc,ion,ISM is the escape fraction of ionizing photons from the ISM. Our reference value is fesc,ion,ISM = 0.02 (Gnedin, Kravtsov, & Chen 2008). In Section 6.1 we discuss more extensively this choice. The calculation of the Lyα photon rateṄLyα is more complicated. There are three different components consid- NLyα =ṄLyα, +Ṅ Lyα, neb +ṄLyα, ISM. (4) NLyα, andṄ Lyα, neb are calculated from the spectrum by averaging over t at 1215.67Å with a line width of 2Å, which is the typical intrinsic line width of LAEs (Partridge & Peebles 1967). In STARBURST99 spectra, the nebular continuum is calculated by converting all the ionizing photons into nebular lines except for Lyα. Therefore,Ṅ Lyα, neb gives the upper limit of the contribution from other nebular lines to this wavelength range.
Our value ofṄion is a factor of two lower than the numbers quoted by e.g. Dayal et al. (2009)   in our simulations. Thus, their values would be comparable to ours for fesc,ion,ISM = 0.04 and fesc,Lyα,ISM = 0.6.

CRASHα Input and Output
The inputs for a CRASHα run are the gridded density, velocity and temperature fields (Sec. 4.1) and a source distribution with photon rate and spectrum (Sec. 4.2). In our reference simulations, we adopt uniform initial fields of ionization, xion = nHII/nH = 0, and temperature, T = 100 K.
A discussion on the different choices for these values is given in Section 6.3. The duration of the RT simulation is set to 2 × 10 8 years. The source locations and luminosity were calculated as explained in the previous Section. Other simulation parameters, Np = 10 6 , Np,Lyα = 10 4 and N em,l = 500, were determined using resolution tests.
Throughout the simulation, CRASHα collects photon packets exiting the box, recording their frequency, time, position and directional information, which can be used to quantify observed properties of the source such as SB maps and spectra. Each cell on a side of the simulation cube emits photon packets in different directions. Just like the photon counts, the true distribution of angular directions of the photons exiting from the cells are also sampled by the photon packets. Each direction points towards a different observer thus making each cell (3D) act as a pixel (2D) in the plane perpendicular to the direction of the observer. Here we define the direction of the observer as the one perpendicular to the plane of the side of the cube. By this definition we get six different observer directions for each source. Both SB maps and spectra can be calculated for each observer direction.
The SB SBem of a pixel of size ps due to a source of luminosity L at a distance d (in Euclidian geometry) is given by: The photons, after travelling through the expanding universe, arrive at the observer at z = 0. In a cosmological context, the observed SB SB obs of a source of luminosity L at redshift z = zs observed in a pixel of proper length ps is (Tolman 1934): where DL is the luminosity distance of the source and DA is its angular diameter distance. Therefore, to calculate the observed SB in a pixel for a specific observer, we need the number of photons in each cell travelling in the direction of the observer. The sampling of the angular distribution of photons depends on the number of photon packets in each cell. But the average number of photon packets in each cell of our simulation is too low (∼ 12) for a smooth sampling of the underlying angular directional distribution. Thus we calculate the angular directional distribution using all the photon packets on a side of the box, by assuming that all cells on that side sample the same angular directional distribution with respect to the true north for each cell, which is defined as the radial direction from the source to the cell. Figure 2 shows the sketch of the method. The star represents the cell with the stellar source. The dotted lines represent the radial direction connecting each cell on the edge of the box to the source (i.e. true north for each cell). The solid lines with arrows represent the directions of photon packets exiting such cells. The dashed lines from each cell to the observer represent the observer direction, which for zs = 7.7 is perpendicular to the side. The packets which are not scattered by the IGM follow the radial direction from the source (i.e. true north), while the scattered photon packets sample random directions. To calculate the SB obs we need to estimate the number of photons reaching the observer at z = 0, i.e. traveling along the dashed lines. To calculate the underlying angular direction distribution, we grid the directions of all photon packets (short solid arrows) in each cell, correcting for the different true north direction of each cell (dotted arrow). The gridded distribution of angles with respect to their true north in all the cells on a side are added and normalized using the total number of photon packets exiting the side. This normalized gridded angular distribution is then imposed on the photon count in each cell to obtain the number of photons in each pixel going in the direction of the observer.
The total energy of photons travelling in a given direction within a pixel in each second, L pixel , is given by: which is calculated by adding up the energy of all the photons in the pixel in the direction θ ≤ ps/d. d is the proper distance of the cell to the source (at the centre of the simulation box). L pixel is converted to the units of SB, i.e. erg s −1 cm −2 rad −2 , using the formula: This is then converted to the SB obs dividing by (1 + z) 4 as in equation (8). The SB obs is then converted to flux, F , as: where ps/DA is the pixel size of the observer in radians and DA is the angular diameter distance of the source. Adding up the flux in all pixels, we can obtain the observed source luminosity L obs as: where DL is the luminosity distance of the source. An important technical caveat is that our method for computing the SB differs from the often used 'next-event estimator' approach which has been adopted in other studies (e.g. Tasitsiomi 2006;Laursen, Sommer-Larsen, & Andersen 2009;Zheng et al. 2010;Kollmeier et al. 2010;Barnes et al. 2011). Idealized test cases though show that the method perfectly reproduces the SB profiles of the Rybicki-Loeb halo (Loeb & Rybicki 1999), and also those produced in cases where a spherical symmetric ionized bubble is carved out of the neutral expanding IGM (as in Dijkstra & Wyithe 2010). This approach though is less appropriate for cases where scattering occurs in only a few (neutral) clumps embedded in an otherwise fully transparent medium. However, the method works well in test cases similar to the configurations we are studying in our simulations, and thus we consider the main results of this study to be robust.

RESULTS
Using the pipeline outlined in the previous Section, we simulate ∼ 130 LAEs from L05-30 to evenly sample the DM mass range 10 9−11 M . Table 2 contains the properties of the haloes used in this work. The simulations were run with our default parameters (see Sec. 4.3). Dependence of the results on the choice of the parameters is discussed in Section 6. Before trying to understand the properties and statistical trends of 130 galaxies, we focus on understanding  how the IGM close to the objects affects the appearance of a single LAE. We investigate the ionization structure in the IGM and how it affects the SB maps. Having understood the effects on an individual galaxy, we study the behaviour of the whole sample and the trends shown in observability and Lyα escape fractions due to the IGM.

Behaviour of an Individual Galaxy
To understand how the density field in the IGM around the source affects the shape of the ionized region, which in turn will reflect on the propagation of the Lyα photons, we look at the ionization structure around the most massive dark matter halo in L05 (same galaxy as in Fig. 1) of DM mass of 2.6 × 10 10 M and stellar mass of 3.6 × 10 8 M . Figure 3 shows the mid-plane of the ionization structure at the end of the RT simulation. We can see that the ionized region is not spherical and its edges are shaped by the high density regions (e.g. filaments) in the IGM. Lyα photons propagating through the ionized gas will travel a large distance unscattered while being redshifted away from the Lyα rest wavelength, thus improving the chances of escape from the neutral IGM beyond the ionized region. Lyα photons encountering a high density filament instead will undergo several scatterings before being able to eventually escape Figure 6. Stellar mass of the objects plotted against their DM mass. The objects from different simulations are marked as -L05 (red crosses), L10 (green squares), L20 (blue circles) and L30 (pink triangles). The best fit line (dashed) is obtained using only haloes from L05 and L10, and it has a log-log slope of 1.12 with a standard deviation σ=0.13 dex. The σ region is shaded in light grey.
the simulation box. Because of the different paths followed by the photons, we expect different SB distributions and spectra, depending on the viewing direction. This can be clearly seen in Figure 4, where the SB maps obtained from the six sides of the simulation box are shown for the same object.
The differences in the maps reflect the structure in the IGM. In general, the lesser the photons are scattered before reaching the edge of the box, the more concentrate the SB profile is. In this case, the central pixels also have a much higher value of SB compared to cases in which the photons are scattered by high density neutral gas. Note that, as a result of the very inhomogeneous structure of the IGM, the maximum SB value in the image differs by more than an order of magnitude for different lines-of-sight. This is more evident in Figure 5, which shows a cross section of the SB maps for the three directions shown in the top row of Figure 4. Plotted are the SB values in the map for an x axis value of 0 against the y axis of the map. As noted earlier, the difference in the value of the central pixels in the maps is more than an order of magnitude for different viewing directions of the same object.
The detectability of the objects simulated in this work depends on the SB thresholds SB th set by different observational campaigns. Since the lowest SB value of a pixel in the maps is SBmin = 10 −23 erg s −1 cm −2 arcsec −2 , only for such a small (or lower) value of SB th would it be possible to detect all the flux from this object. At this SB th level the observed luminosities, calculated using Equations 11 and 12, for the source shown in the previous Figures are in the range L obs = 4 − 12 × 10 41 erg s −1 . For reference, the input luminosity is L esc Lyα = 9.2 × 10 41 erg s −1 . This would lead to an escape fraction fesc,Lyα,IGM = L obs /L esc Lyα in the range 0.43 − 1.3 (see below for a discussion). If instead e.g. SB th = 3 × 10 −21 erg s −1 cm −2 arcsec −2 , only two maps would result in a detection, with luminosities L obs (> SB th ) = 8×10 40 erg s −1 and 7.9×10 39 erg s −1 , com-    get an additional contribution of photons which have been scattered off from other lines-of-sight. A Lyα photon encountering the surface of a high density region (like a filament close to the source) has a high probability of scattering towards the lower density gas around it. Once in a line-of-sight through a void, the photon travels longer distances without scattering while redshifting out of resonance, improving the probability of escape towards the observer. Thus linesof-sight through voids are preferred over the ones through high density regions and get Lyα photon contributions from other lines-of-sight, leading to a more than average photon escape and an occasional effective escape fraction > 1. Note that while in the case of absorption photons removed from a line-of-sight do not contribute to any other line-of-sight leading to an escape fraction always ≤ 1, in the case of scattering photons removed from a line-of-sight appear in another one leading to escape fractions which can be also > 1. This makes a 3D treatment of the Lyα photon scattering of particular relevance for a proper LAE modelling.

Statistical Trends
Next we discuss how the distribution in SB affects the statistical properties of a large sample of LAEs. The difference in SB profiles in fact leads also to a spread in the Lyα luminosity of the same object observed from different directions, as seen earlier. The objects in our sample have an intrinsic correlation between the stellar and the DM mass, with a loglog slope of 1.12 and a standard deviation σ=0.13 dex (see Fig. 6). The best fit line (dashed) has been obtained using only haloes from L05 and L10. In fact, due to low particle resolution in the larger boxes, some of the haloes selected in L20 (blue circles) and L30 (pink triangles) have a stel-lar mass which lies below the expected trend. A correlation similar to the one shown in Figure 6 is also expected to exist between the input Lyα luminosity of the source, L esc Lyα , and the DM mass. Any additional scatter in the observed Lyα luminosity L obs is instead due to the effect of the IGM. Figure 7 shows L obs of the full sample calculated for each of the six sides of the box plotted against the DM mass for each of the simulated LAEs. As in Figure 6, there is a strong correlation with a best fit log-log slope of 1.12, which has been calculated using only haloes from L05 and L10. The scatter now though is 0.17 dex, larger than the intrinsic one, and it is induced by the structure in the IGM.
An alternative way to look at the scatter due to the IGM is to plot the escape fraction fesc,Lyα,IGM against the DM mass of the objects (Fig. 8). This removes the scatter present in the intrinsic luminosity L esc Lyα . As we can see, there is no strong correlation between fesc,Lyα,IGM and the DM mass because the scatter due to the IGM structure dominates. The average escape fraction has a nearly constant value of ∼ 0.73 with a σ scatter of 0.18. Some lines-of-sight have escape fractions greater than 1, as explained earlier.
But the statistical mean of the six observed IGM escape fractions for all the objects in the sample is 0.73. We would get an average escape fraction of 1 only if we observed all the directions for the same object and calculated the mean of all the lines-of-sight. The average of the six viewing directions is < 1, showing that scattering removes a large fraction of the flux for most lines-of-sight. Most of the flux from the object escapes along a small fraction of lines-of-sight through large voids, showing again the importance of IGM structure close to the source for the observability of LAEs. The spread in the correlation decreases for higher DM masses. This could be due to two reasons -a more uniform environ-  ment for haloes of higher masses leading to lesser scatter; or the resolution effects caused by the larger physical size of the grid cells (∝ r200) for higher mass haloes compared to low mass ones. But the general consistency of fesc,Lyα,IGM values in L20 and L30 with L05 seems to indicate that resolution issues might be less important. Higher resolution RT simulations of larger comoving volumes are needed to confirm this effect, which is beyond the scope of this study.
Next we investigate the detectability of this dataset by plotting in Figure 9 the observed luminosity of the objects calculated from the six sides of the simulation boxes for different SB thresholds, SB th . Only one object is detected for SB th ∼ 5 × 10 −20 erg s −1 cm −2 arcsec −2 (pink triangle). Note that this is not the most massive object neither in the full sample nor in the simulation box L05. Nevertheless, it looks as the brightest due to the density and velocity fields in the surrounding IGM, which lead to a very concentrated SB profile and ease the detectability of the source. From this example it is clear that the structure in the IGM plays an important role in the detectability of objects by shaping their SB profile and only a full 3D RT approach can properly capture its effect. As we lower the SB threshold, more objects are detected and a large fraction of the luminosity is accounted for. By SB th = 10 −20 erg s −1 cm −2 arcsec −2 (blue circles), seven objects are detected covering a mass range of two orders of magnitude. Even though there is a correlation of luminosity with DM mass apparent in the data, a large scatter is present at small masses. Decreasing the threshold by another dex, leads to the detection of almost all objects with a scatter in luminosity spanning three dex in each mass range. A SB level of SB th = 10 −25 erg s −1 cm −2 arcsec −2 is needed to observe the total flux from all the sources (the same as Fig. 7). The SB values calculated here depend on the choice of a number of different parameters, which will be further discussed in the next Section.
We caution the reader that current observational SB thresholds are significantly higher than the 'typical' SB levels of Lyα radiation scattered in the IGM as found in our simulations. This leads to an effective escape fraction due to the IGM (also see Laursen et al. 2011) which is very low (a few %). In addition, the observed SB is SB obs ∝ (1 + z) −4 (see Eq. 8), making deep detections increasingly difficult at higher redshifts. Current observational thresholds for narrow band surveys are 2.64 × 10 −18 erg s −1 cm −2 arcsec −2 at z = 5.7 in SXDS (Subaru/XMM-Newton Deep Survey; M. Ouchi; private communication; Zheng et al. 2010). Steidel et al. (2011) achieved deep SB limits of 2.5 × 10 −19 erg s −1 cm −2 arcsec −2 for Lyman Break Galaxies (LBGs) at z ∼ 3 which was obtained through stacking of 92 images. Their total observed Lyα flux translates to an escape fraction of Lyα of ∼ 0.25 − 0.5 (see Dijkstra & Hultman Kramer 2012). This is in agreement with Figure 9. Obtaining these SB thresholds at higher redshifts is a very difficult task. JWST, with its predicted 3 SB th ∼ 10 −19 erg s −1 cm −2 arcsec −2 , would be an important step for the deep detections of LAEs.

PARAMETER STUDY
In this Section we study the dependence of our results on the different parameters adopted. To this aim we use the biggest object in L05 with a dark matter mass of 2.5 × 10 10 M . The reference values are those described in Section 4.3.
A minor contribution to the observed flux comes from the photons emitted by the recombining gas in the IGM, which has been estimated to be ∼ 3%. This is not significant compared to the flux from the central source, but it should be kept in mind that this flux is not affected by the absorption in the ISM and thus has a higher escape probability. Since it is not concentrated as a point source, its contribution to the detectability of the object is negligible; rather, it might reduce the strength of the point source by spreading the flux to the scattered component of the angular directional distribution. We do not include the re-emission photons in our parameter study.

Escape Fraction
As already discussed, we need to specify the value of both ionizing, fesc,ion,ISM, and Lyα, fesc,Lyα,ISM, escape fractions from the ISM. The amount of ionizing photons reaching the IGM is controlled by fesc,ion,ISM, thus determining the extent of the ionized region through which the Lyα photons propagate. The value of the escape fraction also affects the production of Lyα photons by recombinations in the ISM. In our reference simulations, the box around the object has a size of ∼ 35 r200 comoving. Because this is not large enough to contain the HII region produced by fesc,ion,ISM = 0.99, for this parameter study we use an object with DM mass of 2.29 × 10 10 M from L10. We extract a cube with a side of ∼ 117 r200 comoving, which is gridded to a 256 3 grid.
To quantify the effect of changing the value of fesc,ion,ISM on the spectrum, in Figure 10 we show the Lyα spectrum of the photons exiting the box for simulations with fesc,ion,ISM = 0.02, 0.5, 0.99, while fesc,Lyα,ISM is set to the default value of 0.3. Note that the photon packets are always blueshifted back to the source position from the edge of the box. Also note that the spectrum is for all the photons exiting the box throughout the simulation time. The Lyα photon count increases with decreasing fesc,ion,ISM (see Eq. 5), as recombination in the ISM is the dominant contribution to Lyα photon production. It also affects the shape of the spectrum by moving its peak closer to 1215.67Å, as the photons escape with less scatter. This can also be seen from the SB maps in Figure 11. The higher the fesc,ion,ISM is, the more concentrated is the SB profile. But due to the lower production of Lyα photons in the ISM, the SB values in each pixel decrease for higher fesc,ion,ISM. Thus, increasing fesc,ion,ISM reduces the overall detectability of the source at a specific SB th , but improves its detectability at lower SB th with respect to a dimmer galaxy with lower fesc,ion,ISM giving similar total Lyα fluxes.
We have then tested the effect of the escape fraction of Lyα photons due to the ISM, fesc,Lyα,ISM, which controls the fraction of Lyα photons reaching the IGM, as detailed in Equation 6. We use the same object and box of the previous tests, with fesc,Lyα,ISM = 0.01, 0.3 and 0.99. Figure 12 shows the spectrum of the photons exiting the cube, whose intensity increases with increasing fesc,Lyα,ISM as more photons reach the IGM. On the other hand, the shape of the spectrum is not altered. Similarly, the morphology of the SB profile is not affected.
The escape fraction of Lyα photons from the ISM has been estimated in a large number of studies (e.g. Le Delliou et al. 2005Delliou et al. , 2006Davé, Finlator, & Oppenheimer 2006;Nagamine et al. 2010;Ouchi et al. 2008;Dayal et al. 2009;Kobayashi, Totani, & Nagashima 2007;Laursen, Sommer-Larsen, & Andersen 2009;Steidel et al. 2011;Schaerer et al. 2011;Forero-Romero et al. 2011;). Since the amount and type of dust in the ISM of a galaxy and its correlation with other halo properties is highly unknown, we choose a reference value of 0.3. This adopted value is similar to the value that was required by Dayal et al. (2009) in order to reproduce the observed Lyα luminosity functions at z = 5.7 (note that this requirement assumes that the IGM transmits ∼ 50% of all photons: for a more transparent/opaque IGM, the required escape fraction would be lower/higher).

Effect of Input Lyα Spectrum
The RT of Lyα photons through the IGM is affected by their wavelength as they escape the ISM. In our default runs, we use a monochromatic line at 1215.67Å, but we can expect the spectral shape of the photons emitted by the stellar sources to be distorted by their interaction with the ISM before entering the IGM. In fact, observations of LAEs at z ∼ 2 − 3 show a non-monochromatic spectrum (Verhamme et al. 2008). Other authors (e.g. Dayal et al. 2009;Zheng et al. 2010;Dijkstra & Wyithe 2010) use gaussian profiles with a width given by a fraction of the circular velocity at the virial radius of the DM halo, with wider profiles leading to an easier escape of the Lyα photons from the surrounding medium. A lot of effort has gone into modelling more accurately the Lyα line profile for a range of parameters (Schaerer et al. 2011), but the exact shape of the spectrum and its connection to the galaxy properties are still unclear. Because of all these uncertainties, we choose to investigate the effect of the input Lyα spectrum simply by shifting a monochromatic spectrum.  Figure 13 shows the spectrum of all the Lyα photons exiting the box for different values of the input Lyα monochromatic photons, i.e. with velocity shifts of 0 (reference case), 200, 400 and 800 km s −1 . Since the scattering cross section is ∝ ∆ν −2 , the larger the shift, the lesser the scatter suffered by the photons. The velocity field in the IGM also plays a role in determining the scattering probability of these photons. Since the velocity distribution of the IGM around the object has a complex structure, it is not possible to uniquely predict the velocity shift at which the photons stop being scattered. For our test object, as the photons are shifted by 200 km s −1 , the scattering is reduced by a factor of two, which in turn is reflected in a reduction of photons at 1215.67Å. For a shift of 400 km s −1 , the scattering is reduced by an order of magnitude. Monochromatic spectra shifted by ≥ 800 km s −1 get hardly scattered. This also affects the SB profiles, as shown in Figure 14, where the SB in the four cases is plotted for one of the sides of the box (top-left map in Figure 4). As the velocity shift increases, the probability of scatter of photons decreases, leading to a more concentrated SB profile. In addition, only a small fraction of the photons are scattered, forming a low SB halo. This can be seen more clearly in Figure 15 where the cross section of the maps is shown. The larger the velocity shift is, the higher is the SB value in the central pixel, improving the detectability of the source. Also, a larger fraction of the flux escapes from this pixel. For example, 0.004, 0.1, 1.5 and 4 per cent of the total flux escapes through this pixel for the cases of no velocity shift, a velocity shift of 200, 400 and 800 km s −1 , respectively. Therefore, the initial spectra of Lyα photons (i.e. the shape induced by the processing of the Lyα photons in the ISM) has a very important effect in determining the SB profile of the object and on the observability of LAEs. Quantifying this in more details is beyond the scope of this study and will be investigated elsewhere.
From Lyα spectrum models of Schaerer et al. (2011), the peak of the spectra of photons escaping from the galaxy is at wavelengths where the velocity shift is equivalent to ∼ 1.5 − 2 times the outflow velocities. The typical outflow velocities estimated from observations by Verhamme et al. (2008) is ∼ 150 − 200 km s −1 , thus showing a Lyα spectrum peak at 200 − 400 km s −1 . With a simplifying assumption that almost all the flux enters the IGM at these wavelength, we can estimate the expected changes in observed luminosity 4 at different SB th due to this velocity shift. At SB th = 10 −21 erg s −1 cm −2 arcsec −2 , for a velocity shift of 200 (400) km s −1 , f eff esc,Lyα,IGM = 0.01 (0.1), compared to 7 × 10 −4 in our reference case with no velocity shift. Thus detectability of a LAE at a specific SB th can improve by orders of magnitude for realistic velocity shifts of 200 − 400 km s −1 even with xion = 0 (see also Dijkstra & Wyithe 2010). Therefore, understanding the effect of the ISM of the galaxies along with their outflows is a crucial input for better modelling of LAEs. This topic will be addressed elsewhere in further detail. 4 Note that a specific SB threshold might select more than 1 pixel.

Effect of IGM Ionization
Another important factor determining the observability of LAEs is the ionization level of the IGM outside the fully ionized region created by the source. The initial ionization fraction of our simulations is set to a uniform default value of xion = 0 for convenience, but simulations of the reionization process indicate that most of the IGM is substantially ionized at the redshifts studied in this paper (e.g. Ciardi et al. 2011). The IGM ionization level is determined by a combination of a background radiation and local sources. Here we study how these two contributions change the observability of our test source.

Uniform Ionization Distribution
To understand the effect of the IGM ionization level on the Lyα RT, we repeat the reference simulation with an initial volume averaged ionization fraction of xion = 0.5 and 0.89. The first value corresponds to an epoch when only half of the IGM is ionized, while the second is the value obtained from the reionization model E1.2-α1.8 described in Ciardi et al. (2011), which is designed to be consistent with existing observational constraints such as the Thomson scattering optical depth and the photo-ionization rate at z < 6. Figure 16 shows a slice cut through the ionization structure at the end of the simulation for the three cases described above. As expected, the higher xion is, the larger is the fully ionized region produced by the source, although the differences are not dramatic because the high density filaments surrounding the source effectively confine the fully ionized region in all cases and determine its shape. Note that while the edges of the ionized region are controlled by the high density filaments in the IGM, the inclusion of ionizing flux from galaxies in those high density clumps could lead to a different ionization structure (as shown in Sec. 6.3.2).
As a larger ionized region results in less scatter for the Lyα photons, we expect an easier escape in this case. This is shown in Figure 17, where the spectrum of all the Lyα  photons exiting the simulation box is plotted. Due to the lower number of scatterings undergone by the Lyα photons for higher values of xion, the peak shifts to bluer wavelengths and the photon count at the peak increases while reducing the width of the line profile. The SB maps for one of the sides of the box for the different xion values is shown in Figure 18. As xion increases, the SB profile becomes more concentrated while improving detectability. From the maps, we can infer that observations at SB th = 3 × 10 −21 erg s −1 cm −2 arcsec −2 would appear as a point source to the observer. For this viewing direction, the luminosity (effective escape fraction) of the point source detected for xion = 0.5 and 0.89 is 8.4 × 10 38 erg s −1 (0.001) and 3.7 × 10 40 erg s −1 (0.04), respectively, compared to a non detection at xion = 0. More specifically, if we were to consider the maps from all six faces of the cube, at xion = 0.89 the object would be detected as a point source from all the six directions with f eff esc,Lyα,IGM in the range 0.004 − 0.6, compared to only two detections as point source in six maps at xion = 0 with f eff esc,Lyα,IGM = 0.008 and 0.09. Thus at high ionization fractions most of the objects can be detected as point sources at lower SB thresholds compared to our reference case. Our results appear similar to Zheng et al. (2010) who find that even at a mean neutral fraction of 10 −4 , only a small fraction of the photons (8-33 % of the intrinsic luminosity) escape as a point source.
Because this paper is primarily focused on the effect of the IGM in the immediate surroundings of a source of Lyα photons, we do not follow the propagation of the photons to the observer. To do so, we would need much larger boxes, losing the resolution necessary to properly account for the scales we are mainly interested in here. Nevertheless, the IGM outside our simulation box is bound to scatter the photons in an amount dependent on their frequency when they escape the box. The optical depth τ due to scattering away from the line-of-sight to the observer (at z = 0) for a LAE at z = zs is given in Equations 1 and 2.
For xHI = 1 at zs = 7.7, the optical depth τ ∼ 1 for ∆v ∼ 1400 km s −1 , corresponding to a comoving distance of 8.6 h −1 Mpc from the source which is much larger than our cube sizes. Beyond this distance from the source the Lyα photons have redshifted far enough from resonance such that the neutral IGM does not further affect the Lyα SB profiles. However, in some cases analyzed in this work -in particular low mass objects for which smaller boxes are extractedthe photons emitted by the source do not redshift far from the line resonance when they reach the edge of the box. The neutral gas outside the box can then further modify the predicted Lyα SB profile.
As a simple estimate of the scattering in neutral gas beyond the boundaries of the cube we can use Equation 2. The effect of this scattering on the spectra of photons escaping the box for our test object is shown in Figure 19. For xion = 0, the τ integrated over the whole spectrum is 2.9, but as the ionization increases, the τ decreases to 1.8 at xion = 0.5 and to 0.55 at xion = 0.89, leading to more photons being observed. The scattering beyond the cube can therefore be significant for small xion, resulting in lower SB values per pixel extending over a larger angular size. Even deeper SB cuts would then be required to observe the total luminosity of the object. In contrast, for high xion scattering beyond the cube can be safely ignored, and the SB profile would retain the point source like appearance seen in Figure 18. The qualitative results described in the previous sections (e.g. a spread and increment in f ef f esc,Lyα,IGM for lower SB th ) would still hold. Another factor which reduces scattering outside the simulation box is the velocity shift discussed in Section 6.2. A shift of 200 (400) km s −1 reduces the effective τ to 2.8 (2.3) even for xion = 0. This implies that the predicted steepening of the SB profiles with xion (see Sec. 7 below) and velocity off-set are (conservatively) underestimated by our calculations.
We point out that constraints on the reionization history that are provided by current observations of the Lyα forest and the cosmic microwave background, as well as theoretical investigations, suggest that the ionization fraction at z = 7.7 was mostly complete (i.e. well in excess of xion ∼ 0.50; see e.g. Fig. 11 of Pritchard et al. 2010or Fig. 4 of Ciardi et al. 2011. This suggest that for realistic constraints on xion at z = 7.7 the scatter beyond the cube is not significant, and so we do not include this extra scatter in the rest of our analysis.

Clustering of Nearby Sources
While a roughly uniform ionization fraction could be induced by an ionizing background, photons from neighboring haloes clustered around a source would induce further fluctuations in the gas ionization structure. To investigate the impact this has on our reference simulation, we propagate ionizing photons from all the stellar sources present in the cube (referred to as 'case C' for clustering), rather than only from the central source. Figure 20 shows a slice through the ionization structure obtained in case C (middle panel) along with that from the reference simulation (left panel). As expected, the ionization bubble around the object is much larger and complex when the effect of all the sources is con-sidered. In case C, the final value for the volume averaged ionization fraction is xion = 0.340, compared to 0.047 of the reference case.
To understand how the ionization structure due to clustered sources affects the propagation of Lyα photons compared to a uniform ionization, we do another test run (referred to as 'case U' for uniform) producing the same final volume averaged ionization fraction of case C, i.e. 0.340, but this is induced by photons emitted only from the central source embedded in a gas with a uniform initial ionization fraction of 0.293. Case U is shown in the right panel of Figure 20 and exhibits a very different ionization structure compared to case C. While the region with fully ionized gas is smaller in case U, the situation is reversed in the region beyond the bubble with some fully neutral gas still present in case C. Note that in all cases the Lyα photons are emitted only by the central source. Figure 21 shows the SB maps for the three cases discussed above (for the top right panel in Fig. 4). Both case C and case U show a more concentrated SB profile compared to the reference run as expected for a lower xion, but the details in the SB maps differ due to the distribution of neutral gas. In particular, all the maps (from the six faces of the cube) of case U have a more concentrated structure compared to case C due to the higher ionization level of the gas outside the fully ionized bubble. In fact, since the ionization bubble is not large enough for the Lyα photons to be redshifted out of resonance when they reach its border, the neutral gas distribution beyond it is particularly important for the determination of the corresponding SB maps. Because the concentration in the SB improves detectability, case U is expected to favor observations of LAEs. At SB th = 10 −21 erg s −1 cm −2 arcsec −2 , case C has f eff esc,Lyα,IGM ∼ 0.001, while case U has ∼ 0.004, compared to the reference case of ∼ 0.0007.
It might seem that clustering of neighbouring sources does not substantially change the observability of the LAE, but because the flux is redistributed in different directions due to the different ionization structure in the IGM, f eff esc,Lyα,IGM can vary by orders of magnitude depending on the observing direction. For the same SB th , another direction (top right panel in Fig. 4) has f eff esc,Lyα,IGM = 1.5, 0.7 and 0.6 for case C, case U and the reference case, respectively. For reference, the input luminosity is L esc Lyα = 9.2 × 10 41 erg s −1 . Thus it should be kept in mind that even though both neighbouring sources and uniform background flux improves detectability of LAEs by increasing the mean ionization fraction of IGM in the region, the detailed ionization structure around a source plays a crucial role in determining the SB profiles in different directions.

EFFECTS ON ESTIMATES OF THE REIONIZATION HISTORY
Keeping in mind the discussions in the previous Sections, we now turn to study how structure in the IGM affects the methods which use the SB profiles of LAEs to estimate the mean ionization fraction of the universe. One of the most popular methods is to calculate the luminosity function of LAEs and determine its redshift evolution (e.g. Dayal et al. 2009;Kashikawa et al. 2011). Us- Figure 22. Luminosity functions calculated for the 45 most massive objects in L05 for different SB cuts. The runs are for initial uniform volume averaged ionization fraction x ion =0 (thick red), x ion = 0.5 (medium thick green) and x ion = 0.89 (thin blue). The SB thresholds are SB th = 10 −25 erg s −1 cm −2 arcsec −2 (dotted lines), 10 −21 erg s −1 cm −2 arcsec −2 (triple dot dashed lines), 10 −20 erg s −1 cm −2 arcsec −2 (dot dashed lines) and 5 × 10 −20 erg s −1 cm −2 arcsec −2 (dashed lines).
ing analytic models to describe the reionization history, the above authors estimate the evolution of τ (employing equations similar to Eq. 2). The intrinsic LF of LAEs is then calculated from halo mass functions, assuming that the luminosity of an object is proportional to its mass and that it is attenuated by e −τ due IGM absorption/scattering. A comparison with the observed LF provides a constraint on the actual τ and thus on the gas neutral fraction in the IGM. But these methods rely on the assumption that observations of LAEs obtain the total flux from the source detected as a point source. Some recent observations though have detected faint extended SB profiles (Steidel et al. 2011), as produced in our simulations. As we have seen, the probability of detection of a LAE strongly depends on SB th , especially for low xion. Here we investigate how a different choice of SB th affects the LAE LFs for different values of xion, aiming at a qualitative assessment, while a more detailed and quantitative comparison to observations is part of a future study.
We model the 45 most massive objects in L05 using xion=0, 0.5 and 0.89 as described in Section 6.3.1. For each object, we calculate the six SB profiles and their respective observed luminosities L obs (> SB th ), which yields a total of 270 data points (LAEs) at each SB th = 5 × 10 −20 , 10 −20 , 10 −21 and 10 −25 erg s −1 cm −2 arcsec −2 . We consider the six sides of a simulation cube as part of six separate observational fields with comoving side of 5 h −1 Mpc each. Structure in the IGM leads to differences in the SB profiles for the six sides (as seen in Fig. 4), which in turn leads to a scatter in the observed LFs calculated from the different fields. Therefore, LFs are calculated separately for each of the observational fields and then averaged in each luminosity bin to obtain the mean luminosity functions, which are shown in Figure 22.
First we focus on xion = 0, to understand the effects of SB thresholds on detections (as seen in Figure 9), and determine the luminosity functions. At SB th = 5 × 10 −20 erg s −1 cm −2 arcsec −2 only one object is detected with an observed luminosity of L obs (> SB th ) = 3 × 10 38 erg s −1 , while SB th = 10 −20 erg s −1 cm −2 arcsec −2 leads to the detection of 5 LAEs with observed luminosities larger than in the previous case. The LF is thus shifted to higher values in both axis. Going one order of magnitude deeper in SB leads to a huge increase in the number of detections (∼ 130 LAEs). This gives rise to a shift in the LF of about two orders of magnitude in both luminosity and number density. At SB th = 10 −25 erg s −1 cm −2 arcsec −2 all the flux in all the LAEs is observed, returning the intrinsic LF of the simulated LAE sample. It should be kept in mind that the value of SB th needed to observe the full flux depends on the intrinsic luminosity of the sources and neutral IGM structure around it.
In the real universe, both the mean ionization fraction and the intrinsic properties of galaxies change. This leads to a complex evolution of the LF of LAEs with redshift. In addition, the different SB th from real observations add another level of complexity to the redshift evolution of observed LFs. To help in constraining xion with observations of LFs, here we try to understand how the LF depends on the ionization fraction of the IGM for different choices of SB th . For SB th = 10 −25 erg s −1 cm −2 arcsec −2 , the luminosity functions for different values of xion are equivalent, because all the flux from the sources is detected. For lower SB cuts instead, the number of detections changes with xion, leading to an increasing number density of LAEs at a given luminosity for higher xion. The difference between LFs at different xion is higher for smaller values of SB th . It should be noted that shallow detections (e.g. 10 −20 erg s −1 cm −2 arcsec −2 ) at higher xion (e.g. xion = 0.89) have similar LFs as deeper but incomplete detections (e.g. 10 −21 erg s −1 cm −2 arcsec −2 ) at lower xion (e.g. xion = 0). Thus, when comparing luminosity functions at different redshifts, the SB threshold of the observations should be taken into account to obtain a realistic estimate of the decrease in the number density of LAEs. Another point to note is that discrepancies between LFs for different values of SB th are smaller at higher ionization fractions (e.g. xion = 0.89). Thus only mild differences in LFs for orders of magnitude in SB th at a specific redshift can be used as an evidence for high xion. At high xion, absorption outside the edge of the cube (as discussed in Section 6.3.1), which has not been included in this test, would lead to an even larger gap between LFs for different SB th . Nevertheless one has to keep in mind that velocity shifts discussed in Section 6.2 (e.g. due to outflows from the galaxy) could also lead to milder changes in luminosity functions with changing SB th .
Due to the difficulties in obtaining deep observations of LAEs (refer to Section 5), we look into the possibility of gathering additional information from stacked profiles. This was recently done by Steidel et al. (2011) at z = 3 to estimate the low SB extended emission present in individual non-detections. In addition, also Zheng et al. (2011) looked at estimating total Lyα luminosities from stacked profiles for different mass and luminosity bins. In our study, we investigate the shape of the stacked SB profiles without binning for halo properties, to understand how the profiles are changed by xion. The data set used is the same one discussed above. In general, for xion = 0 the individual profiles show a flatter SB profile with substantial structure due to scattering through a mostly neutral but structured IGM. For higher initial ionization fractions, the strength of the point source increases as more photons escape with less scattering resulting in a steeper SB profile. These characteristics can be exploited to study the mean ionization fraction of the universe from stacked LAE SB profiles. Figure 23 shows the stacked SB maps from the 270 simulated maps of the 45 galaxies described previously. The top (bottom) panels refer to the mean (median) values in each pixel of the stack. The mean value is generally affected by the outliers when the distribution is not gaussian. In our stacks, the higher mass haloes act as outliers being brighter and more point-source-like due to larger ionized bubbles around them, but are fewer in number compared to the majority of SB maps of low mass haloes. This leads to the mean profile being steeper and brighter than the median profile at each xion. But it is important to note that in both mean and median stacks, the trend of the stacked SB profiles with xion is the same.
As we have seen in Section 6.3.1, the higher xion is, the steeper is the SB profile of a LAE. This trend is more important for haloes of lower masses, as the ionization bubble around these objects are smaller due to the lower ionizing photon flux, leading to a more diffuse SB profile at low xion. As xion increases, a larger fraction of the photons from these haloes escapes unscattered, leading to a more concentrated SB profile. For high mass haloes, even at low xion = 0, the photons escape with lesser scatter due to larger ionized regions around them, leading to a steeper SB profile. Due to the shape of the halo mass function though, there is a much larger number of low mass haloes in a comoving volume of the universe. Therefore, when stacking SB profiles of LAEs in a volume, the average of the stack is dominated by the low mass haloes. Thus the stacked SB profile at low xion = 0 is flat. The stacked profile has no structural details, because the profiles in individual SB maps have a random distribution, thus averaging out in the stack. At ionization fractions of xion = 0.5, 0.89, the photons escape easier in all directions leading to a steeper SB profile in all directions. This also appears when stacking SB profiles, with the stacked SB profile becoming more concentrated with increasing xion. Thus this method of stacking can be used to estimate the steepness of the profile and in turn the volume averaged ionization fraction of the universe at that redshift. As already discussed, the velocity shifts could also lead to the steepening of the individual SB profile even at xion = 0, thus reducing the difference in concentration between stacked profiles for different xion. But absorption due to IGM outside the cube at xion = 0 would lead to flattening of individual SB profiles. This would help in differentiating between xion in stacked SB profiles.
In short, a lot of information about the galaxies and their surrounding IGM can be obtained from the deep observations of LAE SB profiles. A more detailed modelling and comparison with observations is the work of a future study.

DISCUSSION AND CONCLUSIONS
The Epoch of Reionization (EoR) is an interesting and important event in the history of the universe, the time being when H in the gas changed from a mainly neutral to a highly ionized state. The details of the process are still unclear and are the topic of a large number of theoretical and observational studies. Lyman Alpha Emitters (LAEs), due to their strong Lyα emission line which scatters for even tiny amounts of neutral hydrogen, are used as one of the tools to probe the EoR. Modelling LAEs is a complex task with a large number of parameters and scales being important. Different studies either focus on a detailed understanding of individual aspects of LAEs or use numerical/semi-analytic approaches to model the LAEs and compare them to observations. In this paper, we focus on the Inter-Galactic Medium (IGM) close to the ionizing/Lyα source and its effect on the observability of LAEs due to Lyα RT. The IGM close to the source is important: the density, temperature and ionization structure controls the Lyα radiative transfer (RT), determining the intensity and surface brightness (SB) distribution along different lines-of-sight.
To achieve this, we simulated a sample of more than 100 LAEs, using galaxies and surrounding IGM extracted from simulation boxes of 5-30 h −1 Mpc at z = 7.7. Coupled RT of ionizing and Lyα photons was performed using CRASHα, to determine the ionization structure carved in the IGM by the ionizing photons exiting from the galaxies and the spectrum of the Lyα photons scattering through the remaining neutral hydrogen along different lines-of-sight. The outputs were also used to produce Lyα SB maps for the lines-ofsight perpendicular to the six sides of each simulation cube. Analyses were done to study individual objects as well as statistical trends. A parameter study was also undertaken to understand the dependence of the results on the different, currently uncertain factors involved in the simulations.
Our main results are: (i) Inhomogeneities in the IGM affects Lyα RT, leading to structure in the simulated SB maps of LAEs. The anisotropic escape of Lyα photons through the IGM causes large variations in the total flux in SB maps along different lines-of-sight for the same object. This leads to ∼ 30% more scatter in the observed luminosity-mass relation than the intrinsic one for our sample of simulated LAEs. This also leads to a spread in the values for the escape fraction of Lyα photons from the IGM: fesc,Lyα,IGM = 0.73 ± 0.18 (1σ). Note that some lines-of-sight, especially through voids, can have fesc,Lyα,IGM > 1 due to Lyα photon contributions from other lines-of-sight by scattering and higher probability of Lyα photon escape towards the observer through voids.
(ii) Observational campaigns have SB thresholds, SB th , which limit the fraction of the flux observed from the SB distribution. At SB th (e.g. 10 −20 erg s −1 cm −2 arcsec −2 ), the observability of the LAE strongly depends on the IGM ionization structure and velocity field around the objects. Therefore, in observational campaigns at very low SB th , a single detection within the observed comoving volume needs not be that of the most massive LAE in the region but could be a lower mass object seen through a biased line-of-sight. At even fainter SB th , more detections are possible, but the observed luminosity values can vary by orders of magnitude for a single object depending on the line-of-sight. To obtain the total flux from the object along each line-of-sight, one needs unrealistically deep observations at SB th ∼ 10 −25 erg s −1 cm −2 arcsec −2 . Thus we find that the impact of SB thresholds on estimates of observed luminosity is very important and needs to be taken into account in the calculation of luminosity functions of LAEs from observational campaigns (also see Zheng et al. 2010a).
(iii) One of the main factors which affect the Lyα RT through the IGM is the wavelength of the Lyα photons when they escape from the Inter-Stellar Medium (ISM) into the IGM. Outflows/winds in the ISM would redshift the photons reducing their probability of scattering in the IGM and aiding the escape of Lyα photons towards the observer (Dijkstra & Wyithe 2010). We find that the higher the redshifting of the Lyα photons before entering the IGM is, the more concentrated is the SB profile, making the detection easier at low SB th and improving the observability of LAEs. Therefore, outflows from the galaxies could weaken the imprint of reionization on the statistics of LAEs than thought previously. But due to the uncertainty in the link between the galaxy and outflows, this would add a challenge to estimating the ionization fraction of the IGM at high redshifts.
(iv) The ionization structure of the IGM also plays a very important role in Lyα RT thus determining the SB profiles. We estimate the effects of an ionizing background using RT simulations with an initial non-zero uniform ionization level. At low levels of mean ionization in the IGM (xion ∼ 0), the Lyα photons undergo scattered diffusion through the IGM and the SB profiles thus produced are more extended and faint, making it harder to observe the total flux from the object. At high levels of mean IGM ionization (xion > 0.5), the photons scatter less, leading to a more point source-like SB profile, making it easier to detect more flux for low SB th . Clustering of sources also affects the ionization structure around the object by making the region more ionized, but with a complex structure for the ionized bubble due to the distribution of neighbouring ionizing sources and IGM gas distribution. The qualitative effect on the SB profiles due to clustering is similar to that of a uniform ionization case with the same mean ionization level. But the detailed distribution of flux along different lines-of-sight varies due to the differences in the IGM ionization structure. Thus, this shows that proper treatment of the ionization structure in the region is important for better modelling LAEs.
(v) Due to the difficulty in achieving very high SB th in observational campaigns, stacking of SB maps can be used to extract more information from the current LAE samples mapped upto relatively lower SB th (e.g. Steidel et al. 2011;Zheng et al. 2011). We find that the mean/median stacked SB profile of LAEs becomes steeper at higher mean IGM ionization levels thus giving an additional technique to use LAEs to study EoR. But outflow and halo properties among others could also lead to steepening of SB profiles and should to be considered when using this technique.

ACKNOWLEDGMENTS
Many thanks to the anonymous referee for helpful comments, to Rudiger Pakmor for providing us the Gadget-toGrid code, to Ouchi Masami for sensitivity estimates of LAEs' surveys and to Andrea Ferrara for stimulating comments. AJD acknowledges support from and participation in the International Max Planck Research School on Astrophysics at the Ludwig-Maximilians University. AM acknowledges the support of the DFG Priority Program 1177.