Chasing the beginning of reionization in the JWST era
Christopher Cain Garett Lopez Anson D'Aloisio Julian B. Munoz Rolf A. Jansen Rogier A. Windhorst Nakul Gangolli





Christopher Cain
clcain3@asu.edu

0000-0001-9420-7384]Christopher Cain
School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-6004, USA


Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA



Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA

0000-0002-8984-0465]Julian B. Muñoz
Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA

0000-0003-1268-5230]Rolf A. Jansen
School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-6004, USA

0000-0001-8156-6281]Rogier A. Windhorst 
School of Earth and Space Exploration, Arizona State University,
Tempe, AZ 85287-6004, USA


Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA


§ ABSTRACT 

Recent JWST observations at z > 6 may imply galactic ionizing photon production in excess of prior expectations.  Under observationally motivated assumptions about escape fractions, these suggest a z ∼ 8-9 end to reionization, in strong tension with the z <  6 end required by the Lyα forest.  In this work, we use radiative transfer simulations to understand what different observations tell us about when reionization ended and when it started.  We consider a model that ends too early (at z ≈ 8) alongside two more realistic scenarios that end late at z ≈ 5: one that starts late (z ∼ 9) and another that starts early (z ∼ 13).  We find that the latter requires up to an order-of-magnitude evolution in galaxy ionizing properties at 6 < z < 12, perhaps in tension with recent measurements of ξ_ ion by JWST, which indicate little evolution.  We also study how these models compare to recent measurements of the Lyα forest opacity, mean free path, IGM thermal history, visibility of z > 8 Lyα emitters, and the patchy kSZ signal from the CMB.  We find that neither of the late-ending scenarios is conclusively disfavored by any single data set. However, a majority of these observables, spanning several distinct types of observations, prefer a late start.  Not all probes agree with this conclusion, hinting at a possible lack of concordance between observables.  Observations by multiple experiments (including JWST, Roman, and CMB-S4) in the coming years will either establish a concordance picture of reionization's early stages or reveal systematics in data and/or theoretical modeling.  









§ INTRODUCTION
 

Despite an explosion of new data in the past decade probing cosmic reionization, little is known about how and when the process began.  The ending of reionization, believed to occur at 5 < z < 6, has been constrained by observations of the Lyα forest of high-redshift QSOs <cit.>.  Direct (and indirect) measurements of the mean free path from QSO spectra <cit.>, and measurements of the IGM thermal history at z ≥ 5.5 <cit.> have further corroborated this picture.  The electron scattering optical depth to the CMB <cit.> constrains the midpoint of reionization to be z ∼ 7.5.  Evidence of damping wings in high-redshift z ∼ 7-8 QSOs <cit.> and measurements of the neutral fraction based on (non)detections of Lyα emitters <cit.> indicate that the IGM was partially neutral at these redshifts.  

Understanding the timeline of reionization is crucial for revealing the nature of the ionizing sources that drove it.  Explaining how reionization could be driven by galaxies alone has been historically challenging thanks to the high early measurements of τ_ es <cit.> which required very early and/or extended reionization histories.  This tension eased as the measured value of τ_ es steadily decreased.  Several recent works <cit.> have showed that galaxies can complete reionization by z ≈ 6 under physically reasonable assumptions about their ionizing properties.  Recent evidence for an end later than z = 6 further relaxed demands on galaxy ionizing output (although see ).  Concurrent efforts demonstrated that AGN are unlikely to have contributed the majority of the ionizing budget responsible for reionization (e.g. , although see ).  

These findings paint a relatively simple, consistent picture of reionization: it ended at 5 < z < 6, was in progress at z ∼ 7-8, and was likely driven by galaxies.  Prior to JWST, the precise timing of reionization's midpoint and especially its early stages were not tightly constrained.  Constraints on reionization's midpoint from τ_ es <cit.> spanned a range of ± 0.75 in redshift (at 1σ), and few direct constraints on the first half of reionization existed.  The space of models proposed by the aforementioned works span a wide range of possibilities[The scenarios proposed by  (early-starting, gradually ending reionization), and  (late-starting, rapidly ending) reionization, roughly bracket the proposed possibilities.  ] without contradicting observations.  Indeed, one important goal for JWST is to probe the properties of galaxies at z > 7-8 in hopes of learning more about reionization's early stages. 

However, the first JWST results may be complicating, as much as clarifying, our understanding of reionization.  JWST has allowed for measurements of the UV luminosity function (UVLF) at much higher redshifts than HST, up to z ∼ 14 <cit.>.  It has also allowed us to measure the ionizing efficiency of galaxies, ξ_ ion, above z = 6 <cit.>.  Recently,  pointed out that a face-value interpretation of recent UVLF and ξ_ ion measurements from JWST () combined with observationally motivated assumptions about ionizing escape fractions (f_ esc) suggests reionization ended around z ∼ 8-9, inconsistent at > 2σ with the Planck τ_ es measurement[Although the tension is slightly smaller with the recent re-measurement of τ_ es from Planck data by  - see Figure <ref>.  ].  This result represents a stark reversal from the historical problem of galaxies producing too few ionizing photons to complete reionization on time <cit.>.  Several bright Lyα emitters (LAEs) at z > 8 <cit.> have also been observed, with the highest-redshift detection to date at z = 10.6 <cit.> by JWST.  These observations may be surprising if the IGM is mostly neutral at these redshifts, since observing Lyα emission requires some level of ionization around galaxies <cit.>.  

The top panel of Figure <ref> shows the volume-averaged ionized fraction (x_ HII^ V) for the three reionization models we will study in this work.  The early start/early end model is motivated by the aforementioned findings of , and has a midpoint (endpoint) of z_ mid = 8.5 (z_ end = 8).  The late start/late end model is motivated by Lyα forest observations at 5 < z < 6 and the Planck τ_ es measurement, has z_ mid = 6.5 and z_ end = 5.  The early start/late end model also has z_ end = 5, but an earlier midpoint (z_ mid = 7.5) and a start at z ∼ 13.  These three models broadly represent the three types, or categories, of possible reionization histories.  The bottom panel shows τ_ es for each, compared with the  measurement (black) and the recent re-analysis of Planck data from , with shaded regions denoting ± 1σ uncertainties.  The late start/late end model is consistent within 1 σ with , and the early start/late end case is similarly consistent with .  



In this work, we study the observational properties of the three reionization models shown in Figure <ref> using radiative transfer (RT) simulations.  We will demonstrate, in accord with previous works, that observations from the z ≲ 6 Lyα forest strongly disfavor the early start/early end model, in agreement with τ_ es.  We will also compare the two late-ending models to a wide range of observations with the goal of understanding whether existing data supports a late or an early start to reionization.  This work is organized as follows. <ref> describes how we calibrate our three reionization models and discusses their implications for galaxy properties in light of recent JWST observations.  In <ref>, we describe our methods for running RT simulations of reionization and forward-modeling various observables.  We compare our models to several sets of complementary observations in <ref>, discuss the implications of our findings in <ref>, and conclude in <ref>.  After discussing an observable, we will bold-face the name of the reionization model (if any) preferred by that observable.  Throughout, we assume the following cosmological parameters: Ω_m = 0.305, Ω_Λ = 1 - Ω_m, Ω_b = 0.048, h = 0.68, n_s = 0.9667 and σ_8 = 0.82, consistent with  results. All distances are in co-moving units unless otherwise specified.  



§ IMPLICATIONS OF JWST GALAXY OBSERVATIONS
 



 §.§ Reionization Model Calibration



The main input to our RT simulations (described in <ref>) is the globally averaged ionizing photon emissivity of sources verses redshift, Ṅ_γ(z).  In this section we describe how we use a combination of JWST observations and Lyα forest data to construct (or “calibrate”) Ṅ_γ(z) for our three models.  

Our starting point is the amount of non-ionizing UV light produced by galaxies, which is quantified by the UV luminosity function (UVLF).  This has been measured up to z ∼ 14 by JWST using both photometry and spectroscopy <cit.>.  The top two panels of Figure <ref> show two sets of UVLFs that we use in our analysis.  In both panels, the dashed curves show the z < 8 UVLFs measured by  with HST. In the left panel, the solid curves denote the double-power-law (DPL) fits from  at 8 ≤ z ≤ 12, and the dotted curve is the measurement from  at z = 14.5.  In the right panel, the dotted curves show results for the redshift-dependent UVLF parameters given in Eq. 3-6 of  (which is a best-fit to measurements from ).  At z > 8, the UVLFs in the left panel are a factor of ∼ 3 lower than those in the right panel.  We will use these sets to roughly bracket observational uncertainties on the UVLF.  

The lower left panel shows the integrated UV luminosity density, ρ_ UV, vs. redshift.  The curves show the logarithmic average of ρ_ UV calculated using the two sets of UVLFs in the top panels, and the shaded regions show the spread between them.  For illustration, we integrate the UVLF down to two limiting magnitudes - a bright cutoff of M_ UV^ cut = -17 (black) and a fainter M_ UV^ cut = -13 (magenta).  The background shading denotes the redshift ranges covered by HST and JWST data.  Note that at z > 8, the redshift evolution of ρ_ UV(z) is fairly insensitive to M_ UV^ cut, the main difference being normalization.  This is because both sets of UVLFs have faint end slopes close to α = -2.1 at z ≥ 8, with little evolution across that redshift range <cit.>.  This is in contrast to some pre-JWST expectations, which predicted a steepening of the faint-end slope with redshift and a corresponding shallow evolution in ρ_ UV for faint M_ UV^ cut <cit.>.  Instead, ρ_ UV evolves quickly with redshift, with a factor of ∼ 10 evolution between z = 8 and 13.  



The lower right panel quantifies the redshift evolution of galaxy ionizing properties in our models using

    ⟨ f_ escξ_ ion⟩_L_ UV≡Ṅ_γ(z)/ρ_ UV(z)

This is the UV-luminosity (L_ UV)-weighted average of the product of f_ esc and ξ_ ion for the galaxy population.  The black, red, and blue curves show this quantity for our three reionization models (see legend) assuming the average ρ_ UV curve for M_ UV^ cut = -13 (lower left panel).  The shaded regions show the spread in this quantity (for fixed Ṅ_γ(z)) arising from observational uncertainty in ρ_ UV(z), as shown in the lower left.  We calibrate Ṅ_γ(z) for the early start/early end model such that its ⟨ f_ escξ_ ion⟩_L_ UV agrees with the observationally motivated model from  (which assumes M_ UV^ cut = -13), shown by the faded gray dashed curve.  Consistent with the findings of , this model ends reionization early at z ≈ 8 (see  Figure <ref>).  

Our late start/late end scenario is motivated by the possibility that the tension with τ_ es and the Lyα forest in the early start/early end model can be solved with a simple redshift-independent re-scaling of Ṅ_γ(z).  This could be the case if the measurements of ξ_ ion assumed in  (from ) are systematically biased high, and/or if the same is true of the f_ esc values they inferred from the results of .  We find that scaling Ṅ_γ(z) down by a factor of 5 from the early start/early end case brings the end of reionization to z ∼ 5.  We then further adjust Ṅ_γ(z) at z < 7 at the few-percent level until we achieve good agreement with the Lyα forest at z ≤ 6 (as we will show in <ref>).  The cyan-dashed curve shows the observationally motivated model from , multiplied 3.5, which also has redshift evolution similar to this scenario.   

However, the late start/late end scenario is not the only possibility allowed by τ_ es and the Lyα forest.  It could also be that Ṅ_γ(z) is consistent with the  model at the highest redshifts (z ≳ 10), but declines at lower redshifts in such a way that reionization ends at z < 6, as required by the Lyα forest.  This scenario is represented by our early start/late end model (red dashed curve).  This model also ends reionization at z = 5, and by adjusting its Ṅ_γ(z) at z < 10, we can also achieve agreement with the Lyα forest.  This model requires a factor of ≈ 10 decline in ⟨ f_ escξ_ ion⟩_L_ UV between z = 10 and 6, a decrease much steeper than suggested by the results of of  and .  This could be achieved if ξ_ ion, f_ esc, M_ UV^ cut, or some combination of these evolves significantly across this redshift range (see <ref>).  



 §.§ Is the early start/late end model plausible?



Given the discrepancy between the evolution of ⟨ f_ escξ_ ion⟩_L_ UV in the early start/late end model and in the  and  models, it is natural to ask whether this scenario is plausible.  The top panel of Figure <ref> shows ⟨ f_ escξ_ ion⟩_L_ UV for this model, alongside observationally and theoretically motivated scenarios that make various assumptions about the redshift evolution of ξ_ ion and/or f_ esc.  The dotted curve is the  model scaled down by 0.2, which agrees with the late start/late end case.  The dashed curve is the same, except that we extrapolate ξ_ ion (Eq. 4 in ) outside the range of M_ UV and redshift within which it was fit to data.  For the dot-dashed curve, we further replace the observationally motivated f_ esc prescription assumed in  with the global f_ esc(z) from the flagship THESAN simulation <cit.>.  We have re-scaled each of the gray curves by different constants to bring them as close as possible to the early start/late end model.  We obtain the best agreement for the dot-dashed curve, which boosts the redshift evolution of both ξ_ ion and f_ esc relative to .  This shows that the early start/late end model is plausible given evolution in f_ esc and/or ξ_ ion, but only under “favorable” assumptions about both.    



The bottom panel of Figure <ref> illustrates another mechanism that could work in the direction of making reionization start earlier: evolution in M_ UV^ cut.  The red curve assumes that the M_ UV^ cut evolves from -10 at z = 12 to -19 at z = 6 (M_ UV^ cut(z) = -19 + 3/2(z - 6)), which causes ρ_ UV to decline much less rapidly with redshift than it does in the bottom left panel of Figure <ref>.  This allows ⟨ f_ escξ_ ion⟩_L_ UV in the early start/late end scenario to evolve less quickly, in better agreement with the  model (dotted curve).  This type of behavior in the galaxy population could arise from decreasing dust obscuration <cit.> and/or feedback from the IGM reducing or shutting off star formation in low-mass halos at lower redshifts <cit.>.  The evolution in M_ UV^ cut assumed in this illustrative example is extreme, and likely ruled out by existing observations <cit.>, but serves to show how an evolving M_ UV^ cut could support the early start/late end model.  

These comparisons suggest that a factor of a few of evolution in each of ξ_ ion and f_ esc, and perhaps some in M_ UV^ cut, could explain the early start/late end scenario.  Indeed, there is some observational and theoretical support for the idea that ξ_ ion could increase with redshift and/or be larger in fainter galaxies <cit.>, However, some works find that ξ_ ion is closer to constant with M_ UV, and perhaps redshift <cit.>.  Unfortunately, f_ esc at these redshifts cannot be directly measured, and there is no consensus in the literature from either indirect observational estimates <cit.> or cosmological simulations <cit.>.  Thus, galaxy observations alone cannot rule out the early start/late end case, although such a model does require more significant evolution in galaxy ionizing properties than suggested by current observations. As such, we conclude that the late start/late end model seems more likely based on existing data.  



 §.§ Ionizing photon budget





We show the Ṅ_γ(z) for our models in the top panel of Figure <ref>, alongside several others from the literature that also match the z < 6 Lyα forest.  The cyan-solid and magenta-dashed faded curves show the “fiducial” and “early” models from , and the dot-dashed gray curve shows the reference model[We actually show the reference model + “extra un-resolved sinks” (their Fig. 6) since it assumes the same sub-grid sinks opacity model we use in this work.  ] from .  Our late start/late end model drops off faster than the others at z > 9, and our early start/late end model closely matches the “early” model from .  The bottom panel of Figure <ref> shows the integrated ionizing photon output of galaxies per H atom, with the vertical lines denoting the end of reionization and the horizontal lines the number of photons per H atom needed to complete reionization (the so-called “ionizing photon budget”).  We find a budget of ≈ 2.6, 2.8, and 2.7 photons per H atom for the late start/late end, early start/late end, and early start/early end models, respectively. 

All our models are mildly “absorption-dominated” - meaning N_γ/ H is (slightly) more than twice the number of baryons per H atom in the universe (1.082, counting helium).  The photon budget depends mildly on the reionization history itself, but is determined mainly by the recombination rate predicted by our IGM opacity sub-grid model.  Recently,  made a first attempt to observationally measure the “clumping factor”, which quantifies the recombination rate in the in-homogeneous IGM <cit.>.  They found C ∼ 12 at z > 5, higher than the commonly assumed C = 3 <cit.>.  For a reionization history similar to our late start/late end model, they find this higher C increases the photon budget from ≈ 1.5 to ≈ 3.  Assuming that the number of recombinations, N_γ/ H - 1.082, is ∝ C, the photon budget in the late start/late end model implies an effective C ∼ 9.5, slightly below but consistent at 1σ with the  measurements[Caveats include (1) that the real effective clumping factor predicted by our model is likely not constant in redshift <cit.>, and (2) the integrated number of photons absorbed by the IGM is slightly smaller than the number produced, an effect not accounted for in .  ].  Using this crude approximation, a recombination rate high enough to delay the end of the early start/early end model to z ≈ 5 would require C ∼ 85, an unrealistically high number even compared to the most extreme simulations[For z_ end = 5.5 (6), this estimation gives C ∼ 70 (55).  ].  



§ NUMERICAL METHODS
 

In this section, we discuss some of the relevant technical details of our RT simulations, and our methods for forward-modeling the observables discussed in the next section.  The reader interested only in the results of this work may safely skip this section and pick up at <ref>.  



 §.§ Radiative Transfer Simulations



We ran RT simulations of reionization using FlexRT (Cain & D'Aloisio in prep.).  FlexRT is an adaptive ray tracing code that post-processes a time series of cosmological density fields to simulate the growth of ionized regions, the ionizing background, and the IGM thermal history.  The code uses a sub-grid model for opacity to ionizing photons in ionized gas, which captures the effects of ∼ kpc-scale structure that cannot be directly resolved.  Our sub-grid model is based on a suite of high-resolution, coupled hydro/RT simulations that can resolve the Jeans scale of cold, pre-ionized gas, run with the setup of <cit.>.  We refer the reader to <cit.>, and a forthcoming code paper (Cain & D'Aloisio in prep.) for further details about FlexRT.  

Ionizing sources for the RT calculation are halos taken from a dark matter (DM)-only N-body simulation run with the particle-particle-particle-mesh (P^3M) code of .  This run has a box size of L_ box = 200 h^-1Mpc and N = 3600^3 DM particles, and uses a spherical over-density halo finder to identify halos.  We find that our halo mass function matches the halo mass function of  to within ≈ 5-10% (HMF) for M_ halo > 3 ×  10^9 h^-1M_⊙, and is ≈ 50% incomplete at M_ halo = 1 ×  10^9 h^-1M_⊙.  To ensure we have a significant number of halos at the highest redshifts we simulate (z > 15), we adopt a mass cutoff for our sources of M_ halo > 1 ×  10^9 h^-1M_⊙, even though our HMF is incomplete there.  The density fields used in our RT calculation are taken from a high-resolution hydrodynamics simulation with the same large-scale initial conditions as the N-body run, which are re-binned to N_ RT = 200^3 for the RT.  Our simulations start at z = 18 and end at z = 4.8.  

We assign UV luminosities (L_ UV) to halos by abundance-matching to observed UV luminosity functions (UVLFs).  We use the measurements of  at z < 8, those of  at 8 < z < 14, and that of  at z = 14.5 - those shown in the upper left panel of Figure <ref>.  The ionizing emissivity is distributed between halos following ṅ_γ∝ L_ UV[The amounts to assuming that the product of the escape fraction and ionizing efficiency of galaxies, f_ escξ_ ion, is constant over the ionizing source population at a fixed redshift.   Note that although this is inconsistent with the model for f_ esc and ξ_ ion assumed in <ref>, the spatial distribution of photons between sources is both highly uncertain and is a higher-order effect for the purposes of most of this work.  We will comment whenever it becomes relevant in subsequent sections.  ].  The volume-averaged ionizing emissivity of sources, Ṅ_γ = 1/L_ box^3∑_i ∈ halosṅ_γ^i, is set by hand at each redshift and used to normalize the ṅ_γ proportionality.  As explained in <ref>, we calibrated Ṅ_γ(z) for each scenario using a combination of JWST data (<ref>) and measurements from the z < 6 Lyα forest (<ref>, <ref>).  We bin halos by the RT cell they occupy and treat cells containing one or more halos as sources.  We use a single ionizing frequency[See e.g.  for discussions of the effect of multi-frequency RT on the Lyα forest.  ] (E_γ = 19 eV), chosen to reproduce the same frequency-averaged HI cross-section, ⟨σ_ HI⟩, as a power law spectrum of the form J_ν∝ν^-1.5 between 1 and 4 Rydbergs.  Cells are assigned post ionization-front (I-front) temperatures using the flux-based method prescribed in .  The subsequent thermal history is calculated using their Eq. 6.  



 §.§ Modeling the Lyα forest



We model the Lyα forest in post-processing using the density fields from our aforementioned high-resolution hydrodynamics simulation, which has N = 2048^3 gas cells.  Ionized fractions, photo-ionization rates (Γ_ HI) and temperatures are mapped from the FlexRT simulations onto these density fields, and the residual neutral fraction in ionized regions is calculated assuming photo-ionization equilibrium and the case A recombination rate[Case A is appropriate for the under-dense gas that set the forest transmission at these redshifts.  ].  The native spatial resolution of our density field is Δ x_ cell = 97.6 h^-1kpc, too coarse to fully resolve the low-density voids that set the transmission at 5 < z < 6 <cit.>.  We apply the multiplicative correction prescribed in Appendix A of  to our residual neutral fractions to get an effective resolution of Δ x_ cell = 12.2 h^-1kpc.  The Lyα voigt profile is approximated using the analytic fit from .  

Since our gas temperatures are calculated on the coarse RT grid, we do not capture the temperature-density relation on scales smaller than the RT cell size.  As temperature (usually) correlates positively with density in the IGM, low (high)-density hydro cells embedded in larger RT cells will be assigned temperatures that are too high (low).  To correct for this, we assign a local temperature-density relation to each cell using the procedure described in  (see also their Appendix E), which uses the IGM temperature model of .  We find this lowers the mean transmission by ≈ 10-15%, since the correction cools the under-dense cells, which affect the mean transmission the most.  

We compute Lyα forest statistics by casting 4000 sightlines from random locations and in random directions of length 50 h^-1Mpc, for a total path length of 200 h^-1Gpc.  We calculate transmission statistics from z = 4.8 to z = 6 in Δ z = 0.2 increments.  Since our native resolution of 97.6 h^-1kpc is only 2-3× narrower than the typical width of Lyα line profile (≈ 15 km/s vs. 30-50 km/s), we do the integration over the line at a velocity resolution 4 × higher than that of the hydro sim.  We find this reduces the mean transmission by a few percent at most.   



 §.§ Lyα transmission around galaxies



We have also modeled Lyα transmission on the red side of line center around halos that could host Lyα emitting galaxies (LAEs). 
This allows us to assess the statistics of LAE visibility in our models at z > 7.  LAEs typically emit Lyα red-shifted to both the red and blue sides of line center <cit.>.  Although any emission on the blue side would be absorbed by even the ionized part of the IGM at these redshifts, attenuating the red side requires damping wing absorption from the fully neutral IGM <cit.>.  This makes LAEs a potentially powerful probe of the IGM neutral fraction.  However, interpreting observed red-side Lyα emission (or lack thereof) is complicated by uncertainties in modeling the intrinsic line profile of the LAEs themselves, and other factors such as surrounding inflows/outflows and proximate self-shielding systems <cit.>.   

In this work, we will avoid modeling the intrinsic LAE line profile, instead focusing on the IGM transmission, T_ IGM.  To calculate this, we trace 50 randomly oriented sightlines away from each halo and compute the Lyα transmission profiles along each sightline at ± 5 Å from line center (≈± 800 km/s).  We use the same N = 2048^3 high-resolution hydro simulation for this calculation as for the Lyα forest, and we calculate the ionization state of the ionized gas in the same way[Our results are largely insensitive to how the highly ionized IGM is modeled, however, since the red side transmission is set by the damping wing from the neutral IGM.  ].  We begin integrating the Lyα opacity 500 h^-1kpc (5 hydro cells) away from the location of the halo to avoid gas within the halo itself contributing to T_ IGM.  Gas velocities relative to the halo are computed by subtracting the halo velocity measured from the N-body simulation.  

The gas around massive halos can have sharp line-of-sight velocity gradients owing to inflows near the halo.  Sightlines pointing away from these halos see positive velocity gradients, which narrows the Lyα line in redshift space.  We find that the resulting sharp jumps in velocity between adjacent cells can produce artifacts in our transmission spectra.  To mitigate this, we linearly interpolate the gas velocities (and all other relevant quantities) onto a grid with 4 × higher resolution than that of the simulation when calculating T_ IGM (see  for a description of a similar procedure).  We find the transmission profiles to be well-converged at this resolution.  



 §.§ Modeling the Patchy kSZ signal



CMB photons can scatter off free electrons during reionization.  If the electrons are moving relative the the CMB rest frame, this results in a Doppler shift of the photons.  This can shift the blackbody spectrum and result in additional temperature anisotropies in the CMB <cit.>. This is known as the Sunyaev-Zel'dovich Effect, and its resultant temperature deviation along a line of sight is given by:


    Δ T/T = - σ_T n_e,0∫ e^-τ_ esγ̂·q⃗/cds/a^2


where γ̂ is the line of sight direction, q⃗ = (1 + δ) χv⃗ is the ionized momentum field, σ_T is the Thompson scattering cross section and τ_ es is the CMB optical depth to the last scattering surface, and n_e,0  is the mean electron density at z=0.  The integral can be broken up into a post-reionization, homogeneous kinetic Sunyaev-Zel'dovich Effect (hkSZ), and a high-z patchy kinetic Sunyaev-Zel'dovich Effect (pkSZ) while reionization is still occurring.  Untangling these components is difficult, so observations use templates to subtract the hkSZ component from the measured total kSZ power.  For the purposes of this work, we take any contributions from z ≳ 5 as part of the pkSZ, even in simulations where reionization ends at z > 5.  This keeps the different scenarios directly comparable to each other.  Note that actual measurements (such as the one from SPT by , see <ref>) must assume a fixed z_ end in their analysis.  

We calculate the signal using a method similar to that first suggested by .  The kSZ angular power spectrum can be calculated from the 3D power spectrum of q⃗ by: 


    C_ℓ = ( σ_T n_e,0/c)^2 ∫ds/s^2 a^4 e^-2τP_q_⊥(k = ℓ/s, s)/2

where (2π )^3 P_q_⊥(k) δ^D(k⃗-k⃗'⃗)=   ⟨q_⊥(k⃗) 
  ·q_⊥^* (k⃗'⃗)⟩ 
is the power spectrum of the specific ionized momentum modes perpendicular to the Fourier wave vector, q_⊥ = q - k̂ (q·k̂ ), δ^D(k-k') is the Dirac-δ function and
q = ∫q⃗ e^-i k⃗·q⃗ d^3 x⃗ denotes the Fourier transform of q⃗.  Only the q_⊥ mode contributes significantly to the kSZ signal, due to q_|| contributions canceling out when integrating over the line of sight.  We refer the reader to  and appendix A in  for details.  In general, it is common to write the angular power spectrum in the dimensionless form: D_ℓ = C_ℓ ( ℓ + 1 ) ℓ / (2 π).

Running RT simulations imposes practical limitations on the volume of our box. This limited box size means we miss large scale velocity modes that add to the kSZ power.  To compensate for this, we use the same correction applied in  eq. (B1), and take advantage of linear theory to generate a correction term:


    P_q_⊥^ miss(k,z) = ∫_k<k_ boxd^3 k'/(2π)^3 (1-μ^2) P_χ(1+δ)( | k⃗ - k⃗'⃗ | ) P_vv^ lin ( k')

where k_box = 2 π / L_box is the box-scale wavemode, μ = k̂·k̂'̂, P_χ(1+δ),χ(1+δ)(k) is the ionized matter power spectrum and P_vv^lin(k) = (ȧf/k) P_δδ^lin(k) is the velocity power spectrum in the linear approximation taken from the public code CAMB <cit.>.  

It was found by <cit.> that this kind of approach can still underestimate the D^pkSZ_3000 by ∼10-20% due to the irreducible or connected component of the ⟨δ_χ v ·δ_χv⟩ term being non-negligible because of the non-Gaussianity of reionization.  So our calculations could be seen as conservative estimates of the power.  As we will see in <ref>, an increase in power could strengthen our conclusions.  As such, we do not expect this missing term to affect our qualitative results.  



§ IMPLICATIONS OF OTHER REIONIZATION OBSERVABLES
 

In the rest of this work, we study three other observational windows into reionization: measurements from the spectra of high-redshift QSOs at z ≤ 6.5 (<ref>), observations of Lyα-emitting galaxies at z ≥ 8 (<ref>), and the patchy kSZ effect from reionization (<ref>).  Our goal will be to see if these observations, together with aforementioned JWST data, reveal a consistent picture about when reionization started and when it ended.  



 §.§ QSO Observations at 5 < z < 6
 



  §.§.§ The Lyα Forest
 



The 5 < z < 6 Lyα forest is perhaps the most compelling indicator that reionization ended at z < 6.  This conclusion has emerged from studies of the mean transmission of the Lyα forest and the scatter in Lyα opacities <cit.>.  We will study both in this section.  



Figure <ref> shows the mean transmission of the Lyα forest, ⟨ F_ Lyα⟩, at 5 ≤ z ≤ 6 in our models compared with recent measurements from .  The late start/late end and early start/late end models both agree well with the data - this is by design, since both Ṅ_γ histories were calibrated so that the simulations would match the  measurements of ⟨ F_ Lyα⟩.  As such, these two models cannot be distinguished using the mean forest transmission alone.  We see that the early start/early end model misses the measurements severely, predicting a mean transmission near unity at all redshifts.  Indeed, the average HI photo-ionization rate, Γ_ HI, in ionized gas is ≈ 2.7 × 10^-11 s^-1 at z = 6, ≈ 2 orders of magnitude higher than measurements at that redshift <cit.>.  Thus, the forest transmission measurements strongly disfavor a z ∼ 8 end to reionization.  Because of this, we will omit the early start/early end model from many of our subsequent comparisons, and focus instead on whether observations can distinguish the other two scenarios.  

Figure <ref> shows the cumulative distribution function (CDF) of forest effective optical depths measured over intervals of 50 h^-1Mpc, P(< τ_ eff^50).  The τ_ eff^50 distribution contains more information than ⟨ F_ Lyα⟩, since it is sensitive to the spatial fluctuations in the IGM ionization state.  This shape is sensitive to the IGM neutral fraction, since neutral islands produce high-τ_ eff^50 sightlines at z < 6 in late reionization scenarios <cit.>.  We show P(< τ_ eff^50) at z = 5, 5.2, 5.6, and 6 for each of our models, alongside measurements from .  To show what P(< τ_ eff^50) would look like in a hypothetical early-ending model that is compatible with the forest, we re-scaled Γ_ HI in the early start/early end[The re-scaling factor used here was ∼ 10^2 at all redshifts!  ] model such that ⟨ F_Lyα⟩ agrees with the measurements shown in Figure <ref>.  This allows us to cleanly compare how the shape of P(< τ_ eff^50) is affected by the IGM neutral fraction at fixed ⟨ F_ Lyα⟩.  

At z = 6 (lower right),  the early start/late end model best matches the  data.  The blue dotted curve shows that a z > 6 end to reionization results in a CDF that is too narrow, implying too little scatter in τ_ eff^50, echoing the conclusions of .  Conversely, the late start/late end model predicts too much scatter (the CDF is too wide).  This is because it has a neutral fraction of ≈ 30% at z = 6, as compared to only 15% in the early start/late end case.  At z = 5.6 (lower left), none of the models match the observations very well.  Both of the late-ending models produce too much scatter in τ_ eff^50, and the early-ending one produces too little.  This suggests that a model with a non-zero neutral fraction smaller than that in the early start/late end case (7.5%) would match the data best.  At z = 5.2 (upper right), the early-ending scenario fits the data very well, and the other two models have too much scatter in P(< τ_ eff^50).  This indicates that both late-ending models end reionization slightly too late.  At z = 5 (upper left), when the neutral fraction is < 1% in all three models, they agree well with the observations.  

At face value, observations of P(< τ_ eff^50) at 5 ≤ z ≤ 6 seem to prefer the early start/late end model.  This scenario has a lower neutral fraction at z ≤ 6 than the late start/late end case, and as such better matches the observed scatter in P(< τ_ eff^50).  However, none of the scenarios in Figure <ref> match the observed P(< τ_ eff^50) at all redshifts.  This is likely because our late-ending models finish reionization too late, as evidenced by the z = 5.2 comparison (upper right panel).  The findings of , based on these same measurements, suggest that reionization should be complete by z = 5.3 (vs. z ≈ 5 in our models), a shift of Δ z ≈ 0.3 from our models.  In Appendix <ref>, we estimate what P(< τ_ eff^50) would look like if both late-ending models finished reionization at z = 5.3.  We use the FlexRT outputs at z' = z - 0.3 to calculate P(< τ_ eff^50), then re-scale Γ_ HI in ionized gas until ⟨ F_ Lyα⟩ matches measurements.  We show that this procedure brings our simulations into better agreement with the measurements, and that the early start/late end model remains preferred.  

In , we pointed out several factors that can affect the precise timing of reionization's end in models matched to measurements of ⟨ F_ Lyα⟩.  First, lack of spatial resolution in the forest can lead to an under-estimate of the mean transmission at fixed x_ HI <cit.>, resulting in a spuriously early end to reionization (by Δ z ≈ 0.2, see Fig. 11 of ) when calibrating to measurements.  Our forest calculations include the resolution correction prescribed Appendix A of , so in principle they account for this effect.  However, those corrections were derived in a 25 h^-1Mpc box, and it is unclear how they might change in a larger box.  It is therefore possible that we have over-corrected for resolution.  We also found in  that harder ionizing spectra and less clustered ionizing sources result in an earlier end to reionization at fixed ⟨ F_ Lyα⟩.  Moreover, the uncertain clustering of ionizing sources also affects large-scale fluctuations in the ionizing background, which could be particularly intense if quasars played a large role in the end-stages of reionization <cit.>.  Any of these effects could affect when reionization needs to end to match ⟨ F_ Lyα⟩ measurements, and the shape of P(< τ_ eff^50) at fixed neutral fraction, potentially affecting our interpretation of P(< τ_ eff^50) measurements.  



  §.§.§ The Mean Free Path
 



Next, we will study the mean free path to ionizing photons (MFP, λ_ mfp), the average distance an ionizing photon travels through the IGM before being absorbed.  The MFP is sensitive to the distribution of neutral gas in the IGM and small-scale clumping in the ionized IGM <cit.>.   We calculate the Lyman limit MFP in our simulations using the definition in Appendix C of , 

    λ_ mfp^912 = ⟨∫ x df ⟩/⟨∫ df ⟩ = - ⟨∫_1^0 x df ⟩

where x is the position along a randomly oriented sightline, f(x) is the transmission of 912Å photons, and the angle brackets denote an average over many sightlines.  found that this definition matches well with forward-modeled direct MFP measurements from QSO spectra, even in a partially neutral IGM.  We caution, however, that different ways of estimating the MFP from simulations can give modestly different results <cit.>.  

The left panel of Figure <ref> shows the MFP in our late-ending models, compared to measurements from .  Both scenarios are in broad agreement with the measurements.  However, the direct measurements using QSO Lyman Continuum (LyC) spectra (all but the  points) display a preference for the late start/late end model.  This is largely due to the short z = 6 direct measurements, which prefer the rapid neutral fraction-driven decline in λ_ mfp^912 in the late start/late end case.  By contrast, the early start/late end model is ≈ 2σ away from the central values of the z = 6 measurements from .  Both models are within 1σ of the indirect, Lyα forest-based measurements from  (see also <cit.>). 

This result is consistent with previous findings that ongoing reionization at z = 6 is needed to explain the direct QSO measurements <cit.>.  Another effect at play is that the ionized IGM at z = 6 in the late start/late end was more recently ionized, and thus clumpier <cit.>.  These factors result in direct MFP measurements preferring the late start/late end scenario.  Our earlier result that P(< τ_ eff^50) measurements prefer the early start/late end model hints at a possible tension between the Lyα forest and direct MFP measurements.  This is consistent with the fact that the indirect MFP measurements from  (based on P(< τ_ eff^50) itself) at z = 6 are a factor of ≈ 2 above the direct measurements.  We emphasize that this tension is mild, since the 1 σ error bars of these measurements overlap.   



  §.§.§ IGM Thermal History
 

The IGM temperature at mean density, T_0, is shown in the right panel of Figure <ref>, alongside measurements from .  Both late-ending models display a “bump” in temperature at z ≈ 5.3 due to the end of reionization.  Heating by I-fronts <cit.> increases T_0 until near reionization's end, after which cooling from the expansion of the universe and Compton scattering off the CMB set the evolution of T_0 <cit.>.  The peak in T_0 is higher in the late start/late end model because a larger fraction of the IGM is re-ionized at z < 6.  The redshift of the bump suggested by the  measurements is closer to z ∼ 5.6 - this is consistent with our earlier finding (based on P(< τ_ eff^50)) that reionization may end Δ z ∼ 0.3 too late in our models.  

At face value, the early start/late end model agrees best with T_0 measurements.  Although both models are consistent with the  points, the late start/late end is too hot at z ≥ 5 for the measurements there.  In the early start/late end case, a larger fraction of the IGM has re-ionized at higher redshift, giving it more time to cool by z = 5.  That model is also agrees well with the reionization history in the best-fitting model of , which fits a broad range of IGM temperature measurements down to z = 2.  That model has ionized fractions of ≈ 35% (≈ 15%) at z = 8 (10), similar to our early start/late end scenario (which has 1 - x_ HI≈ 40%, (20%)).  An important caveat is that the thermal history is sensitive at the 20-30% level to the spectrum of the ionizing radiation, through a combination of the post I-front temperature (T_ reion) and photo-heating in ionized gas afterwards <cit.>.  For example, a much softer ionizing spectrum could shift the T_0 histories significantly lower at fixed reionization history (see e.g. the bottom middle panel of Fig. 3 in ).  This could bring the late start/late end model into agreement with z ≤ 5 T_0 measurements.  However, this would also require a later reionization history at fixed ⟨ F_ Lyα⟩ <cit.>, which would worsen the disagreement with the observed P(< τ_ eff^50) in Figure <ref>.  As such, we conclude that measurements of T_0 mildly prefer the early start/late end model.  



  §.§.§ Neutral fraction constraints at z ≤ 6.5



Finally, we compare our reionization models to observational constraints on x_ HI at z ≤ 6.5 obtained using QSO spectra.  Figure <ref> compares our models with constraints from Lyα forest dark pixels <cit.>, dark gaps <cit.>, QSO damping wings <cit.>, P(< τ_ eff^50) <cit.>, and Lyα forest damping wings <cit.>.  The bold lines show the reionization histories in our two late-ending models, while the faded lines show these shifted to the right by Δ z = 0.3, consistent with the discussion surrounding P(< τ_ eff^50) in <ref>.  

Most constraints are upper (lower) limits on the neutral (ionized) fraction.  Several of these are in mild tension with the late start/late end model, and most are consistent with the early start/late end case.  The z = 5.5 dark gap constraint from  and the recent QSO damping wing limits from  disfavor the late start/late end model.  The recent lower limit on the neutral fraction from , derived from Lyα forest damping wings at z = 5.8, disfavors an end to reionization early than this.  The same can be said of the  forest damping wing measurement at z = 5.6, although their measurement actually prefers the late start/late end case.  An important caveat is that if reionization ends earlier by Δ z = 0.3 (as hinted by P(< τ_ eff^50) measurements), the tension with the late start/late end model disappears.  In fact, the neutral fraction in the shifted early start/late end model cannot be much lower without being in tension with the  damping wing limit.  We conclude that neutral fraction constraints at z < 6.5 mildly prefer the early start/late end model, but that relatively small, realistic shifts in the reionization history could change this conclusion.  





 §.§ Lyα Emitters at z > 8
 

In this section, we study the Lyα transmission properties at z ≥ 8 around massive halos that could host bright LAEs, such as GN-z11 <cit.>.  Our goal is to determine whether these observations prefer an early or late start to reionization.  




  §.§.§ Examples of IGM transmission at z = 8



In Figure <ref>, we illustrate how Lyα transmission surrounding bright galaxies differs in the late start/late end and early start/late end models at z = 8.  The solid curves show the IGM transmission (T_ IGM) vs. velocity offset (v_ off) on the red side of systemic averaged over halos with M_ UV < -17.  The vertical magenta line denotes systemic Lyα, v_ off = 0.  T_ IGM goes to 0 at systemic, and at v_ off > 0 displays a shape similar to the characteristic damping-wing profile.  We see much higher Lyα transmission in the early start/late end case, owing to its much higher ionized fraction (indicated in the legend).  At v_ off < 500 km/s, T_ IGM is a factor of 2 or more above the late start/late end model.  The thin lines show individual transmission profiles for 20 sightlines surrounding the brightest galaxy in the box.  These are much higher than the average at v_ off≳ 200-400 km/s, and drop below the mean at smaller v_ off due to fast in-flowing gas around this object (see discussion of inflows in <ref>). The higher transmission at large v_ off owes to this object occupying a larger ionized region than the average galaxy.



The higher T_ IGM in the early start/late end model would seem to naturally explain the detection of bright galaxies hosting Lyα emission at z ≥ 8.  However, it is important to note that even the late start/late end model displays some transmission, even with its 16% ionized fraction, and this can be fairly high around the most biased objects (as the thin black lines show).  Indeed, even around this single halo there is significant sightline-to-sightline scatter.  This suggests that a statistical sample of observations is required to judge conclusively which model is preferred <cit.>.  It also suggests that some LAE detections at z ≥ 8 could be explainable even if reionization starts relatively late.  




  §.§.§ Visibility of LAEs



For an LAE with an intrinsic equivalent width EW_ int, and an average IGM transmission over the emitted line, ⟨ T_ IGM⟩_ line, the observed equivalent width EW_ obs is

    EW_ obs = ⟨ T_ IGM⟩_ line EW_ int

An object is detectable if EW_ obs is greater than some threshold EW_ obs^min.  This condition can be expressed as

    EW_ obs/ EW_ int  = ⟨ T_ IGM⟩_ line >  EW_ obs^min/ EW_ int≡ T_ thresh

where we have defined T_ thresh as the minimum IGM transmission that would make the LAE detectable.  To avoid assumptions about the intrinsic properties of the LAE population, we will parameterize our visibility calculations in terms of T_ thresh.  We will also adopt the common simplification that ⟨ T_ IGM⟩_ line can be approximated by T_ IGM at the v_ off of the line's emission peak[This assumption is not true in general because T_ IGM can vary significantly over the width of the emission line.  ].  This allows us to parameterize LAE visibility in the (T_ thresh, v_ off) parameter space.  In this section, we calculate LAE visibility statistics at z = 8, 9, 10, and 11.  



Recently,  studied the distribution of ionized bubble sizes in reionization models similar to ours.  Generally, it is expected that galaxies must inhabit an ionized bubble of radius ≳ 1 pMpc to guarantee a high level of Lyα transmission on the red side of systemic <cit.>.  They found that a model similar to our early start/late end scenario is required to produce a significant number of such ionized bubbles 8 ≤ z ≤ 10.  In Figure <ref>, we perform a similar analysis using our visibility calculations.  We show the fraction of LAEs with T_ IGM > T_ thresh vs. redshift for several choices of T_ thresh, v_ off, and UV magnitude range (faint vs. bright galaxies).  The caption gives these parameter combinations for each curve as a brightness (M_ UV) and a combination of T_ thresh and v_ off.  Visibility fractions increase with decreasing T_ thresh and increasing v_ off.  The latter is true because T_ IGM increases with v_ off as the damping wing opacity decreases (Figure <ref>).  Visibility is also higher for brighter galaxies, which inhabit the largest ionized bubbles.  The left and right panels show results for the late start/late end and early start/late end models, respectively.  

The solid curves show visibility for bright (M_ UV < -21) LAEs with large velocity offsets (400 km/s) and low visibility thresholds (T_ thresh = 0.2).  Such objects are visible nearly 100% of the time in the early start/late end case, and 40% of the time in the late start/late end model even at z = 11.  Reducing v_ off to 200 km/s (dashed curves) deceases visibility, especially in the late start/late end case, but even then 20-40% of LAEs are visible at 8 < z < 11.  In the early start/late end case, the visibility fraction counter-intuitively increases with redshift.  This is because the evolution in visibility is not being driven by the neutral fraction, but by inflows surrounding massive halos.  At higher redshifts, brighter objects are found in less massive halos, which are surrounded by smaller inflows.  This leads to increased transmission at v_ off = 200 km/s <cit.>.  

The dotted curves show that increasing T_ thresh from 0.2 to 0.5 (for M_ UV < -21 and v_ off = 200 km/s) has a substantial effect on visibility.  In the late start/late end model, < 20% of LAEs are visible at z = 8, and this drops to near-0 at z ≥ 9.  However, in the early start/late end case, 40-50% of such objects are visible across this redshift range.  The dot dashed curve considers faint (-19 < M_ UV < -17) galaxies with high v_ off = 400 km/s and low T_ thresh = 0.2.  These objects are visible 50-90% of the time in the early start/late end model, but < 20% of the time in the late start/late end case.  Finally, the double-dot dashed curves also show faint galaxies, but with v_ off = 200 km/s and T_ thresh = 0.5, a parameter combination that minimizes LAE visibility.  In the late start/late end case, fewer than 10% of such objects are visible at any redshift, while in the early start/late end model, 17% (8%) of such objects are visible at z = 10 (11), and over half are visible at z = 8.  

In the left panel, we show recent measurements of the fraction of galaxies hosting Lyα emitters, (the Lyα fraction, X_ Lyα^ obs), at z = 8.7 and z = 10 from .  The purple points show X_ Lyα^ obs as measured from that work.  Comparing these points directly with the fraction of LAEs that are visible assumes that the intrinsic fraction of galaxies hosting LAEs is unity - that is, X_ Lyα^ int = 1.  For this reason, we display these points as lower limits on X_ Lyα^ obs/X_ Lyα^ int.  The green points show X_ Lyα^ obs/X_ Lyα^ int assuming X_ Lyα^ int = 0.3, the Lyα fraction measured at z = 5 by , which we do not show as limits.  The curves in the left panel show wide spread that is broadly consistent with the measured visibilities.  However, in the right panel, all the curves are on the high end of the measurements (except for the double-dot dashed curve).  At face value, these findings indicate that the observed visibility of LAEs is too low for the early start/late end scenario, instead preferring the late start/late end model.  

For the global LAE visibility fraction to evolve consistently with measurements from  in the early start/late end model, the double-dot dashed curve in Figure <ref> would have to characterize the bulk of the population.  These are faint LAEs with emission at low v_ off that require high IGM transmission to observe.  While it is true that faint LAEs tend to have low v_ off <cit.>, they also tend to have fairly high EW_ intr <cit.>.  A majority of the LAEs observed at z = 5 by  in that M_ UV range have EW_ intr > 50Å, and about half have EW_ intr > 100Å.  With the visibility threshold of EW_ obs^min = 25Å used in , this would imply T_ thresh < 0.5 for the majority of faint objects, and T_ thresh < 0.25 for half of them.  Moreover, a significant fraction of the faint objects observed at z = 5-6 in  have v_ off > 200 km/s.   Brighter galaxies in their sample generally have smaller EW_ intr, but these also tend to have higher v_ off <cit.> and inhabit larger bubbles, such that they would remain visible in the early start/late end model even if they required higher T_ thresh to detect.  By contrast, in the late start/late end model, LAEs with a wide range of properties have visibilities consistent with the  measurements.   



  §.§.§ Does GN-z11 require an early start?



GN-z11 is the highest-redshift LAE detected to date, at z = 10.6.  It has a broad Lyα emission feature centered at v_ off≈ 550 km/s, a full-width at half maximum (FWHM) of Δ v ≈ 400 km/s, and an observed EW of 18Å.  Using a Bayesian analysis based on reionization simulations and an empirically derived model for the intrinsic EW distribution of LAEs from ,  inferred that the IGM must be at least 12% ionized at 2 σ confidence (yellow point in Figure <ref>).  This constraint, at face value, clearly favors the early start/late end model (see Figure <ref> in the next section).  Here, we consider whether the observed properties of GN-z11 require reionization to start early. 

We can estimate the EW_ intr required to produce the observed GN-z11 emission line as follows.  First, we model the intrinsic line as a Gaussian with some central velocity v_ off^ intr, FWHM Δ v^ intr, and amplitude A.  Then, the observed emission profile for a given sightline is given by the intrinsic profile multiplied by T_ IGM(v_ off).  We also model the observed line as a Gaussian, with parameters given in the previous paragraph, and the continuum and normalization chosen to give the observed EW.  For a sample of ∼ 2000 sightlines surrounding -22 < M_ UV < -21 galaxies, we fit for the parameters of the intrinsic line that, after attenuation by the IGM, gives a best fit to the observed line.  The distribution of EW_ intr recovered with this procedure, P( EW_ int |  EW_ obs), at z = 10.6 is shown in Figure <ref> for our late-ending models.  The shaded region denotes the range of EW observed in similarly bright galaxies at lower redshifts , the highest of which is ≈ 60Å (Fig. 1 of ). 



We see a stark contrast between P( EW_ int |  EW_ obs) in our models.  Nearly all the sightlines in the early start/late end model require EW_ intr < 60Å, and about 20% require EW_ intr < 30Å.   This suggests that bright LAEs with EWs on the high end of the observed distribution will produce a GN-z11-like observation most of the time in this scenario.  In the late start/late end model, the distribution is much wider and shifted to much higher EW_ intr.  Only 12% of sightlines allow for EW_ intr < 60Å, and none of them allow < 30Å.  So, although objects such as GN-z11 are expected to be fairly rare in the late end/late start scenario, they would not be impossible to find.  Note that the two M_ UV < -21 LAEs observed in  with the highest EWs (≈ 30 and 60Å) were both detected in Hα.  Based off a clear detection of Hγ emission in GN-z11,  estimated that it should have strong Hα emission.  As such, we can conclude that the detection of GN-z11 does not rule out the late start/late end scenario.  However, if forthcoming observations reveal similar objects to be ubiquitous at z > 10, it would be strong evidence in favor of something similar to the early start/late end model.    



  §.§.§ Constraints on x_ HI with galaxies at z ≥ 6.5



To conclude our discussion of LAEs, we look at measurements of x_ HI from observations of galaxies at z ≥ 6.5.  These include constraints from the statistics of LAE detections, and those based on Lyα damping wing absorption in galaxy spectra.  Figure <ref> shows a collection of these measurements and limits compared to our late-ending reionization models in the same format as Figure <ref> (with references in the caption).  These constraints are all model-dependent to some degree, so showing them on the same plot may not constitute a fair comparison.  Our goal here is to illustrate the diversity of constraints obtained across multiple observations and inference techniques.  



Unlike in Figure <ref>, we see no clear preference for either scenario.  Indeed, at z < 8, several constraints prefer each of the models, while some have error bars too large to distinguish them.  The only consensus that these constraints give collectively is that reionization is in progress at 7 < z < 8.  At z > 10, all the constraints are based on damping wings except that of , which is based on the detection of GN-z11.  There is no clear consensus between these constraints either.  There is a notable dearth of constraints at 8.5 < z < 10, the redshift range where the two models differ the most.  The only exception is the  point at z ∼ 8.7, which falls exactly between our models but has large error bars.  It seems clear from this comparison that constraints on x_ HI from high-redshift galaxies do not, at present, display a clear preference for either a late or early start to reionization. 



 §.§ Patchy kSZ from reionization
 



In this section, we will turn again to the CMB to help distinguish our reionization models.  We display the patchy kSZ power spectra for all three models
(see <ref>) in the left panel of Fig. <ref>, along with the 1σ error bar and 95% confidence upper limits from .  We see that the late start/late end model alone lies within the 1σ of the SPT measurement.  The other two both fall outside this range but still within the 95% confidence upper limit, with the early start/late end case coming closest to the upper limit[We note that any measurement of the pkSZ involves assumptions about the late time, homogeneous kSZ.  The methods in  assume an end to reionization at z_ end = 5.5, whereas our simulations (and our pkSZ contributions) continue until z_ end∼ 5.  Correcting for this would bring the measurement up slightly (∼ 0.1 μK^2), but not enough to qualitatively affect our conclusions.  ].  We also include the revised 2 σ upper limit from , which is somewhat lower than the SPT result and favors the late start/late end model even more.  

To gain intuition for the origin of the differences in pkSZ power, we plot the differential contribution to D_ℓ=3000 per z in the right panel.  Both early-starting models begin contributing power as soon as reionization starts at z = 18, as ionized bubbles form and grow to sufficient scales.  The two begin diverging at z ∼ 10, as the early start/early end case finishes reionization, causing the pkSZ power at ℓ=3000 to drop abruptly at z = 8 (see annotation).  This fall-off in power corresponds to the disappearance of large-scale ionization fluctuations, at which point the features in the kSZ signal on these scales are set by fluctuations in density and velocity only.  In contrast, ionization fluctuations persist longer in the early start/late end model, and so continue to contribute power to the pkSZ signal at ℓ = 3000 at z < 8.  

In the late start/late end case, reionization begins much later but still ends at z = 5, which makes the peak in dD_3000^ pkSZ/dz narrower.  The shaded red region shows that nearly half the power in the early start/late end case arises at z > 10 when the ionized fraction is < 20% in that model.  In the late start/late end case, reionization is just starting around z = 10.  Thus, we see that the pkSZ is highly sensitive to reionization's duration, and particularly its early stages <cit.>.  We see also that although the  measurement does not rule out the early start/late end model at 2σ, it clearly prefers the late start/late end case.  This finding is consistent with the recent limits on the duration of reionization from  using data from SPT and the Herschel-SPIRE experiment.  They found that Δ z_50, the difference between redshifts at 25% and 75% ionized fractions, is < 4.5 at 95% confidence[One caveat is that their constraint assumes that reionization ends by z = 6, whereas in our models it completes at z ≈ 5.  Their constraint would loosen by dz ∼ 1 if z = 5 were assumed.  ] - our early start/late end model has Δ z_50 = 3.1.  



§ DISCUSSION
 



 §.§ “Face-value” interpretations of the data


 

In the previous sections, we studied how the properties of our three models compare to a broad range of observables.  These include measurements of the UVLF and ξ_ ion from JWST (<ref>), inferences from the spectra of high-redshift QSOs (<ref>), Lyα transmission from z > 8 galaxies (<ref>), and constraints from the CMB (<ref>).  We concluded in <ref> that measurements of the Lyα forest at z ≤ 6 strongly disfavor the early start/early end scenario.  However, the evolution of the ⟨ F_ Lyα⟩ alone could not distinguish between an early and late start to reionization.  Most of the other observables we studied individually displayed a preference for one or the other, but none could conclusively rule out either[We note that these preferences were determined primarily qualitative - based on “chi-by-eye” comparisons of models and observations, and as such they represent somewhat qualitative results.  Still, they are useful for gauging the direction that each data set will likely to push constraints on reionization if included in a quantitative analysis <cit.>.  ].  This motivates the key question in this work: when taken together, what story do these observables tell about reionization's early stages?  Indeed, a major goal of the field is to synergize the constraining power of many observations to constrain reionization, and qualitative analyses like the one presented here can help chart the path for more detailed studies.  

We summarize our findings in Table <ref>.  The left-most column lists the “categories” of observables that we studied - the CMB (blue), high-redshift galaxies (red), and z < 6.5 QSOs (magenta).  The second column from the left lists each of the observables, and the remaining two columns denote whether each observable prefers a late or early start to reionization.  We find that τ_ es, ⟨ F_ Lyα⟩, and measurements of x_ HI at z > 6.5 display no preference for either case.  JWST observations of the UVLF/ξ_ ion, the MFP, LAE visibility at z > 8, and the SPT pkSZ measurement prefer the late start/late end case.  By contrast, the Lyα forest P(<τ_ eff^50), the IGM thermal history, and x_ HI measurements at z < 6.5 prefer the early start/late end model.  

A key finding is that not all these observables prefer the same scenario.  This suggests a possible lack of concordance between different data sets with respect to reionization's early stages.  We note, however, that all three observables that favor the early start/late end model are based QSO spectra at z ≤ 6.5.  Indeed, nearly all the data associated with these three observables (with the exception of the QSO damping wings) arises, directly or indirectly, from the Lyα forest, and thus cannot be treated as fully independent[Indeed, much of the data comes from the same QSO survey, XQR-30 <cit.>.  ].  As we explained in <ref>, the conclusions we drew from these probes are sensitive to our modeling choices, which are necessary to link the late stages of reionization (which the forest probes directly) to its early stages.  By contrast, the observables that support a late start are derived from different data sets using vastly different techniques, and at least one probe in every category prefers a late start.  As such, we judge that the consensus of these probes from a wide range of data sets indicates that the late start/late end model is mildly preferred (overall) by observations.  

The lack of consensus between different observables has several possible resolutions.  Perhaps the most straightforward is that existing observations lack the accuracy and/or precision to achieve a unanimous consensus about the early stages of reionization.  Indeed, nearly all the observables considered here still have large uncertainties.  Theoretical modeling uncertainties, required to interpret the data, can similarly affect these conclusions.  As mentioned earlier, the QSO-based observables that seem to support an early start probe only reionization's end stages, requiring a model to infer the early history.  These issues will continue to improve with time as more (and better) observational data is acquired and reionization models become faster and more accurate.  

However, a more concerning possibility is that forthcoming observations and rigorous theoretical analysis will reveal a statistically significant tension between different observables with respect to reionization's early stages.  In this case, the fault must lie with observational systematics and/or hidden deficiencies in theoretical modeling.  Either possibility presents a potential pitfall for efforts to constrain reionization with multiple data sets.  Such constraints may be artificially tight if “tensions” between data sets exist.  This could lead to the pre-mature conclusion that the reionization history is known to high precision.  Forthcoming efforts should be aware of this potential pitfall, and take care to understand the effects of individual data sets on joint constraints.  



 §.§ Forthcoming observational prospects



In this section, we will briefly discuss prospects for future observations that would help strengthen the constraining power of some of the probes discussed here.  The first is to continue improving constraints on the UVLF and ξ_ ion, particularly for faint galaxies.   demonstrated that ξ_ ion could be measured reliably for faint (-17 < M_ UV < -15), lensed galaxies during reionization.  Such studies, together with efforts to directly constrain the faint end of the UVLF, will be crucial for determining the redshift evolution of these quantities and whether there is a fall-off in ionizing output for the faintest galaxies (see bottom panel of Figure <ref>).  Continued efforts to understand how f_ esc correlates with galaxy properties at low redshift, such as the Low-redshift Lyman Continuum Survey (LzLCS, ) will also be crucial for placing reasonable limits on the evolution of f_ esc (see also e.g. ).  

There is also further progress to be made with QSO-based observations at z ∼ 6.  Improved constraints on the mean free path and IGM thermal history may help distinguish an early vs. late start, as Figure <ref> shows.  Forthcoming observations with Euclid <cit.> will dramatically increase the number of known quasars at these redshifts, allowing for spectroscopic follow-up that will improve statistical uncertainties on both sets of measurements.  Efforts to measure the relationship between Lyα forest opacity and galaxy density <cit.> may also help tighten constraints on the reionization history at z < 6.5 <cit.>.  

Further observations with JWST will improve the statistics of z > 8 galaxies displaying significant Lyα emission.  They will also yield constraints on x_ HI from Lyα damping wing absorption at redshifts where very few bright quasars are available.  Forthcoming observations with the Nancy Grace Roman telescope <cit.> will also reveal bright LAEs over a much wider area than JWST, enabling improved constraints on the early reionization history <cit.>. 

Forthcoming improvements on CMB constraints from multiple experiments, including the Atacama Cosmology Telescope <cit.>, SPT <cit.>, Simons Observatory <cit.>, and CMB-S4 <cit.> will improve constraints on τ_ es and pkSZ.  They will may also detect new signals that probe reionization, such as patchy τ <cit.> and higher-order statistics and cross-correlations with other signals <cit.>.  These will help constrain the early stages of reionization because of their sensitivity to its duration and morphology <cit.>.  



§ CONCLUSIONS
 

In this work, we have studied the observational properties of three representative reionization histories. In the first, reionization starts early and ends at z ∼ 8, earlier than suggested by the Lyα forest and τ_ es.   This scenario is motivated by recent JWST observations of the UVLF and ξ_ ion at z > 6, which, when combined with observationally motivated assumptions about f_ esc, suggest copious ionizing photon output by high-redshift galaxies.  We have investigated the observational properties of this model, alongside two others in which reionization ends much later at z ∼ 5, in agreement with the Lyα forest and τ_ es.  One model starts reionization relatively late at z ∼ 9, and the other starts early at z ∼ 13.  



    
  * We find, consistent with previous work, that the early start/early end scenario severely violates high-redshift QSO observations, most notably the Lyα forest.  These observations require reionization to end at 5 < z < 6, or at least not much sooner.  This is consistent with recent measurements of τ_ es from .  Unfortunately, neither the mean transmission of the Lyα forest nor τ_ es display a clear preference for whether reionization started late or early.  The former measures only the global average transmission of the IGM at reionization's end, and the latter is only an integrated constraint, and thus does not uniquely constrain reionization's early stages.   

    
  * Observations of the UVLF and ξ_ ion by JWST, direct measurements of the MFP from QSO spectra, the visibility of z ≥ 8 LAEs, and the recent SPT measurement of the patchy kSZ all prefer a late start to reionization.  In light of the latest UVLF measurements at z ≥ 8 from JWST, our early-starting model requires an order of magnitude of evolution in galaxy ionizing properties (quantified by ⟨ f_ escξ_ ion⟩_L_ UV) between z = 6 and 12.  This is less compatible with observations than the flat evolution in our late-starting model.  Direct measurements of the MFP from QSO spectra also prefer a late start, mainly because of the high neutral fraction needed to match direct measurements from QSO spectra at z = 6.  The steep drop-off in LAE visibility at z > 8 observed by  is more consistent with a late than an early start.  Finally, the low central value of the SPT pkSZ measurement prefers a late start, and disfavors our early-starting model at almost 2 σ.  
    
    
  * By contrast, the distribution of Lyα forest opacities, the thermal history of the IGM, and measurements of x_ HI at z < 6.5 prefer an early start.  The forest P(< τ_ eff) is too wide for the observations in our late start/late end model, preferring instead the lower neutral fraction in the early start/late end case.  Constraints on x_ HI at z ≤ 6.5 from a variety of QSO-based inferences suggest a similar conclusion.  The cooler IGM in the early-starting case at z ≤ 5 is also in better agreement with observations.  

    
  * Our findings suggest that no single probe can conclusively rule out either the late start/late end or early start/late end model in favor of the other.  However, we do find that observations across multiple independent datasets - JWST observations of galaxy properties, LAE detections, the CMB, and QSO absorption spectra - prefer a late start to reionization.  The observables that prefer an early start are all derived from the same type of observations (QSO spectra) and only probe the tail end of reionization.  As such, they are not fully independent probes, and they require a model to derive inferences about reionization's early stages.  As such, we conclude that overall, existing observational data displays a mild preference for the late start/late end scenario.   

    
  * The face-value disagreement we find between different probes suggests that (1) present observations, and the models used to interpret them, are insufficiently precise/accurate to paint a consensus picture of reionization's early stages, and/or (2) there are systematic effects (in observations and/or theoretical modeling) leading to the appearance of tension.  The second possibility motivates care in interpreting the results of analyses using multiple data sets.  Joint analyses using many observables could lead to artificially tight constraints on the reionization history and other quantities if tensions arising from systematics are not carefully understood.   



Forthcoming work, on the observational and theoretical side, should continue working to synergize the information available from many observables.  Our work motivates further efforts targeting the early stages of reionization, which will yield key insights into the evolution of galaxy properties across the reionization epoch and into the cosmic dawn era.  This work is also a cautionary tale that motivates careful understanding of potential systematics in both observations and theoretical modeling.  Such systematics, if not studied carefully, could lead to premature conclusions about reionization.  



CC acknowledges helpful conversations with Seth Cohen, Timothy Carleton, Evan Scannapieco, Frederick Davies, and Kevin Croker, and support from the Beus Center for Cosmic Foundations.  A.D.'s group was supported by grants NSF AST-2045600 and
JWSTAR-02608.001-A.  RAW acknowledges support from NASA JWST Interdisciplinary Scientist grants
NAG5-12460, NNX14AN10G and 80NSSC18K0200 from GSFC. JBM acknowledges support from NSF Grants AST-2307354 and AST-2408637, and thanks the Kavli Institute for Theoretical Physics for their hospitality.  Computations were made possible by NSF ACCESS allocation TG-PHY230158. 




aasjournal





§ P(< Τ_ EFF^50) FOR SHIFTED REIONIZATION HISTORIES



In Figure <ref>, we show P(< τ_ eff^50) for the late start/late end and early start/late end models with their reionization histories shifted earlier by Δ z = 0.3, as described in <ref>.  We show 6 redshifts between z = 5 and 6 in intervals of 0.2.  This exercise shows roughly what P(< τ_ eff^50) would look like if these models ended reionization at z = 5.3 instead of 5.  We see that at z = 5, 5.2, and 5.4, there is little difference between the models.  At z = 5.6 and 5.8,  P(< τ_ eff^50) is slightly narrower in the early start/late end case, as it is in Figure <ref>.  The difference is that the early start/late end model now seems to agree well with the observations, whereas in the late start/late end case, P(< τ_ eff^50) is still too wide.  At z = 6, it is difficult to tell which model is a better fit to the data due to the large number of non-detections in the data, which set the width of the green shaded region.  This test demonstrates that even if reionization ends significantly earlier than it does in our models, the early start/late end scenario remains preferred by P(<τ_ eff^50). 





§ LYΑ VISIBILITY IN THE FULL (T_ THRESH, V_ OFF) PARAMETER SPACE


In this appendix, we show complete results for our LAE visibility analysis in terms of v_ off and T_ thresh presented in <ref>.  Figure <ref> shows the fraction of visible LAEs at T_ thresh < 0.5 and 0 < v_ off < 800 km/s, at the four redshifts (columns) and three magnitude ranges (rows) we considered in the late start/late end model.  Red (blue) regions denote high (low) LAE visibility fractions.  The white contour lines denote visibility fractions of 50%, 25%, and 10% (see annotation in the left-most panel of the middle row).  The hatched white box in the upper left corner of each panel denotes the region where T_ thresh > 0.2 and v_ off < 500 km/s.  We expect a majority of z ≥ 8 LAEs to inhabit this region.  Most LAEs in M_ UV < -17 galaxies observed at slightly lower redshifts (z ∼ 5-6, when T_ IGM is close to unity) have v_ off < 500 km/s and EW_ int < 500 Å <cit.>.  For T_ thresh = 0.2, the latter would correspond to EW_ obs^min < 100 Å.  



At z = 8, a significant portion of the hatched region displays a high visibility fraction, especially for the brightest galaxies.  For M_ UV < -21 galaxies, LAEs with v_ off > 200 km/s have a significant chance of being visible, provided they are bright enough to be observed with a factor of 2 IGM attenuation.  Fainter galaxies require somewhat fainter detection thresholds, since on average they inhabit smaller ionized bubbles and are more sensitive to IGM attenuation.  In the faintest M_ UV bin, v_ off > 400 kms/ and T_ thresh < 0.25 is required for half of LAEs to be visible.  Still, we should expect to observe some LAEs at z = 8 in the late start/late end model, especially the brightest ones.  

At z > 8, the transmission of Lyα declines rapidly, especially for fainter galaxies.  At these redshifts, in all but the brightest M_ UV bin, the 50% contour line does not intersect the hatched region.  As we showed in Figure <ref>, this drop-off in visibility is consistent with the observed decline in X_ Lyα observed by  at z > 8.  The brightest objects remain likely to be observed if v_ off > 400 km/s and T_ thresh < 0.25 all the way to z = 11.  Notably, GN-z11 (M_ UV = -21.5, z = 10.6, ) has v_ off = 550 km/s, meeting this condition.  The recently observed JADES-GS-z9-0 (M_ UV = -20.43, z = 9.4, ) has v_ off = 450 km/s.  These objects would be likely visible in the late start/late end model if their intrinsic EWs were in the neighborhood of 100 Å, ≈ 4× larger than their observed EWs (see <ref>).  These results are consistent with a universe in which most LAEs at z > 8 are obscured by the IGM, but a small number of objects that are relatively bright, have high v_ off, and/or high intrinsic EWs remain visible.    

Figure <ref> is the same as Figure <ref>, but for the early start/late end model.  In contrast to the late start/late end case, a significant fraction of the parameter space displays high visibility, even at z = 11.  Objects with M_ UV < -21 are likely to be visible up to z = 11 as long as T_ thresh≤ 0.5 and v_ off > 200 km/s.  Even the faintest galaxies have a significant chance of being observed at z = 11.  It should then be expected that in such a scenario, a significant fraction of LAEs - even faint ones - with typical physical properties should remain visible up to z = 11.  At face value, the observed sharp decline in LAE visibility across the population of LAEs up to this redshift does not prefer this scenario.  












