Chronicling the reionization history at 6 ≲ z ≲ 7 with emergent quasar damping wings

The spectra of high-redshift ( z ≳ 6) quasars contain valuable information on the progression of the Epoch of Reionization (EoR). At redshifts z < 6, the observed Lyman-series forest shows that the intergalactic medium (IGM) is nearly ionized, while at z > 7 the observed quasar damping wings indicate high neutral gas fractions. However, there remains a gap in neutral gas fraction constraints at 6 ≲ z ≲ 7 where the Lyman series forest becomes saturated but damping wings have yet to fully emerge. In this work, we use a sample of 18 quasar spectra at redshifts 6 . 0 < z < 7 . 1 to close this gap. We apply neural networks to reconstruct the quasars’ continuum emission around the partially absorbed Lyman α line to normalize their spectra, and stack these continuum-normalized spectra in three redshift bins. To increase the robustness of our results, we compare the stacks to a grid of models from two hydrodynamical simulations, ATON and CROC, and we measure the volume-averaged neutral gas fraction, ¯ x HI , while jointly fitting for the mean quasar lifetime, t Q , for each stacked spectrum. We chronicle the evolution of neutral gas fraction using the ATON (CROC) models as follows: ¯ x HI = 0 . 21 +0 . 17 − 0 . 07 (¯ x HI = 0 . 10 +0 . 73 < 10 − 4 ) at ⟨ z ⟩ = 6 . 10, ¯ x HI = 0 . 21 +0 . 33 − 0 . 07 (¯ x HI = 0 . 57 +0 . 26 − 0 . 47 ) at ⟨ z ⟩ = 6 . 46, and ¯ x HI = 0 . 37 +0 . 17 − 0 . 17 (¯ x HI = 0 . 57 +0 . 26 − 0 . 21 ) at ⟨ z ⟩ = 6 . 87. At the same time we constrain the average quasar lifetime to be t Q ≲ 7 Myr across all redshift bins, in good agreement with previous studies.


INTRODUCTION
The Epoch of Reionization (EoR) represents the last major phase transition of our Universe, where the hydrogen in the intergalactic medium (IGM) transitions from a completely neutral state to the ionized state we observe today.The past decade has brought crucial insights into the ionizing sources (e.g.Robertson (2022)), as well as the morphology and timing (e.g.Bosman et al. (2022)) of the reionization process, all of which are important to understand open questions about galaxy formation, cosmology, and the evolution of the earliest supermassive black holes.Despite this progress, many details of the EoR remain unknown.
Corresponding author: Dominika Ďurovčíková dominika@mit.edu The timing of the EoR in particular has been studied extensively.For instance, the reionization optical depth inferred from the cosmic microwave background (CMB) favours a late and fast reionization and yields a midpoint of z re ≳ 7 (Planck Collaboration et al. 2020).At redshifts z ≲ 6, the observed Lyman α and β forests in quasar spectra show that the IGM is already close to fully ionized (Fan et al. 2006;McGreer et al. 2015;Bosman et al. 2018;Eilers et al. 2018aEilers et al. , 2019;;Zhu et al. 2022;Bosman et al. 2022).However, the Lyα absorption saturates already at relatively low volume averaged neutral gas fractions of xHI ∼ 10 −4 , leading to complete flux absorption and the formation of a Gunn-Peterson through (Gunn & Peterson 1965).On the other hand, high neutral gas fractions at redshifts around z ≳ 7 can be probed via off-resonant (damped) absorption features in quasar spectra, i.e. damping wings (Miralda-Escudé 1998) imprinted on the red side of the Lyα emission line.
The observed strengths of the damping wings indicate high neutral gas fractions of xHI ≳ 30% but are limited to a few, high-redshift quasar sightlines (Bañados et al. 2018;Greig et al. 2017Greig et al. , 2019Greig et al. , 2022;;Davies et al. 2018a;Ďurovčíková et al. 2020;Yang et al. 2020a;Wang et al. 2020;Mesinger & Haiman 2007;Schroeder et al. 2013), as uncertainties in the predictions of quasar continua prevent us from securely modeling weaker damping wings along individual sightlines.This leaves us with a range of neutral gas fraction spanning a few orders of magnitude that is challenging to probe with either the Lyman-series forest or the damping wing signature.Thus, there remains a gap in xHI constraints at intermediate neutral gas fractions (corresponding roughly to 6 ≲ z ≲ 7) where most quasars do not show strong damping wings.
Current upper bounds on the neutral hydrogen fraction at 6 ≲ z ≲ 7 come from the Lyα and Lyβ forest dark pixel counts (Jin et al. 2023;Zhu et al. 2023) and are consistent with constraints from the CMB.In addition, galaxy-based constraints at z ≳ 6.5 are beginning to emerge in the literature.Lyman Break galaxies (LBGs) (Mason et al. 2018(Mason et al. , 2019)), Lyman α emitters (LAEs) (Ouchi et al. 2010;Sobacchi & Mesinger 2015) as well as stacked galaxy damping wings (Umeda et al. 2023) all point towards a significant neutral gas fraction, xHI ≳ 40%, at z ≳ 6.5.
This work aims to fill this gap in the history of the EoR.We do so by tackling two main challenges.First, reionization is known to be patchy.This, combined with the insensitivity to many orders of magnitude in xHI , renders single-sightline quasar studies challenging near the end of the EoR, where cosmic variance is simply too large to provide any useful information.Second, the quasar's radiation field affects the surrounding IGM.Quasar activity is accompanied by the release of radiation that ionizes the surrounding IGM over time, creating a region known as the proximity zone, thereby modifying the observed spectrum in the vicinity of the Lyα line.
Not only does the proximity effect lead to a degeneracy in neutral gas fraction constraints, the observed proximity zones also contain a wealth of information about the quasar's past activity.Proximity zone sizes have proved particularly useful in inferring the time period over which a given quasar has been UV luminous, known as the quasar lifetime.Measurements of quasar lifetimes from single-epoch spectra of highredshift quasars have yielded an average quasar lifetime of log 10 (t Q /yr) = 5.7 +0.5  −0.3 for a sample of 15 quasars at 5.8 ≤ z ≤ 6.6 (Morey et al. 2021), which, along with other lifetime measurements of individual objects (Eil-ers et al. 2018b(Eil-ers et al. , 2021;;Davies et al. 2019Davies et al. , 2020;;Andika et al. 2020), has challenged our models of supermassive black hole (SMBH) growth.These measurements have also revealed a population of surprisingly young quasars with billion-solar-mass black holes at z ≳ 6 (Eilers et al. 2021), which pose an even greater strain on our models of black hole growth, suggesting that either super-Eddington accretion or UV-obscured growth phases need to be evoked (Eilers et al. 2018b;Davies et al. 2019;Satyavolu et al. 2023).This imminent challenge calls for further investigation into population-level quasar lifetimes that is symbiotic with constraining the progress of reionization.
To this end, we use a sample of 18 spectra of quasars from the Epoch of Reionization to constrain the evolution of the neutral hydrogen gas fraction, xHI , at redshifts of 6.0 < z < 7.1.We overcome the challenges of cosmic variance and lack of sensitivity by stacking our quasar sample in discrete redshift intervals, thus fitting an average quasar spectrum in contrast to existing studies in the literature which focus on fitting single-sightline spectra.Additionally, we take into account the degeneracy between the neutral gas fraction and quasar radiation imprints by jointly fitting for an average quasar lifetime in each redshift bin.Our study showcases the emergence of damping wings near the end of the EoR and thus presents an important step towards measuring sightline-averaged constraints of the volumeaveraged neutral fraction, as well as average quasar lifetimes across cosmic time.
This paper is structured as follows.In § 2, we describe our observations, as well as the reconstruction method of the intrinsic emission around the Lyα line of our observed quasars, in order to obtain normalized flux transmission profiles.In § 3, we describe the simulations that we use to model our observations, and we provide the details about their inference in § 4. Subsequently, we present the neutral fraction and quasar lifetime constraints in § 5, and provide a short summary in § 6.

Quasar sample
This study uses a sample of 18 quasars at redshifts 6.0 < z < 7.1, shown in Fig. 1.The spectra for the majority of our target quasars were taken with the Magellan/FIRE spectrograph (Simcoe et al. 2013).For three of the quasars in our sample, i.e.J1030+0524, J159-02, and J1120+0641, we use combined Magellan/FIRE and VLT/X-Shooter (Vernet et al. 2011) spectra, while one other spectrum in our sample, i.e.J1148+5251, was observed with Keck/MOSFIRE (McLean et al. 2010(McLean et al. , 2012) ) and Keck/ESI (Sheinis et al. 2002).The total integration time along with other information and references is included in Table 1.Note that the highestredshift quasar in our sample is ULAS J1120+0641, which has been studied extensively (e.g.Simcoe et al. 2012;Bosman et al. 2017;Schindler et al. 2020) as it was the first z > 7 quasar ever discovered (Mortlock et al. 2011).The spectrum of ULAS J1120+0641 that we are using in this paper has been observed with Magellan/FIRE and VLT/X-Shooter, as compared to the VLT/FORS and Gemini/GNIRS spectrum that has been originally published for this object.All spectroscopic data were reduced using the PypeIt package1 (Prochaska et al. 2020;Prochaska et al. 2020).
In particular, we subtract flat fields that were taken during the same observing run and perform sky subtraction by differencing dithered A-B exposure pairs with a further sky line residual elimination according to Bochanski et al. (2009).The subtracted 2D images then undergo optimal extraction (Horne 1986) to obtain the 1D spectra.The wavelength solution is derived from the night sky OH lines.Subsequently, all extracted 1D spectra are flux calibrated by using sensitivity functions from standard stars.For each quasar, we coadd the flux-calibrated 1D spectra across exposures and correct for telluric absorption by jointly fitting an atmospheric model and a quasar model.
A bias in the inferred neutral gas fraction occurs if sightlines to the quasars in our sample intercept dense neutral hydrogen clouds, called proximate damped Lyα systems (pDLAs).Such high neutral hydrogen column density can cause damping wings to appear even in a fully ionized Universe, thus mimicking the effect of a highly neutral IGM.All quasars in this sample were hence chosen to be uncontaminated by known pDLAs -for example, J183+05 had to be removed from our sample (Bañados et al. 2019).
In addition, we imposed a luminosity cut of M 1450 = −26.5 to ensure good quality spectra.The distribution of our quasars in the magnitude-redshift space is shown in Fig. 2, along with a comparison of our sample to a recently published sample (Fan et al. 2023).Quasars in our sample are color-coded to show the three discrete redshift bins, 6.0 ≤ z < 6.3 (blue), 6.3 ≤ z < 6.7 (green), and 6.7 ≤ z ≤ 7.1 (red), that correspond to our final constraints.

Lyα continuum reconstruction
In order to constrain the volume-averaged neutral fraction and the quasar lifetime, we need to first estimate the intrinsic quasar continua around the Lyα line.To this end, we implement the neural network model called QSANNdRA, which was first introduced by Ďurovčíková et al. (2020).Around Lyα , this model has been shown to produce a minimal bias, 0.3%, and a scatter below 10% at 1σ, competitive with the best continuum reconstruction techniques in the literature (Bosman et al. 2021).
In this work, we use QSANNdRA to predict the intrinsic quasar continuum in the wavelength range 1170 Å − 1290 Å (the "blue side") based on the "red-side" spectrum at 1290 Å − 2900 Å, which is agnostic to absorption by the neutral hydrogen in the intervening IGM along the line-of-sight.
We train and test this algorithm on low-redshift quasar spectra from the sixteenth data release version of the SDSS Quasar Catalog (Lyke et al. 2020), and perform the same set of training data cuts, the same spectral smoothing procedure as well as the same network architecture and training hyperparameters as in the original work ( Ďurovčíková et al. 2020).In summary, we utilize the quasar catalog flags to exclude all broad-absorption line (BAL) quasars and quasars with known proximate damped Lyman absorbers (pDLAs) in their sightlines, and we impose a median SNR cut of ≥ 3.0.We choose the low-redshift spectra to have redshifts of 2.09 < z < 2.51 such that both the Lyα and Table 1.Quasar sample used in this study.F01 -Fan et al. (2001), F03 -Fan et al. (2003)  the Mg II emission lines fall into the spectral range of the BOSS and SDSS spectrographs.In order to recover smooth spectral continua without absorption lines, we smooth each spectrum with an updated version of the QSmooth algorithm2 ( Ďurovčíková et al. 2020) that implements piecewise spline fitting in the last smoothing step.We further clean up our training set with the Random Forest (Breiman 2001) rejection method described in Ďurovčíková et al. (2020).The aforementioned steps are implemented using the Scikit-Learn Python package (Pedregosa et al. 2011).
The remaining training set consists of 28,443 quasars, of which 80% are used for training and the remaining 20% are used for testing.QSANNdRA is constructed and trained using the Keras Python package (Chollet et al. 2015).The trained model achieves a mean absolute error of ≲ 14% across the whole prediction wavelength range.Moreover, we checked that this performance is nearly identical for the train and test data sets, which suggests good generalizability of this model to previously unseen spectroscopic data.
The trained model is subsequently applied to our highredshift quasar sample.First, we rebin the high-z spectra to a 100km/s resolution, which roughly matches the spectral resolution of the low-redshift SDSS sample, and we smooth each spectrum using QSmooth as described earlier.Since their rest-frame ultraviolet spectra are redshifted to the near-infrared, our high-z spectra all contain regions of enhanced telluric absorption due to the Earth's atmosphere which we mask before smoothing.Specifically, we mask out the following regions: 13, 000 Å ≲ λ obs ≲ 14, 500 Å, 18, 000 Å ≲ λ obs ≲ 19, 500 Å, and λ obs ≳ 22, 800 Å.We visually inspected each smoothed red-side spectrum, and for a number of quasars we manually adjusted smoothing parameters, the telluric mask, and masked additional absorption features that seemed to be biasing the QSmooth continuum.
Due to QSANNdRA's inability to work with missing input data, we fill in the telluric regions in each quasar spectrum with a mean spectral stack of its 100 nearest neighbours in the SDSS training sample.To  identify these nearest neighbours, we perform a simple Nearest Neighbor search using the Scikit-Learn package (Pedregosa et al. 2011) in a compressed, 10-component PCA space of their red-side spectrum.
After filling in the telluric gaps, we apply QSANNdRA to all high-z quasars in our sample and display the resultant continuum predictions in Fig. 3. Based on the input red-side spectrum (shown in red), the model outputs 100 individual samples of the continuum (shown in faint blue), as well as the mean prediction for each quasar (thick blue).The SDSS nearest-neighbor stack is shown in orange for reference, but is not used in the remainder of this work.Note that differences between the SDSS composites and the predictions from QSANNdRA (and other continuum reconstruction techniques) have been seen before (e.g.Bosman et al. (2021)) and are not surprising -the biggest differences in the Lyα continuum are seen for quasars whose CIV line at λ rest = 1545.86Å does not match the SDSS neighbor stack very well.This mismatch is thus a manifestation of the mild redshift evolution in the CIV line properties, as has been noted in Fan et al. (2023) (see also Meyer et al. (2019); Schindler et al. (2020)).Also note that we ultimately forward model the redshift uncertainty in our later analysis and thus the exact redshift that is used to bring the quasars into their rest frame should not affect our results.
We discuss and display these continuum predictions in more detail in Appendix A.

Stacking Quasar Spectra
With the continuum predictions at hand, we proceed by dividing our quasar sample into discrete redshift intervals and stacking all continuum-normalized spectra within each bin to compute an average transmitted quasar spectrum.Stacking has two main advantages: Firstly, it allows us to account for the large cosmic variance near the end of reionization and thus overcome the insensitivity to the intermediate range of neutral gas fractions, xHI .Note that all damping wing based neutral fraction constraints to date (Wang et al. 2021;Greig et al. 2017Greig et al. , 2019;;Ďurovčíková et al. 2020;Davies et al. 2018a;Bañados et al. 2018;Yang et al. 2020a) only sample single sightlines in the sky, which impedes their usability near the end of the EoR when the neutral gas fraction is very patchy.Secondly, stacking and averaging over multiple density fields helps with "gaussianizing" the noise in the spectrum, which allows us to approximate the likelihood function used for inference as a multivariate Gaussian distribution (see § 4).
Due to the aforementioned degeneracy between the neutral fraction and quasar radiation effects, we also obtain a sample-average constraint on the quasars' life- times.This stacking procedure is inspired by Morey et al. (2021), who previously used stacking to determine the mean quasar lifetime in a sample of 15 quasars at 5.8 ≤ z ≤ 6.6.Note that the interpretation of the inferred mean quasar lifetimes depends on the state of the IGM; for a fully ionized IGM, proximity zone sizes are sensitive to the episodic lifetime (i.e. the most recent "on" phase), whereas a mostly neutral IGM enables a measurement of the integrated quasar lifetime (Eilers et al. 2017).Therefore, due to stacking, we expect to be more sensitive to the episodic lifetime around z ∼ 6, and the integrated lifetime towards z ∼ 7, although the exact interpretation is complicated.
The resultant stacked spectra are shown in Fig. 4. We show spectral stacks in three redshift bins, 6.0 ≤ z < 6.3 (blue), 6.3 ≤ z < 6.7 (green), and 6.7 ≤ z ≤ 7.1 (red).These redshift intervals were chosen such that each bin contains approximately the same number of objects while avoiding bins of vastly unequal width.This is particularly relevant at the lower end of our redshift range, where we expect damping wings to have not yet fully emerged.Crucially, we use histogram binning with a resolution of 600 km/s to avoid correlating the uncertainties in neighbouring pixels and to smooth over the density fields in the vicinity of the quasars.
It should be noted that our analysis is still limited to a small sample of quasars along different sightlines.How-ever, by performing this stacking procedure we highlight the possibility of using damping wings of quasars also near the end of the EoR.

HYDRODYNAMICAL SIMULATIONS
In order to model both the effect of the neutral gas in the surrounding IGM and the effect of the quasar radiation field, we employ hydrodynamical simulations with radiative transfer.We use two different reionization simulations, namely the P-GADGET-3 simulations post-processed with the ATON code for radiative transfer (henceforth referred to as ATON, Satyavolu et al. ( 2023)), and the CROC simulations with fully coupled Optically Thin Variable Eddington Tensor (OTVET) radiative transfer (henceforth referred to as CROC, Chen & Gnedin (2021)).In both cases, we insert a mock quasar into the most massive halos in the simulation snapshots and draw multiple sightlines in order to simulate cosmic variance of the line-of-sight density and ionization fields.The following sections provide further details about the two simulations and a brief discussion of their main differences.Our use of two simulations is motivated by checking for consistency and increasing the robustness of our results, and we stress that the aim of this paper is not to assess the validity of either of them.

P-GADGET-3 simulations with ATON radiative transfer
The first suite of models comes from post-processing the P-GADGET-3 simulation (modified version of GADGET-2, Springel (2005)) using the ATON code (Aubert & Teyssier 2008, 2010) for three-dimensional radiative transfer (Kulkarni et al. 2019).ATON solves the radiative transfer using the M1 approximation (Aubert & Teyssier 2008;Gnedin & Abel 2001) and obtains the gas ionized fraction and temperature self-consistently.All our ionizing sources are placed in halos of mass ≳ 10 9 M ⊙ .The radiative transfer is run with a single photon frequency to reduce computational cost.The box size of this simulation is 160 cMpc/h with 20483 gas and dark matter particles and a resolution of 78.125 ckpc/h.The simulation is run between z = 99 and z = 4, and physical quantities such as the gas densities are saved in intervals of 40 Myr.As discussed in Kulkarni et al. (2019) and Keating et al. (2020), these simulations are calibrated to match the observed mean Lyα flux at z > 5 (Bosman et al. 2018;Becker et al. 2015) and are consistent with a number of high-redshift observations (Planck Collaboration et al. 2020;Keating et al. 2020;Becker et al. 2015;Greig et al. 2017Greig et al. , 2019;;Davies et al. 2018a;Wang et al. 2020;Weinberger et al. 2018Weinberger et al. , 2019;;Gaikwad et al. 2020).The end-point of reionization in this simulation is at z = 5.3, with the process half complete at z = 7.This picture remains consistent with the latest observations of the effective Lyα opacity by Bosman et al. (2022).
The creation of synthetic quasar spectra is described in Satyavolu et al. (2023).To save computational costs, we only use the simulation snapshots with neutral fractions of xHI = {3.7 × 10 −5 , 0.07, 0.13, 0.21, 0.37, 0.54, 0.75, 0.89} in this work.In each snapshot, quasars are placed in halos having masses in the range 10 11 M ⊙ ≲ M halo ≲ 10 12 M ⊙ 3 .The quasar light curve is assumed to be a 'light bulb'.For each of the quasars in our sample, the magnitude at 1450 Å (M 1450 ) is converted into the total number of ionizing photons Ṅtot by assuming the quasar to have a broken power-law spectral profile (Lusso et al. 2015) The simulated quasar lifetimes have been chosen to have the following values: log t Q = {4.0,5.0, 6.0, 7.0, 8.0} Myr.
For each quasar in our sample, we rescale the physical length scales in each simulation snapshot by a factor of (1+z)/(1+z sim ) and the densities along the line of sight by a factor of (1 + z) 3 /(1 + z sim ) 3 , where z is the quasar redshift and z sim is the redshift of the simulation snapshot.We then perform a line of sight radiative transfer using the method described in Satyavolu et al. (2023) and draw between 500 and 1000 sightlines in each snapshot.This method involves solving the thermochemistry equations to determine the abundances of ionized hydrogen and helium, alongside tracking the evolution of temperature.We use the post-processed neutral gas fraction, temperature, along with the rescaled gas densities and peculiar velocities from the underlying cosmological simulation to compute the Lyα optical depth by assuming a Voigt profile (García 2006).

CROC simulations with OTVET radiative transfer
The second suite of models is generated by postprocessing the Cosmic Reionization on Computers (CROC) simulations (Gnedin 2014).We utilize all six 40 cMpc/h CROC simulation runs to sample a wide range of large-scale structures.All these simulations are run with the same physics with the sole difference being the random seed utilized to generate the initial conditions.Therefore, they should be treated as six typical random places in the universe.These CROC simulations are run with the Adaptive Refinement Tree (ART) code (Kravtsov 1999;Kravtsov et al. 2002;Rudd et al. 2008) to reach high spatial resolution using the adaptive mesh refinement approach.The base grid is 39 ckpc/h in size, and the peak resolution is ≈ 100 pc (in physical units).CROC simulations include relevant physics such as gas cooling, heating, star formation and stellar feedback.In the simulations, individual star particles are the main radiation sources and the radiative transfer is done using the Optically Thin Variable Eddington Tensor (OTVET) method (Gnedin & Abel 2001), which is fully coupled to gas dynamics.
To create synthetic spectra for quasars, we postprocess the sightlines with 1D radiative transfer code described in Chen & Gnedin (2021).The sightlines are drawn from the full simulation snapshots.The full snapshots contain complete information of the simulation such as neutral fraction of hydrogen and helium, density, temperature and peculiar velocity of the gas.Because of the limited storage, full simulation snapshots are saved sparsely, and thus we only have 4 full snapshots with volume-averaged neutral fraction xHI between 0.99 and 10 −4 for each run.Because of this, we use all six boxes to create a more regular grid of xHI by binning the snapshots into groups around the following values: xHI = {4 × 10 −4 , 0.10, 0.36, 0.57, 0.83, 0.94, 0.98}, such that each xHI bin contains 2-3 snapshots.In addition to making the neutral fraction grid more regular, we do this to increase the cosmic variance captured by the CROC models at a given neutral gas fraction, while keeping the variance comparable across the different neutral fraction models.
From each simulation snapshot, we select the 20 most massive halos and draw 50 random sightlines centered on them.The masses of these halos at z = 6.8 are between 1.1 × 10 11 M ⊙ < M h < 1.1 × 10 12 M ⊙ .To simulate each quasar, we keep the neutral fraction and temperature of each cell unchanged while we scale the physical length and density to the redshift of the quasar by the expansion factor (same as in the ATON simulation).The total number of ionizing photons is calculated from M 1450 by adopting the average broken power-law spectral shape measured by Lusso et al. (2015), see Eq. (1) (again, same as in ATON).We post-process the sightlines in each snapshot for the following range of quasar lifetimes: log t Q = {3.0,3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 7.8, 8.0} Myr.Note that for the highest-z bin, all the 50 sightlines per halo are post-processed with quasar parameters according to the observed values, however, due to computational limitations only 5 sightlines per halo are post-processed for most of the quasars in the two lower redshift bins.We checked that changing the exact number of sightlines does not significantly alter our results.

Comparison of CROC and ATON
This section summarizes the main similarities and differences of the two hydrodynamical simulations employed in this work.We emphasize that the aim of this paper is not to perform an in-depth comparison of reionization simulations, but rather to use the two simulations to estimate systematic uncertainties and increase the robustness of our constraints.
In both simulations, we choose the most massive halos to host mock quasars, each of which is magnitude and redshift matched to quasars in our sample.Note that redshift-matching amounts to rescaling the physical scales and the density fields to the redshift of the quasar across snapshots at different simulation redshifts -the simulation redshift is only a proxy to the volumeaveraged neutral gas fraction.Furthermore, in both simulations we use the same broken power law spectral template (Lusso et al. 2015) to convert the magnitudes of the quasars in our sample to the number of ionizing photons.
The main differences between the two simulations are how radiative transfer is treated, the simulated volumes, and the resolution.The CROC simulation implements a 3D radiative transfer that is fully coupled to gas dynamics, whereas in ATON the radiative transfer is only done in postprocessing.Satyavolu et al. (2023) tested how using their 1D radiative transfer impacts the hydrogen and helium ionization fractions and the gas temperature as compared to a full 3D radiative transfer, and found the difference minimal.
The size of the simulated box is four times larger in ATON (160 cMpc/h) than CROC (40 cMpc/h), which can have an impact on the sizes of ionized bubbled that are modelled in each case.Iliev et al. (2014) has shown that ionized bubbles require box sizes of > 100 cMpc/h to converge.Smaller simulation volumes can thus contain systematically smaller ionized bubbles and hence produce stronger damping wings in the mock sightlines (see Keating et al. 2023).
Lastly, the high-spatial resolution of the CROC simulation produces DLAs which can also bias the simulated sightlines towards stronger damping wings.To reduce this bias, we exclude all simulated sightlines with neutral hydrogen column densities of N HI > 5 × 10 19 cm −2 (motivated by Chen & Gnedin (2021)).Note that we also checked the impact of an increased spatial resolution on the resultant sightlines in the lower-resolution ATON simulation, but found the differences negligible.We forward model the uncertainties by adding the smooth continuum prediction noise (second column), random measurement noise (third column), and finally the corresponding redshift error along the x-axis (fourth column; σz = 0.01 for illustration purposes).We repeat this sightline processing for all quasars within a given redshift bin (NQSO), and then stack NQSO processed sightlines to form a sample stack (rightmost column).The whole procedure is repeated 10000 times for the same combination of xHI and log tQ in order to calculate the final model and noise statistics that are used for inference.
Overall, performing two independent analyses using two very different simulations helps us mitigate systematics associated with modelling, thus increasing the robustness of the resultant constraints.At the same time, it serves as a consistency check for reionization simulations, however, a more detailed comparison is beyond the scope of this paper.

INFERRING xHI AND log t
In this section, we explain the construction of models from the two hydrodynamical simulations and how they are used to infer xHI and log t Q for each stacked quasar spectrum.

Forward modeling of uncertainties and model construction
In order to construct realistic models for a grid of xHI and log t Q values, we forward-model four different sources of uncertainties present in our observed spectral transmission profiles: the uncertainty from the spectral measurement (coming from the reduction pipeline), the uncertainty from the continuum prediction (originating from QSANNdRA), the uncertainty in quasar redshifts, and cosmic variance (due to the limited number of sightlines in our analysis).
In order to account for cosmic variance, we stack the simulation sightlines corresponding to different quasars in the same way as we do with the observed spectra -we illustrate this forward modeling in Fig. 5. Specifically, we sample one simulation sightline per quasar in the corresponding redshift bin (leftmost panel of Fig. 5).For this simulation sightline, we randomly draw a continuum prediction sample from QSANNdRA (second panel from the left in Fig. 5) as well as a noise vector from a Gaussian distribution centered at 0 and with a standard deviation defined by the measurement noise at each spectral pixel (third panel of Fig. 5).This combination of errors ensures that we model the randomness of the measurement error while still accounting for the fact that the error in continuum prediction is correlated between neighbouring pixels (i.e. it is smoothly varying over wavelength).Afterwards, we sample a new quasar redshift from a Gaussian distribution defined by the quasar redshift z and the corresponding redshift uncertainty σ z as given in Table 1, and we shift the sightline wavelength accordingly (fourth panel of Fig. 5).Having added the first three sources of uncertainty to each quasar sightline, we stack the resultant sightlines to form a sample  stacked sightline for each redshift bin (rightmost panel of Fig. 5.At this stage, we also simultaneously rebin the model sightline to the same resolution as our spectral stacks, i.e. 600 km/s.We repeat this procedure 10000 times to create 10000 sample spectral stacks for each redshift bin and for each combination of xHI and log t Q provided by the two simulations.The mean of these sample stacks constitutes our model, m(x HI , log t Q , ⟨z⟩).
Note the dependence on ⟨z⟩ -each redshift bin has a separate set of models.We show examples of models con-structed this way in Fig. 6 for both ATON and CROC models (top and bottom panels, respectively).

Maximum likelihood inference
With models, m(x HI , log t Q , ⟨z⟩), at hand, we use maximum likelihood estimation to find a combination of xHI and log t Q that best fits each of our stacked spectra.We use a multivariate Gaussian log-likelihood function, i.e.where y = y(⟨z⟩) is the observed stacked spectrum at the mean redshift ⟨z⟩, m = m(x HI , log t Q , ⟨z⟩) is the model stacked spectrum, and C = C(x HI , log t Q , ⟨z⟩) is the covariance matrix computed from the 10000 stacked spectra samples for each model.The use of covariance matrices allows us to obtain unbiased constraints on xHI and log t Q as the spectral flux is invariably correlated across neighboring pixels.Note that we compute the covariance matrices by forward-modeling noise in the models, as each stacked spectrum contains too few sightlines to calculate covariances on the data.We show an example of a covariance matrix at xHI = 0.37(0.38)and log t Q = 6.0 for both ATON and CROC models in Fig. 7.
To infer the neutral fraction and quasar lifetime constraints, we evaluate the log-likelihood function across the parameter grids corresponding to the two simulation models.We calculate the joint probability p(x HI , log t Q ) by normalizing the likelihood such that xHI,log tQ L(x HI , log t Q ) = xHI,log tQ p(x HI , log t Q ) = 1, and we plot this probability normalized to its peak in Fig. 8 for each of the three aforementioned redshiftstacked spectra using ATON and CROC models, respectively.The peak of p(x HI , log t Q ) is marked by a star and corresponds to the best fit model shown in the bottom panels of Fig. 8.We also tested the performance of this method on mock spectral stacks with models and covariances from the respective simulations in § B.
Further, we calculate the posteriors p(x HI ) and p(log t Q ) as p(x HI ) ∝ Σ log tQ L(x HI , log t Q ), (3) such that xHI p(x HI ) = 1 and log tQ p(log t Q ) = 1, from which we infer the 1σ uncertainties on our constraints.Due to the coarseness of our parameter space combined with the degeneracy between xHI and log t Q , the posterior distribution and hence the uncertainties may be asymmetric.
For comparison, we also include 1D inferences of the neutral gas fraction in Appendix C for both ATON and CROC simulations.There, we assume a fixed mean quasar lifetime of log t Q = 6.0 in each redshift bin, a value motivated by the measurement of the average quasar lifetime for a sample of 15 quasars at z ∼ 6 by Morey et al. (2021).Such inference is simpler as it bypasses the problem of the neutral fraction/quasar lifetime degeneracy and thus serves as a consistency check.
Note that not all quasars in each redshift bin are expected to have the same lifetime -in fact it has been shown that quasar lifetimes at z ∼ 3 − 4 can be modeled well with a lognormal distribution (Khrykin et al. 2021).Our model construction draws only sightlines corresponding to the same quasar lifetime for a particular value of log t Q and hence potentially reduces the underlying variance.We tested our inference pipeline against this concern in § D, where we show that our models perform considerably well even on mock spectral stacks with an underlying distribution of lifetimes.In order to avoid additional computational costs, we decided not to relax this assumption.

RESULTS
In this section we present our constraints on the volume-averaged neutral gas fraction, as well as the sample-averaged quasar lifetimes for each redshift bin.The results are summarized in Table 2. Figure 8. Maximum likelihood estimation with ATON models (top two rows) and CROC models (bottom two rows) for all three redshift bins.In the top panel, we show the 2D joint probability distribution evaluated for each combination of log tQ and xHI, as well as the 1σ and 2σ contours as white dashed and dotted curves, respectively.Note that the parameter space where we do not have ATON models is hatched.In the bottom panel, we plot stacked spectrum for each bin along with the model corresponding to the maximum likelihood point on the 2D grid in the top panel (marked by a star).
The maximum likelihood fits of the two higher redshift bins in ATON (green and red in the top of Fig. 8) seem to underestimate the strength of the observed damping wings.This is simply due to the coarseness of our parameter grid, whereby a neighboring model with a higher neutral gas fraction simply does not provide a better fit than the one shown in Fig. 8. Also note that our constraints do not directly follow the differences in the damping wings in Fig. 4 -for example, the ⟨z⟩ = 6.10 and the ⟨z⟩ = 6.46 ATON constraints both yield the same maximum likelihood neutral gas fraction despite the differences in the transmitted flux redward of Lyα .This is again due to the coarseness of our grid as well as due to the fact that we are jointly fitting for xHI as well as the quasar lifetime.This limitation could be overcome by only fitting the red-side damping wings of our spectral stacks, as suggested by Chen (2023), which are agnostic to the stochasticity of the density fields along the individual quasar sightlines -this is beyond the scope of this paper and will be explored in a separate work.
In some cases, particularly in the case of the lowest two redshift bins (blue and green), the error bars on the CROC constraints are much larger than the error bars on the ATON constraints.In fact, the 1σ of the CROC posterior reaches the lowest neutral gas fraction in our grid in the ⟨z⟩ = 6.10 case, which is why we indicate these constraints as limits at xHI ∼ 10 −4 .We hypothesize that the difference in constraining power could be due to the limited number of simulated sightlines (up to 100 for each snapshot compared with > 500 per snapshot in ATON), as we do not see such large difference in the highest redshift bin (red).If this were true, the models as well as the covariances would be dominated by the cosmic variance in the simulation and too noisy to provide a precise constraint.The situation is further complicated by the insensitivity to varying neutral gas fractions for the longest simulated lifetimes, which arises from the degeneracy between xHI and log t Q and is difficult to mitigate.
It is also worth noting that the ⟨z⟩ = 6.46 likelihood map as well as the posteriors exhibit a bimodality that is consistent between ATON and CROC -a good fit is either provided by a longer lifetime quasar embedded in a more neutral IGM or a shorter lifetime quasar in a less neutral IGM.This degeneracy causes the inferred error bars in Fig. 9 to be asymmetric, and, even more importantly, explains the apparent tension between the two constraints at ⟨z⟩ = 6.46 -the maximum likelihood models of CROC and ATON, respectively, fall at the two different peaks of this bimodality.The bimodal posterior could be either caused by having two different populations of quasars in our intermediate redshift bin (therefore seeing a combination of small and large proximity zones in the stacked spectrum), or it could arise simply due to large cosmic variance -in the latter case, the increased transmission around λ rest ≲ 1210 Å is due to the sampled density fields rather than the quasar lifetimes.Either way, an increased number of quasars in this redshift bin should provide a more conclusive neutral fraction constraint.
Note that the ⟨z⟩ = 6.87 bin also exhibits a multimodal likelihood map for the CROC models -in this case the apparent multimodality is more likely a consequence of the grid coarseness rather than cosmic variance.
In Fig. 9, we show a comparison of our neutral gas fraction constraints to existing constraints in the literature.In the lower end of our redshift range, z ∼ 6, our constraints are consistent with limits from Lyα and Lyβ dark pixels McGreer et al. (2015); Jin et al. (2023) and dark gaps (Zhu et al. 2022), as well as with limits on the Lyα forest opacity (Fan et al. 2006;Yang et al. 2020b;Bosman et al. 2021).Towards redshifts above z ≳ 7, we are consistent with constraints from single-sightline quasar damping wings (Wang et al. 2021;Greig et al. 2017Greig et al. , 2019;;Ďurovčíková et al. 2020;Davies et al. 2018a;Bañados et al. 2018;Yang et al. 2020a), as well as with galaxy limits (Sobacchi & Mesinger 2015;Mason et al. 2018;Umeda et al. 2023).Furthermore, our constraints are also consistent with the dark pixel limits from Jin et al. (2023) at 6 ≲ z ≲ 7. Overall, the neutral gas fraction constraints in this work bridge the gap between redshifts z ∼ 6 and z ∼ 7 and complete the chronicle of the Epoch of Reionization.
Having measurements where previously only limits were present, we now discuss the impact our results have on reionization models.Although there is a growing consensus that the EoR was driven by star-forming galaxies as opposed to quasars (Robertson et al. 2015;Robertson 2022), the population of galaxies that drove reionization remains disagreed upon.Models driven primarily by the most massive galaxies tend to yield late and rapid reionization (Naidu et al. 2020), whereas a slower, ).On the lower redshift end, squares denote constraints from the Lyα forest optical depth (Fan et al. 2006;Yang et al. 2020b;Bosman et al. 2021), and downward and upward triangles label constraints from Lyα and Lyβ dark pixels (McGreer et al. 2015;Jin et al. 2023) and dark gaps (Zhu et al. 2022), respectively.At higher redshifts, single-sightline quasar damping wing results are shown as diamonds (Wang et al. 2021;Greig et al. 2017Greig et al. , 2019;;Ďurovčíková et al. 2020;Davies et al. 2018a;Bañados et al. 2018;Yang et al. 2020a), and galaxy constraints from Lyman break galaxies (LBGs), Lyα emitters (LAEs) and the recent galaxy damping wing measurements are shown as circles (Mason et al. 2018;Ouchi et al. 2010;Sobacchi & Mesinger 2015;Mason et al. 2019;Umeda et al. 2023).We also show two examples of reionization models: one that reionizes late and rapidly (e.g.Naidu et al. 2020), and one with an earlier and slower reionization (e.g.Finkelstein et al. 2019).Note that the ATON and CROC constraints are slightly offset in redshift for better visibility.
earlier reionization has been demonstrated if the much more abundant, fainter galaxies dominate the ionizing photon budget (Finkelstein et al. 2019;Rosdahl et al. 2022).However, the implication of late vs. early reionization on the population of ionizing sources rests on many assumptions.For instance, the dependence of the Lyα escape fraction on the galaxy mass varies between models and can lead to late reionization that is not necessarily driven by the most massive galaxies (Keating et al. 2020;Kulkarni et al. 2019).Our measurements, as shown in Fig. 9, are consistent with both of these alternatives, but slightly favour a late reionization history.
At present, we are not able to rule out the early model due to the large uncertainties on our constraints, but the extension of this analysis to a larger sample of quasars should be able to yield a more conclusive answer.

Quasar lifetime constraints
We further marginalize the joint probability over the grid of neutral gas fractions and obtain the following average quasar lifetime constraints in the three redshift bins (Table 2): log t Q = 6.0 +1.0 −1.0 at ⟨z⟩ = 6.10, log t Q = 5.0 +2.0 −1.0 at ⟨z⟩ = 6.46, and log t Q = 5.0 +1.0 −1.0 at ⟨z⟩ = 6.87 (ATON, top two rows in Fig. 8), and log t Q = 6.5 +1.3 −0.5 at ⟨z⟩ = 6.10, log t Q = 7.0 +0.5 −0.5 at The measurements from this work are shown in red and blue for ATON and CROC models, respectively.The figure is adapted from (Eilers et al. 2021), but here we only include datapoints that represent population averages for clarity (all in gray).Constraints from different methods are encoded by symbol: hydrogen and helium proximity zones are shown as diamonds (Morey et al. 2021) and squares (Khrykin et al. 2021), respectively, pentagons denote lifetimes measured from quasar clustering (Laurent et al. 2017;Shankar et al. 2010;Shen et al. 2007;White et al. 2012), constraints from the extents of Lyα nebulae are shown as upward triangles (Trainor & Steidel 2013;Borisova et al. 2016), downward triangles encode measurements of the transverse proximity effect (Kirkman & Tytler 2008;Schmidt et al. 2017;Oppenheimer et al. 2018), and host population duty cycle studies are shown as circles (Yu & Tremaine 2002;Chen & Gnedin 2018).The horizontal red-shaded band indicates the range of lifetimes one would expect to measure if the black holes underwent a continuous, Eddington-limited accretion.Note that the ATON and CROC constraints are slightly offset in redshift for better visibility.⟨z⟩ = 6.46, and log t Q = 6.5 +0.5 −0.5 at ⟨z⟩ = 6.87 (CROC, bottom two rows in Fig. 8).
Again, the lifetime constraints resulting from the two different simulations are consistent within error bars, with the CROC constraints indicating systematically longer quasar lifetimes.We again observe the effect of the aforementioned bimodality in the ⟨z⟩ = 6.46 bin, where the ATON and CROC results deviate more than in the other two bins.As explained in § 5.1, this either suggests two populations of quasars with different lifetimes or that the cosmic variance is too large in this redshift range -in either case, a larger quasar sample will provide a clarification.
We also note that most of the error bars on lifetime constraints are limited by the coarseness of our model grids.This is particularly true for the ATON models, where we only sample order-of-magnitude differencesfor example, the confidence intervals on the ⟨z⟩ = 6.10 and ⟨z⟩ = 6.89 constraints are quite conservative for this reason.
These quasar lifetime measurements are consistent with other constraints from the literature, as shown in Fig. 10 (note that only population averages instead of individual lifetime measurements are shown for the sake of clarity).Of particular importance is the diamond data point at z ∼ 6 which is most comparable to our method and represents the only other HI proximity zone population average in the literature to date (Morey et al. 2021).Additionally, our quasar sample overlaps with the quasar sample used by Morey et al. (2021), and so it is crucial to see our constraints be consistent with this measurement.Furthermore, single-sightline damping wing analyses of quasars at z ≳ 7 by Davies et al. (2018a); Wang et al. (2020); Yang et al. (2020a) all find quasar lifetimes log t Q ≲ 6.5, consistent with both of our highest-redshift bin measurements.
Technically speaking, the exact interpretation of the lifetimes inferred from HI proximity zones depends on the ionization state of the IGM.This is because quasars are unlikely to accrete at a constant rate in a single episode -in fact, quasars are believed to exhibit a flickering light curve whereby the black hole experiences episodes of high and low/no accretion (e.g.Zhou et al. 2023;Hopkins & Hernquist 2009).For a highly ionized IGM, the proximity zone sizes are sensitive to the episodic lifetime, i.e. the most recent episode of quasar activity, as the imprints of any previous accretion phases would have been washed out by the recombination of hydrogen within the proximity zone during the quasar's "off"-phase.On the other hand, a highly neutral IGM enables the measurement of the integrated lifetime, as the abundance of neutral hydrogen impedes a complete recombination within the proximity zone during the quasar's "off" phases (Davies et al. 2019).Hence, on the lower redshift end of our quasar sample, we expect to be sensitive to the episodic lifetime, whereas the highest redshift quasar spectra are likely beginning to be sensitive to the integrated lifetime.However, stacking multiple quasar sightlines of varying neutral gas fraction renders distinguishing these two scenarios tricky and so we adopt a light-bulb model (where the episodic and the integrated lifetimes are the same) in order to interpret our results.Table 2. Constraints on the neutral gas fraction, xHI, and the mean quasar lifetime, log tQ, for each of the three redshift bins.The quoted uncertainties are the 16th and 84th percentiles of the marginalized posterior.Due to the discreteness of our parameter grid, we quote the distance to the nearest grid value in cases where either percentile coincides with the constraint value -in this case the uncertainties are not properly resolved by our grid and the 1σ constraints are therefore conservative.
Under the assumption of a light-bulb light curve, our measurements provide further support for the existence of a young quasar population at high redshifts (Eilers et al. 2018b(Eilers et al. , 2021;;Morey et al. 2021) and pose a challenge to the standard models of SMBH growth.The masses of most quasars in our sample have been measured and show a mean mass of 1.5 × 10 9 M ⊙ (Bigwood et al. in prep, also Farina et al. 2022;Mazzucchelli et al. 2023)).If a continuous, Eddington-limited accretion with a radiative efficiency of ϵ ∼ 0.1 (Soltan 1982;Yu & Tremaine 2002) is assumed, it takes almost a billion years to grow a black hole mass of 10 9 M ⊙ from a 100 M ⊙ seed (Inayoshi et al. 2020), i.e. several orders of magnitude longer than our measurements suggest.
Continuous Eddington-limited accretion is unlikely to occur over epochs of billions of years.As has been recently shown in a self-consistent manner (Zhou et al. 2023), quasars are expected to undergo flickering light curves whereby the black hole accretion rate and thus the number of ionizing photons that enter the IGM varies.As such, growing the supermassive black holes at the centers of these quasars might not be so problematic towards the end of Reionization, as small proximity zones simply indicate only the last "on" phase.For higher neutral gas fractions, however, the story of SMBH growth is still quite unclear.Hence, other growth pathways, such as quasar obscuration or radiatively inefficient accretion need to be invoked, details of which are discussed in previous works (Eilers et al. 2017(Eilers et al. , 2018b;;Davies et al. 2019;Satyavolu et al. 2023).

SUMMARY
In summary, we have used a sample of 18 quasars at redshifts 6.0 ≤ z ≤ 7.1 and computed continuumnormalized spectral stacks to constrain the neutral gas fraction as well as the mean quasar lifetimes in three redshift bins.
• For the first time, we showcase the emergence of damping wings between redshifts 6 ≲ z ≲ 7 in stacked quasar spectra (Fig. 4).
• Our measurements of the neutral gas fraction and quasar lifetimes are consistent between the two simulations used in this work, ATON and CROC, despite their differences.This serves as an important consistency check for simulation based modeling and increases the robustness of our results.

DATA RELEASE
The quasar spectra used in this work will be made public upon the acceptance of this paper.We publish the reduced spectral files along with the smoothed spectral fits and the mean continuum predictions computed in this work (see § 2.2 for details).See Table 3 for a description of the published files.We would also like to thank Carlos Contreras, Matías Díaz, Carla Fuentes, Mauricio Martínez, Alberto Pastén, Roger Leiton, Hugo Rivera, and Gabriel Prieto for their help and support during the Magellan/FIRE observations.This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.
Based on observations collected at the European Southern Observatory under ESO programmes 084.A-0360, 086.A-0162, 087.A-0607, 098.B-0537, 286.A-5025, 089.A-0814, and 093.A-0707.Some of the data presented herein were obtained at Keck Observatory, which is a private 501(c)3 non-profit organization operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration.The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the Native Hawaiian community.We are most fortunate to have the opportunity to conduct observations from this mountain.
Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S.  transmitted flux around the N V line.The interesting thing to note is that in both of these cases, the SDSS nearest neighbor composite shown in Fig. 3 also underpredicts the flux in this region -this means that either the C III line that is masked by the region of telluric absorption is unusual in these two quasars, or this discrepancy shows a limitation of the training set of the neutral network.The former explanation would mean that the C III line profile is important in an accurate reconstruction of the Lyα profile.On the other hand, J1148+5251 is predicted to have an exceptionally strong damping wing, in agreement with the SDSS nearest neighbor composite.Upon closer inspection, one can notice that the C IV line is partially masked due to atmospheric absorption, which renders this prediction less reliable due to our inability to make out the full profile of this emission line.

B. INFERENCE PIPELINE TESTING
In order to ensure that the inference pipeline described in § 4 works as expected, we tested that it works on mock stacked spectra from simulations.To this end, we perform the same procedure as outlined in § 4.1 for each simulation to create mock spectral samples for each redshift bin -more specifically, we draw the corresponding number of sightlines for each quasar, add the measurement, continuum prediction and redshift uncertainties, and stack these mock sightlines to produce a sample stacked spectrum at a known neutral fraction, xHI , and quasar lifetime, log t Q .Subsequently, we run our inference pipeline on this mock stacked spectrum to test whether we are able to fit the true xHI and log t Q .We repeated this procedure 100 times for each combination of the neutral fraction and quasar lifetime, and we plot the resultant performance as confusion matrices in Fig. 12.
In the case of both simulations, ATON (top) and CROC (bottom), we see that the resultant confusion matrices are mostly diagonal, with the best accuracy achieved for the highest neutral gas fractions and the longest quasar lifetimes.It is understandable that we see some off-diagonal terms -these are consequences of the different sources of noise that we are forward modeling as well as the degeneracy between the neutral fraction and quasar lifetimes.Overall, Fig. 12 indicates that our inference pipeline works well.

C. FITTING xHI AT A FIXED LIFETIME
For completeness, here we show the maximum likelihood inferences of the neutral gas fraction assuming a fixed quasar lifetime, log t Q = 6.0.The choice of this quasar lifetime is motivated by the first measurement of the mean quasar lifetime in a sample of z ∼ 6 quasars (Morey et al. 2021).
The resultant inferences are shown in Fig. 13 for both the ATON simulation (top two rows) and the CROC simulation (bottom two rows).By assuming a prior on the mean lifetime in each redshift bin, the inference on the neutral gas fraction becomes simpler -in particular, we do not observe the posterior bimodality in the ⟨z⟩ = 6.46 that arises in the joint inference in Fig. 8.
The fits shown in Fig. 13 showcase a monotonical neutral gas fraction evolution with redshift, as one would expect based on the visual inspection of the damping wing profiles Fig. 4.However, only estimating the neutral gas fraction does not capture the full picture due to the degeneracy between xHI and log t Q , which is why we only include this analysis in the appendix and focus on the joint analysis in the main text.

D. MEAN LIFETIME VS DISTRIBUTION OF LIFETIMES
In § 4.1, where we described the forward modeling of uncertainties and the construction of models, we noted that the simulation sightlines for a model at a fixed log t Q were all drawn from runs corresponding to the same quasar lifetime.This is equivalent to assuming that all quasars in a particular redshift bin have the same lifetime, equal to the mean lifetime that we infer alongside the neutral gas fraction.In this section, we test our models against relaxing this assumption.Specifically, we test whether a mock simulation stack that constains sightlines at a range of lifetimes for a particular xHI can be accurately fit with our models for its true neutral gas fraction and its true mean quasar lifetime.
To this end, we follow the same procedure as in § B but for a given xHI , we stack sightlines corresponding to N QSO draws of log t Q from a lognormal distribution centered on the mean quasar lifetime.Inspired by Khrykin et al. (2021), we assume the variance of the lifetime distributions to be 1 dex.We plot the resultant confusion matrices in Fig. 14 for both ATON (top) and CROC (bottom).Note that we keep the same forward modeling of uncertainties as in § B.
The right panels in Fig. 14 now show a more "fuzzy" diagonal structure -it is now harder for our models to infer the correct mean quasar lifetime in the particular mock stack.Notwithstanding, the resultant confusion matrices for .Confusion matrices displaying the performance of our inference pipeline on mock stacked spectra for both ATON (top panels) and CROC (bottom panels) simulations.In all four panels, we plot the true value of xHI (left) and log tQ (right) on the y-axis against the respective best fit values on the x-axis.All four confusion matrices are mostly diagonal, with the off-diagonal terms being a consequence of the forward-modelled uncertainties as well as the aforementioned degeneracy.
the neutral gas fractions (the left panels) are still fairly diagonal, which showcases that our method recovers the true neutral gas fraction considerably well even if the quasars in our stacked spectra (Fig. 4) have different lifetimes.  .Confusion matrices displaying the performance of our inference pipeline on mock stacked spectra at a range of quasar lifetimes for both ATON (top panels) and CROC (bottom panels) simulations.In all four panels, we plot the true value of xHI (left) and the true mean quasar lifetime of the mock stack, log tQ, (right) on the y-axis against the respective best fit values on the x-axis.This test demonstrates that our inference of the neutral gas fraction is robust against quasars in our stacked spectra having different lifetimes.

]Figure 1 .
Figure1.The 18 quasar spectra used in this study sorted and color-coded by redshift.The error vector is shown as a thin gray line.

Figure 2 .
Figure2.The magnitude vs. redshift distribution of quasars in our sample (colors corresponding to the three different redshift bins used in our analysis) compared to the high-redshift quasar database published byFan et al. (2023) shown as gray contours.

Figure 3 .
Figure3.Quasar continuum reconstruction for all objects in our sample predicted by QSANNdRA.We show the observed data in black and the measurement noise vector in light gray.QSANNdRA outputs a range of predictions (in light blue) as well as a mean prediction (dark blue) based on the red spectral fit above 1290 Å (red curves).The flux in the shaded telluric regions is determined from the nearest neighbour SDSS stack shown in orange.

Figure 4 .
Figure4.Stacked spectra of quasars in our sample in three distinct redshift bins show the emergence of damping wings between z ∼ 6 and z ∼ 7 (thick colored curves, the location of Lyα is marked by the dashed black line).The number of quasars within a given redshift bin is noted in brackets, and the shaded region around the stacked spectra shows the variance in the mean flux.The thin colored lines show a higherresolution (∆v = 100 km/s) version of the three stacked spectra.

Figure 5 .
Figure5.Illustration of the forward modeling of uncertainties and model construction for both simulations (ATON in the top panel, CROC in the bottom panel).For a particular quasar and a combination of xHI and log tQ values, we draw one sightline at random from the corresponding simulation (first column from the left).We forward model the uncertainties by adding the smooth continuum prediction noise (second column), random measurement noise (third column), and finally the corresponding redshift error along the x-axis (fourth column; σz = 0.01 for illustration purposes).We repeat this sightline processing for all quasars within a given redshift bin (NQSO), and then stack NQSO processed sightlines to form a sample stack (rightmost column).The whole procedure is repeated 10000 times for the same combination of xHI and log tQ in order to calculate the final model and noise statistics that are used for inference.

Figure 6 .
Figure6.Example models from both ATON (top) and CROC (bottom) for a fixed neutral gas fraction (left) and a fixed quasar lifetime (right).The vertical black dashed line marks the position of Lyα .To ease visual comparison, we have highlighted the most comparable model (xHI = 0.37 (ATON) and xHI = 0.36 (CROC) with log tQ = 6.0) in a thicker line in all four panels.

Figure 7 .
Figure 7. Example covariance matrix for each simulation.Here we show the covariance corresponding to xHI = 0.37 and xHI = 0.38, and log tQ = 6.0 for ATON (top) and CROC (bottom), respectively.Despite originating from different simulations, the covariance matrices are remarkably similar.
p(x HI , log t Q ))

Figure 9 .
Figure 9. Constraints on the neutral fraction evolution across cosmic time.The constraints from this work are shown in red and blue for ATON and CROC models, respectively.All existing constraints are shown in gray and are overplotted on constraints from the Cosmic Microwave Background (CMB; Planck Collaboration et al. 2020).On the lower redshift end, squares denote constraints from the Lyα forest optical depth(Fan et al. 2006; Yang et al. 2020b;Bosman et al. 2021), and downward and upward triangles label constraints from Lyα and Lyβ dark pixels(McGreer et al. 2015;Jin et al. 2023) and dark gaps(Zhu et al. 2022), respectively.At higher redshifts, single-sightline quasar damping wing results are shown as diamonds(Wang et al. 2021;Greig et al. 2017Greig et al. , 2019;;Ďurovčíková et al. 2020;Davies et al. 2018a;Bañados et al. 2018; Yang et al. 2020a), and galaxy constraints from Lyman break galaxies (LBGs), Lyα emitters (LAEs) and the recent galaxy damping wing measurements are shown as circles(Mason et al. 2018;Ouchi et al. 2010;Sobacchi & Mesinger 2015;Mason et al. 2019;Umeda et al. 2023).We also show two examples of reionization models: one that reionizes late and rapidly (e.g.Naidu et al. 2020), and one with an earlier and slower reionization (e.g.Finkelstein et al. 2019).Note that the ATON and CROC constraints are slightly offset in redshift for better visibility.

Figure 10 .
Figure10.Constraints on quasar lifetimes across redshifts.The measurements from this work are shown in red and blue for ATON and CROC models, respectively.The figure is adapted from(Eilers et al. 2021), but here we only include datapoints that represent population averages for clarity (all in gray).Constraints from different methods are encoded by symbol: hydrogen and helium proximity zones are shown as diamonds(Morey et al. 2021) and squares(Khrykin et al. 2021), respectively, pentagons denote lifetimes measured from quasar clustering(Laurent et al. 2017;Shankar et al. 2010;Shen et al. 2007;White et al. 2012), constraints from the extents of Lyα nebulae are shown as upward triangles(Trainor & Steidel 2013;Borisova et al. 2016), downward triangles encode measurements of the transverse proximity effect(Kirkman & Tytler 2008;Schmidt et al. 2017;Oppenheimer et al. 2018), and host population duty cycle studies are shown as circles(Yu & Tremaine 2002;Chen & Gnedin 2018).The horizontal red-shaded band indicates the range of lifetimes one would expect to measure if the black holes underwent a continuous, Eddington-limited accretion.Note that the ATON and CROC constraints are slightly offset in redshift for better visibility.

Figure 11 .
Figure11.Predicted continua from QSANNdRA (blue) for all quasars in our sample, zoomed in on the wavelength range around the Lyα emission that is relevant for the inference of xHI and log tQ.The three columns represent the different redshift bins, from the highest redshift bin on the left to the lowest redshift bin on the right.
Figure12.Confusion matrices displaying the performance of our inference pipeline on mock stacked spectra for both ATON (top panels) and CROC (bottom panels) simulations.In all four panels, we plot the true value of xHI (left) and log tQ (right) on the y-axis against the respective best fit values on the x-axis.All four confusion matrices are mostly diagonal, with the off-diagonal terms being a consequence of the forward-modelled uncertainties as well as the aforementioned degeneracy.

Figure 13 .
Figure13.Maximum likelihood inference of the neutral gas fraction, xHI, for a fixed quasar lifetime, log tQ = 6.0.The choice of this lifetime value was inspired by the mean quasar lifetime measured for a sample of z ∼ 6 quasars(Morey et al. 2021).The top two rows correspond to the ATON simulation, whereas the bottom two rows correspond to the CROC simulation.The top panels in each case show the log-space normalized posterior probability p(xHI), and the stars mark the position of the best fit plotted in the bottom panels.
Figure14.Confusion matrices displaying the performance of our inference pipeline on mock stacked spectra at a range of quasar lifetimes for both ATON (top panels) and CROC (bottom panels) simulations.In all four panels, we plot the true value of xHI (left) and the true mean quasar lifetime of the mock stack, log tQ, (right) on the y-axis against the respective best fit values on the x-axis.This test demonstrates that our inference of the neutral gas fraction is robust against quasars in our stacked spectra having different lifetimes.

Table 3 .
Content description of published FITS files for each quasar.
Department of Energy Office of Science, and the Participating Institutions.SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah.The SDSS website is www.sdss4.org.SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science,