Recent JWST observations at z gt 6 may imply galactic ionizing photon production in excess of prior expectations. Under observationally motivated assumptions about escape fractions these suggest a z sim 89 end to reionization in strong tension with the z lt 6 end required by the Lyalpha 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 approx 8 alongside two more realistic scenarios that end late at z approx 5 one that starts late z sim 9 and another that starts early z sim 13. We find that the latter requires up to an orderofmagnitude evolution in galaxy ionizing properties at 6 lt z lt 12 perhaps in tension with recent measurements of xirm ion by JWST which indicate little evolution. We also study how these models compare to recent measurements of the Lyalpha forest opacity mean free path IGM thermal history visibility of z gt 8 Lyalpha emitters and the patchy kSZ signal from the CMB. We find that neither of the lateending 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 CMBS4 in the coming years will either establish a concordance picture of reionizations early stages or reveal systematics in data andor theoretical modeling.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Chasing the beginning of reionization in the JWST era
{'Chasing the beginning of reionization in the JWST era': "Christopher Cain, 1 Garett Lopez, 2 Anson D'Aloisio, 2 Julian B. Mu˜noz, 3 Rolf A. Jansen, 1 Rogier A. Windhorst, 1 and Nakul Gangolli 2  \n1  \nSchool of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287-6004, USA 2 Department of Physics and Astronomy, University of California, Riverside, CA 92521, USA Department of Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA  \n3", '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.", '1. 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 highredshift QSOs (Becker et al. 2015; Kulkarni et al. 2019; Keating et al. 2020; Nasir & D'Aloisio 2020; Qin et al. 2021; Bosman et al. 2022). Direct (and indirect) measurements of the mean free path from QSO spectra (Worseck et al. 2014; Becker et al. 2021; Zhu et al. 2023; Gaikwad et al. 2023; Davies et al. 2024b), and measurements of the IGM thermal history at z ≥ 5 . 5 (Gaikwad et al. 2020) have further corroborated this picture. The electron scattering optical depth to the CMB ( τ es , Planck Collaboration et al. 2020) constrains the midpoint of reionization to be z ∼ 7 . 5. Ev-  \nCorresponding author: Christopher Cain  \nclcain3@asu.edu  \nidence of damping wings in high-redshift z ∼ 7 -8 QSOs (Davies et al. 2018; Wang et al. 2020; Yang et al. 2020) and measurements of the neutral fraction based on (non)detections of Ly α emitters (LAEs, e.g. Mason et al. 2018; Mason et al. 2019; Whitler et al. 2020; Jung et al. 2020) indicate that the IGM was partially neutral at these redshifts.  \nUnderstanding 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 (Dunkley et al. 2009; Komatsu et al. 2011) which required very early and/or extended reionization histories. This tension eased as the measured value of τ es steadily decreased. Several recent works (e.g. Bouwens et al. 2015; Robertson et al. 2015; Finkelstein et al. 2019; Matthee et al. 2022) 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 Davies et al. 2021, 2024a).  \nConcurrent efforts demonstrated that AGN are unlikely to have contributed the majority of the ionizing budget responsible for reionization (e.g. Dayal et al. 2020; Trebitsch et al. 2023, although see Madau & Haardt 2015; Madau et al. 2024).  \nThese 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 (Planck Collaboration et al. 2020) 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 1 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.  \nHowever, 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 (Finkelstein et al. 2024; Adams et al. 2024; Donnan et al. 2024). It has also allowed us to measure the ionizing efficiency of galaxies, ξ ion , above z = 6 (Endsley et al. 2024; Simmonds et al. 2024; Pahl et al. 2024). Recently, Mu˜noz et al. 2024 pointed out that a face-value interpretation of recent UVLF and ξ ion measurements from JWST (Donnan et al. 2024; Simmonds et al. 2024) 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 2 . This result represents a stark reversal from the historical problem of galaxies producing too few ionizing photons to complete reionization on time (see Robertson et al. 2015, and references therein). Several bright Ly α emitters (LAEs) at z > 8 (Zitrin et al. 2015; Larson et al. 2022; Bunker et al. 2023; Curti et al. 2024; Tang et al. 2024) have also been observed, with the highest-redshift detection to date at z = 10 . 6 (Bunker et al. 2023) by JWST. These observations may be surprising if the IGM is mostly neutral at these redshifts, since observing Ly α  \nFigure 1. Summary of the reionization scenarios studied in this work. Broadly speaking, they represent the three types of possible reionization histories. Top: volume-averaged ionized fraction x V HII vs. redshift for the late start/late end , early start/late end , and early start/early end models. These models have z end ≈ 5, 5, 8 and z mid ≈ 6 . 5, 7 . 5, 8 . 5, respectively. Bottom: CMB electron scattering optical depth τ es for each model, compared to the results of Planck Collaboration et al. 2020; de Belsunce et al. 2021 (shaded regions denote ± 1 σ ). The late start/late end model is within 1 σ of Planck, while early start/late end is within 1 σ of the re-analysis by de Belsunce et al. 2021.  \n<!-- image -->  \nemission requires some level of ionization around galaxies (Mason & Gronke 2020).  \nThe top panel of Figure 1 shows the volume-averaged ionized fraction ( x V HII ) for the three reionization models we will study in this work. The early start/early end model is motivated by the aforementioned findings of Mu˜noz et al. 2024, 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 Planck Collaboration et al. 2020 measurement (black) and the recent re-analysis of Planck data from de Belsunce et al. 2021, with shaded regions denoting ± 1 σ uncertainties. The late start/late end model is consistent within 1 σ with Planck Collaboration et al. 2020, and the early start/late end case is similarly consistent with de Belsunce et al. 2021.  \nIn this work, we study the observational properties of the three reionization models shown in Figure 1 us-  \ning 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. § 2 describes how we calibrate our three reionization models and discusses their implications for galaxy properties in light of recent JWST observations. In § 3, 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 § 4, discuss the implications of our findings in § 5, and conclude in § 6. 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 Planck Collaboration et al. 2020 results. All distances are in co-moving units unless otherwise specified.", '2.1. Reionization Model Calibration': "The main input to our RT simulations (described in § 3) is the globally averaged ionizing photon emissivity of sources verses redshift, ˙ N γ ( z ). In this section we describe how we use a combination of JWST observations and Ly α forest data to construct (or 'calibrate') ˙ N γ ( z ) for our three models.  \nOur 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 (e.g. Harikane et al. 2024; Finkelstein et al. 2024; Adams et al. 2024; Donnan et al. 2024). The top two panels of Figure 2 show two sets of UVLFs that we use in our analysis. In both panels, the dashed curves show the z < 8 UVLFs measured by Bouwens et al. 2021 with HST. In the left panel, the solid curves denote the double-power-law (DPL) fits from Adams et al. 2024 at 8 ≤ z ≤ 12, and the dotted curve is the measurement from Donnan et al. 2024 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 Donnan et al. 2024 (which is a best-fit to measurements from Bowler et al. 2016, 2020; Donnan et al. 2022; McLeod et al. 2023). At z > 8, the UVLFs in the left panel are a factor of ∼ 3 lower than those in the right  \npanel. We will use these sets to roughly bracket observational uncertainties on the UVLF.  \nThe 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 cut UV = -17 (black) and a fainter M cut UV = -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 cut UV , 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 (Kravtsov & Belokurov 2024). 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 cut UV (e.g. Finkelstein et al. 2019). Instead, ρ UV evolves quickly with redshift, with a factor of ∼ 10 evolution between z = 8 and 13.  \nThe lower right panel quantifies the redshift evolution of galaxy ionizing properties in our models using  \n⟨ f esc ξ ion ⟩ L UV ≡ ˙ N γ ( z ) ρ UV ( z ) (1)  \nThis 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 cut UV = -13 (lower left panel). The shaded regions show the spread in this quantity (for fixed ˙ N γ ( z )) arising from observational uncertainty in ρ UV ( z ), as shown in the lower left. We calibrate ˙ N γ ( z ) for the early start/early end model such that its ⟨ f esc ξ ion ⟩ L UV agrees with the observationally motivated model from Mu˜noz et al. 2024 (which assumes M cut UV = -13), shown by the faded gray dashed curve. Consistent with the findings of Mu˜noz et al. 2024, this model ends reionization early at z ≈ 8 (see Figure 1).  \nOur 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 ˙ N γ ( z ). This could be the case if the measurements of ξ ion assumed in Mu˜noz et al. 2024 (from Simmonds et al. 2024) are systematically biased high, and/or if the same is true of the f esc values they inferred from the results of Chisholm et al. 2022; Zhao & Furlanetto 2024. We find that scaling ˙ N γ ( z ) down by a factor of 5 from the early start/early end case brings the end of reion-  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 2. Evolution of galaxy properties and calibration of ˙ N γ ( z ) for our reionization scenarios. Top Left: measured UVLFs from Bouwens et al. 2021 at z < 8, Adams et al. 2024 at 8 < z < 12 . 5, and Donnan et al. 2024 at z = 14 . 5. Top Right: The same, but using UVLF parameters from Eq. 3-6 of Donnan et al. 2024 at z ≥ 8. These two sets are different by a factor of ∼ 3 at z > 8, roughly bracketing observational uncertainty in the UVLF. Lower left: integrated ρ UV ( z ) for cutoff magnitudes M cut UV = -17 and -13. The shaded regions denote the spread between the two sets of UVLFs in the top panels, and the curves denote the logarithmic averages. Lower Right: product of f esc and ξ ion , weighted by L UV and averaged over the galaxy population (Eq. 1). The colored lines show this quantity for our models, and the shaded regions show the spread resulting from observational uncertainty in the UVLF (at fixed ˙ N γ ( z )). The faded gray dashed curve shows one of the observationally motivated models assumed in Mu˜noz et al. 2024 (with M cut UV = -13) and the faded cyan dashed curve shows the model of Pahl et al. 2024 scaled up by a factor of 3 . 5. Our early start/early end scenario (by construction) approximately matches the Mu˜noz et al. 2024 model. The late start/late end case shares the same redshift evolution, but is scaled down by a factor of 5 to recover agreement with the z ≤ 6 Ly α forest. The early start/late end model is calibrated to agree with the Mu˜noz et al. 2024 result at z ≥ 10, but it also required to agree with the Ly α forest. This model demands a factor of ≈ 10 decline in ⟨ f esc ξ ion ⟩ L UV between z = 12 and 6 - much steeper evolution than suggested by the Mu˜noz et al. 2024 or Pahl et al. 2024 findings.  \n<!-- image -->  \nto z ∼ 5. We then further adjust ˙ N γ ( 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 § 4.1.1). The cyan-dashed curve shows the observationally motivated model from Pahl et al. 2024, multiplied 3 . 5, which also has redshift evolution similar to this scenario.  \nHowever, the late start/late end scenario is not the only possibility allowed by τ es and the Ly α forest. It could also be that ˙ N γ ( z ) is consistent with the Mu˜noz et al. 2024 model at the highest redshifts ( z ≳ 10), but  \ndeclines 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 ˙ N γ ( 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 Mu˜noz et al. 2024 and Pahl et al. 2024. This could be achieved if ξ ion , f esc , M cut UV ,  \nor some combination of these evolves significantly across this redshift range (see § 2.2).", '2.2. 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 Mu˜noz et al. 2024 and Pahl et al. 2024 models, it is natural to ask whether this scenario is plausible. The top panel of Figure 3 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 Mu˜noz et al. 2024 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 Mu˜noz et al. 2024) 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 Mu˜noz et al. 2024 with the global f esc ( z ) from the flagship THESAN simulation (Yeh et al. 2023). 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 Mu˜noz et al. 2024. 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.  \nThe bottom panel of Figure 3 illustrates another mechanism that could work in the direction of making reionization start earlier: evolution in M cut UV . The red curve assumes that the M cut UV evolves from -10 at z = 12 to -19 at z = 6 ( M cut UV ( 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 2. This allows ⟨ f esc ξ ion ⟩ L UV in the early start/late end scenario to evolve less quickly, in better agreement with the Mu˜noz et al. 2024 model (dotted curve). This type of behavior in the galaxy population could arise from decreasing dust obscuration (Topping et al. 2024) and/or feedback from the IGM reducing or shutting off star formation in low-mass halos at lower redshifts (Wu et al. 2019; Ocvirk et al. 2021). The evolution in M cut UV assumed in this illustrative example is extreme, and likely ruled out by existing observations (Atek et al. 2018), but serves to show how an evolving M cut UV could support the early start/late end model.  \nThese comparisons suggest that a factor of a few of evolution in each of ξ ion and f esc , and perhaps some in M cut UV , could explain the early start/late end scenario. Indeed, there is some observational and theoret-  \nFigure 3. Summary of mechanisms that could allow for the early start/late end model. Top: ⟨ f esc ξ ion ⟩ L UV (as in Figure 2), alongside several observationally and theoretically motivated scenarios. The grey dotted curve shows the Mu˜noz et al. 2024 model multiplied by 0 . 2, which matches the late start/late end model . The gray dashed curve is the same model, but with caps on extrapolation of ξ ion with M UV and redshift imposed by Mu˜noz et al. 2024 removed (see text). The dot-dashed curve further replaces the model for f esc used in Mu˜noz et al. 2024 with the average f esc ( z ) measured in the flagship THESAN simulation (Yeh et al. 2023). This last curve reproduces evolution qualitatively consistent with the early start/late end model. Bottom: the same, but assuming M cut UV evolves linearly from -19 at z = 6 to -10 at z = 12 (see annotation). A trends towards fainter M cut UV at higher z causes ρ UV to decline less steeply, requiring shallower redshift evolution of ionizing properties. See text for further details.  \n<!-- image -->  \ncal support for the idea that ξ ion could increase with redshift and/or be larger in fainter galaxies (Atek et al. 2024a; Cameron et al. 2024; Simmonds et al. 2024), However, some works find that ξ ion is closer to constant with M UV , and perhaps redshift (Matthee et al. 2022; Pahl et al. 2024). Unfortunately, f esc at these redshifts cannot be directly measured, and there is no consensus in the literature from either indirect observational estimates (e.g. Naidu et al. 2022; Citro et al. 2024) or cosmological simulations (e.g. Trebitsch et al. 2018; Rosdahl et al. 2022; Kostyuk et al. 2023). Thus, galaxy observations alone cannot rule out the early  \n<!-- image -->  \nFigure 4. Top: ionizing photon emissivity for our models, alongside several others from the literature that are calibrated to match the Ly α forest observations at z ≤ 6. Our late start/late end model has the steepest evolution at z > 9, and our early start/late end model agrees well with the 'early' model from Asthana et al. 2024. Bottom: ionizing photon production per H atom by galaxies. The vertical dashed lines denote the end of reionization in each model, and the horizontal lines indicate the number of photons per H atom required to complete reionization.  \n<!-- image -->  \nstart/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 ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model seems more likely based on existing data.", '2.3. Ionizing photon budget': "We show the ˙ N γ ( z ) for our models in the top panel of Figure 4, alongside several others from the literature that also match the z < 6 Ly α forest. The cyansolid and magenta-dashed faded curves show the 'fiducial' and 'early' models from Asthana et al. 2024, and the dot-dashed gray curve shows the reference model 3 from Cain et al. 2024. Our late start/late end model drops off faster than the others at z > 9, and  \nour early start/late end model closely matches the 'early' model from Asthana et al. 2024. The bottom panel of Figure 4 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.  \nAll 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, Davies et al. 2024a made a first attempt to observationally measure the 'clumping factor', which quantifies the recombination rate in the in-homogeneous IGM (Gnedin & Ostriker 1997; Pawlik et al. 2010). They found C ∼ 12 at z > 5, higher than the commonly assumed C = 3 (e.g. Finlator et al. 2011; McQuinn et al. 2011). 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 Davies et al. 2024a measurements 4 . 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 5 .", '3. 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 § 4.', '3.1. Radiative Transfer Simulations': "We ran RT simulations of reionization using FlexRT (Cain & D'Aloisio in prep.). FlexRT is an adaptive ray  \ntracing 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 (D'Aloisio et al. 2020; Nasir et al. 2021). We refer the reader to (Cain et al. 2021, 2023), and a forthcoming code paper (Cain & D'Aloisio in prep.) for further details about FlexRT.  \nIonizing 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 3 M) code of Trac et al. 2015. This run has a box size of L box = 200 h -1 Mpc 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 Trac et al. 2015 to within ≈ 5 -10% (HMF) for M halo > 3 × 10 9 h -1 M ⊙ , and is ≈ 50% incomplete at M halo = 1 × 10 9 h -1 M ⊙ . 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 -1 M ⊙ , 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.  \nWe assign UV luminosities ( L UV ) to halos by abundance-matching to observed UV luminosity functions (UVLFs). We use the measurements of Bouwens et al. 2021 at z < 8, those of Adams et al. 2024 at 8 < z < 14, and that of Donnan et al. 2024 at z = 14 . 5 - those shown in the upper left panel of Figure 2. The ionizing emissivity is distributed between halos following ˙ n γ ∝ L UV 6 . The volume-averaged ionizing emissivity of sources, ˙ N γ = 1 L 3 box ∑ i ∈ halos ˙ n i γ , is set by hand at each redshift and used to normalize the ˙ n γ proportionality. As explained in § 2.1, we calibrated ˙ N γ ( z ) for each scenario using a combination of JWST data ( § 2.1) and  \nmeasurements from the z < 6 Ly α forest ( § 3.2, § 4.1.1). 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 7 ( 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 D'Aloisio et al. 2019. The subsequent thermal history is calculated using their Eq. 6.", '3.2. 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 8 . The native spatial resolution of our density field is ∆ x cell = 97 . 6 h -1 kpc, too coarse to fully resolve the low-density voids that set the transmission at 5 < z < 6 (Doughty et al. 2023). We apply the multiplicative correction prescribed in Appendix A of D'Aloisio et al. 2018 to our residual neutral fractions to get an effective resolution of ∆ x cell = 12 . 2 h -1 kpc. The Ly α voigt profile is approximated using the analytic fit from Tepper-Garc'ıa 2006.  \nSince our gas temperatures are calculated on the coarse RT grid, we do not capture the temperaturedensity 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 Cain et al. 2024 (see also their Appendix E), which uses the IGM temperature model of McQuinn & Upton Sanderbeck 2016. 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.  \nWe compute Ly α forest statistics by casting 4000 sightlines from random locations and in random directions of length 50 h -1 Mpc, for a total path length of 200 h -1 Gpc. We calculate transmission statistics from  \nz = 4 . 8 to z = 6 in ∆ z = 0 . 2 increments. Since our native resolution of 97 . 6 h -1 kpc 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.", '3.3. 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 (Verhamme, A. et al. 2006). 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 (see Mason & Gronke 2020, and references therein). 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 (Park et al. 2021; Smith et al. 2022).  \nIn 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 ˚ A 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 9 . We begin integrating the Ly α opacity 500 h -1 kpc (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.  \nThe gas around massive halos can have sharp line-ofsight 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 in-  \nerpolate 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 Gangolli et al. 2024 for a description of a similar procedure). We find the transmission profiles to be well-converged at this resolution.', '3.4. 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 (Sunyaev & Zeldovich 1980). This is known as the Sunyaev-Zel'dovich Effect, and its resultant temperature deviation along a line of sight is given by:  \n∆ T T = -σ T n e, 0 ∫ e -τ es ˆ γ · ⃗q c ds a 2 (2)  \nwhere ˆ γ 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 SunyaevZel'dovich Effect (hkSZ), and a highz 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 Reichardt et al. 2021, see § 4.3) must assume a fixed z end in their analysis.  \nWe calculate the signal using a method similar to that first suggested by Park et al. 2013. The kSZ angular power spectrum can be calculated from the 3D power spectrum of ⃗q by:  \nC ℓ = ( σ T n e, 0 c ) 2 ∫ ds s 2 a 4 e -2 τ P q ⊥ ( k = ℓ/s, s ) 2 (3)  \nwhere (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 = ∫ ⃗qe -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  \nwhen integrating over the line of sight. We refer the reader to Ma & Fry 2002 and appendix A in Park et al. 2013 for details. In general, it is common to write the angular power spectrum in the dimensionless form: D ℓ = C ℓ ( ℓ +1) ℓ/ (2 π ).  \nRunning 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 Park et al. 2013 eq. (B1), and take advantage of linear theory to generate a correction term:  \nP miss q ⊥ ( k, z ) = ∫ k<k box d 3 k ' (2 π ) 3 (1 -µ 2 ) P χ (1+ δ ) ( | ⃗ k -⃗ k ' | ) P lin vv ( k ' ) (4)  \nwhere k box = 2 π/L box is the box-scale wavemode, µ = ˆ k · ˆ k ' , P χ (1+ δ ) ,χ (1+ δ ) ( k ) is the ionized matter power spectrum and P lin vv ( k ) = (˙ af/k ) P lin δδ ( k ) is the velocity power spectrum in the linear approximation taken from the public code CAMB (Lewis et al. 2000).  \nIt was found by Alvarez (2016) 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 nonGaussianity of reionization. So our calculations could be seen as conservative estimates of the power. As we will see in § 4.3, an increase in power could strengthen our conclusions. As such, we do not expect this missing term to affect our qualitative results.", '4. 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 ( § 4.1), observations of Ly α -emitting galaxies at z ≥ 8 ( § 4.2), and the patchy kSZ effect from reionization ( § 4.3). 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.', '4.1. QSO Observations at 5 < z < 6 4.1.1. 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 (Kulkarni et al. 2019; Keating et al. 2020; Nasir & D'Aloisio 2020; Bosman et al. 2022). We will study both in this section.  \nFigure 5 shows the mean transmission of the Ly α forest, ⟨ F Ly α ⟩ , at 5 ≤ z ≤ 6 in our models compared with  \nFigure 5. Ly α forest mean transmission, compared with recent observations (see text). The early start/early end model is in severe disagreement with measurements, predicting transmission near unity at all redshifts. The other two models agree with the measurements of Bosman et al. 2022 by construction, since their ˙ N γ ( z ) histories were calibrated to match those measurements (see discussion in § 2.1). The fact that such a calibration is possible for both late-ending models indicates that ⟨ F Ly α ⟩ alone cannot distinguish between them.  \n<!-- image -->  \nrecent measurements from Becker & Bolton 2013; Eilers et al. 2018; Bosman et al. 2018; Yang et al. 2020; Bosman et al. 2022. The late start/late end and early start/late end models both agree well with the data this is by design, since both ˙ N γ histories were calibrated so that the simulations would match the Bosman et al. 2022 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 (D'Aloisio et al. 2018; Gaikwad et al. 2023). 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.  \nFigure 6 shows the cumulative distribution function (CDF) of forest effective optical depths measured over intervals of 50 h -1 Mpc, P ( < τ 50 eff ). The τ 50 eff 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τ 50 eff sight-  \nFigure 6. Distribution of effective Ly α optical depths over 50 h -1 Mpc segments of the Ly α forest, P ( < τ 50 eff ), compared to measurements from Bosman et al. 2022. We show results and measurements at z = 5, 5 . 2, 5 . 6, and 6 (see panel titles). For the early start/early end model, we have re-scaled the Ly α opacities along each all sightlines by a constant value such that the mean transmission matches the Bosman et al. 2022 measurements. This allows us to make a clean comparison of how the IGM neutral fraction affects the shape of P ( < τ 50 eff ) at fixed ⟨ F Ly α ⟩ . At z = 6, the re-scaled early start/early end model produces too little scatter in τ 50 eff , the late start/late end case produces too much, and the early start/late end model is a good match. At z = 5 . 6, none of the models match the observations particularly well. The observations lie between the blue dotted and red dashed curves, suggesting a non-zero neutral fraction smaller than that of the early start/late end model (7 . 5%). At z = 5 . 2, the re-scaled early start/early end matches the data best, suggesting that reionization should end slightly earlier than it does in both late-ending models. All models match the data at z = 5, when reionization is complete.  \n<!-- image -->  \nz < 6 in late reionization scenarios (Kulkarni et al. 2019; Nasir & D'Aloisio 2020; Qin et al. 2021). We show P ( < τ 50 eff ) at z = 5, 5 . 2, 5 . 6, and 6 for each of our models, alongside measurements from Bosman et al. 2022. To show what P ( < τ 50 eff ) 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 10 model such that ⟨ F Lyα ⟩ agrees with the measurements shown in Figure 5. This allows us to cleanly compare how the shape of P ( < τ 50 eff ) is affected by the IGM neutral fraction at fixed ⟨ F Ly α ⟩ .  \nAt z = 6 (lower right), the early start/late end model best matches the Bosman et al. 2022 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 τ 50 eff , echoing the conclusions of Bosman et al. 2022. 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 lateending models produce too much scatter in τ 50 eff , 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  \nmatch 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 ( < τ 50 eff ). 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.  \nAt face value, observations of P ( < τ 50 eff ) 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 ( < τ 50 eff ). However, none of the scenarios in Figure 6 match the observed P ( < τ 50 eff ) 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 Bosman et al. 2022, 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 A, we estimate what P ( < τ 50 eff ) 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 ( < τ 50 eff ), 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 ❡❛/a114❧② /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model remains preferred.  \nIn Cain et al. 2024, 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 underestimate of the mean transmission at fixed x HI (Doughty et al. 2023), resulting in a spuriously early end to reionization (by ∆ z ≈ 0 . 2, see Fig. 11 of Cain et al. 2024) when calibrating to measurements. Our forest calculations include the resolution correction prescribed Appendix A of D'Aloisio et al. 2018, so in principle they account for this effect. However, those corrections were derived in a 25 h -1 Mpc 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 Cain et al. 2024 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 (Chardin et al. 2015; Madau et al. 2024). Any of these effects could affect when reionization needs to end to match ⟨ F Ly α ⟩ measurements, and the shape of P ( < τ 50 eff ) at fixed neutral fraction, poten-  \ninterpretation of P ( < τ 50 eff ) measurements.", '4.1.2. 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 (Emberson et al. 2013; Park et al. 2016; D'Aloisio et al. 2020; Chan et al. 2024). We calculate the Lyman limit MFP in our simulations using the definition in Appendix C of Chardin et al. 2015,  \nλ 912 mfp = ⟨ ∫ xdf ⟩ ⟨ ∫ df ⟩ = -〈∫ 0 1 xdf 〉 (5)  \nwhere x is the position along a randomly oriented sightline, f ( x ) is the transmission of 912 ˚ A photons, and the angle brackets denote an average over many sightlines. Roth et al. 2024 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 (Lewis et al. 2022).  \nThe left panel of Figure 7 shows the MFP in our late-ending models, compared to measurements from Worseck et al. 2014; Becker et al. 2021; Zhu et al. 2023; Gaikwad et al. 2023. Both scenarios are in broad agreement with the measurements. However, the direct measurements using QSO Lyman Continuum (LyC) spectra (all but the Gaikwad et al. 2023 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 λ 912 mfp 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 Becker et al. 2021; Zhu et al. 2023. Both models are within 1 σ of the indirect, Ly α forest-based measurements from Gaikwad et al. 2023 (see also Davies et al. (2024b)).  \nThis result is consistent with previous findings that ongoing reionization at z = 6 is needed to explain the direct QSO measurements (Cain et al. 2021; Lewis et al. 2022; Garaldi et al. 2022; Lewis et al. 2022; Satyavolu et al. 2024; Roth et al. 2024). 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 (Park et al. 2016; D'Aloisio et al. 2020). These factors result in direct MFP measurements preferring the ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ scenario. Our earlier result that  \nFigure 7. Other (non-Ly α forest) QSO-based measurements of global IGM properties at 5 < z < 6 in our late-ending models. Left: the mean free path to ionizing photons at 912 ˚ A compared to recent measurements. We show direct measurements from the Lyman Continuum spectra of QSOs (Worseck et al. 2014; Becker et al. 2021; Zhu et al. 2023), and the recent indirect measurements of Gaikwad et al. 2023 based on the Ly α forest. Both models are consistent with the forest-based measurements from Gaikwad et al. 2023, but the direct QSO-based measurements (particularly at z = 6) prefer the late start/late end scenario. The faster redshift evolution of λ mfp in that model, driven by the evolving neutral fraction, is in better agreement with the z = 6 direct measurements. Right: Average temperature at mean density ( T 0 ), compared to measurements from Becker et al. 2011; Boera et al. 2019; Walther et al. 2019; Gaikwad et al. 2020. Both models display reionization-driven temperature peaks at z ∼ 5 . 3, consistent with the Gaikwad et al. 2020 measurements. The late start/late end model peaks at a higher temperature, since a larger fraction of the IGM has been recently heated by fast-moving I-fronts (see annotation). The early start/late end model is more consistent with measurements at z ≥ 5, indicating a mild preference for that scenario (but with some caveats - see text).  \n<!-- image -->  \nP ( < τ 50 eff ) 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 Gaikwad et al. 2023 (based on P ( < τ 50 eff ) 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.", '4.1.3. IGM Thermal History': "The IGM temperature at mean density, T 0 , is shown in the right panel of Figure 7, alongside measurements from Becker et al. 2011; Boera et al. 2019; Walther et al. 2019; Gaikwad et al. 2020. Both late-ending models display a 'bump' in temperature at z ≈ 5 . 3 due to the end of reionization. Heating by I-fronts (D'Aloisio et al. 2019; Zeng & Hirata 2021) 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 (McQuinn & Upton Sanderbeck 2016). The peak in T 0 is higher in the late start/late end model because a larger fraction of the IGM is reionized at z < 6. The redshift of the bump suggested by the Gaikwad et al. 2020 measurements is closer to z ∼ 5 . 6 - this is consistent with our earlier finding (based  \non P ( < τ 50 eff )) that reionization may end ∆ z ∼ 0 . 3 too late in our models.  \nAt face value, the early start/late end model agrees best with T 0 measurements. Although both models are consistent with the Gaikwad et al. 2020 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 Villasenor et al. 2022, 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 Ifront temperature ( T reion ) and photo-heating in ionized gas afterwards (D'Aloisio et al. 2019). 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 Asthana et al. 2024). This could bring the late start/late  \nend model into agreement with z ≤ 5 T 0 measurements. However, this would also require a later reionization history at fixed ⟨ F Ly α ⟩ (Figure 5 of Cain et al. 2024), which would worsen the disagreement with the observed P ( < τ 50 eff ) in Figure 6. As such, we conclude that measurements of T 0 mildly prefer the ❡❛/a114❧② /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model.", '4.1.4. 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 8 compares our models with constraints from Ly α forest dark pixels (McGreer et al. 2015; Jin et al. 2023), dark gaps (Zhu et al. 2022), QSO damping wings (Greig et al. 2024), P ( < τ 50 eff ) (Choudhury et al. 2021; Gaikwad et al. 2023), and Ly α forest damping wings (Zhu et al. 2024; Spina et al. 2024). The bold lines show the reionization histories in our two lateending models, while the faded lines show these shifted to the right by ∆ z = 0 . 3, consistent with the discussion surrounding P ( < τ 50 eff ) in § 4.1.1.  \nMost 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 Zhu et al. 2022 and the recent QSO damping wing limits from Greig et al. 2024 disfavor the late start/late end model. The recent lower limit on the neutral fraction from Zhu et al. 2024, 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 Spina et al. 2024 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 ( < τ 50 eff ) 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 Zhu et al. 2024 damping wing limit. We conclude that neutral fraction constraints at z < 6 . 5 mildly prefer the ❡❛/a114❧② /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model, but that relatively small, realistic shifts in the reionization history could change this conclusion.', '4.2. 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 (Bunker et al. 2023). Our goal is to determine whether these observations prefer an early or late start to reionization.  \nFigure 8. Constraints on the neutral fraction at z ≤ 6 . 5 from dark pixels, dark gaps, damping wings, P ( < τ 50 eff ), and forest damping wings (see text for references). The bold curves show our late-ending reionization histories, and the faded curves show these shifted to the right by ∆ z = 0 . 3 (see § 4.1.1). The upper limit at z = 5 . 8 from Zhu et al. 2024 from Ly α forest damping wings is inconsistent with an earlier ending to reionization, and the same can be said of the z = 5 . 6 measurement from Spina et al. 2024. Some of the lower limits on the ionized fraction are mildly inconsistent with the late start/late end model. However, shifting the histories earlier by ∆ z = 0 . 3 removes these tensions.  \n<!-- image -->', '4.2.1. Examples of IGM transmission at z = 8': 'In Figure 9, 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 § 3.3). The higher transmission at large v off owes to this object occupying a larger ionized region than the average galaxy.  \nThe 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,  \nFigure 9. Mean IGM Ly α transmission vs. velocity offset around galaxies with M UV < -17 at z = 8 in our lateending models. The thick curves show the mean transmission profiles, which go to 0 near systemic Ly α ( v off = 0, magenta line), and rise at v off > 0, displaying a shape similar to the damping wing absorption profile. The thin lines denote 20 profiles for individual sightlines surrounding the brightest galaxy (most massive halo) in our box. The early start/late end model has much more transmission on average than the late start/late end case, by a factor of 2 or more at v off < 500 km/s. This owes to its higher ionized fraction (see legend). The sightlines surrounding the brightest galaxy display more transmission at v off ≳ 200 -400 km/s than the average, and less at smaller v off owing to large-scale inflows surrounding that object. The higher transmission owes to the large, biased ionized bubble inhabited by the most massive halo in the box. This illustrates that even in models with low average Ly α transmission, individual sightlines can still be transmissive, especially towards the brightest galaxies.  \n<!-- image -->  \nit 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 (Smith et al. 2022; Perez et al. 2023). It also suggests that some LAE detections at z ≥ 8 could be explainable even if reionization starts relatively late.', '4.2.2. 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  \nEW obs = ⟨ T IGM ⟩ line EW int (6)  \nAn object is detectable if EW obs is greater than some threshold EW min obs . This condition can be expressed as  \nEW obs EW int = ⟨ T IGM ⟩ line > EW min obs EW int ≡ T thresh (7)  \nwhere 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 11 . 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.  \nRecently, Asthana et al. 2024 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 (Weinberger et al. 2018; Mason & Gronke 2020). 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 10, 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 9). 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.  \nThe 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 be-  \nFigure 10. Evolution of LAE visibility with redshift for several combinations of v off (in km/s), T thresh , and M UV range - denoted as ( v off , T thresh , M UV range) in the legend. The left and right panels show the late start/late end and early start/late end models, respectively. The solid, dot-dashed, and dashed curves show visibilities for bright galaxies with M UV < -21 - these three curves have ( T thresh , v off ) = (400 km/s , 0 . 2), (200 km/s , 0 . 2), and (200 km/s , 0 . 5), respectively. The dotted and double-dot dashed curves show faint galaxies ( -19 < M UV < -17), with ( T thresh , v off ) = (400 km/s , 0 . 2) and (200 km/s , 0 . 5), respectively. Most choices of these parameters yield ≈ 50 -100% visibility in the early start/late end model, the only exception being faint LAEs with low v off and high T thresh (double-dot dashed). By contrast, the late start/late end model displays a wide range of expected LAE visibility statistics for different types of LAEs. In the left panel, we show recent measurements of the LAE fraction from Tang et al. 2024. The magenta points show the fraction of galaxies observed to have Ly α EW > 25 ˚ A. Interpreting this as a fraction of LAEs that are visible means assuming X int Ly α = 1, so these points are shown as lower limits. The green points assume X int Ly α = 0 . 3, as measured at z = 5 by Tang et al. 2024. The low measured visibility fractions are consistent with range of visibility statistics found in the late start/late end case. However, in the early start/late end model, only the most 'restrictive' parameter combination that we show T thresh = 0 . 5, v off = 200 km/s, and -17 < M UV < -19 displays similar statistics.  \n<!-- image -->  \ncause 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 (Park et al. 2021).  \nThe 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 doubledot 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 ob-  \ncts 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.  \nIn the left panel, we show recent measurements of the fraction of galaxies hosting Ly α emitters, (the Ly α fraction, X obs Ly α ), at z = 8 . 7 and z = 10 from Tang et al. 2024. The purple points show X obs Ly α 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 int Ly α = 1. For this reason, we display these points as lower limits on X obs Ly α /X int Ly α . The green points show X obs Ly α /X int Ly α assuming X int Ly α = 0 . 3, the Ly α fraction measured at z = 5 by Tang et al. 2024, 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 ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model.  \nFor the global LAE visibility fraction to evolve consistently with measurements from Tang et al. 2024 in the early start/late end model, the double-dot dashed curve in Figure 10 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 (Mason et al. 2018), they also tend to have fairly high EW intr (Dijkstra & Wyithe 2012; Tang et al. 2024). A majority of the LAEs observed at z = 5 by Tang et al. 2024 in that M UV range have EW intr > 50 ˚ A, and about half have EW intr > 100 ˚ A. With the visibility threshold of EW min obs = 25 ˚ A used in Tang et al. 2024, 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 Tang et al. 2024 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 (see also Mason et al. 2018) 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 Tang et al. 2024 measurements.", '4.2.3. 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 ˚ A. Using a Bayesian analysis based on reionization simulations and an empirically derived model for the intrinsic EW distribution of LAEs from Mason et al. 2018, Bruton et al. 2023 inferred that the IGM must be at least 12% ionized at 2 σ confidence (yellow point in Figure 12). This constraint, at face value, clearly favors the early start/late end model (see Figure 12 in the next section). Here, we consider whether the observed properties of GN-z11 require reionization to start early.  \nWe 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 intr off , 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,  \nFigure 11. Distribution of EW intr required to produce observed an LAE with the properties of GN-z11 at z = 10 . 6 (see text for details). The shaded region denotes EW intr < 60 ˚ A, the range observed for similarly bright galaxies at low redshift (see text for references). We see that nearly all sightlines to bright galaxies in the early start/late end model have EW intr within this shaded region, and with 20% (95%) of sightlines requiring EW intr < 30 ˚ A (60 ˚ A). By contrast, the late start/late end model has a wider distribution centered at higher EW intr . About 12% of sightlines in this model require EW intr < 60 ˚ A though, suggesting that a particularly bright LAE with relatively high EW viewed along a sightline with properties that occur ≈ 10% of the time could produce an observation like GN-z11 in the late start/late end scenario.  \n<!-- image -->  \ngives 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 11 for our late-ending models. The shaded region denotes the range of EW observed in similarly bright galaxies at lower redshifts Endsley et al. 2022; Tang et al. 2023; Saxena et al. 2024; Tang et al. 2024, the highest of which is ≈ 60 ˚ A (Fig. 1 of Tang et al. 2024).  \nWe 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 ˚ A, and about 20% require EW intr < 30 ˚ A. 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 ˚ A, and none of them allow < 30 ˚ A. 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 Tang et al. 2024 with the highest EWs ( ≈ 30 and 60 ˚ A) were both detected in H α . Based off a  \nFigure 12. Summary of neutral fraction measurements in the literature based on either Ly α emission or Ly α damping wings in galaxy spectra. Most of the measurements are at 7 < z < 8, and taken together they do not strongly prefer either late-ending model over the other. There is a dearth of constraints at 8 . 5 < z < 10 (the exception being the Tang et al. 2024 point at z ∼ 8 . 7), and those at z > 10 do not show a clear consensus either. The only measurement that shows a strong preference for the early start/late end model is the constraint from Bruton et al. 2023, based on GN-z11 (see § 4.2.3). Measurements are from Mason et al. 2018; Mason et al. 2019; Hoag et al. 2019; Whitler et al. 2020; Jung et al. 2020; Morales et al. 2021; Bolan et al. 2022; Wold et al. 2022; Morishita et al. 2023; Nakane et al. 2024; Bruton et al. 2023; Hsiao et al. 2023; Curtis-Lake et al. 2023; Tang et al. 2024.  \n<!-- image -->  \nclear detection of H γ emission in GN-z11, Bunker et al. 2023 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.', '4.2.4. 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 12 shows a collection of these measurements and limits compared to our lateending reionization models in the same format as Figure 8 (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.  \nUnlike in Figure 8, 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 Bruton et al. 2023, 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 Tang et al. 2024 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.', '4.3. 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 § 3.4) in the left panel of Fig. 13, along with the 1 σ error bar and 95% confidence upper limits from Reichardt et al. 2021. 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 12 . We also include the revised 2 σ upper limit from Gorce et al. 2022, which is somewhat lower than the SPT result and favors the late start/late end model even more.  \nTo 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  \nFigure 13. Patchy kSZ from reionization for all three of our model. Left: patchy kSZ power vs. ℓ , compared to the recent measurement and 2 σ upper limit at ℓ = 3000 from Reichardt et al. 2021. We also include the 2 σ upper limit from the re-analysis of the SPT data by Gorce et al. 2022 (offset for readability), which disfavors the early start/late end model. Only the late start/late end model falls within the 1 σ errors, and the early start/late end case has the most power and lies close to the 2 σ upper limit from Reichardt et al. 2021. Right: differential contribution to the total power at ℓ = 3000 as a function of redshift. The shaded red region shows the range where the universe is less than 20% ionized in the early start/late end model, from which nearly half the ℓ = 3000 power originates. The early start/early end model has slightly less power than the early start/late end case because it ends reionization earlier and thus has a shorter duration (see annotation).  \n<!-- image -->  \nvelocity 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.  \nIn the late start/late end case, reionization begins much later but still ends at z = 5, which makes the peak in dD pkSZ 3000 /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 (Battaglia et al. 2013; Chen et al. 2023). We see also that although the Reichardt et al. 2021 measurement does not rule out the early start/late end model at 2 σ , it clearly prefers the ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ case. This finding is consistent with the recent limits on the duration of reionization from Raghunathan et al. 2024 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  \n< 4 . 5 at 95% confidence 13 - our early start/late end model has ∆ z 50 = 3 . 1.", "5.1. '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 ( § 2), inferences from the spectra of high-redshift QSOs ( § 4.1), Ly α transmission from z > 8 galaxies ( § 4.2), and constraints from the CMB ( § 4.3). We concluded in § 4.1 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 14 . This motivates the key question in this work:  \nTable 1. Summary of the 'face value' preferences of each observable we studied for a late or early start to reionization. The left-most column gives the category for each observable - whether it derives from the CMB (blue), high-redshift galaxies (red), or z < 6 . 5 QSOs (magenta). The second column lists the observables, grouped by category. The right two columns give the preference of each probe for a late or early start to reionization. Both columns list 'no preference' if an observable does not clearly lean towards either scenario. Importantly, neither scenario is 'confirmed' or 'ruled out' by any observables. We see that not all the probes prefer the same scenario. Three of the observables derived from QSO spectra ( P ( < τ 50 eff ), T 0 , and x HI ( z < 6 . 5), see § 4.1) prefer an early start, while the remaining probes all either prefer a late start or cannot distinguish the two scenarios. All three categories, which use largely independent data sets with vastly different analysis techniques, have at least one probe that prefers a late start. However, all the probes that prefer an early start derive from the same type of data (QSO spectra). These in particular are challenging to interpret due to uncertainties in modeling the z < 6 . 5 reionization history and connecting this to what happens at higher redshifts. As such, we conclude that observations display a mild preference for the ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ scenario.  \nwhen 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.  \nWe summarize our findings in Table 1. The leftmost 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 ( < τ 50 eff ), the IGM thermal history, and x HI measurements at z < 6 . 5 prefer the early start/late end model.  \nA 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 15 . As we explained in § 4.1, 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  \na late start. As such, we judge that the consensus of these probes from a wide range of data sets indicates that the ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ model is mildly preferred (overall) by observations.  \nThe 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.  \nHowever, 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.", '5.2. 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. Atek et al. 2024a 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 3). Continued efforts to understand how f esc correlates with galaxy properties at low redshift, such as the Low-redshift Lyman Continuum Survey (LzLCS, Chisholm et al. 2022; Flury et al. 2022; Jaskot et al. 2024a,b) will also be crucial for placing rea-  \nnable limits on the evolution of f esc (see also e.g. Smith et al. 2020; Pahl et al. 2021; Wang et al. 2023).  \nThere is also further progress to be made with QSObased 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 7 shows. Forthcoming observations with Euclid (Atek et al. 2024b) 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 (Christenson et al. 2021; Ishimoto et al. 2022; Christenson et al. 2023; Kashino et al. 2023) may also help tighten constraints on the reionization history at z < 6 . 5 (Garaldi et al. 2022; Gangolli et al. 2024).  \nFurther 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 (Wold et al. 2024) will also reveal bright LAEs over a much wider area than JWST, enabling improved constraints on the early reionization history (Perez et al. 2023).  \nForthcoming improvements on CMB constraints from multiple experiments, including the Atacama Cosmology Telescope (ACT, Hlozek et al. 2012), SPT (Raghunathan et al. 2024), Simons Observatory (Bhimani et al. 2024), and CMB-S4 (Alvarez et al. 2021) will improve constraints on τ es and pkSZ. They will may also detect new signals that probe reionization, such as patchy τ (Coulton et al. 2024) and higher-order statistics and cross-correlations with other signals (e.g. La Plante et al. 2020). These will help constrain the early stages of reionization because of their sensitivity to its duration and morphology (Chen et al. 2023).', '6. 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.  \n- · 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 Planck Collaboration et al. 2020; de Belsunce et al. 2021. 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.\n- · 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 earlystarting 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 Tang et al. 2024 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 σ .\n- · 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.\n- · 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  \nacross 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 ❧❛/a116❡ /a115/a116❛/a114/a116✴❧❛/a116❡ ❡♥❞ scenario.  \n- · 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.  \nForthcoming 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.  \nCC 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 JWSTAR02608.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 TGPHY230158.", 'REFERENCES': "- Sunyaev, R. A., & Zeldovich, Y. B. 1980, Mon. Not. R. Astron. Soc., 190, 413\n- Tang, M., Stark, D. P., Topping, M. W., Mason, C., & Ellis, R. S. 2024, arXiv e-prints, arXiv:2408.01507, doi: 10.48550/arXiv.2408.01507\n- Tang, M., Stark, D. P., Chen, Z., et al. 2023, MNRAS, 526, 1657, doi: 10.1093/mnras/stad2763\n- Tang, M., Stark, D. P., Ellis, R. S., et al. 2024, Monthly Notices of the Royal Astronomical Society, 531, 2701, doi: 10.1093/mnras/stae1338\n- Tepper-Garc'ıa, T. 2006, MNRAS, 369, 2025, doi: 10.1111/j.1365-2966.2006.10450.x\n- Topping, M. W., Stark, D. P., Endsley, R., et al. 2024, MNRAS, 529, 4087, doi: 10.1093/mnras/stae800\n- Trac, H., Cen, R., & Mansfield, P. 2015, ApJ, 813, 54, doi: 10.1088/0004-637X/813/1/54\n- Trebitsch, M., Hutter, A., Dayal, P., et al. 2023, MNRAS, 518, 3576, doi: 10.1093/mnras/stac2138\n- Trebitsch, M., Volonteri, M., Dubois, Y., & Madau, P. 2018, MNRAS, 478, 5607, doi: 10.1093/mnras/sty1406\n- Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397, doi: 10.1051/0004-6361:20065554\n- Villasenor, B., Robertson, B., Madau, P., & Schneider, E. 2022, ApJ, 933, 59, doi: 10.3847/1538-4357/ac704e\n- Walther, M., O˜norbe, J., Hennawi, J. F., & Luki'c, Z. 2019, ApJ, 872, 13, doi: 10.3847/1538-4357/aafad1\n- Wang, F., Davies, F. B., Yang, J., et al. 2020, ApJ, 896, 23, doi: 10.3847/1538-4357/ab8c45\n- Wang, X., Teplitz, H. I., Smith, B. M., et al. 2023, arXiv e-prints, arXiv:2308.09064,  \ndoi: 10.48550/arXiv.2308.09064  \n- Weinberger, L. H., Kulkarni, G., Haehnelt, M. G., Choudhury, T. R., & Puchwein, E. 2018, Monthly Notices of the Royal Astronomical Society, 479, 2564, doi: 10.1093/mnras/sty1563  \nWhitler, L. R., Mason, C. A., Ren, K., et al. 2020,  \nMNRAS, 495, 3602, doi: 10.1093/mnras/staa1178  \n- Wold, I. G. B., Malhotra, S., Rhoads, J. E., Tilvi, V., & Gabrielpillai, A. 2024, AJ, 167, 157,  \ndoi: 10.3847/1538-3881/ad2adf  \n- Wold, I. G. B., Malhotra, S., Rhoads, J., et al. 2022, ApJ, 927, 36, doi: 10.3847/1538-4357/ac4997\n- Worseck, G., Prochaska, J. X., O'Meara, J. M., et al. 2014, MNRAS, 445, 1745, doi: 10.1093/mnras/stu1827  \nWu, X., Kannan, R., Marinacci, F., Vogelsberger, M., &  \nHernquist, L. 2019, MNRAS, 488, 419,  \ndoi: 10.1093/mnras/stz1726  \n- Yang, J., Wang, F., Fan, X., et al. 2020, ApJL, 897, L14, doi: 10.3847/2041-8213/ab9c26\n- Yeh, J. Y. C., Smith, A., Kannan, R., et al. 2023, MNRAS, 520, 2757, doi: 10.1093/mnras/stad210\n- Zeng, C., & Hirata, C. M. 2021, ApJ, 906, 124, doi: 10.3847/1538-4357/abca38\n- Zhao, J., & Furlanetto, S. R. 2024, arXiv e-prints, arXiv:2401.07893, doi: 10.48550/arXiv.2401.07893\n- Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2022, ApJ, 932, 76, doi: 10.3847/1538-4357/ac6e60\n- Zhu, Y., Becker, G. D., Christenson, H. M., et al. 2023, ApJ, 955, 115, doi: 10.3847/1538-4357/aceef4\n- Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2024, MNRAS, doi: 10.1093/mnrasl/slae061\n- Zitrin, A., Labb'e, I., Belli, S., et al. 2015, ApJL, 810, L12, doi: 10.1088/2041-8205/810/1/L12  \nFigure 14. Ly α forest P ( < τ 50 eff ) with the reionization histories in our late-ending models shifted earlier by ∆ z = 0 . 3, as described in § 4.1.1. This test estimates what these models would look like if reionization ended at z = 5 . 3 (as suggested by the analysis of Bosman et al. 2022) instead of at z = 5. At z = 5, 5 . 2, and 5 . 4, there is little difference between the two models, as reionization is complete (or nearly so) in both cases. At z = 5 . 6 and 5 . 8, P ( < τ 50 eff ) is narrower in the early start/late end model and agrees well with the observations, whereas P ( < τ 50 eff ) is still too wide in the late start/late end case. At z = 6, it is difficult to tell which model agrees best with the data.  \n<!-- image -->", 'A. P ( < τ 50 EFF ) FOR SHIFTED REIONIZATION HISTORIES': 'In Figure 14, we show P ( < τ 50 eff ) 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 § 4.1.1. We show 6 redshifts between z = 5 and 6 in intervals of 0 . 2. This exercise shows roughly what P ( < τ 50 eff ) 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 ( < τ 50 eff ) is slightly narrower in the early start/late end case, as it is in Figure 6. 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 ( < τ 50 eff ) 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 ( < τ 50 eff ).', 'B. 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 § 4.2.2. Figure 15 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%  \nFigure 15. Summary of LAE transmission statistics in the late start/late end model. The y axis is T thresh (Eq. 7) and the x axis is the v off at line center. The color plot shows the fraction of visible LAEs (with T IGM > T thresh ) across this parameter space. The panels denote different redshifts (columns) and ranges of M UV (rows). The white contour lines denote 'visibility fractions' of 10%, 25%, and 50%. The hatched region in the upper left of each panel denotes T thresh > 0 . 2 and v off < 500 km/s.  \n<!-- image -->  \n(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 ˚ A (Goovaerts et al. 2023; Tang et al. 2024). For T thresh = 0 . 2, the latter would correspond to EW min obs < 100 ˚ A.  \nAt 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.  \nAt 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 10, this drop-off in visibility is consistent with the observed decline in X Ly α observed by Tang et al. 2024 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, Bunker et al. 2023) has v off = 550 km/s, meeting this condition. The recently observed JADES-GS-z9-0 ( M UV = -20 . 43, z = 9 . 4, Curti et al. 2024) 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 ˚ A, ≈ 4 ×  \nFigure 16. Same as Figure 15, but for the early start/late end model. A significant fraction of LAEs with typical properties (even faint ones) should remain visible in this scenario even at z = 11, in contrast with the observations Tang et al. 2024.  \n<!-- image -->  \nlarger than their observed EWs (see § 4.2.3). 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.  \nFigure 16 is the same as Figure 15, 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."}