The Universe is at Most 88% Neutral at z=10.6

Recent observations of GN-z11 with JWST have revealed a Ly$\alpha$ emission line with an equivalent width of 18$\pm 2$ angstroms. At z=10.6, this galaxy is expected to lie in the heart of reionization. We use a series of inhomogeneous reionization simulations to derive the distribution of the Ly$\alpha$ EW after traveling through the neutral intergalactic medium with varying average neutral gas fraction, $x_{HI}$. We use these distribution to place an upper limit of $x_{HI}<$ 0.88 at z=10.6 at 95% confidence level. We compare our upper limit to different reionization history models, which include the recently identified enhancement at the bright end of the luminosity function at z>8. We find that models in which faint galaxies have higher escape fraction compared to bright galaxies are favored by the new data.


INTRODUCTION
In less than a year, JW ST has revolutionized our knowledge of the very-high-redshift universe, spectroscopically confirming and photometrically identifying galaxy candidates as early as only a few hundred million years after the Big Bang (Schaerer et al. 2022;Brinchmann 2022;Finkelstein et al. 2022;Curtis-Lake et al. 2023;Cameron et al. 2023;Tang et al. 2023). Recently, Bunker et al. (2023) presented the spectroscopic confirmation of the Lyman-break galaxy candidate GN-z11. Using multiple rest-frame optical and UV emission lines, the redshift of GN-z11 was measured to be z = 10.603 ± 0.0013, currently the highest known from emission lines for a likely starforming galaxy (see discussion in Bunker et al. (2023)).
At redshift well above 10, GN-z11 is expected to be embedded in a largely neutral universe, as inferred from Corresponding author: Sean Bruton bruto012@umn.edu a number of available data at z 8 (e.g., Fan et al. 2006;McGreer et al. 2011McGreer et al. , 2015Ono et al. 2012;Schroeder et al. 2013;Choudhury et al. 2015;Greig et al. 2016;Mason et al. 2018aMason et al. , 2019Greig et al. 2019;Hoag et al. 2019;Jung et al. 2020). Inferring from lower-redshift observational studies suggests that the universe volume weighted neutral fraction (x HI ) is well above 90% at z 10, although no observational constraints currently exist at such early epochs. The spectrum of GN-z11 offers the unique opportunity to understand the state of the diffuse hydrogen at z > 10. Indeed, the spectrum reveals the presence of a strong Lyα emission line, with a rest-frame equivalent width (EW) of 18 ± 2Å, and a velocity shift of the Lyα with respect to systemic velocity of ∆v Lyα ≈ 550km s −1 . This implies that GN-z11 resides in a an ionized bubble, perhaps not surprising given that it is one of the brightest galaxies known at z > 10.
Lyα emission has been used extensively to place constraints on reionization. Lyα is a resonant transition in neutral hydrogen and is expected to completely disappear from the spectra of high-redshift galaxies when these are embedded in a fully neutral intergalactic medium (IGM; Dijkstra 2014Dijkstra , 2017. The evolution of the absolute number density of Lyα emitters has been used to constrain the timeline of reionization (e.g., Rhoads & Malhotra 2001;Malhotra & Rhoads 2004Hu et al. 2019;Wold et al. 2021;Morales et al. 2021;Ning et al. 2022). However, because Lyα is the byproduct of star-formation, its disappearance cannot be uniquely interpreted as signaling a completely neutral IGM. Studies of the reionization history of the universe are now often based on relative properties, e.g., on the relative fraction of Lyα emitters among star-forming galaxies (Stark et al. 2010;Pentericci et al. 2011;Jung et al. 2018), on the Lyα equivalent width distribution (EW Lyα ; Mason et al. 2018b;Hoag et al. 2019;Whitler et al. 2020), or on the clustering of Lyα emitters Sobacchi & Mesinger 2015;Ouchi et al. 2018;Yoshioka et al. 2022). Additionally, the environment in which a galaxy resides has an impact on the visibility of Lyα. A galaxy located within a region of the universe previously ionized by either itself, or a previous generation of stars, i.e., residing in what are commonly referred to as "ionized bubbles," would be more likely to show Lyα in emission because Lyα photons would have time to redshift out of resonance before encountering the neutral IGM. The probability of observing Lyα emission then depends on the size of the ionized bubbles and the Lyα velocity shifts with respect to the IGM.
In this paper, we use the observations of GN-z11 to constrain x HI at z = 10.6. Specifically, we use a set of simulations to forward model the evolution of the Lyα EW distribution in the presence of an increasingly neutral IGM and use these distributions to construct the likelihood of observing a galaxy with the observed EW Lyα and M UV properties as GN-z11.
The structure of the paper is as follows. In Section 2, we postprocess 21cmFast simulations to create a model of Lyman-α emitters (LAEs). We then present our inference model. In Section 3, we explore the implication of our result for the reionization history and discuss how the result depends on the assumptions. Finally we conclude in Section 4. We assume a ΛCDM cosmology with H 0 = 67.66 km s −1 Mpc −1 (Planck-Collaboration et al. 2020).

MODELING
The end goal of our modeling is to make an inference on the global neutral fraction of the universe, x HI , at z = 10.6, in light of the fact that we observe one LAE with EW Lyα = 18 ± 2Å. To estimate x HI , we will use a hierarchical Bayesian approach, similar to the one used in Bruton et al. (2023) and Mason et al. (2018b). We will populate a series of large inhomogenous reionization simulations with LAEs and calculate the distribution of the EW Lyα after transmission through the IGM of different ionization fractions.

Populating a Reionization Simulation with LAEs
Our model is built on the backbone of 21cmFastv2 (Mesinger & Furlanetto 2007;Mesinger et al. 2011Mesinger et al. , 2016, a seminumerical code that combines excursion set formalism and perturbation theory to create simulations of the universe throughout reionization. We use nine simulations with a volume average neutral fraction varying between x HI =0.01 to x HI =0.92, each 4.1Gpc 3 in volume. When modeling the evolution of x HI , 21cmFast takes into account recombinations, photoheating star-formation suppression, supernova feedback, and radiation. For each dark matter halo of mass, M h , 21cmFast computes the integrated Lyα optical depth as a function of the velocity offset from line center. The simulations naturally account for the evolving presence of ionized bubbles around galaxies of different masses and residing in different environments. In what follows, we define the emergent Lyα (Lyα emer ), as the Lyα line luminosity escaping from the interstellar and circumgalactic medium of a galaxy, while the Lyα observed after passing through the IGM is referred to as transmitted Lyα (Lyα tran ). To assign an emergent Lyα to a dark matter halo, we follow Bruton et al. (2023). Briefly, after using the M h -M UV relation from Mason et al. (2015) to assign an M UV magnitude to each halo, we use the EW Lyα -M UV probability distribution from Mason et al. (2018a), calibrated with De Barros et al. (2017) observations, to randomly draw the emergent EW Lyα for each halo. This is discussed in more detail in the following paragraph. The emergent Lyα line profile shape is assumed to be a truncated Gaussian, with velocity offset and Full-Width-Half-Maximum (FWHM) dependent on Lyα luminosity. The velocity offset scales linearly with the log of the Lyα luminosity (see Bruton et al. (2023) for more details), and we set FWHM equal to velocity offset, which is consistent with the findings in Verhamme et al. (2018), wherein they fit the relation between observed low-and high-z LAEs' velocity offsets and FWHMs. They find the empirical relation to be consistent with the one-toone relation predicted from radiation transfer modeling. We also note that Hayes & Scarlata (2023) find, through a Bayesian hierarchical inference model, that GN-z11's intrinsic red wing Lyα emission has a velocity offset of 400 km s −1 and a FWHM of 433 km s −1 , consistent with the one-to-one. Finally, we use the line profiles, together with the velocity dependent optical depths of Lyα, to calculate the Lyα transmission, T = Lyαtran Lyαemer , and the transmitted EW Lyα as EW tran Lyα = T × EW emer Lyα . Our analysis rests on the assumption that the distribution of emergent EW Lyα at M UV = −21.5 does not change with redshift. For a generic M UV , this distribution takes the form, from Mason et al. (2018a), is the Heaviside step function and δ(EW Lyα ) is the Dirac delta function. The parameter A accounts for the fraction of galaxies that do not emit Lyα and W c determines the exponential decline of the probability distribution function toward larger EW Lyα . Both parameters are functions of M UV ; the general behavior is that brighter galaxies are more likely to be nonemitters and have a stronger exponential cutoff so that UV bright galaxies are less likely to have large EW Lyα . In Figure 1 we show the MUSE Hubble Ultra Deep Field (HUDF) Survey measurements of the average emergent EW Lyα at different redshifts as a function of M UV (Hashimoto et al. 2017), z ∼ 6 observations from De Barros et al. (2017), and the prediction of our model. We apply an observed flux lower limit > 8 × 10 −19 erg s −1 for Lyα to match the selection limit in the MUSE survey (the flux limit in the De Barros et al. (2017) data is similar at 2.2 × 10 −18 erg s −1 across the entire wavelength range, deeper in wavelength ranges without skylines).
We limit to data to z 6, to limit the impact of the transmission through the IGM and isolate the effects of the interstellar and circumgalactic medium. This permits a comparison between EW Lyα distributions at different redshifts to see if emergent EW Lyα is independent of redshift. The comparison in Figure 1 confirms that there is no evidence of evolution in the median emergent EW Lyα from z = 3.6 to z = 4.9, and weak evidence of evolution from z = 4.9 to z = 6.0, a redshift range corresponding to 800 million years. We note that a partially neutral IGM at z = 6 would remove LAEs with weak EW Lyα from the observations and thus bias the distribution to be higher, and that the universe may not be fully ionized by z = 6 (Qin et al. 2021).
Still, we test the impact of the intrinsic EW Lyα distribution by adopting a different exponential cutoff strength, W c , for our intrinsic EW Lyα distribution and redo our inference. For M UV = −21.5, W c = 19, as fit from the De Barros et al. (2017) data. As a test, we weaken the exponential cutoff, setting W c = 43, which, allowing galaxies to take on EW Lyα values greater than the GN-z11 value of 18Å more easily, decreases the probability of having a galaxy with an EW Lyα of 18Å by ≈20%. When redoing the inference with this new distribution, we infer the same 95% upper limit on x HI . This is an extremely different intrinsic EW Lyα distribution, but the impact on the inferred value of x HI is negligible. This shows that even if there is strong evolution in the intrinsic EW Lyα distribution toward z = 10.6 allowing galaxies to have larger EW Lyα , its impact on our result will be small.
The Lyα transmission through an inhomogenous IGM, particularly in the case of a highly neutral IGM, changes the EW Lyα distribution, which is no longer well described by Equation 1. This is demonstrated in Figure  2, where we show p(EW tran Lyα | M UV = −21.5, x HI ) for two values of x HI . For the highly neutral universe with x HI = 0.9, the probability distribution function (PDF) of the transmitted EW Lyα can no longer be modeled with the first term in Eq. 1. The observed excess of high-EW LAEs compared to the prediction from the ionized universe is a direct result of the inhomogeneity of the reionization process. Some galaxies reside in ionized bubbles, and their Lyα emer is not very attenuated by the IGM. This results in a distribution of EW Lyα that does not follow an exponential decline. This also highlights an important aspect of our simulations-the LAEs reside in a variety of local neutral fractions, which may differ from the global neutral fraction. However, when these LAEs are taken in conglomerate, the probability of having a given EW Lyα is dependent on the global neutral fraction, not the local neutral fraction. Thus, Figure 2. The probability of EW(Lyα) for varying values of xHI, given that the galaxy has Lyα emission (i.e. excluding nonemitters). The solid lines are empirically taken from the simulation after IGM attenuation-the dashed lines are an attempt to fit the first term of Eq. 1 to the distribution. The fit fails when the universe is partially neutral because there is excess probability of having high EWLyα compared to an exponential cutoff.
the inference arising from this probability yields the global value, rather than a local value. Seeing as we have only one object to feed into our inference, we expect the posterior to be weakly constrained; more objects, probing more lines of sight, would be much more constraining. ) is the PDF of the transmitted EW computed in the previous section for M UV = −21.5. Lyα optical depths vary smoothly with x HI , so we interpolate between our nine simulations to build a function for the EW Lyα distributions dependent on x HI and M UV . We assume that the likelihood of EW obs Lyα is a normal, with known standard deviation that we set equal to the measurement uncertainty provided by Bunker et al. (2023).

Hierarchical Bayesian
We use a Metropolis-Hastings algorithm to sample the posterior, accounting for measurement errors on EW Lyα .

RESULTS AND DISCUSSION
The main result of our analysis is presented in Figure 3, where we show the x HI posterior in the left panel. Not surprisingly, with only one galaxy, the value of x HI is poorly constrained. However, it is clear that x HI is smaller than 0.88 at 95% confidence level (the 95% credibility region is indicated by the vertical line in Figure 3). In the right panel of Figure 3, we show the constraint on x HI derived from GN-z11 on the reionization history timeline, and compare it with constraints on the neutral fraction from other observations: Lyα EW of galaxies (Mason et al. 2018a(Mason et al. , 2019Hoag et al. 2019); the clustering of Lyα emitter galaxies Greig et al. 2016); Lyα and Lyβ dark fraction (McGreer et al. 2015); QSO damping wings (Davies et al. 2018). While the GN-z11 observation does not constrain the absolute value on x HI (one can just as easily say that x HI is > 0.04 with 95%), the fact that the observation is very far back in the expected timeline of reionization offers a lever, allowing one to distinguish between reionization models.
Predicting the redshift evolution of the volume averaged neutral fraction depends on the number of ionizing sources present at each time, their ionizing spectrum, and the ability of ionizing radiation to escape into the IGM. Once these are known, x H i can be computed by solving: (Madau et al. 1999;Robertson et al. 2013;Ishigaki et al. 2018) whereṄ ion is the ionizing photon production rate, n H is the comoving gas number density, and t rec is the recombination time scale. n H and t rec are defined as where X p , Y p are the primordial mass fraction of hydrogen and helium, Ω b is the baryon energy density fraction, ρ c is the critical density, (Ω b ρ c = 4.2 × 10 −31 [g/cm 3 ]), C H ii ≡ n H ii 2 / n H ii 2 is the clumping factor, and α B (T ) is the case B recombination coefficient. Here we assume C H ii = 3, and α B = 2.6 × 10 −13 cm 3 s −1 , for an electron temperature of 10 4 K. We solve Equation 2 iteratively, assuming the boundary condition that x H i = 1.0 at z = 18, i.e., we assume that the first sources of ionizing photons appear at this redshift.  (Bouwens et al. 2021) and double power-law LF (Harikane et al. 2023). The red and blue lines are the models where bright and faint galaxies have higher escape fraction, respectively. The observational constraints are shown in white markers, and the constraint in constraint inferred by GN-z11 is marked as yellow star.
The major uncertainties in the calculation of x HI above are introduced in the ionizing photon production rate,Ṅ ion , that can be expressed as the product of three components:Ṅ where f esc is the absolute LyC escape fraction, ξ ion is the ionizing photon production efficiency, and ρ U V is the UV luminosity density, i.e., the integral of the galaxy UV luminosity function (LF).
The UV LF has changed very quickly after JWST observations. Early JW ST results indicate that bright galaxies at z > 8 are more numerous than previously expected (Finkelstein et al. 2022), and a double powerlaw LF (as opposed of a classical Schechter LF) was proposed to account for the newly discovered population of bright galaxies (Bouwens et al. 2021;Harikane et al. 2023). Are these new bright sources responsible for the observed neutral fraction at z = 10.6? The relative contribution of bright and faint galaxies to the reionization budget is unconstrained. Which objects prevail has an impact on the reionization timeline: reionization starts late but completes rapidly if bright galaxies dominate, (Sharma et al. 2016;Naidu et al. 2020), while it starts early but proceeds slowly if faint galaxies played the major role (Finkelstein et al. 2019).
In what follows we consider four simple reionization history models that include the recent developments in the LF studies and vary the contribution of bright and faint galaxies to the ionizing budget. For all models, we assume log(ξ ion /[Hz erg −1 ])=25.7, characteristic of high-redshift star-forming galaxies (Ning et al. 2023;Tang et al. 2023).
First, we consider the reionization models that assume different luminosity functions at z > 8: a Schechter LF (dashed line, Bouwens et al. 2021), and a double powerlaw LF (solid line; Harikane et al. 2023). At z ≤ 8, both models are calculated with the Schechter LF in Bouwens et al. (2021). We adopt a constant f esc =0.12 at all redshifts in both models. We find that although there are more UV bright (M U V < −22) galaxies when the double power-law LF is used, the early (z > 9) reionization histories are roughly equivalent in the two models. Neither model, however, is able to reionize the universe early enough to account for the x HI upper limit at z = 10.6, as in both models x HI reaches 90% only by z 90 ≈ 9.4. Next, we consider the models where bright (M UV < −18) or faint galaxies (M UV ≥ −18) dominate the ionizing photon budget. To achieve this goal we vary the escape fraction as a function of both M UV and redshift: with flattening at f esc =0.5. We refer to these models as bright and faint models. In both models, we adopt the DP LF from Harikane et al. (2023) at z > 8 and the Schechter LF from Bouwens et al. (2021) at z ≤ 8.
The bright model (red curve in Figure 3) shows a late and rapid reionization similar to the models with constant f esc , despite the fact that the dominant contributors to the ionizing photon budget are different (in the constant f esc models, faint galaxies dominate because of their high space density compared to bright galaxies). In the faint model (blue curve in Figure 3), the ionizing photon contribution of galaxies with M UV > −18 is enhanced with respect to the constant f esc models, causing reionization to start earlier. Specifically, we find that x H i reaches 90% by z 90 ≈ 10.9 and progresses more slowly compared to models where brighter galaxies dominate. The faint model is the preferred scenario for the early stage of reionization (z > 10) in light of the x HI < 0.88 constraint inferred from GN-z11.
An independent constraint on the reionization history comes from observations of the cosmic microwave background. Recently, Planck (Planck-Collaboration et al. 2020) measured the (integrated) optical depth to Thomson scattering to be τ = 0.056 ± 0.007, ruling out the need of a large contribution from galaxies at z > 10, and favoring a fast and late reionization process. The redshift evolution of the free electron fraction, x e (z), constrained by Planck, however, depends on the specific form used to model the reionization history, as τ is only sensitive to the integral of the electron density over time, and not its detailed shape. The FlexKnot model assumed in PlanckCollaboration VI limits the contribution to the optical depth of z > 15 galaxies to less than 1%, but leaves space to the possibility that the universe was not fully neutral at z = 10.6, with a free electron fraction of ≈ 8−10%, consistent with the upper limit presented here.

CONCLUSIONS
At a record distance of z = 10.6, GN-z11 is the highest-redshift Lyα-emitting galaxy known, well within the heart of the reionization epoch. We use an inhomogeneous reionization simulation to derive the probability distribution of the transmitted EW Lyα through the IGM as a function of M UV and the average neutral gas fraction, x HI . We use these distributions to estimate the posterior distribution function on x HI at z = 10.6. With data for only one galaxy, we place an upper limit on the global neutral fraction at z = 10.6, i.e., x HI 0.88 with a probability of 95%. With this constraint we are able to exclude reionization histories dominated by bright galaxies; a scenario wherein faint galaxies have a higher escape fraction of ionizing photons, and so drive reionization, is favored.