The Impact of Cosmic Variance on Inferences of Global Neutral Fraction Derived from Lyα Luminosity Functions during Reionization

We investigate the impact of field-to-field variation, deriving from cosmic variance, in measured Lyα emitter (LAE) luminosity functions (LFs) and this variation’s impact on inferences of the neutral fraction of the intergalactic medium (IGM) during reionization. We post-process a z = 7 IGM simulation to populate the dark matter halos with LAEs. These LAEs have realistic UV magnitudes, Lyα fluxes, and Lyα line profiles. We calculate the attenuation of Lyα emission in universes with varying IGM neutral fraction, x¯HI . In a x¯HI = 0.3 simulation, we perform 100 realizations of a mock 2 deg2 survey with a redshift window Δz = 0.5 and flux limit f Lyα > 1 × 10−17 erg s−1 cm−2; such a survey is typical in depth and volume of the largest LAE surveys conducted today. For each realization, we compute the LAE LF and use it to recover the input x¯HI . Comparing the inferred values of x¯HI across the ensemble of the surveys, we find that cosmic variance, deriving from large-scale structure and variation in the neutral gas along the sightline, imposes a floor in the uncertainty of Δx¯HI∼0.2 when x¯HI = 0.3. We explore mitigation strategies to decrease this uncertainty, such as increasing the volume, decreasing the flux limit, or probing the volume with many independent fields. Increasing the area and/or depth of the survey does not mitigate the uncertainty, but composing a survey with many independent fields is effective. This finding highlights the best strategy for LAE surveys aiming at constraining the x¯HI of the universe during reionization.


INTRODUCTION
Reionization is the most recent phase change of the universe, in which the hydrogen in the intergalactic medium (IGM) went from being virtually fully neutral to nearly completely ionized.It is thought to have been complete within about the first billion years of the universe.It is expected to have been a patchy process, as ionized bubbles formed and grew around the first sources of ionizing radiation in the universe.Generally, it is thought that stars in the first galaxies provided the bulk of the necessary ionizing radiation to drive reionization (Finkelstein et al. 2015;Robertson et al. 2015;Bouwens et al. 2015a), though which galaxies dominate the ionizing photon budget is still being investigated.
Over the past several years, observations have been made which put constraints on the timeline of reionization using complementary methods sensitive to ionized and neutral hydrogen.Some examples: Planck-Collaboration et al. (2020) use the scattering signature of electrons freed during reionization on the cosmic microwave background to constrain the midpoint on instantaneous reionization to redshift z re = 7.67 ± 0.73.The fraction of dark pixels in the Lyα and Lyβ forests constrains reionization to be very nearly complete at z ∼ 6 (Becker et al. 2021;Qin et al. 2021;Bosman et al. 2022).Quasar spectra at z > 7 offer insight into the global neutral fraction (x HI ) of the IGM throughout reionization based on the IGM's Lyα damping wing imprints (Greig et al. 2019;Davies et al. 2018).
Galaxies, and in particular Lyα-emitters (LAEs) are also used to constrain the IGM hydrogen neutral fraction, and its evolution with time.Lyα is a resonant transition in hydrogen and Lyα photons scatter multiple times even in small amounts of neutral gas (Dijkstra 2014(Dijkstra , 2017)).Thus, observations of Lyα are particularly useful to probe the H i spatial distribution, column density and velocity fields.Observations of Lyα freely propagating through the IGM is indicative of a largely ionized IGM or a significantly redshifted Lyα line.
The LAE fraction, the fraction of galaxies which are found to emit Lyα, has been found to rapidly decline at z > 6 (Stark et al. 2010;Pentericci et al. 2011;Jung et al. 2018), and this decline is interpreted as a sign of increasing neutral hydrogen fraction in the IGM (Treu et al. 2013;Mesinger et al. 2015;Mason et al. 2018a).There are complications, however, in that the observed decreasing LAE fraction could also result from evolving galaxy properties, such as the escape fraction of Lyα photons (Dijkstra 2014), an increase in the number of Lyman Limit Systems (LLS) (Bolton & Haehnelt 2013), or other considerations (see Finkelstein (2016) for a review).Still, Mesinger et al. (2015) and Mason et al. (2018a) argue that the increasing neutral fraction of the IGM is the most plausible explanation, and in this context the rapid disappearance of Lyα at z > 6 paints a picture of late and rapid reionization (Ouchi et al. 2018;Hoag et al. 2019;Mason et al. 2019;Yoshioka et al. 2022).Such a timeline for reionization may favor UV bright, massive galaxies being the drivers of reionization (Naidu et al. 2020), as opposed to an undetected large population of faint objects (Finkelstein et al. 2019).
The moderate tension between these inferences on xHI can perhaps be attributed to cosmic variance.There is stochasticity in the number of LAEs observed in a given survey arising from both large scale structure and the inhomogeneity of reionization.
It is noteworthy that LAEs at z > 6 are fairly rare objects-dedicated narrowband surveys at z ∼ 7 typically detect tens of objects in volumes ∼ 1 × 10 6 cMpc (Ota et al. 2017;Itoh et al. 2018;Hu et al. 2019;Wold et al. 2021); these small number statistics could easily lead to field-to-field variations in the measurement of the z = 7 LAE LF.This effect has been observed in features like the LAE LF's "bright end bump" (Zheng et al. 2017;Hu et al. 2019;Wold et al. 2021), the result of very bright LAEs falling within the field of observation.The "bright end bump" is not observed in all LAE surveys at z = 7 (or even between different fields of the same survey), presumably either because of these objects inherent rarity or their concealment behind a neutral IGM.This is the crux of the issue we will explore in this paper: there is a degeneracy between the stochasticity in the number of LAEs observed and the inferred global neutral fraction of the universe.Indeed, Mesinger & Furlanetto (2008a) showed that, owing to the inhomogeneity of reionization, there is intrinsic scatter in the inferred global neutral fraction of the universe from observations of the Lyα damping wing in high-z quasars or gamma-ray bursts (GRBs).McQuinn et al. (2008) similarly found that a single high-z GRB could not place a constraint on the global neutral fraction with an uncertainty better than ∆x HI ∼ 0.3.Mason et al. (2018b) demonstrated the same principle with regards to Lyα equivalent widths; inference on the global neutral fraction is inherently stochastic when using observations of the equivalent widths (EWs) of Lyα in UV bright galaxies.Further, Mesinger & Furlanetto (2008a), Mc-Quinn et al. (2008), andMason et al. (2018b) all showed that the uncertainty on xHI is a function of xHI itself, ∆x HI (x HI ).Generally, the uncertainty tends to decrease as xHI increases because the size distribution of ionized bubbles gets smaller, since they have not yet had time to grow and merge in early reionization (Mason et al. 2018b).We will demonstrate that there is a similar effect, a floor in the uncertainty on inferring xHI , when using the evolving Lya LF and explore mitigation strategies to reduce this uncertainty.
We post-process a cosmological inhomogeneous reionization simulation (Mesinger & Furlanetto 2007;Mesinger et al. 2011Mesinger et al. , 2016)), populating the dark matter halos with simulated galaxies from empirical relations.We first verify our simulation against the intrinsic UV and Lyα luminosity functions (LFs) at z ∼ 6, then model the Lyα LFs after IGM attenuation at z ∼ 7. We perform mock observations on our simulation, reproducing the probed volumes of surveys which are carried out today.Across each of the many realizations of these mock observational programs, we measure the LAE LF and infer xHI .We investigate how the inferred z = 7 value of xHI changes as a result of statistical variance in the observed LAE LFs.
In §2 we give an overview of the simulation and the post processing we have performed.In §3 we explore the inherent uncertainty on inferred global neutral fraction at z ∼ 7 deriving from cosmic variance.We also explore if this uncertainty decreases by significantly increasing the survey area, depth, or using many independent fields.We conclude in §4.Throughout the paper we use the AB magnitude system and the Planck-Collaboration et al. ( 2016) cosmology with (Ω Λ , Ω m , H 0 ) = (0.69, 0.31, 68 km s −1 M pc −1 ),

THE REIONIZATION SIMULATION AND ITS POST-PROCESSING
In this section, we first summarize the simulation of dark matter halo masses, positions, and their Lyα optical depths.We then present an overview of our methods for assigning galaxy properties to the halos.

The 21cmFAST Simulation
We use nine custom simulations produced with the 21cmFASTv2 software (Mesinger & Furlanetto 2007;Mesinger et al. 2011Mesinger et al. , 2016)).21cmFASTv2 calculates the evolution of the hydrogen neutral fraction in the early universe using a semi-numerical approach.Inside the simulation cube, 1.6 cGpc on each side with a 1024 3 cell resolution, the code tracks dark matter and the phase of hydrogen in each cell while accounting for recombinations, photoheating star-formation suppression, supernova feedback, and radiation.Dark matter halos are identified from the 3072 3 initial conditions, and mapped to Eulerian positions at z=7 using perturbation theory (Mesinger & Furlanetto 2007).More detailed information can be found in Mesinger et al. (2016).Each simulation contains the comoving Cartesian coordinate positions of halos with masses ranging from 10 10.25 − 10 12 M .The halos start with only 14 discrete halo masses, evenly spread in log space.To avoid affects resulting from this discretized mass function, we redistribute the masses to produce a smooth distribution in log space, preserving each halos positions while approximating the original mass.The simulation also computes Lyα optical depths as a function of velocity offset from Lyα line center, τ (∆v), up to 525km s −1 (∆λ = 2.1 Å) in steps of 75km s −1 (∆λ ∼ 0.3 Å).These optical depths are calculated by integrating the neutral hydrogen along a 300 comoving Mpc path; 300 cMpc was chosen experimentally, as it ensures convergence of the optical depth (Mesinger & Furlanetto 2008b).Intrinsically, the optical depths are smooth functions of velocity offset from the line center, but we only have coarsely sampled optical depths at particular velocity offsets.To mitigate any effect from this coarse sampling, we interpolate the optical depths across velocity offsets when calculating the IGM attenuation on the Lyα line profile.
In total, we have nine z = 7 21cmFASTv2 simulations with different global neutral fractions, xHI , in the range 0.01 < xHI < 0.92.Each simulation takes the same z = 7 dark matter halo distribution and overlays on a different IGM ionization map, which corresponds to changing the ionization efficiency (McQuinn et al. 2007;Mason et al. 2018a).To get an intuitive picture, one may imagine that in a scenario with large xHI (most IGM hydrogen neutral), any given halo will have a small ionized bubble around it, if any.In a low xHI scenario (most IGM hydrogen ionized), those same halos will have larger ionized bubbles around them, which may have begun to overlap with neighboring halos' bubbles.This means that each individual halo's optical depth, at a given velocity offset, smoothly increases as neutral fraction increases, which allows us to interpolate between the neutral fraction files, obtaining a finer grid in neutral fraction than our original nine simulations.

Propagating Halo Growth
In order to make predictions at z = 7, we need to calibrate the Lyα LF using observations at lower redshifts (z ∼ 6), where there is evidence that reionization is very nearly complete (Qin et al. 2021).Though reionization may not be complete at z ∼ 6, residual xHI (i.e. if xHI ∼ 0.1) has a very small effect on the optical depth of Lyα and so the LAE LF will not be significantly changed (Mason et al. 2018a).As redshift decreases, halos grow because they accrete matter through mergers.Thus, we need to forward propagate the halo mass distribution from z = 7 to lower redshifts for proper comparison to observations.
We use the median halo growth trajectories from Behroozi et al. (2019) to grow our halos.We interpolate between the trajectories to create a fine grid in halo mass, then assign each of the z = 7 halo masses to a particular growth curve.To assign a growth trajectory to a halo mass, we find the curve which has z = 7 halo mass closest to the given halo mass in the simulation.Then, moving along the matched growth curves to the lower redshift of interest, we get the new masses.These new masses are then redistributed into a smooth log-space distribution.

Galaxy Properties
To achieve our goal of simulating LAEs, we need to assign intrinsic Lyα luminosities and Lyα line properties to the dark matter halos.Our road map will be: 2. Assign a Lyα EW (EW Lyα ) to each halo, using a UV dependent EW probability distribution function (PDF) ( §2.3.2) 3. Assign Lyα line profiles, with line velocity offset and width dependent on the halo masses ( §2.3.3)Item 3 is required to compute the opacity to Lyα photons from the neutral IGM.In what follows, we first describe the methodology used to assign M UV and Lyα EW to each halo and how this methodology was validated with real data.We then discuss the new prescription to assign the Lyα line profile and the comparison with observations.

The M h -MUV Relation
We use the M h -M UV relation from Mason et al. (2015) to assign UV magnitudes to our halos.This relation has an in-built redshift dependence, since halo distributions and their growth rates change with redshift (Mason et al. 2015).We use relations tabulated at z = 6 and z = 6.8 in this paper.Considering our small range of redshifts of interest (5.7 < z < 7) and the relatively small evolution between the z = 6 and z = 6.8 M h -M UV relation, we do not derive a new M h -M UV relation for every redshift we consider.Instead, we apply the z = 6 relation if our simulation redshift is 5.7 < z < 6.4 and z = 6.8 relation if our simulation redshift is 6.4 < z < 7.
The relations provided in Mason et al. (2015) do not include a prescription for the scatter, so we follow methods akin to Ren et al. (2019) and Whitler et al. (2020) and assume a normal scatter in M UV .The specific value of the scatter is chosen so that the resulting UV LF reproduces the measured LFs at z ∼ 6 and z ∼ 7 simultaneously.In Figure 1, we show UV LF observations at z=6 in black and z=7 in gray.At both redshifts, we plot three UV LFs resulting from different scatter prescriptions: 0.1 in green, 0.3 in blue, 0.5 in orange.Solid lines are used for z=6, and dotted lines for z=7.We apply each relation with its scatter 50 times and indicate the range within which 95% of the UV LFs fall with the shaded regions.
Motivated by the agreement between the UV LFs with σ M U V = 0.3 and the observations at z = 6 and z = 7, we fix the value of σ M U V to 0.3.The fact that we are able to reproduce UV LFs at two redshifts simultaneously with only one varying parameter is encouraging; it validates both our halo growth prescription and the M h -M UV relation.

MUV-EWLyα Distribution
The distribution of rest-frame EW Lyα given M UV takes the form Mason et al. (2018a).H(EW Lyα ) is the Heaviside step function and δ(EW Lyα ) is the Dirac delta function.Mason et al. (2018a) fit the parameters A and W c using  2017) were photometrically selected using the Lyman break, so called Lyman break galaxies (LBGs), and it is worth noting that not all LBGs show Lyα emission.The parameter A accounts for the fraction of galaxies which do not emit Lyα and W c determines the exponential decline of the probability distribution function towards larger EW Lyα .Both parameters are allowed to vary as a function of M UV , and we take their value from Mason et al. (2018a).To give a sense of the resulting distribution, the log 2D distribution of the M UV -EW Lyα plane for z = 6.6 galaxies, along with traces for constant Lyα flux limits, are shown in Figure 2. UV faint galaxies have a larger probability of taking on a large EW Lyα value than UV bright galaxies.
We are inherently assuming that the underlying EW Lyα distributions are the same at z ∼ 6, where the De Barros et al. (2017) observations are made, and at z ∼ 7, where we will apply the EW Lyα distribution.Mason et al. (2018a) argue that, while oversimplifying, this assumption is justified by the relatively short time frame; less than 200 Myr pass between z = 6 and z = 7.Further, this means that any evolution in the EW Lyα distribution between z = 6 and z = 7 would be attributed to an evolving xHI in our model, as explored in Mason et al. (2018a).

Lyα Line Profile
Since the Lyα optical depth due to the neutral IGM depends on the wavelength offset from the line center, we need to assign a Lyα line profile that accounts for radia-tive transfer effects in the galaxies' interstellar medium (ISM) and circumgalactic medium (CGM).We will take a relatively simple approach in modelling the Lyα line profiles, modifying the methods in Mason et al. (2018a).
We assume that z ∼ 6 Lyα line profiles are statistically the same as those at z = 7 when ignoring attenuation from a neutral IGM.Under this assumption, we take observed Lyα line properties at z ∼ 6, correlate them with observable galaxy properties, and assign them to our z = 7 halos.The line profiles we will assign, then, have the effects of the ISM and CGM built-in.Attenuation by the IGM could change the observed Lyα line profile statistics at z = 7.
We assume that the line profile modified by the ISM and CGM is a Gaussian with an offset with respect to the rest-frame line center (∆v Lyα ).Unlike Mason et al. (2018a), we assign the velocity offset with respect to the galaxy rest-frame using an empirical relationship from Hayes et al. (2021), which is explained further in the next subsection.
The line's full width half maximum (FWHM) is set equal to ∆v Lyα , as found in Verhamme et al. (2018).The Lyα line profiles are truncated blueward of the Lyα rest-frame line center because anything blueward of Lyα line center will redshift into resonance and be scattered by residual neutral fraction.The Lyα line profile after leaving the ISM and CGM is (2) Note, though, that roughly one quarter of ultraluminous (log(L Lyα /[ergs s −1 ]) > 43.5) LAEs at z > 6 emit flux blueward of line center (Hu et al. 2016;Songaila et al. 2018;Meyer et al. 2020).These LAEs likely reside in large ionized bubbles in the IGM.These bubbles may have been carved out by ionizing radiation from the ultraluminous LAEs themselves, a population of nearby ionizing sources, an obscured quaser, or a combination of contributors (see Bagley et al. (2017) and Mason & Gronke (2020) for a discussion).Due to their rarity, we do not modify our line profiles to account for these ultraluminous galaxies.
Lyα Velocity Offsets -Hayes et al. ( 2021) compiled 229 Lyα selected galaxies at 2.9 < z < 6.6 from observations with the Very Large Telescope's (VLT) Multi Unit Spectroscopic Explorer (MUSE) and 74 galaxies at z < 0.44 from Hubble's Cosmic Origins Spectrograph (COS) and compared the two populations' Lyα line profiles.The VLT galaxies' Lyα velocity offsets relative to their systemic redshifts (i.e. the line profiles' 1st moments) were found by calibrating the observations with low-z data and a model (Runnholm et al. 2021).These velocity offsets, ∆v Lyα , and the Lyα luminosities are correlated.The data are described in more detail in Hayes et al. (2021).Note that these galaxies were selected with the Lyα line, while the De Barros et al. ( 2017) sample which provided the EW Lyα distribution were selected with the Lyman break.While not all LBGs have Lyα emission, we have captured this behavior in the EW Lyα distribution; those that do have Lyα emission may be detected in the Hayes et al. (2021) sample, and so these empirical velocity offsets from Lyα selected galaxies may be applied to LBGs as well.Note, however, that we find that the transmission fraction of Lyα is largely insensitive to the exact velocity offset value in our model, provided that the line has a width and so some of the flux is away from line center.For example, if we set the velocity offsets to 0 km/s, drastically different from the typical values from our relation of 200-300 km/s, we see a maximum of a 6% difference in the Lyα LFs in the range from 42.4 < log 10 (L Lyα ) < 43.7 for xHI =0.3 (this luminosity range was chosen because of the observational range of the Lyα LFs at z = 7, see Figure 5).As we will only be concerned with Lyα luminosities larger than L Lyα 10 42 ergs s −1 , the fact that our relation will predict a negative velocity offset at Lyα luminosities 4 × 10 39 ergs s −1 is not an issue.
With each halo assigned a M UV , EW Lyα , Lyα optical depth (τ IGM ), and Lyα line profile, we have everything we need to compute the transmitted fraction of Lyα flux.Explicitly, the fractional transmission through the IGM is: (4)

Validation
With all the pieces required to simulate Lyα-emitters in place, we proceed to validate the recipes by comparing  Figure 5 shows the LAE LF calculated from the full simulation volume for various xHI values.As xHI increases, attenuation of intrinsic Lyα emission increases and the resulting LAE LF is suppressed; between a fully ionized universe and a universe with xHI ∼ 0.85, the suppression is about an order of magnitude.The observational data at z = 7 agree well with a fully ionized universe, which implies that in our simulation, the evolution of the Lyα LF from z = 5.7 to z = 7 in observations can be fully explained via halo growth and changing intrinsic Lyα emission.Previous works, such as Dijkstra et al. (2014), have argued that at least some of the Lyα LF evolution likely comes from evolution of galaxy properties.

COSMIC VARIANCE IN LAES
With our simulation constructed, we begin analysis to determine the importance of cosmic variance on the inference of global neutral fraction.As laid out in the introduction, Lyα photons are easily scattered by neutral hydrogen-this has the effect of decreasing the number of observed LAEs when neutral fraction of H i increases.This decreases the volume density of observed LAEs.The differential change in the Lyα LF between redshifts has been used to infer the global neutral hydrogen fraction of the universe (Malhotra & Rhoads 2004;Ouchi et al. 2010;Kashikawa et al. 2011;Konno et al. 2014;Zheng et al. 2017;Ota et al. 2017;Itoh et al. 2018;Konno et al. 2018;Hu et al. 2019;Wold et al. 2021; The observational data at z = 7 is consistent with a completely or mostly ionized IGM, as demonstrated with the data's agreement with the blue and orange curves.This, combined with the z = 5.7 Lyα LF agreement, shows that the evolution of the z = 5.7 to z = 7 Lyα LF is explainable with smaller halos and lower intrinsic Lyα emission, rather than an increase in global neutral fraction in our simulation.Morales et al. 2021).Our goal is to quantify how cosmic variance impacts this inference and explore mitigation strategies in survey design with an updated, empirical LAE model.

A Nominal LAE Survey
To gain some footing, we first consider a nominal survey, representative of volumes probed by typical surveys carried out by ground-based telescopes equipped with narrowband filters.We will "observe" our simulation with the specifics of the nominal survey and use Bayesian inference to estimate xHI from the "observed" Lyα LF.
We simulate a circular survey covering 2 square degrees to a Lyα flux limit of f Lyα > 1×10 −17 ergs s −1 cm −2 (or L Lyα 6 × 10 42 ergs s −1 at z = 7).Such a survey pushes faint enough in flux/luminosity to detect LAEs fainter than the knee of the Lyα LF at z=7 (Ota et al. 2017;Itoh et al. 2018;Hu et al. 2019).This depth is comparable to those reached by narrowband surveys, and so provides an interesting baseline.Notably, we will consider a LAE "observed" simply if it falls into the geometric cutout and exceeds the minimum line flux limit; in this way, our survey ignores any effects of contamination or incompleteness, and so any stochasticity on the inferred values of xHI between surveys is the result of stochasticity in the number of observed LAEs alone.We impose a redshift window of ∆z ∼ 0.5 centered at z = 7 by selecting galaxies in a slice of 160 cMpc within the simulated cube.Such a survey has a volume of 7.58 × 10 6 cMpc 3 at z = 7, a volume comparable to z = 7 narrowband surveys which have larger areas but smaller redshift window (Ota et al. 2017;Itoh et al. 2018;Hu et al. 2019;Wold et al. 2021).One should keep in mind that the simulation is a snapshot at z = 7; there is no differential evolution between the "front" end of our survey and the "back" end taken into account.This allows us isolate the affect of cosmic variance by treating the global neutral fraction as a constant over our survey window.We then have a well-defined global neutral fraction that should be inferred in our simulation from the mock observations.Note, however, that this lack of differential evolution is not important with regards to any specific LAE.The cumulative optical depth is dominated by the gas nearby the galaxy, so much so that the neutral patches with optical depth 1 are within 10 Mpc, corresponding to ∆z ≈ 0.02 at z = 7 (Mesinger & Furlanetto 2008a).This implies that the differential evolution is unimportant when calculating the optical depth for a given LAE.
We create a z = 7 simulation with a realistic input global neutral fraction of xHI = 0.3 (Morales et al. 2021) and observe an area in a randomized location within the cube.This randomized location can, in principal, have a different local neutral fraction than the input neutral fraction, but note that we are interested in inferring the latter, rather than the former.This survey results in a total sample of 274 LAEs.The LAE LF for one specific realization is shown in Figure 6.We split the LAEs into bins of width ∆log(L Lyα ) = 0.1.The model Lyα LFs (as measured from the entire simulation volume) for xHI = 0.01, 0.30 (input value), and 0.50 are shown as the blue, black, and green curves, respectively.The error bars associated with the simulated LF measurements come from Poisson statistics on the counts in each bin if N LAE > 20.For smaller numbers of LAEs, we use the estimates tabulated in Gehrels (1986) for small numer Poisson errors.

Inference on xHI from the Nominal Survey
In order to speed up the inference on xHI , we fit a function to the model LAE LFs from Section 2.5 as a function of xHI .More information on this analytical function can be found in Appendix A. Some example LAE LFs for different xHI are visualized in Figure 6.
We use a Bayesian approach (implemented using the MCMC fitting from the pymc3 package) to estimate the best fit hydrogen neutral fraction, xHI , given our model Lyα LFs (Section 2.5) and the simulated LF measurement (Section 3.1).The 2D surface fit in Appendix A is used to create the predicted LAE volume densities, µ, for a given neutral fraction.Then, the likelihood is given by where the subscript i denotes each of the measured luminosity bins.The unknown model parameters to be fit are the precision, τ , of the observed LAE LF data and xHI .We consider a Gamma function prior on τ with α = β = 1 and a flat prior on xHI .The posterior distribution function for xHI is then given by where p(τ ) is the Gamma prior on τ and p(x HI ) is the uniform prior on xHI .
We note that we are assuming here that the precision of all of the LF measurements, τ , is the same.The results of our paper do not change if we instead assign adaptive luminosity bins to ensure that this assumption is true, that the precision between all the measurements is the same, i.e. the same number of LAEs go into each bin.We use fixed bin sizes in log(L Lyα /[ergs s −1 ]), rather than these adaptive bins, for ease of comparison to measurements in literature.
The posterior on xHI for the LAE LF shown in Figure 6 is displayed in Figure 7.The input neutral fraction is denoted with the vertical black dashed line.We can see that the width of the posterior is rather wide and the posterior PDF continues increasing toward low-values of xHI .The posterior's median is indicated with the vertical dashed blue line-it underestimates the real neutral hydrogen fraction.Looking back at the LAE LF in Figure 6, the origin of this offset is apparent: more of the observed LF data points lie above the input xHI than below it, and a higher number density of LAEs leads to an inference of a lower xHI .In other words, the global neutral fraction is underestimated because the number of LAEs in the specific realization of the nominal area is higher than the predicted number, i.e. the randomized survey position fell on a slightly over-dense region of the simulation (with respect to LAEs).
Motivated by this result, in the following sections we explore the range of LAE densities expected at z = 7, quantify how this range of environments propagates to a range on the inferred xHI , and then test methods to mitigate the effect of varying environments on the inference of xHI .

The Effect of Cosmic Variance on xHI Inference
We simulate 10,000 realizations of the nominal survey scattered randomly about the simulation volume, always with a global neutral fraction xHI = 0.3, and consider the total number of LAEs, N LAE , in each realization which meet the detection threshold.The distribution of N LAE is shown in blue in the left panel of Figure 8.For comparison, a Poisson distribution with the same population mean is shown in black; two particular survey realizations, one at the 2.5 and one at the 97.5 percentile of the distribution of N LAE , are noted with green and orange vertical lines, respectively.We will look at these particular realizations to determine a conservative (i.e.covering 95% of surveys) absolute variance in xHI inferred from the different over-and under-dense environments.
We calculate the posteriors on xHI from the LAE LF measured from the over-and under-dense regions and show them in the right panel of Figure 8; orange shows the posterior on xHI inferred from the LAE rich region, and green shows xHI inferred from the LAE sparse region.The difference between the two medians is quite large, about ∆x HI ∼ 0.19.This difference is critically important, as it is impossible to know a priori if a particular survey has landed in a region which is over-dense, under-dense, or something close to the average density, and this has a direct effect on the inferred xHI .This issue is still present for studies which use the differential evolution of the LAE LF to infer xHI , because we cannot know if the comparative surveys are in over-or underdense regions of the universe with respect to LAEs.
This inference was done in a very ideal scenario, where there is no impurity, no incompleteness, and the model LAE LFs we are comparing to our "observations" are from the same simulation, and so we know we have the physics in the modelling precisely correct.In reality, there would be additional uncertainty entering from the fact that comparison models will make assumptions on the underlying physics.Despite this extremely optimistic scenario, ∆x HI ∼ 0.19 at z = 7 when inferring from our nominal survey's LAE LF, a rather large uncertainty.It is important to note that this uncertainty in xHI is a function of xHI itself, ∆x HI (x HI ).We find, in agreement with Mason et al. (2018b) that the uncertainty decreases as xHI increases, owing to a narrower distribution of ionized bubbles sizes at high xHI values.We save a full quantification of ∆x HI (x HI ), akin to the work done for quasars and GRBs in Mesinger & Furlanetto (2008a), for future work, and instead take the uncertainty at xHI = 0.3 as a baseline and investigate mitigation strategies.
It is noteworthy that this nominal survey is already fairly large and pushes to moderate flux depths.Taking into account the considered redshift window, 6.75 < z < 7.25, the survey's volume, 7.58 × 10 6 cMpc 3 , is larger than many of the z = 7 LAE surveys observed to date, and comparable to Wold et al. (2021)  systematic uncertainty, so we now explore how ∆x HI changes when varying the survey strategy.

Mitigating the Systematics
We consider three scenarios to mitigate systematic uncertainty on the inference of xHI : increasing the survey area, increasing the survey depth, and considering a survey composed of many independent fields.

Increasing Area
To determine survey areas which may be effective at decreasing the effect of cosmic variance on the inference of xHI , we consider four increasingly large areas centered on the two positions identified in Figure 8, i.e. an overdense and an under-dense region.We simulate areas of 2, 5, 10, and 20 square degrees, while keeping the survey line-flux limit fixed to the nominal survey.For each survey, we use the same procedure described in Section 3.2 to infer the neutral hydrogen fraction.The resulting posteriors for both the over-and under-dense regions are shown in Figure 9.
Figure 9 shows that increasing the area by a factor of 5 to 10 has a tendency to modestly decrease the width of the posterior-in the LAE sparse region, the standard deviation of the 2 square degree survey's posterior is 0.15, while the 20 square degree survey's standard deviation is 0.13.
We take the extremely optimistic case, an observationally perfect 20 square degree survey with ∆z = 0.5 and luminosity limit L Lyα > 6×10 42 ergss −1 and investigate it further.
We compute 100 realizations of this survey, randomly positioning them within the simulated volume, and infer xHI from each.We calculate the median for each of these 100 posteriors and plot the distribution in Figure 10.For a baseline comparison, we repeat the same procedure with the nominal 2 square degree survey (which has the same flux limit as this 20 square degree survey); 100 realizations are created and each is used to make an inference on xHI .The medians of the 100 nominal survey posteriors are plotted in Figure 10 as the underlying black distribution.95% of the nominal survey's medians are within a range ∆x HI = 0.25, while the 20 square degree medians have 95% of their realizations within ∆x HI = 0.21.The width of the distribution of these medians gives us a sense of the systematic uncertainty in xHI arising from cosmic variance between the two survey strategies; there is an extremely modest tightening of posteriors for the 20 square degree surveys compared to the 2 square degree surveys, which does not seem to be statistically significant.It is important to emphasize that the narrowing of the posteriors in Figure 9 was for one particular survey only, whereas Figure 10 shows that the uncertainty arising from cosmic variance is not improved by increasing the area.This is demonstrated by the fact that the distributions of the medians of the 100 posteriors for the 2 and 20 square degrees surveys are approximately equivalent.
Thus, even increasing our area by a factor of 10x to 20 square degrees, considering a fairly wide redshift window, and detecting LAEs fainter than the luminosity function knee at z=7 is not enough to significantly re-  duce the uncertainty on xHI .For surveys larger than 20 square degrees, we begin to get to the point where even our extremely large simulation cannot provide a large number of independent fields of view; we thus cannot conclusively say what volume is necessary to drive the uncertainty on xHI below ∼ 0.2 for a survey with flux limit f Lyα > 1 × 10 −17 ergs s −1 cm −2 , we can only say that it is larger than 7.58 × 10 7 cMpc 3 , a significantly larger volume than those which have been probed in LAE surveys so far.Evidently, simply taking a survey with a larger area is not enough to mitigate the effect of cosmic variance.

Increasing Depth
One may wonder what effect additional depth would have on the inference of xHI , motivated by the fact that less luminous LAEs will generally occupy less massive halos, and so be less biased, and thus be less sensitive to cosmic variance.To investigate this idea, we repeat the exercise from Section 3.4.1,but modify the survey-we instead consider a 10 square degree survey with a flux limit of f Lyα > 3 × 10 −18 ergs s −1 cm −2 (or L Lyα 1.8 × 10 42 ergs s −1 at z = 7).This survey's flux limit is about 3.3x deeper than the nominal survey and it is 5x larger, so it would take about 55x as long to execute.Going deeper in flux results in a drastic increase in the number of LAEs observed due to the shape of the luminosity function; going from a flux limit of f Lyα > 1 × 10 −17 ergs s −1 cm −2 to f Lyα > 3 × 10 −18 ergs s −1 cm −2 increases the number of observed LAEs by about a factor of 10.Following the steps from the previous section, we calculate the pos- teriors on xHI from 100 realizations of such a survey; again, the distribution of the medians for the 100 surveys are shown in Figure 11.Again, the distribution for the 100 nominal surveys is shown in black as a baseline to compare against.
The width which encompass 95% of the posterior medians for the deep survey is virtually the same as the nominal survey, ∆x HI = 0.23.Thus, even pushing ∼ 3.3x deeper in flux and increasing the area by 5x over the nominal survey, requiring ∼55x more exposure time, has not decreased the systematic uncertainty.

Independent Fields
Lastly, we consider a survey composed of many independent fields, the motivation being that each field will probe a different environment for the LAEs, averaging out high-and low-density environments.
Further motivation can be seen examining Figure 12.The left panel shows the distributions of the number of LAEs per square degree observed by 1000 surveys of areas 2 and 20 square degrees (and to depth L Lyα 6 × 10 42 ergs s −1 ).The blue dashed line indicates a 2 square degree survey realization centered on an under-dense region (at the 2.5 percentile) and the red dashed line shows a 20 square degree survey centered on that same position7.As expected, the distributions of the number of LAEs per square degree converges towards the overall simulation mean as the survey area increases; the square degree surveys' distribution (red) is far narrower than the two square degree surveys' (blue).
However, the imprint of the under-dense region identified with a 2 square degree survey is still dominating the number of LAEs observed in the 20 square degree survey when compared to the parent distribution, as shown in the right panel.The 20 square degree survey centered on the under-dense region (dashed red vertical line) is extremely close to the parent distribution's 2.5 percentile (solid black vertical line).This is indicative of the fact that adding area to a given survey has a tendency to add the average number of LAEs per square degree, by definition.This will have no effect on moving a given realization towards the center of the distribution of N LAE for the parent distribution of that survey size; only adding an underdense or over-dense region changes its position relative to the parent distribution.Then, a survey centered on an under-dense (over-dense) region will continue to exhibit low (high) LAE counts relative to other similar surveys, even as the area increases drastically.
Instead, we simulate a survey of 100 independent 0.2 square degree fields, totalling 20 square degrees, to a depth of f Lyα > 3 × 10 −18 ergs s −1 cm −2 , the same total area as our largest survey and as deep as our deepest simulated survey.All of these surveys are simulated in a simulation with global neutral fraction xHI =0.3.Each of the 100 0.2 square degree fields has the inference done on them separately, and the resulting 100 posteriors are multiplied together to form a joint posterior on xHI .
Running a survey in this particular manner is computationally expensive, and so we limit ourselves to 90 realizations of such a survey; we show the distribution of the posterior medians in Figure 13.The width which encompasses 95% of the surveys has decreased-∆x HI = 0.05.Probing many independent fields is the only effective way to decrease the systematic uncertainty introduced by cosmic variance, decreasing the 95% width by about a factor of 4. Note that the inferred global neutral fraction is biased slightly high by about 0.03.This is because the analytical fit to the model Lyα LFs is only accurate to 4%.In practice, it would be possible to remove this bias by bypassing the step of fitting the LFs analytically; however, for our analysis we need to compare many observing strategies and do hundreds of instances of inference for each.The analytic function speeds up the inference process and is the only way to complete our tests in a feasible amount of time.
This method has an additional advantage; Mesinger & Furlanetto (2008b) showed that counts-in-cell statistics can be an effective way to use the enhanced clustering of LAEs during reionization to constrain xHI .While the counts-in-cells method does not strictly require that the fields be discontiguous, a survey strategy such as ours   has the additional benefit of overcoming any effect of cosmic variance, which could impact a clustering analysis as well.Thus, such a survey would have two methods of constraining xHI , both of which are robust against cosmic variance, which could be checked against each other for consistency.

The Implication of the Mitigation Techniques
We are left to conclude that the only effective way to break through the floor in the uncertainty in xHI arising from the cosmic variance of LAEs is with surveys composed of many independent fields.Drastically increasing the survey area and pushing to much deeper fluxes alone proved ineffective.Of note, the survey composed of independent fields, with a flux limit of f Lyα > 3 × 10 −18 ergs s −1 cm −2 and an area of 20 square degrees would take ∼ 110 times as long to carry out as the nominal survey, which is currently representative of the largest dedicated LAE surveys to date-thus, this serves more as a proof of concept than a realistic goal for a survey design in the immediate future.
It is worth keeping in mind that an uncertainty of ∆x HI = 0.2, as we have demonstrated is characteristic when the universe has xHI = 0.3, is still a competitive constraint.Particularly, if such measurements are made at a range of redshifts, we may be able to construct a fairly constraining timeline of reionization; thus, the LAE LF will continue to be a useful tool in constraining reionization, in addition to other techniques, such as damping wings in quasars and GRBs, LAE clustering, LAE fraction, and Lyα EW distributions.

CONCLUSIONS
We have post-processed a z = 7 comoving 1.6Gpc 3 simulation to create a mock LAE catalog using empirical relations.Our simulation simultaneously reproduces the z = 6 and z = 7 UV LFs and the z = 5.7 Lyα LF.
We produce z = 7 Lyα LF models, ranging from 42 < log(L Lyα /[ergs s −1 ]) < 44, as a function of xHI .We created mock surveys with different observation strategies and compare the measured Lyα LFs to these models, inferring posteriors on xHI .The observed LAE LF varies, deriving from the large scale structure and stochasticity in the amount of neutral gas along the sightline between the LAEs and the observer, and the shape and position of the peak of the posterior distribution function on xHI changes in turn.
We show that the precision of xHI estimates drawn from an LAE LF derived from a single 2 square degree field with f Lyα > 1×10 −17 ergs s −1 cm −2 is limited by the cosmic variance of LAEs.We find an uncertainty floor of ∆x HI ∼ 0.2 when xHI = 0.3.Note that a 2 square degree survey with a redshift window of ∆z = 0.5 has a volume of 7.58 × 10 6 cMpc 3 , about the same size as existing narrowband LAE surveys.
We investigated three methods to push this floor in the uncertainty down.(1) Increasing the area covered by a factor of 10, (2) pushing to ∼ 3.3x deeper flux limits, and (3) breaking up the survey into independent fields.For the added area and deeper survey strategies, the variance in the medians of the posteriors on xHI remains virtually unchanged from the nominal 2 square degree survey.This demonstrates that the xHI inferred from observed LFs in contiguous fields are largely at the whim of what cosmic environment one's survey happened to land in.
However, a large, deep survey composed of smaller independent fields proved effective in reducing the systematic uncertainty on xHI .Under this strategy, we see width of 95% of the posterior medians of xHI decrease from ∆x HI ∼ 0.2 to ∆x HI ∼ 0.05.Thus, probing multiple independent fields is critically important in constraining xHI .
This result demonstrates that surveys of LAEs aiming at constraining xHI should adopt a strategy similar to the one proposed in Section 3.4.3 in order to minimize the systematic uncertainty on xHI .The Parallel Application of Slitless Spectroscopy to Analyze Galaxy Evolution (PASSAGE, PI M. Malkan) JWST survey, as a pure parallel survey, follow this strategy, though JWST's small field-of-view implies that the total area will be much smaller than the survey proposed in Section 3.4.3.PASSAGE is being carried out during Cycle 1 of JWST observations and its observations of high-z LAES will allow the inference of new information about the timeline of reionization.The Nancy Grace Roman Space Telescope, with its large field of view, is a very good candidate to run a survey analogous to the one discussed in Section 3.4.3,though with a higher flux limit.

APPENDIX
A. LUMINOSITY FUNCTION ANALYTICAL FIT We assume the following 2D functional form, which reproduces the logs of the modelled Lyα LFs to 4% accuracy with a median -0.1 % difference across all LAE LF bins and neutral fractions:   -2436.6 115.7 -2563.8 -1.4 -1.4 3.3 6146.7 -286.6 120.0

Figure 1 .
Figure1.Three z = 6 UV LFs, resulting from three different M h -MUV scatter values, are plotted in solid green (σM U V = 0.1), blue (σM U V = 0.3), and orange (σM U V = 0.5).We calculate the UV LF for each σM U V value across 50 different realizations of the simulation and show the region within which 95% of UV LFs fall with the shaded areas.The same color scheme is used for the z=7 UV LFs, which are plotted as dotted lines.The z=6 data are plotted as white markers with black borders; the squares areBouwens et al. (2021) and the circles are Ono et al. (2018).The z=7 data are shown in gray; the squares are Finkelstein et al. (2015), stars are Atek et al. (2015), circles are Bouwens et al. (2011), pentagons are Bouwens et al. (2015b), crosses are Bowler et al. (2014), diamonds are Castellano et al. (2010), down triangles are Livermore et al. (2017), up triangles are McLure et al.(2013), left pointing triangles areOuchi et al. (2009), right pointing triangles areSchenker et al. (2013), and tri-down areTilvi et al. (2013).The M h -MUV relation with scatter of σM U V = 0.3 reproduces the z=6 and z=7 data simultaneously.The bright end provides the demarcation between models with different σM U V values.

Figure 2 .
Figure 2. The log 2D distribution of galaxies at z=6.6 in the MUV − EWLyα plane.The colored lines denoted in the legend indicate traces of constant limiting Lyα line flux in ergs s −1 cm −2 .A survey with a given line flux limit could detect LAEs which lie above the trace corresponding to its flux limit.Note that there is an abundance of galaxies at EWLyα=0 along the very bottom of the distribution, corresponding to the non-Lyα emitting population.
For xHI =0.6, this percent difference grows to 20%.We use emcee to fit a linear relation to the observed high-z LAEs' line profiles' ∆v Lyα and the log of their Lyα luminosities.The best fit linear relation is given by ∆v Lyα = −3156 +882 −895 + 80 +21 −21 log 10 (L Lyα ) (3) with a Gaussian scatter of 13.5 +6.1 −3.5 km/s.The best fit relation, with the MUSE data, is shown in Figure 3.The z < 3.5 MUSE data are shown in orange, and the z > 5.5 data are shown in green.Both redshift ranges are consistent with the same linear relation, indicating little to no redshift evolution.This linear relation is used to center the Gaussian before truncation; random Gaussian scatter of 13.5 km/s is added to the offsets.

Figure 3 .
Figure 3.We show the 1st moment of the Lyα lines versus the log of the Lyα luminosity for the MUSE data at z < 3.5 (green) and z > 5.5 (red).Our linear fit is shown with the blue line, while the black hashed lines show the 1 − σ Gaussian scatter about the relation that we fit with emcee.

Figure 4 .
Figure 4.The z = 5.7 Lyα LF is shown in blue with constraints from Konno et al. (2018) and Ouchi et al. (2010).We see very good agreement between our simulation and the observational constraints.

Figure 5 .
Figure 5.The curves show the LAE LF after IGM attenuation for various global neutral fractions, xHI, at z = 7.The shaded regions show the range in which 95% of 30 simulation realizations reside.Observational data from Ota et al. (2017); Itoh et al. (2018); Hu et al. (2019) are shown with various markers.The data have been offset by 0.01 with respect to each other on the abscissa to aid with visualization.The observational data at z = 7 is consistent with a completely or mostly ionized IGM, as demonstrated with the data's agreement with the blue and orange curves.This, combined with the z = 5.7 Lyα LF agreement, shows that the evolution of the z = 5.7 to z = 7 Lyα LF is explainable with smaller halos and lower intrinsic Lyα emission, rather than an increase in global neutral fraction in our simulation.

Figure 6 .
Figure 6.The black points show a LAE LF measured from a 2 square degree survey with perfect observations to a depth of fLyα > 1 × 10 −17 [ergs s −1 cm −2 ] and a redshift window of 6.75 < z < 7.25.The model LFs for xHI = [0.01,0.3, 0.5] are shown as the blue, black, and green curves, respectively.The error bars are assumed to be Poissonian for N > 20, otherwise they are tabulated from Gehrels (1986).The intrinsic scatter of the black points about the 0.3 model LF, the input neutral fraction, comes from cosmic variance.

Figure 7 .
Figure 7.The posterior on xHI inferred from the LAE LF measured in Figure 6 is shown in blue.The simulation's input xHI, 0.3, is indicated with the vertical dashed black line.The posterior favors a slightly lower xHI than the true value because the measured region of space happened to have a slightly over-dense number of LAEs.The large width of posterior indicates that xHI is not very well constrained; 95% of the posterior is contained in the range 0.01 < xHI < 0.45.

Figure 8 .
Figure 8.The left panel shows the distribution of the number of LAEs observed, NLAE, in 10,000 realizations of the nominal 2 square degree survey in a xHI=0.3 at z = 7 universe in blue.We investigate one over-dense region (orange vertical line) and one under-dense region (green vertical line) further.A Poisson distribution with the same sample mean is shown in black-the excess variance over the Poisson distribution is due to the clumpy nature of the large-scale-structure and the rarity of LAEs.The right panel shows the posteriors on xHI for the extrema of the 95th percentile range of LAE densities in orange (high LAE density) and green (low LAE density).The simulation's input xHI is shown with the vertical black dashed line.The difference between the medians of these distributions is ∆ ∼ 0.19.

Figure 9 .
Figure9.We increase the area of the nominal surveys on the previously considered density extrema locations while keeping the luminosity limit and redshift window constant.The posteriors in the low NLAE region flatten out at a non-zero value because of a single bright LAE which happened to fall within the survey.Increasing the area to 20 square degrees centers moves the peak of the posterior towards the true value of xHI, so we will further investigate 20 square degree surveys to determine the intrinsic spread in inferred xHI from cosmic variance.

Figure 10 .
Figure10.The distribution of 100 posterior medians for the 20 square degree survey (blue) and nominal 2 square degree survey (black) for MCMC analyses.Both of these surveys have ∆z = 0.5 and luminosity limit LLyα > 6 × 10 42 [ergs s −1 ].The area of the 20 square degree survey was chosen because it seemed to remove the offset in inferred xHI with respect to the true value seen in Figure9.The scatter in the posteriors medians is consistent between the two survey strategies, indicating that even increasing the area of the survey by a factor of 10 is not enough to overcome the fact that the uncertainty on xHI is dominated by cosmic variance and imposes an uncertainty ∆xHI ∼ 0.2.

Figure 11 .
Figure 11.The distribution of 100 posterior medians for the deep (LLyα > 1.8 × 10 42 [ergs s −1 ]) 10 square degree surveys (blue) and the posterior medians of the 100 nominal surveys (black).For both survey strategies ∆z = 0.5.The widths of the distributions of both survey strategies is about the same; 95% of the nominal survey's medians are within a range ∆xHI= 0.25, while the deep 10 square degree medians have 95% of their realizations within ∆xHI= 0.23.This indicates that pushing to deeper flux limits is not a viable strategy in mitigating the effect of cosmic variance.

Figure 12 .
Figure12.The left panel shows the distribution of the number of LAEs per square degree observed by 1000 surveys of areas 2 and 20 square degrees (and to depth LLyα 6 × 10 42 ergs s −1 ) in blue and red, respectively.Surveys for each area, centered on the under-dense region identified in the nominal 2 square degree survey, are indicated with the dashed vertical lines of the same colors.The distributions of NLAE per square degree tend to narrow as survey area increases.The right panel shows the distribution of the NLAE observed in 1000 realizations of the 20 square degree.The underlying green distribution shows the distribution of NLAE in 300 surveys of 100 0.2 square degree fields, totalling 20 square degrees.This survey strategy will be discussed more in Section 3.4.3;note that the distribution of NLAE is narrower for this strategy.The dashed red vertical line again indicates the under-dense region selected with the 2 square degree survey, and the distribution's 2.5 percentile is indicated with the solid black vertical line.The imprint of the under-dense region identified with a 2 square degree survey is still dominating the 20 square degree survey's total number of LAEs compared to other 20 square degree surveys, illustrated by the fact that the red dashed line does not move much towards the distribution's mean.

Figure 13 .
Figure 13.The distribution of posterior medians for 90 surveys composed of 100 0.2 square degree fields to a depth of LLyα 1.8 × 10 42 ergs s −1 is shown in blue.The distribution of the 100 nominal survey's medians is shown in black.95% of the nominal survey's medians are within a range ∆xHI= 0.25, while the 90 surveys composed of 100 independent fields have 95% of their realizations within ∆xHI= 0.05.This strategy is effective in reducing the systematic uncertainty on xHI.