The Effects of Lyman-Limit Systems on the Evolution and Observability of the Epoch of Reionization

We present the first large-scale, full radiative transfer simulations of the reionization of the intergalactic medium in the presence of Lyman-limit systems (LLSs). To illustrate the impact of LLS opacity, possibly missed by previous simulations, we add either a uniform or spatially-varying hydrogen bound-free opacity. This opacity, implemented as the mean free path (mfp) of the ionizing photons, extrapolates the observed, post-reionization redshift dependence into the epoch of reionization. In qualitative agreement with previous studies, we find that at late times the presence of LLSs slows down the ionization fronts, and alters the size distribution of H II regions. We quantitatively characterize the size distribution and morphological evolution of H II regions and examine the effects of the LLSs on the redshifted 21-cm signal from the patchy reionization. The presence of LLSs extends the ionization history by $\Delta z \sim 0.8$. The LLS absorbers significantly impede the late-time growth of the H II regions. The position dependent LLS distribution slows reionization further and additionally limits the late growth of the ionized regions. However, there is no"freeze out"of the H II regions and the largest regions grow to the size of the simulation volume. The 21-cm power spectra show that at large scales the power drops by a factor of 2 for 50% and 75% ionization stages (at $k = 0.1$ $\text{h} \, \text{Mpc}^{-1} $) reflecting the limiting effect of the LLSs on the growth of ionized patches. The statistical observables such as the RMS of the brightness temperature fluctuations and the peak amplitudes of the 21-cm power spectra at large-scales ($k = 0.05 - 0.1$ $\text{h} \, \text{Mpc}^{-1} $) are diminished by the presence of LLS.


INTRODUCTION
The first generations of baryonic structures, which formed before the Universe was a billion years old, released sufficient numbers of hydrogen ionizing photons into the intergalactic medium (IGM) to completely ionize it. This process of cosmic reionization is generally assumed to be driven by ionizing photons produced by stellar sources, most of them with energies at or slightly above the Lyman limit. As a consequence, reionization was a very patchy process with ionized (H II) regions expanding around the positions of photon email: hshukla@lbl.gov (affiliate) sources. The optical depth of the IGM initially defined the upper limit of the mean free path (mfp) of these ionizing photons to the distance at the edge of the ionized region within which their sources were located.
Over time, more sources formed creating newer ionized regions that grew and merged into an expanding patchy tapestry that eventually led to the completely ionized IGM at z ∼ 6. In this post-reionization era the mfp for ionizing photons, measured directly from quasar spectra, is found to be shorter than it would have been in a uniform IGM with mean baryon density exposed to the average ionizing UV background -for a review on quasar spectra see Bechtold (2003). The diminished value of mfp suggests the presence of discrete absorbing systems along the line-of-sight which have a high enough H I column density to absorb most of the photons near the Lyman limit, i.e., NH I > 10 17 cm −2 . These systems are known as Lyman Limit Systems (LLSs). The mfp due to these systems evolves with redshift and has been observationally determined over the redshift range 2.3 < z < 6 (Songaila & Cowie 2010;Worseck et al. 2014).
Conceivably, similar optically thick systems should have been present during reionization. This implies that the LLSs would determine the mfp of the ionizing photons when the sizes of the ionized regions grow larger than the mfp. This, in turn, would result in a slowing down of the growth of ionized regions compared to the case where such systems are absent. Proper modeling of the reionization process should therefore include the effects of such LLSs. The mfp evolution for z < 6 shows that the mfp at z ∼ 6 is around 9 proper Mpc (pMpc) and rapidly declines with increasing redshift suggesting that the LLS will have a substantial impact on reionization.
The reionization process, however, is complicated by the fact that LLSs are not the sole absorbers of ionizing photons. During the EoR a significant fraction of photons is absorbed in the ionized IGM due to recombinations. As the density in the IGM fluctuates, so does the recombination rate and a proper estimate of the total effect of recombinations requires resolving the small scale density variations, including corrections for self-shielded systems (in which the highest density parts remain neutral and do not contribute to the recombinations) as well as the effect of photo-heating, which modifies the distribution of small scale density variations.
The distinction between LLS and recombinations in the ionized IGM is in fact somewhat artificial, as the lowest column density LLSs will be high density but still ionized regions whose remaining H I column densities give them optical depths larger than 1. At higher column densities the systems will contain a fully neutral, high density core and can be described as self-shielding. Still, most of the absorption in such systems happens in the ionized layer sitting around the neutral core. In a recent paper, Kaurov & Gnedin (2013) present an analysis of the results of coupled radiationhydrodynamic cosmological simulations in the context of ionized IGM density fluctuations and self-shielded systems in order to establish whether a distinction between their contributions as photon sinks during reionization has any physical basis. Although, they see a continuous variation in properties, they do conclude that such a distinction is physically motivated.
The small scale structures in the IGM cannot yet be resolved in the large volumes required for studying the reionization process (> 100 Mpc) and therefore need to be included using sub-grid recipes. This is far from trivial due to the complexities described above and also because there are no observational measurements of the density and distribution of LLSs nor of the mfp for ionizing photons during the EoR. As a result a range of different approaches has been used. These approaches can be divided into two main categories, namely, (a) imposing a certain mfp for the ionizing photons, and (b) modifying the recombination rates to account for the unresolved density fluctuations. The former is equivalent to introducing discrete absorbers while the latter capture the effects of inhomogeneous recombinations throughout the IGM.
The first category is a necessary ingredient for semi-numerical simulations (Mesinger & Furlanetto 2007;Choudhury 2009;Santos et al. 2010;Mesinger et al. 2011;Alvarez & Abel 2012) and has also been used in some numerical simulations (Iliev et al. 2014). It is equivalent to imposing a hard limit on how far ionizing photons can travel but does not affect the photons for distances smaller than the mfp. Due to the uncertain evolution of the mfp for z > 6 the value adopted is typically a free parameter, see Greig & Mesinger (2015), which may or may not evolve with redshift. Earlier numerical studies focused on so called mini-halos as discrete absorbers of ionizing photons (Ciardi 2006;Mc-Quinn et al. 2007), after detailed radiative-hydrodynamics simulations of mini-halo photoevaporation (Shapiro et al. 2004;Iliev et al. 2005b) provided the evaporation times and photon consumption rates to assign to individual mini-halos overtaken by the global ionization fronts during the EoR. These were assumed to have a certain distribution and life time inside ionized regions and were then added as an additional optical depth for the ionizing photons. The mfp values can be calculated from the assumed distribution of the absorbers. A related semi-analytical approach was taken in which this minihalo photon consumption was accounted for in the velocity of the I-fronts that lead the expansion of intergalactic H II regions during reionization, including a statistical treatment of the biased clustering of minihalos surrounding the source halos (Iliev et al. 2005a;Shapiro et al. 2006). However, mini-halos will all have evaporated by the time reionization ends (Iliev et al. 2005c).
The second approach, i.e. accounting for unresolved recombinations through a sub-grid prescription, has been more commonly used in numerical simulations (Iliev et al. , 2007(Iliev et al. , 2008Kohler & Gnedin 2007) and was recently introduced in semi-numerical simulations by Sobacchi & Mesinger (2014). In this approach the mfp is not imposed but has to be calculated from the simulation results.
Although the details vary, all of the above studies find that the size distribution of ionized regions is affected by the presence of the LLS. Especially at the later stages of reionization, the ionized regions do not grow as fast as when no LLS are included. This typically leads to a reduction of the large scale power in power spectra.
Very few large scale numerical reionization simulations have included discrete absorbers, the exceptions being the ones mentioned above. This means that a number of questions remain unanswered. In this paper we want to address some of these questions using an improved algorithm to limit the mfp of ionizing photons. The questions we focus on are, (i) How a redshift dependent mfp, as suggested by the post-reionization observations, affects the reionization process? We use the homogeneous distribution of the LLS for examining the effects. (ii) Whether a density dependent distribution of absorbers changes their impact? We implement the density (position) dependent LLS to explore the same.
The units in this paper are defined such that all the lengths in the units of h −1 Mpc or h Mpc −1 are comoving unless mentioned otherwise. The cMpc and pMpc are comoving and proper Mpc respectively. This paper is organized in the following sections. In § 2 we discuss the Lyman-limit systems and their observations. The N-body and radiative transfer simulations along with the Lyman-limit systems implementation are presented in § 3. The size distribution of the H II regions and the associated power spectra derived from our simulations are discussed in § 4. In § 5 we present the effects of LLSs on the 21-cm observables. We compare the results of the position dependent implementation of the LLS in § 6 followed by our conclusions in § 7.

LYMAN-LIMIT SYSTEMS
The H I column densities, estimated by the absorption lines of the intervening hydrogen observed in post-reionization quasar spectra, are used to categorize the hydrogen absorption systems into three overlapping states, viz., Lyman-α forest (for H I column densities NHI < 10 17 cm −2 ), Lyman-limit systems (LLSs, 10 17 cm −2 < NHI < 10 20 cm −2 ), and damped Ly-α systems (DLAs, NHI > 10 20 cm −2 ). The Ly-α forest consists of low density and highly ionized structures, in contrast to DLAs, which are high density and partly neutral, thus exhibiting a strong damping wing of the Ly-α line. The studies of these systems have enabled precise measurements of the NHI values leading to high precision constraining of observables such as the primordial power spectrum.
The Lyman-α forest and DLAs do not affect the reionization process significantly due to their low optical depth (forest) and relative rarity (DLAs). In contrast, LLSs have both a relatively high optical depth and abundance, and thus the potential to considerably influence the later stages of the reionization (Alvarez & Abel 2012).
First observed as quasar absorption lines in surveys (Tytler 1982), the Lyman-limit systems appear as abrupt discontinuities in the quasar absorption line spectra at the rest-frame Lyman limit at wavelength λ ∼ 912 A. Prochaska et al. (2010) define LLS as regions with Lyman continuum optical depth of τLLS 2, i.e., NHI 10 17.5 cm −2 . The LLSs are assumed to be located in and around galactic halos. These systems are relatively easily identifiable even with low resolution and poor signal-to-noise. However, unlike the Lyα forest and DLAs, the LLSs are poorly understood largely because NHI estimations require complete spectral coverage of both the Lyα line and the Lyman break. At lower redshifts (z < 2.6) the Lyman limit is shifted into the UV spectrum and thus is unobservable from the ground. High redshift surveys (Songaila & Cowie 2010;Prochaska et al. 2010) present measurements of the number density function of the LLS. Another recent survey (Ribaudo et al. 2011) with Hubble Space Telescope archival data identifies 206 LLSs for z < 2.6.
All these surveys identify the LLSs in the absorption line spectra and estimate the number of LLSs per unit redshift per H I column (NHI) density function -f (NHI, z) ∝ ∂ 2 N /∂z∂NHI. The best fit power-law index, β, to this function constraints the column densities of the LLSs. In addition, earlier simulation studies (Kohler & Gnedin 2007;McQuinn et al. 2011;Sobacchi & Mesinger 2014) of LLS abundances and the mfp of ionizing photons agree reasonably well with the observations. Using the parametric values of the density function from the observations ( The general effect of the LLSs is to restrict the mfp of ionizing photons and consequently impede the evolution and merging of H II regions during reionization. Early on during reionization, when the ionized regions are still small, the mfp is determined by the size of these regions. Towards the end of the reionization, the H II regions reach sizes larger than the measured mfp due to LLS at z ∼ 6, of approximately 60 comoving Mpc and thus the mfp should be set by the LLS. This is the stage where the LLSs begin to regulate the ionization history. At late times, there are many groups of local ionizing sources that add to the ionization fronts further complicating the morphological evolution of the ionized regions. The simulations discussed herein attempt to quantify the effects of the LLSs on the reionization history, with the caveat that the lack of high redshift observational data leads to an ad hoc implementation of the LLSs at early times. The extrapolation of mfp from observed low redshifts to early times as far as z = 20 clearly yields unrealistic values, see Figure 1. However, as discussed below, LLSs only start to affect the simulations much later after z = 15.96 (LLS1) and z = 13.30 (LLS2) with mfp value of 0.1 Mpc in proper units. More accurate effects of LLS at higher redshifts could only be modeled based upon the actual distribution of the high redshift LLSs. However, the very first numerical simulations of a large volume and higher dynamic range presented here, elucidates various useful insights beneficial for future research in the field.

SIMULATIONS
The evolving matter density fields in a given comoving volume for the desired redshift ranges are generated with the Nbody code entitled CubeP 3 M (Harnois-Déraps et al. 2013). The evolution of the density fields is based upon the initial conditions specified by the standard Zel'dovich approximation and primordial power spectrum transfer function derived by CAMB 1 (Lewis et al. 2000), originally based on CMBFAST (Seljak & Zaldarriaga 1996). The cosmological parameters used are for the flat ΛCDM model of the Universe based on WMAP 5-year data combined with constraints from baryonic acoustic oscillations and high-redshift supernovae, given as, (ΩM = 0.27, ΩΛ = 0.73, h = 0.7, Ω b = 0.044, σ8 = 0.8, ns = 0.96). To ensure against numerical artifacts (Crocce et al. 2006) the initial conditions are generated at sufficiently high redshift (here zi = 300).
The CubeP 3 M code is public domain N-body code designed for simulating large-scale cosmological systems. The code is accurate, efficient, scalable, and parallel across distributed (MPI) and shared (OpenMP) memory systems. The underlying N-body algorithm estimates short-range  Table 3. The filled black squares are the data points from Songaila & Cowie (2010). The blue symbols are from the ionization models from Emberson et al. (2013) for three values of the ultraviolet background (Γ −12 ).
(sub-grid distances) gravitational forces using the particleparticle (P-P) method. While for the long-range forces, a 2-level particle-mesh (PM) method is applied. Computationally, in comparison to the P-P method which is of the order O(N 2 ), the P 3 M method has significantly lower overhead and is of the order O(N log N ), where N is the number of particles.
For the current study, two sets of comoving volumes of sizes, 37 h −1 Mpc and 114 h −1 Mpc, are used. The same N-body simulations were presented in detail in (Iliev et al. 2012). The smaller volume, 37 h −1 Mpc, is used for testing and validation purposes only. The results from the smaller volume are not included in this paper.
The simulations for 114 h −1 Mpc use 3072 3 ≈ 28.9 billion dark matter particles distributed in 6144 3 mesh cells. Each particle has mass of 5×10 6 M . Through the simulation steps as the structures start to form, the halos are identified using a spherical overdensity halo finder with overdensity parameter of ∆ = 178 with respect to the mean density. The halos with more than 20 particles (M > 10 8 M ) are considered resolved. The number density of the halos increases for lower redshifts and the mass function approaches the Sheth-Tormen mass function.  and Watson et al. (2013) provide detailed fits to the high-redshift halo mass function.

C 2 -Ray -Radiative Transfer
The second stage of the simulation performs radiative transfer using the C 2 -Ray (Conservative, Causal Ray-Tracing) code ). The conservative part of the code ensures spatial and temporal photon conservation, while the causal ray-tracing is implemented using the shortcharacteristic method.
The C 2 -Ray implements a discrete spatial and temporal version of the ionization rate equation as (Osterbrock 1989), where, Γ(r) is the ionization rate at distance r from the hydrogen ionizing source, Lν is the spectral energy distribution of the ionizing source at frequency ν, σν is the crosssection for the ionizing photons, and τν is the frequency dependent optical depth of the hydrogen gas for bound-free transitions.
Due to the spatial discretization of the computational domain, we do not know τν as a continuous function of position. As was shown in Abel et al. (1999), the photoionization rate of one cell whose center has a distance r from the source can be calculated as, where, V shell is the volume of the spherical shell with radius r and width ∆r, the shell is filled with neutral hydrogen of number density nHI,Ṅ (r − ∆r/2) is the rate of ionizing photons arriving andṄ (r + ∆r/2) the number of photons leaving this shell. Since,Ṅ Equation 2 implies that the local photo-ionization rate Γ depends on the difference between τin ≡ τ (r − ∆r 2 ) and τout ≡ τ (r + ∆r 2 ). In Section 3.3 we explain how these optical depths are modified to include the effect of LLS. In our reionization model, the ionizing luminosity of the collapsed halos is proportional to their mass M . Each halo produces a number of photons, for every n-body output of ∆t = 11.46 Myr. The efficiency factor fγ is the product fescf N , where, fesc is the ionizing photon escape fraction, f is the star formation efficiency, and N is the number of ionizing photon per stellar atoms, and mp is the proton mass. The parameter N depends on the initial mass function (IMF) of the stellar population producing the ionizing radiation. Its value for a Pop II population (Salpeter IMF) is ∼4000 and for a Pop III population (top-heavy IMF) it can reach ∼100,000. Due to the uncertainties in fesc and f , the value of fγ is not well constrained. In this study we use 10 for high mass halos (HMACH: High Mass Atomicaly-Cooling Halos) and 150 for low mass halos (LMACH: Low Mass Atomicaly-Cooling Halos), see Table 2. The higher value for the LMACHs is motivated either by a larger contribution of metal-free/poor stars or by a larger escape fraction. The LMACHs are also assumed to be susceptible to negative radiative feedback. When the cell in which an LMACH is present is ionized at the start of a new 11.46 Myr time step, the LMACH will not produce any ionizing photons.  Table 2. Simulation parameters for the 114 h −1 Mpc box with LLS. fγ is the star formation efficiency for high and low mass, and RT grid is the coarser grid for radiative transfer ray-tracing. The underlying cosmology uses the WMAP 5-year results.
From previous work we know that the efficiency factors chosen result in reasonable reionization history in accord with the WMAP optical depth value of τes = 0.089 ± 0.014 (Hinshaw et al. 2013). We note that the recent results from the Planck mission favour a lower value of 0.066 ± 0.016 (Planck Collaboration et al. 2015) for this parameter.

Simulating the Effects of Lyman Limit Systems
It is thought that the LLSs correspond to the denser, ionized gas associated with collapsed objects and dense filaments of the cosmic web which, although ionized, have sufficiently high column density of H I to result in an optical depth > 1. This optical depth is usually described in terms of a mean free path for the ionizing photons and the available observational constraints are expressed in the same way. There are two possible approaches to implementing the effects of the LLSs in our calculations -either as an enhanced local recombination rate, or as an additional absorber optical depth along each light ray. The first approach is closer to the physical mechanism responsible for these objects, and recombinations in dense structures below the resolution limit of the code are described with clumping factors. This has for example been the approach of Sobacchi & Mesinger (2014). Introducing a local clumping factor to reflect an enhanced recombination rate does not however directly relate to a value of the mean free path (mfp) due to such absorbers. It is thus difficult to impose such a constraint on the LLS models used in the simulations. Therefore, for this study we decided to work directly with the mfp instead. This approach is similar to prior (semi-numerical) work by Alvarez & Abel (2012), however, these authors implemented the mfp as a hard boundary which photons could not cross, which is not very physical. In our implementation the mfp defines an additional ionizing photon-absorbing component which after one mfp reaches an optical depth of 1. This approach allows us to directly use the published observational mfp expressions. An additional advantage of working with the mfp directly is that it is easier to connect the effect on the sizes of H II regions to the imposed mfp.
The mfp in the ionized medium is not a measured quantity beyond z = 6, however. Therefore, we extrapolate it from the lower redshift results. To the extent that the ionized regions during reionization can be viewed as being locally post-reionization, there is some justification in performing this extrapolation. The additional absorption that the implemented LLS component adds in the fully neutral regions is marginal, and therefore only affects the already ionized regions. To perform the extrapolation of the mfp beyond z = 6, we use the parametrization given by Songaila & Cowie (2010). Based upon the observational data, the number density of the LLSs per unit redshift path dz is parametrized as, where, f (NHI, z = 3.5) is the number density at z = 3.5. Estimating the log-likelihood function for the entire redshift range of 0 < z < 6 in Songaila & Cowie (2010) yields the values for the parameters to be f (NHI, z = 3.5) = 2.8±0.33 and γ = 2.04 +0.29 −0.37 . Furthermore, in the approximation that the column density function f (NHI, z)dNHI ∝ N −β HI dNHI (Petitjean et al. 1993), the mean free path is related to the number density as (Miralda-Escudé 2003), where, f (NHI, z) is measured above the column density NHI corresponding to a value a for the optical depth (here, a ≈ 1), τ is the optical depth at the Lyman limit, and ν0 is ionization edge frequency. Here, Γ is the gamma function and not the ionization rate (Songaila & Cowie 2010).
For simulation L1 (without LLS) at z = 8.39 (overlap), the average Lyman limit optical depth is 0.83 over the volume size of 114 h −1 Mpc. This translates into an mfp of 20.6 pMpc. This value is an order of magnitude larger than expected from an extrapolation of the observed lower redshift values. This is due the fact that our radiative transfer resolution of 0.45 h −1 Mpc does not resolve the sub-grid density inhomogeneities which determine the actual mfp value.
In order to test the impact of different evolutions of the mfp we use two different extrapolations, one using the parameters given by Songaila & Cowie (2010) and one using parameters derived from fitting the curve in the inset of Figure 1 in McQuinn et al. (2011). The simulations for these two choices are labeled as LLS1 and LLS2, respectively and the parameters used are listed in Table 3. The evolution of the mfp for these two sets of parameters is shown in Figure 1. For comparison, we also show the evolution of the mfps for three different photoionization rates in terms of Γ−12measured ∼ 0.3 s −1 for Lyα forest at z ∼ 6 -from the simulations of Emberson et al. (2013). Note that these mfps are based on the Lyman limit optical depth from the radiative transfer calculations of ionization by a fixed radiation background in a 0.5 Mpc sized box.
From Figure 1 it is obvious that both extrapolations result in very small mfps at very high redshifts. Below our cell size it does not make sense to implement the LLS model as we already use an assumed escape fraction for absorptions within the source cell. Furthermore, at very high redshifts the cosmic structures are much less developed, thus fewer LLS systems should exist, implying longer mean free paths due to these systems. It is likely, therefore, that these extrapolations become unreliable at such high redshifts. As we want to concentrate on the later stages of reionization when the ionized regions have reached sizes of 10-20 cMpc and LLSs should have more impact, we chose to only switch on our LLS absorption if the mfp is larger than 5 grid cells (3 cMpc). This occurs around z ∼ 15. We found that this choice did not impact the evolution around this transition redshift.

LLS implementation in C 2 -Ray
To include the effects of LLS due to mfp in the radiative transfer calculation, an additional optical depth term, τLLS, is added. The implementation of the additional optical depth is based upon the notion that after one mean free path, τLLS acquires the typical value of the optical depth of a Lymanlimit system at the Lyman limit. The opacity may be added in two different ways depending upon the assumptions made for the spatial distribution of LLSs. In the first scenario, a uniform distribution of the LLSs is assumed, while in the second, the LLSs are assumed to be concentrated in the higher density regions. The latter case is motivated by the observation that LLSs are likely to be located at the outer regions of the halos, embedded in the filaments of the large scale structure (Furlanetto et al. 2006;Erkal 2014).
In the uniform distribution of LLSs, we assume that each cell contributes equally to τLLS with a value ∆τLLS. In C 2 -Ray there are two optical depth values associated with each cell. The value, τin, is the optical depth between the source and the entry point of the ray into the cell, and ∆τ , is the optical depth of the ray that is contained within the cell. The difference between these two values is used to calculate the ionization rate Γ in the cell, see Equation 2. The sum of the two values is the input optical depth of the adjacent cell along the line of traversal of the ray.
For two cells, n and n + 1, with the same ray crossing their centers, the input optical depth for cell n + 1, τ n+1 in , is equal to the output optical depth of cell n, τ n out ≡ τ n in + ∆τ n . To include the LLS optical depth this equality is changed to This effectively adds the additional optical depth due to LLS in between the two cells which means that the LLS optical depth is not used in the calculation of Γ for a cell. In principle, one could add ∆τLLS to ∆τ but this would then impact Γ which would give incorrect results. The fundamental purpose to implement the sub-grid LLS model is to remove (absorb) photons due to LLS without impacting any other parts of the radiative transfer calculation. For the position dependent implementation of the LLS in C 2 -Ray, the primary difference is in the estimation of the ∆τLLS. The value of ∆τLLS is proportional to the normalized sum of the cross sections of all halos in a cell. The normalization is such that the average of all the cells has the average total geometric cross-section required to reproduce the assumed mean LLS opacity per cell. The cross section is defined as πr 2 vir , where rvir is the virial radius for a halo. The opacity in a cell is proportional to the cross section and the mean LLS opacity (over all cells), which is set by the fit to the redshift evolution of the mfp. The assumption is that larger halos with larger virial radii have more in terms of LLS structures associated with them. All halos are used for this estimation. Even though ionization may suppress star formation in the low mass halos, presumably the density fluctuations will survive better.
Given the assumed redshift evolution of the mfp of LLS (ν0, z), see Equation 7, we estimate the average additional optical depth per cell as, where, ∆x is the cell size. For the homogeneous case, this is the additional optical depth per cell. For the position dependent case, this is the average value per cell, with the actual value for a cell being proportional to the total cross section of halos in that cell. the cells with a larger total cross section for halos getting higher values and the cells without any halos getting none. Equation 9 ensures that the after a distance (ν0, z) all rays will have picked up an additional optical τLLS. We take τLLS = 2.
For the cells whose centers do not lie on the same ray, the short characteristic ray tracing algorithm constructs τin through interpolation of the relevant τout. The above modification works equally well if one replaces each τout by τ n out + ∆τLLS. As explained in Section 3.2, this additional optical depth is only applied when (ν0, z) > 3 cMpc. For smaller values we set ∆τLLS to zero.
In the following sections we use the fiducial homogeneous distribution implementation of the LLSs for various studies. The comparison of the homogeneous and position dependent methods is discussed in § 6.

RESULTS
In this section we summarize the results of the various analysis methods we use to characterize the simulations and compare them between the L1 (no LLS), LLS1 and LLS2 cases. Figure 2 shows the evolution of the globally averaged massweighted xm and volume-weighted xv ionized fractions of the three models as a function of redshift. In both the top and the bottom panels of the figure we note that the models start diverging at z = 14 and that the differences become more pronounced after z ∼ 10 − 11. The shorter mean free paths imposed by the LLSs delay the overall process of the expansion and merging of the H II regions. The simulation L1 (without LLS) reaches a global average volumeweighted ionized fraction of xv = 0.98 at z = 8.34, while for LLS1 and LLS2, the same level of ionization is reached at z = 7.61, and z = 7.71, respectively. The top panel shows the ratio xm / xv while the bottom panel shows only the averaged xm. From the top panel it is evident that the massweighted ionized fraction is always significantly higher than the volume-weighted fraction for all the three models. This means that reionization is inside-out, that is, the dense regions surrounding the sources are preferentially ionized first yielding higher xm averages. The ionization fronts expand outwards, eventually reaching the less dense regions and voids. The ratio xm/xv is the mean over-density of the ionized regions ) as shown below,

Globally Averaged Quantities
where,ρ is the mean density of the Universe. Figure 2 indicates that the over-density of the ionized regions is always larger than one and for the LLS cases is even larger than for L1. This is because less dense regions such as voids do not have ionizing sources and therefore require photons from afar to get ionized. The simulations LLS1 and LLS2 reach a numerical value of 1.0 at redshifts z = 7.76 and z = 7.617 respectively, while for the L1 model this value is reached much earlier at redshift z = 8.397, thereby, delaying reionization by ∆z = 0.78 − 0.64.

Photon Statistics
The reionization period is defined by the complex interaction and evolution of the ionizing photon sources and sinks. In the simulations the sources of the ionizing photons are the halos. The very first HMACHs form around z ∼ 21. For a more indepth study of the mass and the number density distribution of the halos, see Iliev et al. (2012). The clustering of halos defines the photon emanating regions and sets the initial conditions for the formation and evolution of H II regions. Initially, the mean free path is dominated by the size of the ionized regions in the IGM. While the LLSs continuously absorb ionizing photons, it is only at later times (z ∼ 14−10) that the absorption contributions due to the LLSs become significant.
The photon and baryon populations in the simulations are recorded to extract statistical properties of interest. The integrated Thompson electron-scattering optical depth τes for the three cases as a function of redshift as compared to WMAP-9 mean optical depth estimates lie within the 1-σ spread. They also fall within 1-σ range of the Planck estimate (0.066±0.016). The values of the optical depth, τes, in the simulations at z = 25.33 for L1, LLS1, and LLS2 are 0.0819, 0.0796, and 0.0788, respectively. The optical depth in the presence of LLS is thus diminished by about ∼ 0.002. This is expected as there are overall less ionized electrons available for scattering the CMB photons.
The three reionization histories are the result of available ionizing photons from the sources. Figure 3 shows the cumulative number of photons per baryons as a function of redshift. At the end of the ionization, the photons per atoms are twice as many for the LLS case. The LLSs therefore absorb approximately one extra photon before reionization is completed with not much difference between the two LLS cases. Their effect dominates over the recombinations included in the simulation which only consume 0.5 photon per baryon by the end of reionization.

Morphology of HII Regions
The morphology of ionized regions is complex. We use several methods to qualitatively and quantitatively study the  H II region sizes, distribution, and evolution in our set of simulations. These methods, as discussed below, provide complementary information.
One of the expected outcomes of different mean free paths (z) for the ionizing photons, is a changed evolution of the size distribution of the H II regions and the consequent ionization fraction history. Once the H II regions grow larger than the mfp in certain directions, not all sources inside the region can contribute to their growth and therefore they will not grow as fast as in the case without LLS. The H II regions can still grow larger than the mfp because they are driven by many sources some of which lie closer than the mfp to the edge of the region. Figure 4 shows examples of the morphologies and growth of the ionized patches in the set of our simulations. The top panels show the ionization history of the fiducial model L1 (case without LLS) spanning redshifts, z = 16.9-8.5, approximately 370 million years. The panels from left to right show the slow ionization process, which reaches a mass weighted ionized fraction of 1% only at z = 16.9, even though the first halos in the simulation (with M > 10 8 M ) form at z = 21.

Evolution of HII Regions
The bottom three panels of Figure 4 emphasize the morphological and topological difference between the three simulations. The three spatial ionization slices compare the global ionization of xm = 50% for models, L1, LLS1, and LLS2. The most obvious features are the different sizes of the larger H II regions, especially between the non-LLS and the LLS models. This is indicative of the presence of LLSs slowing down the merger process. Between models LLS1 and LLS2 the differences in shapes and sizes of the ionized regions are not severe. However, in the detailed statistical analyses discussed below, some differences emerge. As expected, the slow growth of ionized regions delays the completion of reionization for the LLS simulations.

Size Distribution of HII Regions
In this section we quantify the results seen in Figure 4 by using three different methods to study the size evolution of the H II regions in our numerical simulations. The statistical property measured in this analysis is the probability distribution function of the volumes (radii) of the H II regions. The three approaches employed to estimate these size distributions are the friends-of-friends (FoF, see , the spherical average (SPA, see Zahn et al. 2007) and 3D power spectra methods. All of these algorithms differ in their approach of defining the size of the H II regions as discussed below. However, the different techniques complement each other and together provide a greater insight into the morphologies of the H II regions (for further details see Friedrich et al. 2011).

Friends-of-Friends
The friends-of-friends (FoF) algorithm operates on the ionized fractions and generates a catalog of connected ionized H II regions. For a chosen ionization threshold, x th , the algorithm connects all the ionized neighboring cells and classifies them in a friendship based topology using the 'equivalence class' or 'sameness' method of the Numerical Recipes (Press et al. 1992). The H II regions catalogs based upon the ionization threshold and volume size are thus generated. These catalogs provide detailed insight in the evolution of the number densities and size distributions of the topologically-connected ionized regions. This method was first introduced in . Figure 5 shows the probability distribution function V dp/dV at different friends-of-friends threshold values, x th = 0.1, 0.5, and 0.9, versus the volume of the H II regions for the L1 model at a global ionization fraction of xm = 0.5 (z = 9.457). The figure highlights the effects of the threshold. The lower threshold values reduce the num-ber of smaller H II regions and increase the number of the larger ones. The results of the FoF method therefore depend on the choice of the threshold (as already shown in Friedrich et al. 2011). However, for a given value of the threshold the results of different simulations can still be compared in a meaningful way.
The topological evolution of the H II regions for the L1 and LLS2 models is shown in Figure 6 for a threshold value of x th = 0.5. The difference between the two LLS models was not discernible in the FoF analysis, therefore, only one LLS model is shown in Figure 6. The colors represents the probability distribution V dp/dV . It is evident from Figure 6 that the H II regions grow as the ionization fraction increases up to a point where the distribution of volumes separates into two populations comprising of very large and relatively smaller sized regions. The emergence of the dichotomy is primarily due to the merging of smaller regions into larger volumes as the ionization fronts travel outwards from the higher density areas. As expected, the larger H II regions of the order of ∼ 10 6 Mpc 3 appear slightly earlier in the L1 (z = 13.48) case as compared to the LLS2 (z = 12.31) model. Another difference between the two models is that the largest of the "small" H II regions disappear faster in the L1 case than in the LLS models. For the no LLS case L1, at z = 8.34 ( xm = 0.99), the largest of the "small" H II regions are 50% smaller than the equivalent population in the LLS2 simulation. This is indicative of fewer mergers in the LLS cases as these regions only disappear when they merge with the larger regions. Towards the end of the ionization, the main contributors to the global average of the ionization fraction are the largest regions. This emphasizes the trends we have noticed earlier where the ionizing photons of the shorter mean free paths are absorbed and fail to contribute in the formation of H II regions that grow and merge. From the observational perspective, the large H II regions could be tomographically imaged with SKA-class interferometers. Figure 7. Probability distribution function R dp/dR per radial bins of spherical HII regions as measured by the spherical averaging algorithm for the no LLS (L1, blue), and LLS models (LLS1, red; LLS2, black). The three sets are for a mass-weighted global ionization fraction xm = 50% (solid), 70% (dashed), and 90% (dotted). The threshold used is x th = 0.9. The vertical lines correspond to the mean free path for the two LLS models as listed in Figure 1 for 50% and 90% ionization fractions.
The volume distribution of such regions may help put limits on the mean free path and therefore on the LLS models.

Spherical Averaging
The second algorithm to evaluate the size distribution statistics was developed in Zahn et al. (2007). This spherical averaging technique constructs spheres of varying radii around each cell in the ionization simulation and estimates the enclosed ionization fraction. The largest spheres with ionization fraction greater than the defined threshold, x th , define the spherically averaged H II regions. In contrast to the FoF method, the SPA technique yields a smoother and spherical distribution function, biased towards the shorter axes of triaxial structures.
The SPA analysis highlights a similar behavior in the evolution of the ionized regions as seen earlier with the FoF method. The shorter mean free path for the LLS simulation affects the growth as measured in the radii of the spherical regions. In Figure 7, the solid, dashed, and dotted lines correspond to the 50%, 70%, and 90% global ionization rates respectively. The color motif remains the same throughout the paper with L1 (blue), LLS1 (red), and LLS2 (black). The redshifts for which the distributions are calculated z = 9. 457,8.958,8.636 for L1,z = 9.236,8.636,8.172 for LLS1,and z = 9.164,8.515,8.012 for LLS2. The vertical lines are the mean free path for the two LLS models as shown in Figure 1. These lines are plotted for the redshifts of 50 and 90% ionization.
It is evident in Figure 7 that the number of regions with smaller radii (< 0.2 − 0.3 Mpc) is relatively similar for all models at different stages of the ionization history. However, for the larger radii (< 0.3 Mpc) differences in the proba- Figure 8. The log-log plot of the dimensionless 3D power spectra of the ionized fraction, at xm = 70% and xm = 90% for the three models.
bility distributions between L1 and both the LLS models emerge; especially as the ionization progresses. For example, at the 50% and 70% ionization stages, the radii for the peak probability in the L1 simulation are larger by a factor of ∼ 1.2 compared to the LLS cases. To a lesser extent, the similar trend is seen when comparing the longer mfp LLS1 model to the LLS2. This is indicative of impeded ionization for the LLS models. At 90% ionization most of the H II regions have merged and therefore the distributions of spherical volumes are similar. In addition, at radii values reaching ∼ 80 − 90 Mpc the H II regions are as big as the simulation box (R114 = 81.42 Mpc). The spherical averaging algorithm reaches the limits of its applicability at this stage.
The maximum radii difference between L1 and LLS2 for 50, 70, and 90% ionization varies from roughly 18, 33, to 18%. Another characteristic that is apparent from the SPA analysis is that for early times (10% ionization) smaller (radius < 0.3 Mpc) spherically averaged spheres contribute most to the probability distribution, R dp/dR. As the ionization progresses, the contribution from smaller spheres decreases and from the larger spheres increases indicating growing H II regions and mergers. This trend is same for all the three models with impeded growth of radii for LLS models as evident in Figure 7. The bubbles of larger radii merge earlier for the L1 model. When ionization reaches 50%, the maximum radii of the H II regions are comparable to the mfp, more so for the LLS2 than for the LLS1 model. However, as the ionization progresses, the radii of the H II regions grow beyond the mfp due to mergers.
These results are consistent with ones from FoF method. The unimpeded ionizing photons in the L1 model noticeably differentiate the ionization and bubble growth history from that of the LLS models. There is not much difference between the LLS ionizing histories themselves, with radii less than 10% different at different stages of the ionization rendering them hard to distinguish.

Ionized fraction Power Spectrum
The third method for the volumetric analysis is the power spectrum of the ionized fraction field. The power spectrum measures the contribution from different spatial frequencies and therefore is sensitive to the underlying structures. We calculated the dimensionless power spectrum per comoving wavenumber ∆ 2 xx (k) = k 3 Pxx(k)/2π, where k [h Mpc −1 ] is the wavenumber and Pxx(k) is the spherically averaged square of the absolute value of the 3D Fourier transform of the ionized fraction in coeval volume. Figure 8 shows the power spectra for tow ionization stages, namely, at 70% and 90%. At the 70% ionization stage, the largest difference between the L1 and LLS cases is for wavenumbers smaller than ∼ 0.3 h Mpc −1 where L1 has approximately 1.6 times more power than LLS2 and 1.3 times more for LLS1. This is indicative of the presence of more structures on larger scales for L1. For smaller scales the differences are small. As reionization reaches 90%, L1 has ∼ 1.5 times less power at wavenumbers larger than ∼ 0.3 h Mpc −1 and 1.4 time more power around k ∼ 0.1 h Mpc −1 , implying that the LLS models retain more small scale structure and have somewhat less large scale structure in the ionization field.

OBSERVING REDSHIFTED 21-CM
Measurements of bubble sizes and shapes will require high signal to noise images at scales of a few arcminutes. This will be possible with the future Square Kilometre Array (SKA). The first generation of radio telescopes will instead focus on measuring statistical quantities of the 21-cm signal, such as the power spectrum. The discussion below explores the 21cm signal using the variance of the brightness temperature fluctuations and 21-cm power spectra. In addition, we also show the morphology of the 21cm signal in the so-called light cone slices.

Variance of the 21-cm background
The differential brightness temperature of the redshifted 21cm emission with respect to the CMB is given by the spin temperature, TS, of the neutral hydrogen and its density, ρHI, and in the limit of TS TCMB, is given by (Field 1959), where, z is the redshift, TCMB is the temperature of the CMB radiation at z, and τ is the corresponding 21-cm optical depth.
As seen in the previous analyses, the overall effect of the LLS is to slow down the ionization process and impede the growth of the H II regions. This should manifest itself as two observable properties in the variance of the fluctuations. Firstly, the peak of the brightness temperature fluctuations for the LLS simulations should be delayed and therefore should be visible at relatively higher frequencies.
Secondly, the amplitudes of the peaks should be diminished due to the relatively smaller size of the H II regions in the LLS simulations. Figure 9 shows the evolution of the RMS fluctuations of the mean differential brightness temperature for the three simulations as convolved with LOFAR-like boxcar beam of size 3 at a frequency bandwidth of 0.2 MHz. As depicted in the figure, at lower frequencies (early times) the fluctuations for all the three models are similar. The temperature fluctuations peak at 141 MHz for the L1 model and 147 MHz for the LLS1 and LLS2 models. The peak value of the brightness temperature RMS for the L1 model is 6.06 mK with the brightness temperature of 16.66 mK. The RMS is lower by about 9% and 8.7% for LLS1 and LLS2 models respectively. The RMS of the temperature fluctuations vs. the mass-weighted ionization global average is shown on the righthand side of Figure 9. The brightness temperature fluctuations are again seen following each other very closely at early times. However, as ionization reaches 20% the temperature fluctuations for different models start to diverge and peak at about 65-70% ionization. After this the temperature fluctuations decrease and are indistinguishable as reionization completes.
As mentioned earlier, these differences in brightness temperatures are manifested by the varying distribution and growth of H II regions in the different models as seen in the statistical analyses of previous sections. These fluctuations have been averaged over by a LOFAR-like beam and bandwidth. This is a simple first order estimate. Detailed and more accurate estimates require defining a noise budget including system temperatures, gains and phase errors, along with propagation effects (foregrounds, ionosphere etc.) and telescope based visibility sampling functions, see e.g. Patil et al. (2014). Expectedly, increasing bandwidth reduces the RMS as the fluctuations are smoothed out for a wider bandwidth. Increasing of the resolution of the beam increases and broadens the RMS. This is also expected as a smaller beam is sensitive to small scale fluctuations that are smoothed out by larger beams. Similar to the analyses in the previous sections, the differences are more pronounced between non-LLS and LLS models. However, based solely upon brightness temperature fluctuations it will be non-trivial to distinguish between the LLS models, see Figure 1. Figure 10 shows slices through the simulation cubes along the redshift (frequency) axis, also known as light cone slices. The ionized fraction of the simulation cubes is converted to the 21-cm emission differential brightness temperature, shown in log scale in mK, for the three models, shown from the top, L1, LLS1, and LLS2. For the desired range of redshifts, the data from the cubes is interpolated along the redshift/frequency axis. The interpolation is performed along the plane with an oblique angle of 10 • across the cubes in order to observe different structures along a random line of sight. The neutral regions are shown in red and the H II regions cover the dynamic range through blue as shown by Figure 11. The 3-D spherically averaged 21-cm differential brightness temperature fluctuation power spectra for models L1 (blue), LLS1 (red), and LLS2 (black) at mass weighted global ionization of 25% (solid lines), 50% (dashed lines), and 75% (dotted lines).

Evolution of the patchiness
the color bar of Figure 10. No corrections for redshift space distortions are applied.
From Figure 10 it is evident that at high redshifts the H II regions are small and distributed sparsely. These regions closely trace the ionizing halos. The effect of the different mfps of the the three models on the ionization becomes visually evident at redshift z ∼ 9.8 increasing with lower redshifts. The H II regions for the three models grow and merge at different paces. For this reason by redshift z ∼ 8.3 the L1 model is fully ionized, whereas, for the models with LLS, the mass weighted global average ionization for LLS1 is xm ∼ 83% and LLS2, xm ∼ 78%. The spatial axis of the coeval simulation box at redshift z = 9.457 subtends an angle of ∼ 0.97 • in the sky with each pixel of size 13.76 . Without any astrophysical and instrumentational propagation effects the H II regions at lower redshifts are large enough to be directly observed by arrays with ∼ 1 angular resolution capabilities. The effects due to the synthesized beam smooth the fine reionization structure but the statistical measurement of the signal is still achievable.

Power Spectrum of 21-cm
The power spectrum of the differential brightness temperature distribution, δT b , is defined as, where, δT b is the Fourier transform of the differential brightness temperature, δT b is the complex conjugate, P21 is the spherically averaged power spectrum, and δ D is the three-dimensional Dirac delta function representing the sampling function of the Fourier transformed quantity. The power spectrum is in units of mK 2 and is also used in dimensionless form as, Figure 12. The time evolution of the two k-modes, k = 0.05 and 0.9 h Mpc −1 of the spherically averaged 3-D power spectra of the differential brightness temperature fluctuations for the three models.
There are various schemes for estimating the power spectrum. In this paper, results from the mesh-to-mesh realto-redshift-space mapping (MM-RRM) methodology (Mao et al. 2012) are used. The MM-RRM uses the ionization fraction, density, and velocity data as input to estimate the power spectrum. This methodology takes into account the effects of redshift-space distortions and predicts accurate estimates of the 21-cm background with the caveat that at k (256) N /4 > 1.75 −1 h Mpc −1 the errors in the estimated PS are large.
The measurement of the power spectrum lends itself naturally to the radio interferometric observations since the visibilities of the interferometric measurement are sampling the Fourier transform of the sky and the power spectrum is the Fourier transform of the two-point correlation function. Figure 11 shows the dimensionless 21-cm differential brightness fluctuation power spectra of the three models at three representative stages, namely, at global ionization average of 25%, 50%, and 75%. The most distinct characteristic visible in the figure is that the fluctuations in the brightness temperature at large scales ( k < 0.2 h Mpc −1 ) grow by roughly 1.5 orders of magnitude for the fiducial L1 case and less than an order of magnitude of the LLS cases. This is a signature of the larger H II regions causing larger temperature fluctuations. The smaller fluctuations, on the other hand, flatten out as reionization progresses. Another noteworthy feature is the divergence of the fiducial L1 model from the LLS models at the largest scales traced by the simulations. Figure 12 shows the evolution of the power spectra with redshift for the three models for two k values, 0.05 and 0.9 h Mpc −1 representing large and small scale fluctuations respectively. The features in the 21-cm power spectra are consistent with the prior analyses. The recurring theme in the analysis of size and distribution of the H II regions is that for the LLS models the growth of ionized regions is obstructed and ionization is delayed. This is well captured in the 21-cm power spectra, defining implications for the upcoming experiments. For both the scales (k = 0.05 and 0.9 h Mpc −1 ) the observational frequency range from 140-150 MHz is where signals peak and the models are most differentiable. For large-scales (k = 0.05 h Mpc −1 ) the signal also goes through a minimum in the 123-127 MHz range with the lowest signal occurring for LLS2 followed by LLS1 and L1 models. While at this minimum the LLS models are only 15% lower than L1, a more significant difference occurs at the maximum where LLS1 and LLS2 show a 27% and 31% decrease compared to L1. As can be noted, the difference between the two models is 0.63 mK 2 . For small-scales (k = 0.9 h Mpc −1 ) the spectra peak very early on at z = 17.85 corresponding to a frequency of 75 MHz. This is region of the power spectrum where the models are virtually identical, but where they are also more likely to be affected by fluctuations in the spin temperature. The next peak for the small-scale features tracks the largescale but occurs slightly earlier around 133 MHz with the corresponding dip occurring at 117 MHz.
It is discernible from the power spectra that the contribution at small-and large-scales lag in the ionization process due to the reduction in the mfp caused by the dense LLSs. The delay in the overall ionization is also evident in the figure with the peak for the non LLS case rising at earlier times followed by models LLS1 and LLS2.

SPATIALLY-VARYING LLS OPTICAL DEPTH
In this section we examine the case of spatially-varying, or position dependent LLS distribution for the Songaila model (LLS1), titled LLS1-PD, and compare it with the corresponding homogeneous distribution case. The code implementation for both cases was discussed in §3.3. By construction, the average LLS optical depth and the number of ionizing photons emitted at any one time is the same in the two cases. This yields essentially identical evolutions of the number of photons per atom and very similar mean reionization histories (not shown). However, the H II region morphology is affected by the differing LLS distributions. In the spatially-varying case the LLSs are concentrated in the vicinity of the ionizing sources, acting like a type of a screen around them. This slows the growth of the ionized patches in the LLS1-PD case, resulting in more numerous, but smaller H II regions early on than are found in the homogeneous case (Fig. 13). In the former case there are also more partially or fully neutral gas pockets inside the ionized regions. At late times the larger ionized patches are more fragmented, with more structure.
These trends in the H II region sizes can be quantified using the the spherical averaging (SPA) technique (Figure 14). At xm = 50% for LLS1-PD the whole distribution is shifted noticeably to smaller sizes, with the peak probability distribution of the radii of the H II regions is less than by a factor of 2. In addition, what becomes more evident in the SPA analyses is the growth pattern. At xm = 70%, the probability distributions of LLS1 and LLS1-PD are virtually indistinguishable. This is an artifact of how the spheres are fit to the H II regions in the SPA technique, which among other things erases the small-scale structure. However, as ionization progresses to xm = 90% the LLS1-PD model again shows on average smaller regions, although still a very Figure 14. Probability distribution function R dp/dR per radial bins of spherical HII regions as measured by the spherical averaging algorithm for the LLS models (LLS1, red) and (LLS1-PD, green) The three sets are for the mass-weighted global ionization fraction xm = 50% (solid), 70% (dashed), and 90% (dotted). The threshold used is x th = 0.9. Figure 15. The time evolution of the two k-modes, k = 1.0 and 0.1 h Mpc −1 of the spherically averaged 3-D power spectra of the differential brightness temperature fluctuations for the LLS1 (red) and LLS1-PD (green) models. similar PDF shape. The lags of the peak amplitudes of PDF diminishes for both the models as the ionization progress. This indicates that the most typical size of the H II regions is almost the same, but there are still more small ionized patches in the LLS1-PD case, while the largest ones are smaller than in LLS1.
These changes in the H II region sizes should also be reflected in the 21-cm brightness temperature fluctuations. In Figure 15, we show the evolution of two k-modes of the 21-cm power spectra, k = 1 h Mpc −1 and k = 0.1h Mpc −1 , corresponding to small and large scale fluctuations, respectively. For both modes the peak fluctuations for LLS1-PD are shifted to slightly later times compared to LLS1. The differences are more pronounced for large H II regions (k = 0.1h Mpc −1 ), with the peak lower by ∼ 1.8 mK 2 and shifted by ∆ν = 3.7 MHz.

Photoionization rates
A noteworthy effect of the presence of the spatial variation of the LLS distribution, as shown in Figure 16, is that it reduces the mean photoionization rates in the ionized IGM during the peak and late reionization. This can potentially help alleviate the tension that sometimes exist between the measured relatively low values of the photoionization rate (implying a 'photon-starved' reionization) and the higher values often found in numerical simulations, which are required to complete reionization by z ∼ 6 and reproduce the measured electron-scattering optical depth. The lower values of the ultraviolet background as measured by Calverley et al. (2011) and Wyithe & Bolton (2011) at z ∼ 6, shown in Figure 16, are Γ = −12.84 ± 0.18 and −12.74 ± 0.3 respectively.
The presence of LLS systems increases the opacity in the dense regions and fewer photons can escape into the IGM, which notably reduces the mean photoionization rates regardless of the model details (see Figure 16). This reduction rises towards overlap and reaches a factor of ∼ 3. Interestingly, the position dependent LLS distribution further reduces the mean photoionization rates throughout most of the evolution, roughly doubling the effect compared to the uniform distribution cases (which are very similar to each other). However, the mean photoionization rate for the LLS1-PD case rises fast towards overlap and becomes roughly the same as the uniform LLS cases, while remaining much lower than in the L1 case.
Although the rise in the photoionization rate overshoots the observed values, the introduction of LLS clearly alleviates the discrepancy. A better match can likely be achieved by lowering the source efficiencies.
The spatial variation of the LLS distribution results in non-trivial changes in the reionization morphology and observable signatures. Even for a very similar reionization history the ionized patches become smaller and more fragmented, reducing slightly the 21-cm fluctuations, particularly at larger scales, and shifting the peak value to higher frequencies. The most interesting effect of the LLSs is that they reduce the mean photoionization rate in the IGM, which in effect is further enhanced by the spatial variation in the LLS distribution.

CONCLUSIONS
We have presented the results of the first large-scale (114 h −1 Mpc) numerical simulations of the epoch of reionization with sources of minimum halo mass 10 8 M in the presence of LLSs. The effects of LLS are implemented in C 2 -Ray using values of the mfp due to LLS extrapolated to the redshift range of reionization from lower redshift observational data, namely, Songaila & Cowie (2010) and McQuinn et al. (2011). Instead of imposing a horizon for the distance ionizing photons can travel, the effect of LLS is distributed over the entire IGM, homogeneously or position dependently, in the latter case based on the halo population.
We have analyzed the simulation results with different techniques to explore the underlying physics defining the size distributions, morphologies, and growth rate of the ionized regions in the presence of LLSs and to establish the efficacy of the observable parameters such as the brightness temperature fluctuations and 21-cm power spectra. Our main conclusions are as following, (i) We note that by introducing the dampening effects of the LLSs on the ionizing photons, the ionization process is delayed by ∆z ∼ 0.8 (for 99% ionization) for both the LLS models as compared with the fiducial non-LLS model. (ii) The photon statistics analysis shows that by the time reionization is complete, the LLS cases have used around 2.5 ionizing photons per baryon, whereas the case without LLS used 1.5. The reduction in the mfp due to LLS thus requires an additional photon per baryon to reionize the Universe. (iii) The topological differences between the large H II regions in the three models are visible in the simulation data at xm = 0.5 and indicate a slower merging of ionized regions in the presence of LLS. The friend-of-friend analysis of sizes of H II regions shows for all cases the emergence of two distinct populations of H II regions around xm = 0.1. However, the cases with LLS clearly retain more small H II regions. The results of the spherical averaging technique for characterizing the sizes of ionized bubbles also show that around xm = 0.5-0.7, the radii of the H II regions for the LLS cases are smaller by a factor of ∼1.2. (iv) We also note that while the sizes and mergers of the H II regions are impeded by the presence of LLSs, the "freeze out" of the size as reported in Sobacchi & Mesinger (2014) is not observed. Both the FoF and spherical averaging methods for characterizing bubble sizes show that the largest H II regions converge to the entire simulation box towards the end of reionization for both LLS models.
(v) The 21-cm brightness temperature power spectra show a factor of two less power in fluctuations at large scales (k < 0.1 h Mpc −1 ) for the LLS cases. (vi) The peak value of the RMS of the brightness temperature fluctuations over observable scales (angular scales larger than 3 and frequency intervals larger than 0.2 MHz), shows only a decrease of about 9% when LLS are present. (vii) The position dependent LLS results are similar to the homogeneous ones when one considers global quantities. However, the morphologies of ionized regions are affected both in overall shape and due to the presence of neutral islands. Also the mean photoionization rate in the ionized IGM is further reduced in the position dependent LLS case. We conclude that position dependent effects are likely to be important when considering LLS.
The method used to implement the LLS in this paper is an improvement over the more usual approach of imposing a hard horizon for ionizing photons. However, it remains ad hoc and depends on the assumed mfp evolution. Our method is useful for ensuring that the mfp in simulations remains bounded to some (redshfit dependent) value but does not capture all of the effects of the sinks of reionization. Major progress in this area will come from high resolution cosmological hydrodynamic radiative transfer simulations and, for the regimes for which these are not feasible, sub-grid models developed based on these. As reionization is a process taking place on the largest cosmological scales but is determined by both sources and sinks on the smallest, galactic scales, it will remain a challenging phenomenon to simulate for the many years to come.