The mean free path of ionizing photons at 5

The mean free path of ionizing photons, $\lambda_{\rm mfp}$, is a key factor in the photoionization of the intergalactic medium (IGM). At $z \gtrsim 5$, however, $\lambda_{\rm mfp}$ may be short enough that measurements towards QSOs are biased by the QSO proximity effect. We present new direct measurements of $\lambda_{\rm mfp}$ that address this bias and extend up to $z \sim 6$ for the first time. Our measurements at $z \sim 5$ are based on data from the Giant Gemini GMOS survey and new Keck LRIS observations of low-luminosity QSOs. At $z \sim 6$ we use QSO spectra from Keck ESI and VLT X-Shooter. We measure $\lambda_{\rm mfp} = 9.09^{+1.62}_{-1.28}$ proper Mpc and $0.75^{+0.65}_{-0.45}$ proper Mpc (68% confidence) at $z = 5.1$ and 6.0, respectively. The results at $z = 5.1$ are consistent with existing measurements, suggesting that bias from the proximity effect is minor at this redshift. At $z = 6.0$, however, we find that neglecting the proximity effect biases the result high by a factor of two or more. Our measurement at $z = 6.0$ falls well below extrapolations from lower redshifts, indicating rapid evolution in $\lambda_{\rm mfp}$ over $5<z<6$. This evolution disfavors models in which reionization ended early enough that the IGM had time to fully relax hydrodynamically by $z = 6$, but is qualitatively consistent with models wherein reionization completed at $z = 6$ or even significantly later. Our mean free path results are most consistent with late reionization models wherein the IGM is still 20% neutral at $z=6$, although our measurement at $z = 6.0$ is even lower than these models prefer.


INTRODUCTION
The metagalactic UV background is a fundamental link between the intergalactic medium (IGM) and the sources of ionizing radiation (stars and active galactic nuclei). Much of our knowledge of the IGM comes from observations of the Ly forest, whose opacity depends directly on the hydrogen ionization rate, Γ. For a given ionizing emissivity, , the ionization rate scales roughly as Γ ∝ mfp (e.g., Haardt & Madau 2012), where mfp is the mean free path of ionizing photons. Accurate measurements of mfp are therefore essential for translating the measured properties of the IGM into constraints on the ionizing sources.
The redshift evolution of mfp may also reflect the timing of reionization (e.g., Rahmati & Schaye 2018). A number of observations now suggest that reionization had a midpoint around ∼7-8 and ended near ∼ 6, or even later. These include (i) the electron optical depth to CMB photons (Planck Collaboration et al. 2020), ★ E-mail: george.becker@ucr.edu (ii) the decline in Ly emission from galaxies at > 6 (e.g., Jung et al. 2020; Morales et al. 2021, and references therein), (iii) largescale opacity fluctuations in the Ly forest at < 6 (Fan et al. 2006;Becker et al. 2015;Bosman et al. 2018;Eilers et al. 2018;Yang et al. 2020), (iv) the association of large Ly troughs at ∼ 5.7 with galaxy underdensities Kashino et al. 2020), (v) Ly damping wings seen in the spectra of ∼ 7 QSOs (Mortlock et al. 2011;Greig et al. 2017Greig et al. , 2019Davies et al. 2018;Bañados et al. 2018;Wang et al. 2020), (vi) the thermal history of the IGM at > 5 (Boera et al. 2019;Walther et al. 2019;Gaikwad et al. 2020), and (vii) the evolution in the number density of neutral metal absorbers near ∼ 6 Cooper et al. 2019;Doughty & Finlator 2019). If reionization did end near or below = 6, then the mean free path at < 6 should increase rapidly with time as large H bubbles merge and the last remaining neutral islands are ionized (e.g., Wyithe et al. 2008). Indeed, recent models of late reionization exhibit a rapid evolution in mfp over 5 < < 6 (e.g., Kulkarni et al. 2019;Keating et al. 2020a,b;Cain et al. 2021). Additionally, absorbers in recently reionized gas are photoevaporated or pressure smoothed over a time scale Δ ∼ 100 Myr, contributing further to the rapid evolution in LyC opacity. In constrast, a significantly earlier reionization would give the IGM more time to relax hydrodynamically, producing a more gradual evolution in mfp at < 6 (Park et al. 2016;D'Aloisio et al. 2020;Cain et al. 2021).
Multiple techniques have been used to measure the mean free path. One approach is to calculate the ionizing opacity from the incidence rate of individual H absorbers (e.g., Miralda-Escude & Ostriker 1990;Meiksin & Madau 1993;Haardt & Madau 1996;Faucher-Giguère et al. 2008;Songaila & Cowie 2010;Rudie et al. 2013;Prochaska et al. 2014). Alternatively, one may directly estimate the opacity from the shape of the transmitted flux profile blueward of the Lyman limit in the mean spectra of QSOs (Prochaska et al. 2009). The latter approach has arguably produced the most precise estimates of the mean free path at high redshifts, with results now spanning 2 5 (Prochaska et al. 2009;Fumagalli et al. 2013;O'Meara et al. 2013;Worseck et al. 2014;Lusso et al. 2018). One can also infer the mean free path from the average of free paths along individual QSO lines of sight (Romano et al. 2019).
One challenge at > 5 is that the mean free path may be comparable to or shorter than the typical size of a QSO proximity zone. In that case, the ionizing flux from a QSO will tend to decrease the opacity in its vicinity, leading to mean free path measurements based on QSO spectra that are biased high (e.g., Worseck et al. 2014;D'Aloisio et al. 2018). One possible solution is to use fainter QSOs, for which the impact of the proximity zone will be decreased (see discussion in Worseck et al. 2014). This is observationally challenging, however, particularly given the low levels of transmission expected at 5.5. At ∼ 6 it is currently impractical to obtain enough high-quality spectra of QSOs that are sufficiently faint to meaningfully avoid the proximity effect.
In this work we perform new measurements of the mean free path at > 5, including the first direct measurement at ∼ 6, that address the proximity effect in two ways. First, we modify the direct measurement approach of Prochaska et al. (2009) to include a scaling of the opacity with the local ionization rate. This allows us to account for the decrease in opacity in the vicinity of a QSO. Second, we measure mfp at 5 from two groups of QSOs spanning a factor of five in mean luminosity. This provides additional leverage in separating the background opacity from the impact of the QSOs.
The rest of the paper is organized as follows. We describe the individual QSO spectra and the composites in Section 2. In Section 3 we outline our model formalism, perform tests with mock spectra, and derive measurements of mfp . We then discuss the implications of our results for reionization in Section 4 before summarizing the results in Section 5. Our observational results assume a ΛCDM cosmology with (Ω m , Ω Λ , 0 ) = (0.3, 0.7, 70 km s −1 Mpc −1 ). Distances are quoted in proper Mpc (pMpc) except where noted.

Samples
This work uses spectra from three QSO samples. First, we use a subset of the spectra from the Giant Gemini GMOS (GGG) survey presented by Worseck et al. (2014). Specifically, we use 40 QSOs spanning 5.00 < < 5.42, which have a mean redshift of = 5.16 and an absolute magnitude corresponding to the mean luminosity at rest-frame 1450 Å of 1450 = −26.8. Second, we include a sample of lower-luminosity QSOs at ∼ 5 observed with the Keck For we each sample we list the number of QSOs, the mean redshift, and the absolute magnitude at rest-frame 1450 Å corresponding to the mean luminosity.  Worseck et al. (2014) is shown for reference (dark blue pentagons) with vertical dashed lines marking their redshift bins. In this work we analyze the GGG QSOs at > 5, along with the samples observed with LRIS (light blue circles) and ESI + X-Shooter (green squares). Details of the individual QSOs are given in Table 2. LRIS spectrograph. The LRIS sample includes 23 QSOs spanning 4.93 < < 5.24 with = 5.09 and an absolute magnitude corresponding to the mean luminosity of 1450 = −25.1. This is a factor of five fainter than the GGG sample. Finally, we use a sample of 13 QSOs at ∼ 6 observed with the Keck ESI and VLT X-Shooter spectrographs. This sample spans 5.82 < < 6.08 with = 5.97 and an absolute magnitude corresponding to the mean luminosity of 1450 = −27.0. The samples are summarized in Table 1, and the QSOs included in each sample are listed in Table 2. We plot the rest-frame 1450 Å absolute magnitudes as a function of redshift for all of our QSOs in Fig. 1.

GGG spectra
Our subset of the GGG data includes all objects at > 5 observed in that survey apart from one flagged as a broad absorption line (BAL) QSO and one whose flux was affected by very poor sky subtraction. The 40 QSOs selected are listed in Table 2. Details of the observation and data reduction are given in Worseck et al. (2014). Here we note that the Lyman continuum portion of the spectra were observed with the GMOS B600 grating through a 1" slit, which gives a FWHM resolution of roughly 320 km s −1 . We also note that the spectra contain noticeable variations in the skylevel zero point, as discussed by Worseck et al. (2014). We account for these variations when fitting models to the mean flux profile (see Section 3). Redshifts for GMOS QSOs are adopted from Worseck et al. (2014). Other redshifts quoted to three decimal places are based on the apparent start of the Ly forest. See text for details. 1450 values for GMOS QSOs were calculated from the flux-calibrated spectra published by Worseck et al. (2014). For LRIS QSOs they are adopted from McGreer et al. (2013McGreer et al. ( , 2018. For ESI and X-Shooter QSOs the 1450 values are from Bañados et al. (2016) and references therein.

LRIS observations
We observed 27 faint ( 1450 ∼ −25) ∼ 5 QSOs in March and September 2019 using the Keck Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995). The targets were drawn from the surveys for faint QSOs conducted by McGreer et al. (2013McGreer et al. ( , 2018 in the SDSS Stripe 82 and the CFHT Legacy Survey fields. We used a 1.0 slit with the D680 dichroic. On the blue side we used the 300/5000 grism, which provided the maximum sensitivity near the Lyman limit for the QSOs in our sample (observed wavelengths near 5400-5700 Å). The resolution from this grism is relatively low (FWHM 490 km s −1 , measured from skylines) but sufficient for the mean free path measurement described in Section 3. On the red side we used the 831/8200 grating (FWHM 110 km s −1 ) centered at 7989 Å, which allowed us to identify individual absorption lines near the start of the Ly forest. The spectra were reduced using a custom reduction package similar to the one described in Becker et al. (2012) and Lopez et al. (2016). Individual frames were sky-subtracted using an optimal algorithm based on Kelson (2003). Preliminary one-dimensional spectra were then optimally extracted following Horne (1986) . For each exposure, a telluric absorption model was fit to the red side and then propagated back to the two-dimensional sky-subtracted frames for both the blue and the red side. A final one-dimensional spectrum for each side was then extracted simultaneously from all exposures of a given object. One complication of our chosen setup is that the D680 dichroic combined with the 300/5000 grism allows contamination from second-order light. This is nominally not a problem for our QSOs, which have essentially no flux blueward of ∼5000 Å; however, it does impact the spectra of blue standard stars, which in turn can impact the flux calibration reward of ∼6000 Å. We addressed this problem by using the type dG-K standard star G158-100, whose flux peaks near 5000 Å and declines rapidly towards the blue. Flux calibration derived from this standard produced a good match between the blue-and red-side spectra of our QSOs. The blue (red) side was extracted in wavelength bins of 120 (60) km s −1 .
Out of this sample, 23 QSOs were selected to create the composite described in Section 2.6. These objects are listed in Table 2 and their spectra are plotted in Appendix A. The remaining four QSOs were rejected either due to the presence of BAL features (J2245+0024, J0210+0003, J0218−002) or due to difficulty in measuring the redshift (J0215−0529).

ESI and X-Shooter spectra
Our ∼ 6 sample is drawn from the Keck ESI and VLT X-Shooter spectra used by Becker et al. (2019). A lower redshift bound of > 5.8 was chosen so that the entire spectrum blueward of the Ly emission line down to a rest-frame wavelength of 820 Å falls entirely in the VIS arm of X-Shooter. An upper bound of < 6.1 was chosen so that the Lyman series opacity of the IGM still allows some possibility of measuring flux blueward of the Lyman limit. Due to the high sensitivity required to detect any continuum transmission at these redshifts, we also required a minimum signalto-noise ratio in the continuum near rest-frame 1285 Å of / ≥ 20 per 30 km s −1 interval. After rejecting BALs and objects with strong associated metal absorption and/or associated Ly damping wing absorption (typically with associated narrow metal lines), we selected 13 QSOs. These are listed in Table 2. The mean redshift in this sample is qso = 5.97. As described in Becker et al. (2019), the ESI spectra have a typical resolution of FWHM 45 km s −1 and were extracted in bins of 15 km s −1 , while the X-Shooter spectra have a typical resolution of FWHM 25 km s −1 in the VIS arm and were extracted in bins of 10 km s −1 . Individual spectra are plotted in Appendix A. The rarity of obvious transmitted flux blueward of the Lyman limit highlights the challenge of directly measuring mfp at these redshifts.

QSO redshifts
Following Worseck et al. (2014), we measured QSO redshifts from the apparent start of Ly forest absorption, forest . Five of our ∼ 6 objects also have precise systemic redshifts measured from either [C ] 158 m emission or narrow nebular Ly emission (see references listed in Table 2). An additional six QSOs 1 from Becker et al. (2019) have CO redshifts but were not included in the composite because they were at slightly higher redshifts or their spectra fell below our S/N requirement. We used the combined sample of eleven objects to estimate the error in our forest estimates, finding that forest was lower than the systemic redshift by an average of 180 km s −1 , with a standard deviation of 180 km s −1 . For LRIS, ESI, and X-Shooter QSOs without a systemic redshift measurement we offset the forest measurements by this amount to arrive at an adopted systemic redshift. The results are listed in Table 2. Given the decrease in the opacity of the Ly forest from ∼ 6 to 5 and the somewhat lower resolution of the red-side LRIS spectra versus the X-Shooter and ESI spectra, it is not entirely clear that the same offset should apply to our forest estimates at ∼ 5. On the other hand, 180 km s −1 corresponds to an offset 0.32 pMpc at = 5, which is relatively small compared to the statistical uncertainties in our measurement of mfp at that redshift arising from cosmic variance (see also Worseck et al. 2014). We therefore adopt this correction  to the forest measurements for the LRIS spectra. Redshifts for the GGG sample are adopted from Worseck et al. (2014).

Composite spectra
We created composite spectra from each of our three samples using the following procedure. We first shifted each spectrum to rest-frame wavelengths. We then divided each spectrum by its continuum flux measured over wavelengths where the flux from broad emission lines is minimal. For the GGG spectra we used the continuum flux near 1450 Å, following Worseck et al. (2014), while for LRIS, ESI, and X-Shooter we used the median flux over 1270-1380 Å. The choice of wavelength range for the continuum estimate has little impact on results because the normalization of the Lyman continuum profile is treated as a free parameter. For the LRIS, ESI, and X-Shooter spectra we corrected for residual zero-point errors by subtracting the median flux measured over a wavelength range expected to be free of transmitted flux. These wavelength ranges (750-800 Å rest frame for LRIS and 820-860 Å for ESI and X-Shooter) were verified to lie well blueward of where the fitted profiles reach zero flux (see Section 3.6). For the ∼ 6 data the lower wavelength bound was chosen to avoid the noisy edge of the X-Shooter VIS coverage, as well residuals from the 5577 Å skyline. The zero-point estimates for these spectra were subtracted prior to creating the composites; however, we do not require the corrections to be perfect. For the GGG sample, moreover, the wavelength coverage of the blue-side spectra does not provide a window where the zero point can be estimated safely blueward of the edge of the transmitted flux. In all cases, therefore, we include the zero point as a free parameter when fitting models to a composite.
For an alternate treatment of the zero-point errors in the GGG data see Worseck et al. (2014). Wavelength regions affected by skyline subtraction residuals were identified via peaks in the error arrays and masked. The ESI and X-Shooter spectra were also lightly median filtered using a 3pixel sliding window to reject spurious bad pixels. Mean composite spectra were then computed in bins of 170 km s −1 for GGG (similar to the binning used by Worseck et al. 2014) and 120 km s −1 for the LRIS and ESI + X-Shooter data. The results are shown in Fig. 2.
In the measurements described below we use bootstrap resampling to estimate the uncertainty in mfp . In each realization, we randomly select qso objects from each sample, with replacement, where qso corresponds to the numbers in Table 1. Before creating the new composite we add a random redshift offset to each spectrum (excluding those with [C ] 158 m or nebular Ly redshifts) drawn from a Gaussian distribution with = 180 km s −1 (see Section 2.5). As noted by Worseck et al. (2014), we found that the redshift errors produce an uncertainty in mfp that is small compared to the uncertainty from cosmic variance. We nevertheless include them for completeness. The bootstrap trials are also used to estimate the pixel-to-pixel errors in the flux, which we smooth using a polynomial fit over the wavelength range used to measure mfp . Additional sources of error are described in Section 3.3.
When fitting the composites we include wavelengths down to 826 Å for GGG, which is limited by the wavelength coverage of the data. For LRIS we fit down to 800 Å, while for ESI + X-Shooter we fit down to 820 Å. We note that wavelength range used to fit the composite overlaps with the wavelength range used to measure zero point offsets in the ESI and X-Shooter spectra. We find, however, that this choice does not have a significant impact on our results. The upper bound in wavelength is 910 Å for all composites, a choice we describe in Section 3.1.

Formalism
We measure a mean free path from the composite spectra using an approach based on the method first developed by Prochaska et al. (2009) and adapted by Worseck et al. (2014) to higher redshifts. The major change included here is to allow the ionizing opacity of the IGM to scale with the local photoionization rate. As demonstrated below, this change is necessary for extending the direct measurement method to ∼ 6.
The observed flux, obs , will be the mean intrinsic QSO spectral energy distribution, SED , attenuated by the effective Lyman series opacity of the foreground IGM, Lyman eff , and the Lyman continuum effective optical depth, LyC eff , Here, 0 is a zero-point correction that we include as a free parameter (see Section 2.6). We discuss the foreground Lyman series opacity in Section 3.4. The intrinsic SED blueward of the Lyman limit is modeled as a power law of the form SED = 912 912 Å − ion . The normalization 912 is treated as a free parameter that incorporates the intrinsic QSO SED, any Lyman continuum attenuation directly associated with the QSOs, and any relative flux calibration error between 912 Å 2 and the rest-frame wavelengths at which the individual QSO spectra are normalized. We adopt a nominal powerlaw exponent of ion = 0.5 (see Section 3.3). Given how rapidly Lyman eff and LyC eff evolve with wavelength, we find that our results for mfp are highly insensitive to this choice except as it impacts our calculations for the ionizing luminosity of a QSO (see Section 3.3).
The effective Lyman continuum opacity for a photon emitted at redshift qso that redshifts to 912 Å at redshift 912 will be where 912 ( ) is the Lyman continuum opacity at 912 Å at redshift (Prochaska et al. 2009). The wavelength dependence of the ionizing absorption cross-section is approximated here as ( ) ∝ −2.75 following O'Meara et al. (2013) and Worseck et al. (2014). Previous works at ≥ 3 have held 912 fixed when fitting a single QSO composite spectrum (Prochaska et al. 2009;Fumagalli et al. 2013;Worseck et al. 2014). The difficulty with this approach at > 5, however, is that mfp may become comparable to or smaller than a typical QSO proximity zone. If the ionizing flux from a QSO decreases the opacity of the IGM in its proximity zone then this will lead to a measurement of mfp that is biased high with respect to its value far from the QSO (see discussions in Worseck et al. 2014;D'Aloisio et al. 2018). This effect can be diminished by selecting QSOs that are relatively faint and hence have shorter proximity zones, as we have done for the LRIS sample. The measurement may still be biased, however, depending on the intrinsic value of mfp . At ∼ 6, moreover, mfp is expected to be significantly shorter than the typical proximity zone of any QSO bright enough to obtain a useful spectrum.
We therefore attempt to account for the proximity effect by modeling the impact of ionizing flux from a QSO on the Lyman continuum attenuation in its vicinity. We parametrize the dependence of the opacity on the local H ionization rate, Γ, as a power law of the form where bg 912 is the background opacity and Γ bg is the average background photoionization rate. 3 This form is motivated by analytic models of the IGM opacity (Miralda-Escudé et al. 2000;Furlanetto & Oh 2005), as well as radiative transfer simulations of Lyman limit systems (McQuinn et al. 2011). These studies suggest values of ∼ 2/3 at > 5, which has been adopted in recent models of the Ly forest opacity fluctuations at these redshifts (Davies & Furlanetto 2016;D'Aloisio et al. 2018;Nasir & D'Aloisio 2020). The uniform opacity model used by Worseck et al. (2014) corresponds to = 0. We discuss our priors on further in Section 3.5.
2 Throughout this paper we use 912 Å to represent the Lyman limit wavelength of 911.76 Å. 3 Here, "background" quantities refer to spatially averaged values in the absence of the QSO. We test the case where fluctuations in the UV background are present in Section 3.2. The local photoionization rate will be the sum of the background rate and the contribution from the QSO, which decreases with distance, giving Γ = Γ bg + Γ qso ( ). The Lyman limit opacity will therefore increase with distance from the QSO as Following Calverley et al. (2011), we characterize the ionizing luminosity of a QSO relative to the ionizing background according to the distance from the QSO, eq , at which Γ qso would be equal to Γ bg in the absence of any absorption or redshifting of ionizing photons from the QSO. We note that the actual distance at which Γ qso = Γ bg will tend to be less than eq due to absorption. Nevertheless, eq is a convenient parameter for helping to quantify how 912 is modified near a QSO. 4 For a QSO with luminosity 1450 at 4 We echo the discussion in Calverley et al. (2011) that eq differs from the observational definition of proximity zone size applied elsewhere at 6. eq is calculated directly from a QSO's ionizing spectrum and Γ bg . It is therefore effectively a prediction for the distance to which the ionizing flux from a QSO would dominate over the background in the absence of any attenuation. Observationally, in contrast, the proximity zone "size" at 6 is typically the distance from a QSO out to which the fraction of transmitted Ly flux exceeds 10% (e.g., Fan et al. 2006;Carilli et al. 2010), and is therefore a measure of where the total (QSO + background) ionization rate drops below the level required for the IGM to meet this transmission threshold. A = 6 QSO with 1450 = −27.0 would have eq = 11.4 pMpc for the nominal parameters given in Section 3.3. This is roughly twice the typical proximity zone size measured by Eilers et al. (2017) for QSOs near this luminosity. This suggests, perhaps not surprisingly, that at 6 rest-frame 1450 Å and a broken power-law continuum of the form the luminosity at 912 Å will be 912 = 1450 ( 912 / 1450 ) − UV and this distance will be eq = 912 0 Here, 0 is the H ionization cross-section at 912 Å. We calculate 1450 from the absolute magnitudes listed in Table 2. In Section 3.3 we calculate mean eq values for our samples and discuss constraints on Γ bg , UV , and ion . The ionizing flux from the QSO will be diluted geometrically and attenuated by Lyman continuum absorption, which increases with distance as Γ qso decreases. We therefore solve for Γ qso ( ) and 912 ( ) numerically under the assumption that 912 ( = 0) = 0. Specifically, we divide the line of sight into small steps of distance . For the first step we assume that Γ qso decreases purely geometrically, i.e., the ionizing flux from the QSO can dominate over the background out to distances that are significantly larger than those indicated by the extent of the observed Ly transmission.
In principle, 912 in equation (6) could be modified by an escape fraction, esc (e.g., Cristiani et al. 2016). For simplicity, however, we assume that the QSOs in our sample are roughly bimodal in terms of their escape fraction, having either esc ∼ 0 or 1, with esc independent of luminosity. Other than cases where there is an obvious, strong associated absorber such as DLA (see Section 2.4), we do not wish to bias our results by attempting to exclude QSOs with low esc . Fortunately, QSOs with esc = 0 will have zero flux blueward of 912 Å. Including these objects should therefore only rescale the mean Lyman continuum profile, which will be captured by the normalization parameter, 912 . Redshift errors may cause absorption from associated high-order Lyman series lines to be blended into the composite flux below 912 Å. We mitigate this by restricting our fits to rest < 910 Å, i.e., ∼600 km s −1 blueward of the nominal QSO redshifts.
In total, therefore, our model for the Lyman continuum flux includes five parameters, 912 , 0 , bg 912 , eq , and . The quantity we wish to obtain is the background mean free path that would be expected in the absence of the proximity effect. The mean free path is defined here to be the distance travelled by photons (emitted at a wavelength somewhat shorter than 912 Å) that would be attenuated by a factor of 1/ by Lyman continuum absorption. In order to calculate this quantity with the proximity effect removed, we recompute the effective Lyman continuum opacity by setting 912 = bg 912 in equation (2). Given the relatively short mean free path at these redshifts, we neglect any redshift evolution of . For all models we fix 912 = 1 and 0 = 0. The fiducial models were chosen to be similar to those measured from the data (see Sections 3.3 and 3.6). We then show how the profile varies with mfp , , and eq . Changes in and eq have a wavelength (radial) dependence that is significantly different from mfp because and eq mainly impact the transmission profile within the proximity zone. As expected, the relative importance of the proximity effect is larger at = 6.0, where a change of ±1/3 in or a factor of two change in eq produces a comparable change in the transmission profile as a factor of two change in mfp . Even so, these examples suggest that it is possible to measure mfp at = 6 given reasonable constraints on and eq , even when eq is a factor of ten larger than mfp . Our constraints on and eq are discussed further below.

Tests with mock spectra
Here we investigate how well our analytic model recovers the relevant parameters from mock spectra drawn from simulations. We refer the reader to Section 4 of D' Aloisio et al. (2018) for a description of the simulations. In summary, we assign QSOs with luminosities taken from the GGG, LRIS, and ESI + X-Shooter samples to the most massive halos in a cosmological hydrodynamics simulation with box = 200 ℎ −1 Mpc and gas = dm = 2048 3 gas and dark matter resolution elements. The hydrodynamics simulation was run with a modified version of the code of Trac & Pen (2004). The QSO halos masses range from 1.3 to 8.0 × 10 12 ℎ −1 M . The QSO luminosities at 1450 Å rest-frame are taken from Table 2. We compute the ionizing luminosity of each QSO assuming a broken power-law of the form given by Lusso et al. (2015), which is similar to what we assume for the data (see Section 3.3).
One QSO is populated in the box at a time and we use the attenuation model of Davies & Furlanetto (2016) to compute the Γ and mfp fields in the box. These iterative calculations include galactic sources, spatially varying mfp and the backreaction of Γ on local mfp values. The background ionization rates are Γ bg = 5 × 10 −13 s −1 and 1 × 10 −13 s −1 at = 5.2 and 6, respectively, which are somewhat different than the values we use when fitting the data (see below). For the Γ and mfp computations we use uniform grids with 64 3 cells. We compute 1, 000 transmission profiles along random sight lines emanating from each QSO in a given sample. We then construct 1, 000 mock composite spectra by averaging over the QSOs in the sample.
We test how well our fitting approach recovers the "true" mfp and values by fitting our model to the mock composite spectra. We compute eq for each QSO using the Γ bg values quoted above. For consistency, we use the same QSO SED that was used to compute the ionizing luminosities for the mock sample. The mean eq values are 7.4, 3.5, and 16.9 pMpc for the mock GGG, LRIS, and ESI + X-Shooter samples, respectively, which we adopt when fitting the models. These values are somewhat larger than the values we compute for the data (see Section 3.3), mainly due to the difference in Γ bg . Our fits to the mock composites have three free parameters: 912 , , and 912 . For comparison, we also fit a constant opacity model that ignores the QSO proximity effect ( = 0). We employ a chi-squared approach assuming equal variance in each wavelength bin. The mocks do not include foreground Lyman series absorption or variations due to intrinsic QSO SEDs. They therefore allow us to determine how well the mfp and values are recovered under ideal circumstances.
Fits to the mock transmission profiles are shown in Fig. 4. In each case where we include the proximity effect in the fit we recover the correct mfp to within 17%. This is true even in the "short" ( mfp = 1.6 pMpc) case at = 6, where eq is a factor of ten larger than mfp . We also recover the correct to within ∼0.1 in all cases except the mfp = 9.2 pMpc, = 0.33 case with the mock LRIS composite, where the impact of the proximity effect is weakest. In contrast, ignoring the proximity effect can produce a significant overestimate of the mean free path (and overestimates of the normalization, a fact that may be evident when fitting high-/ composites). For mfp = 4.6 pMpc and = 0.67, the mfp values returned for the LRIS and GGG mocks are too large by factors of 1.5 and 1.9, respectively. This suggests that accounting for the proximity effect may be necessary even for fainter QSOs, depending on the true value of mfp . Errors for the constant opacity model are largest at = 6, with mfp overestimated by up by factors of two to four.
In summary, we find that reasonable estimates of mfp can be obtained even when the mean free path is much shorter than the proximity zone size provided that the proximity effect is taken into account. Fitting a constant opacity model to Lyman continuum profiles at > 5, in contrast, can lead to significant overestimates of the mean free path, even for samples of relatively faint QSOs. In principle, at least, it is also possible to recover the scaling of Lyman continuum opacity with local ionization rate. Directly constraining requires extremely good data, however, a point we return to below.
An important caveat is that the simulations on which we validated our technique for simultaneously fitting mfp and do not include dynamical effects that are especially relevant if reionization ended near = 6. Park et al. (2016) andD'Aloisio et al. (2020) found that impulsive changes to the UVB (e.g. reionization or a QSO turning on suddenly) shape the density structure of the IGM over Δ ∼ 100 Myr through the interplay between self-shielding and hydrodynamic response of the gas to photoheating. One implication raised by D' Aloisio et al. (2020) is that the dependence of mfp on Γ may be more complex than can be captured with a universal power law. The simulations also assume an infinite QSO lifetime. If the QSOs are much younger than the ∼ 100 Myr relaxation timescale of the optically thick absorbers that set mfp , another distinct possibility is that the local mean free paths have not had sufficient time to respond to the enhanced UV intensities. In this case, the proximity effect would be less apparent in the measurements of mfp .

eq values for observed QSOs
The eq estimates for our QSOs are derived from observational constraints on the metagalactic hydrogen ionization rate and the mean SED of high-redshift QSOs. Similar to previous works (e.g., Becker & Bolton 2013), we estimate Γ bg based on the mean intergalactic Ly transmission at these redshifts. Our nominal evolution in the mean Ly transmission, described in Section 3.4, corresponds to T Ly = 0.14 at = 5.1 and T Ly = 0.0072 at = 6.0. These values are based on measurements made from QSO spectra well outside the proximity zone (see below). We use a hydrodynamical simulation to translate these T Ly values into Γ bg estimates by rescaling the simulated UV background such that the mean Ly transmission of the simulation box matches observations. Specifically, we use the 40 ℎ −1 Mpc box with 2 × 2048 3 particles (40-2048) from the Sherwood simulation suite (Bolton et al. 2017), whose IGM temperatures over 5 < < 6 are broadly consistent with existing measurements (Bolton et Gaikwad et al. 2020). This procedure yields Γ bg 7 × 10 −13 s −1 and 3 × 10 −13 s −1 at = 5.1 and 6.0, respectively. The uncertainties affecting Γ bg , including those related to T Ly , the temperature-density relation, and numerical effects, are similar to those in Becker & Bolton (2013). We therefore adopt a similar overall error on our Γ bg estimates, namely ±0.15 dex. For the QSO SED in equation (5) we adopt UV = 0.6 ± 0.1 and ion = 1.5 ± 0.3 ( ion = 0.5 ± 0.3). The choice of UV is taken from fits to composite QSO spectra by Lusso et al. (2015), and is generally consistent with other similar works (Vanden Berk et al. 2001;Shull et al. 2012;Stevans et al. 2014). Here we adopt a larger error than found by Lusso et al. (2015) in order to allow for greater sample variance. Our choice of ion is broadly consistent with fits to composite spectra from Telfer et al. (2002), Stevans et al. (2014), and Lusso et al. (2015) (though see Scott et al. 2004, who find a harder ionizing slope for low-redshift AGN). For the above parameters and the 1450 values listed in Table 2 we calculate mean eq values of 6.4, 3.0, and 11.1 pMpc for the GGG, LRIS, and ESI + X-Shooter samples, respectively.
For each bootstrap composite that is used to estimate the uncertainty in mfp (see Section 2.6) we randomly sample the above error distributions for Γ bg , UV , and ion and propagate these into the estimates of eq for each object. We then recompute the mean eq based on the objects in that bootstrap sample. The same value of ion is used to model the Lyman continuum transmission profile for a given bootstrap trial. When fitting the GGG and LRIS profiles simultaneously, the same random realizations of Γ bg and the QSO spectral indices are applied to both data sets. For reference, the 68% (95%) ranges of the mean eq values from the bootstrap trials are 5.4-7.8 (4.5-9.3) pMpc, 2.5-3.6 (2.1-4.3) pMpc, and 10.3-15.1 (8.6-18.2) pMpc for the GGG, LRIS, and ESI + X-Shooter samples, respectively.

Foreground Lyman series transmission
The Lyman series opacity at an observed wavelength obs < (912 Å) (1 + qso ) will include foreground contributions from all Lyman series lines, Lyman eff where eff ( ) is the effective opacity of transition at redshift , (1 + ) = obs , and is the rest-frame wavelength of transition . We compute Lyman eff using the 40-2048 Sherwood simulation described above. The simulation outputs are spaced in redshift intervals of Δ = 0.1, with 5000 lines of sight drawn from each output. At each simulation redshift we first compute baseline Ly optical depths by rescaling the native simulated Ly optical depths to reproduce the observed mean IGM Ly transmission. We note that we are computing the Lyman series transmission from an optically thin simulation that does not include elements such as galactic outflows and self-shielded gas that may modify the neutral hydrogen density distribution, and hence impact the ratio of eff / for high-order lines. The numerical resolution of the simulations may also have an effect. We tested the numerical resolution using the 40-1024 run from the Sherwood suite, which also uses a 40 ℎ −1 Mpc box but is a factor of eight lower in mass resolution than our fiducial 40-2048 run. Using the lower resolution run increased the total Lyman series transmission over 890−912 Å in the rest frame by 2% (10%) for QSOs at = 5.1 (6.0). We tested the impact of galaxy physics using the 40-1024-ps13 runs from Bolton et al. (2017), which include a subgrid implementation of star formation and galactic outflows from Puchwein & Springel (2013). These decreased the transmission relative to the 40-1024 run by 4% (∼3%) at = 5.1 (6.0). We also tested the impact of self-shielding using a version of the 40-1024-ps13 run in which self-shielding was added in post-processing following Rahmati et al. (2013) at < 5 and Chardin et al. (2018) at > 5. This decreased the mean transmission by a further 3% (2%) at = 5.1 (6.0). Fortunately, in all cases the effect was mainly to rescale the transmission below 912 Å and not to change the shape of the profile in a way that would significantly impact our mfp measurements. These effects may nevertheless need to be considered in future works.
An additional factor here is the QSO proximity effect. We include the proximity effect for each Lyman series line following the same numerical approach used to compute the Lyman continuum opacities. For a given combination of bg 912 and we compute eff as a function of wavelength over a grid in QSO redshift and eq , interpolating between simulation redshifts as needed. For each composite or bootstrap sample we then compute Lyman eff ( obs ) individually for each QSO using equation (9). We then compute the transmission as T Lyman = exp − Lyman eff , and average the transmission over all lines of sight.
In Fig. 5 we plot the Lyman series absorption for different combinations of mfp , , and eq at = 5.1 and 6.0. At = 5.1 the transmission is not strongly affected by mfp or because the decrease in Γ tot with distance from the QSO is mainly driven by geometric dilution. Including the proximity effect increases T Lyman by a factor of ∼1.3 at rest-frame 912 Å for eq = 5 pMpc, similar to the mean value in the GGG sample. It also modifies the shape of the Lyman series transmission with respect to the no proximity effect ( eq = 0) case. At = 6.0 the effect is even larger, with T Lyman increasing at 912 Å by a factor of 2.5 for eq = 5 pMpc, similar to the mean value for the ESI + X-Shooter sample. There is also a greater dependence on mfp and . We find, however, that our final results are not highly dependent on the choice of mfp and used for the Lyman series transmission. When computing T Lyman , therefore, we hold these parameters fixed at the nominal values shown in Fig. 5, which are comparable to our best-fit results.

Priors on
The scaling of 912 with Γ is highly uncertain, especially at the high redshifts that are relevant for this study. From a theoretical viewpoint, the value of is tied to the shape of the gas density distribution function near the self-shielding threshold. Adopting the Miralda-Escudé et al. (2000) model of IGM opacity, and assuming that the density profile of a typical self-shielding absorber is isothermal, it can be shown that 912 ∝ Γ −2/3 , i.e. = 2/3 (Furlanetto & Oh 2005;McQuinn et al. 2011). Indeed, this value has been adopted in recent models of the fluctuating UVB at > 5 (e.g. Davies & Furlanetto 2016 . It should be noted, however, that the radiative transfer in their study was applied in post processing to absorbers extracted from hydrodynamic simulations. This approach misses the effect of the UVB on the density structure of the absorbers. More recently, D'Aloisio et al. (2020) used fully coupled radiation hydrodynamics simulations to study self-shielding systems (see also Park et al. 2016). Their findings suggest a more complex dependence of 912 on Γ owing to the interplay between self-shielding and the hydrodynamic response of the gas to photoheating, which occurs on a time scale of hundreds of Myr. We can nonetheless examine their gas density distribution functions in an attempt to gain insight into (see their Fig. 5). At densities well above self-shielding, the probability distribution of Δ is reasonably approximated by ∝ Δ −1.8 , where Δ is the gas density in units of the cosmic mean. Applying the analytic arguments of Furlanetto & Oh (2005) and McQuinn et al. (2011) yields a milder scaling of ≈ 0.33. This would be the scaling for a short time after a bright source turned on suddenly, before the gas had time to react to the impulse. We note, however, some important caveats which suggest that may be larger than this. First, the (Δ) of D'Aloisio et al. (2020) are generally not well-described by a power law near self-shielding. Indeed, Δ 3 appears to flatten at densities closer to self-shielding, implying a stronger dependence of 912 on Γ. Secondly, the dependence would likely evolve as the density structure of the gas readjusted to the changing UVB. Based on these considerations, we argue here that = 0.33 may serve as an approximate lower limit. On the other hand, = 1 is the scaling for the case of a uniform IGM in photoionization equilibrium. This limit is approached if the opacity is dominated by diffuse gas near the mean density, rather than over-dense peaks. In our fits we adopt a nominal value of = 0.67 and a range = 0.33-1.0 with a flat prior from which we randomly sample when performing bootstrap trials. We also perform fits with fixed to 0.33, 0.67, and 1.0. In principle, one can measure directly from the data. Even with good constraints on eq this is difficult, however, because at = 5.1 the dependence of the transmitted flux on is relatively weak unless the mean free path is short (Fig. 3), while at = 6.0 the data are too noisy to distinguish between variations in mfp and . In a joint fit to the GGG and LRIS data we find = 0.56, consistent with theoretical expectations, but with a 68% (95%) confidence range of 0.20 to 1.20 (-0.06 to 2.28). Much of this parameter space is strongly disfavored on theoretical grounds, as described above. The choice of ultimately has little impact at = 5.1. Setting = 0.33 (1.0) increases (decreases) our nominal result by 8% (6%). The impact of is more significant at = 6.0, where the proximity effect is more pronounced. There, setting = 0.33 (1.0) increases (decreases) our nominal result by 69% (68%). This represents a substantial portion of our error budget at = 6.0. In future works it may be possible to better constrain directly from the data.

Fits to the data
At = 5.1 we fit the GGG and LRIS composites individually as well as jointly. For our nominal results we use = 0.67, as noted above, and hold the mean eq for each composite fixed to the values given in Section 3.3. We also include the foreground Lyman series The histogram in each panel shows the observed flux blueward of the Lyman limit for the sample indicated, normalized as in Fig. 2. Solid lines show the best-fitting model profiles, which for GGG and LRIS is a simultaneous fit to both data sets. Dark and light shaded regions show the 68% and 95% intervals, respectively, spanned by fits to the bootstrap trials.
transmission described in Section 3.4. We then fit for bg 912 , which is used to calculate mfp , along with 912 and 0 separately for each composite. In the bootstrap trials we draw randomly from the flat distribution over [0.3, 1.0] while the mean eq is varied according to the procedure outlined above. From the individual fits we obtain mfp = 8.85 +1.63 −1.31 pMpc (68% confidence intervals assuming a flat prior on ) from the GGG data and 11.64 +4.12 −3.63 pMpc from the LRIS data. The results are thus highly consistent with one another within the errors. From the joint fit we obtain mfp = 9.09 +1.62 −1.28 pMpc, which we adopt as our nominal result at = 5.1. At = 6.0 we measure mfp = 0.75 +0.65 −0.45 pMpc. The nominal fits along with the ranges spanned by bootstrap trials are shown in Fig. 6. The cumulative probability density functions for mfp at the two redshifts are shown in Fig. 7. The main results are summarized in Table 3, where we also give results for fixed values of .
Our value of mfp at = 5.1 is consistent with the results from Worseck et al. (2014). This suggests that mfp at this redshift is large enough that the impact of the QSO proximity effect is relatively modest, even for the brighter GGG sample where eq 0.7 mfp . Indeed, if we neglect the proximity effect by setting = 0 and eq = 0 for both the Lyman continuum opacity and the foreground Lyman series, emulating the approach of Worseck et al. (2014), our result for mfp increases by only 12% for the GGG composite and remains essentially unchanged for the LRIS composite. This is somewhat less than the bias found with mock spectra in Section 3.2 because the errors in the Lyman continuum modeling are partially offset by errors in the Lyman series modeling. The difference increases to 28% for the GGG sample if we include the proximity effect in  the foreground Lyman series transmission but not in the Lyman continuum, a scenario closer to Figure 3, where accurate modeling of the Lyman series is assumed. At = 6.0, in contrast, the mean value of eq is a factor of 15 larger than our value of mfp , making it critical to take the proximity effect into account. In this case, fully neglecting the proximity effect increases our mfp measurement by a factor of 2.9 above our nominal = 0.67 value. If we attempt to emulate the mock trials by including the proximity effect in the Lyman series absorption but not in the Lyman continuum then our result for mfp increases by a factor of 3.6. This is consistent with the bias expected from the mock trials in Section 3.2, and emphasizes the importance of properly accounting for the proximity effect at ∼ 6.
In Fig. 8 we plot our mfp values as a function of redshift, along with measurements from the literature (Prochaska et al. 2009;Fumagalli et al. 2013;O'Meara et al. 2013;Worseck et al. 2014;Lusso et al. 2018). The Lusso et al. (2018) value at = 2.44 is their fit to the data from O' Meara et al. (2013). Lusso et al. (2018) find somewhat lower values of mfp at ∼ 2 towards QSO pairs, potentially due to an increased incidence of optically thick absorbers in pair environments. We note that Romano et al. (2019) measured the mean free path towards QSOs at ∼ 4. They find values that are ∼10-20% higher than those of Prochaska et al. (2009) andWorseck et al. (2014) over the same redshifts. In trials using the two lowerredshift GGG composites from Worseck et al. (2014) we found that this discrepancy is well explained by the lack of foreground Lyman series absorption in the Romano et al. (2019) analysis. Worseck et al. (2014) fit a power law of the form mfp ( ) ∝ (1 + ) −5.4 over 2.44 < < 5.16 (dotted line in Fig. 8). Extrapolating this fit out to = 6 overshoots our nominal ESI + X-Shooter measurement by a factor of six, and is excluded by the data with >99.99% confidence. We therefore find strong evidence that the evolution of mfp ( ) with redshift steepens at 5. This steepening is broadly consistent with the results of Songaila & Cowie (2010) based on their measurements of discrete Lyman limit absorbers towards QSOs over 5 < < 6.

Implications for reionization
Our measurements are consistent with a low value of mfp at = 6 and a rapid increase from = 6 to 5. Taken at face value, perhaps the most interesting possibility is that this evolution is tied to the end of reionization. In Fig. 9 we compare our measurements to predictions for mfp ( ) from different reionization models. We begin with the simplistic models in D' Aloisio et al. (2020), which employ results from a suite of radiative hydrodynamics simulations of the ionizing photon sinks at > 5. The dotted curve shows a model in which reionization ended long before = 6 such that the IGM has had sufficient time to relax hydrodynamically. This model predicts a redshift evolution of mfp ∝ (1 + ) −5.4 and a mfp ( = 6) that is a factor of ∼7 longer than our measurement. It is worth noting that this model assumes only the ΛCDM cosmology and a constant UVB intensity; yet it yields a redshift evolution for mfp that is identical to the empirical fit of Worseck et al. (2014). The fully relaxed model is inconsistent with our = 6.0 measurement at the 99.9% level ( (< mfp ) = 0.999). For comparison, the solid curves show the "rapid" and "gradual" reionization models of D' Aloisio et al. (2020) wherein reionization is 50% complete at = 7.3 and 9.1, respectively, and ends at = 6. Although mfp at = 6 is lower than in the fully relaxed models, the data are still inconsistent at the 98-99% confidence levels.
It is also possible that reionization ended later than = 6, a scenario that has been proposed recently to explain the large scatter in the > 5 Ly forest opacity (Kulkarni et al. 2019;Keating et al. 2020a,b;Nasir & D'Aloisio 2020;Choudhury et al. 2020;Qin et al. 2021). The dashed curves in Fig. 9 show the "Low CMB " and "Hot  (2020). The dotted line shows the expected evolution if the IGM reionized early enough that the absorbers have had time to fully relax hydrodynamically by = 6. Blue (upper) and red (lower)solid lines show their "gradual" and "rapid" reionization models wherein reionization is 50% complete at = 9.1 and 7.3, respectively, and complete by = 6. Dashed lines are from Keating et al. (2020b). The orange (upper) line shows their "High CMB " wherein reionization is 50% complete near 8.4 and ends near 5.3. The cyan (middle at < 5.5) and red (lower at < 5.5) dashed lines show their "Low CMB " models wherein reionization still ends 5.3 but is 50% complete at 6.7.
Low CMB " models of Keating et al. (2020b), wherein reionization is 50% complete at 6.7 and ends at 5.3. In these models the IGM at = 6.0 is still ∼20% neutral. We also plot their "High CMB " model wherein reionization ends at the same redshift but is 50% complete at 8.4. In this model the IGM at = 6.0 is ∼8% neutral. The High CMB model is excluded at the 99% level. The Low CMB models are more consistent with our measurement at = 6.0, although the data still prefer a lower mfp at the 97% confidence level.
We note that mfp evolves rapidly near = 6 in all of these reionization models, and that they therefore become more consistent with the data if they are shifted slightly in redshift. For example, shifting the models by Δ = −0.2 decreases (< mfp ) at = 6.0 to 0.97 for the D'Aloisio et al. (2020) "rapid" model and 0.86 for the Keating et al. (2020b) Low CMB model. The low value of mfp we measure at = 6.0 may therefore suggest that reionization occurs even later than these models propose.
We further note that the tension with existing models may be reduced if is near the low end of our adopted range. Our nominal mfp value at = 6.0 is a factor of two higher for = 0.33 than for = 0.67 (see Table 3), a result that comes from attributing less of the transmission in Figure 6 to the proximity effect. Moreover, (< mfp ) for the Low CMB model at = 6.0 decreases to 0.91 when we hold fixed to 0.33. It is possible, therefore, that reconciling the reionization history with our measurements of mfp may require the ionizing sinks near ∼ 6 to be less sensitive to photoionization effects than some models assume (for further discussion, see D'Aloisio et al. 2020).

Ionizing emissivity
We can use our estimates of Γ bg and mfp to infer the ionizing emissivity at = 5-6. Here we use the local source approximation, which neglects the redshifting of ionizing photons (e.g., Schirber & Bullock 2003;Kuhlen & Faucher-Giguère 2012). This is a reasonable choice given the short mean free path at > 5 (e.g., Becker & Bolton 2013). Under this approximation the comoving ionizing emissivity is given by Here, s is the slope of the ionizing spectrum of the sources ( s ∝ − s ), and bg is the slope of the ionizing background after filtering through the IGM. If the column density distribution of H absorbers producing most of the Lyman continuun opacity is a power law of the form (  The results at = 5.1 are similar to the values of ion over 2.4 < < 4.75 found by Becker & Bolton (2013), and suggest that the ionizing emissivity over 2 < < 5 may change relatively little over this period even as the source populations of star-forming galaxies and AGN evolve considerably. In contrast, although the errors are large and we have ignored possible fluctuations in the ionizing background, the emissivity at = 6.0 is potentially significantly higher. If confirmed, this would suggest that the mean production efficiency and/or escape fraction of ionizing photons is higher for sources at 6 than for sources at lower redshifts. The nominal value of ion at = 6.0 corresponds to ∼17 ionizing photons per hydrogen atom per Gyr, a rate that may help explain how reionization could have been completed in only a few hundred Myr.

Caveats and future work
Finally, we note that this work has some limitations. Our measurement at = 6.0 is based on a relatively small sample of 13 QSOs, and within this sample there are clearly outliers in terms of Lyman continuum transmission. The spectrum of SDSS J0836+0054, for example, shows discrete transmission peaks down to rest-frame 870 Å (i.e., 22 pMpc from the QSO), whereas none of the other ∼ 6 objects shows obvious transmission below 900 Å. While a skewed distribution of free paths along individual lines of sight is expected (e.g., Romano et al. 2019), and while this particular QSO is the brightest (and lowest-redshift) one in our ESI + X-Shooter sample, a larger sample at ∼ 6 would help to characterize the spatial variations in ionizing opacity near the end of reionization. Given the rapid increase in mfp between = 6 to 5, and the different evolutions predicted by the models over this redshift range (e.g., Figure 9), it also clearly of interest to constrain mfp near = 5.5.
In term of the modeling, the uncertain scaling of 912 with Γ has significant implications for mfp at ∼ 6, as discussed above. We also note that some of the formalism we applied herein assumes an ionized IGM. If reionization is incomplete at ∼ 6, then the mfp we measure at that redshift may correspond to the mean opacity only within the ionized phase provided that the ionized regions surrounding bright QSOs are larger than the proximity zone size ( eq ). The tests presented in Section 3.2 suggest that our approach should be robust to the UVB fluctuations expected near the end stages of reionization. Additional trials with more realistic late reionization simulations, however, would help to clarify how well these tools can be applied when the IGM is partly neutral.
Finally, consistent with previous works, we have not attempted to model the foreground Lyman series transmission in a fully selfconsistent way. Although we do not expect this to significantly impact our mfp results, as discussed above, simultaneously fitting the Lyman series and Lyman continuum transmission may provide insight into the properties (e.g., the H column density distribution) of the absorbers that dominate the ionizing opacity at these redshifts.

SUMMARY
In this work we measure the mean free path of ionizing photons at 5-6 using composite QSO spectra. We introduce a fitting approach that accounts for the QSO proximity effect by modeling the change in ionizing opacity with the local photoionization rate. This is also the first work to extend direct measurements of mfp to ∼ 6, where they are sensitive to the ionizing opacity near the end of reionization. At = 5.1 we measure mfp = 9.09 +1.62 −1.28 pMpc (68% errors) from a combination of bright QSOs from the GGG survey and fainter QSOs observed with LRIS. This is consistent with results from the GGG sample alone obtained by Worseck et al. (2014), who did not attempt to account for the proximity effect. This suggests that mfp is sufficiently long at ∼ 5 that the proximity effect does not greatly impact the transmission of Lyman continuum photons in QSO spectra. At = 6.0 we measure mfp = 0.75 +0.65 −0.45 pMpc using spectra from ESI and X-Shooter. In contrast to lower redshifts, we find that neglecting the proximity effect here can bias the result high by a factor of two or more. Our value lies well below extrapolations from lower redshifts, and suggests that the mean free path evolves rapidly over 5 < < 6. A short mean free path at = 6.0 and a rapid increase from = 6 to 5 are qualitatively consistent with models wherein reionization ends at ∼ 6, or even later (e.g., Kulkarni et al. 2019;Keating et al. 2020a,b;Nasir & D'Aloisio 2020), but disfavor models wherein reionization ended early enough that the IGM has had time to fully relax by ∼ 6 (see D'Aloisio et al. 2020).
Models with later and more rapid reionization (i.e., the "rapid" model of D'Aloisio et al. 2020 and the "Low CMB " models of Keating et al. 2020b) fall closest to our mfp measurements, yet our value at = 6.0 lies below even models wherein the IGM at this redshift is still ∼20% neutral (Keating et al. 2020b). This may indicate that the end of reionization occurred even later than previously thought. Alternatively, the models may be missing some of the absorption systems that limit the mean free path near the end of reionization. Further work will help to clarify how strongly the reionization history can be constrained by mean free path measurements such as the ones in this work. and the anonymous referee for helpful comments and suggestions. GDB, HMC, and YZ are supported by the National Science Foundation through grant AST-1751404. HMC is also supported by an NSF GRFP through grant DGE-1326120. AD is supported by HST grant HST-AR15013.005-A and NASA grant 19-ATP19-0191. JSB is supported by STFC consolidated grant ST/T000171/1. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership between the California Institute of Technology and the University of California; it was made possible by the generous 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 indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This research also made use of the Keck Observatory Archive ( Computational models were made possible by NSF XSEDE allocation TG-AST120066. The Sherwood simulations were performed using the Curie supercomputer at the Tre Grand Centre de Calcul (TGCC), and the DiRAC Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.

DATA AVAILABILITY
The raw data underlying this article are available from the Gemini archive at https://archive.gemini.edu, the Keck Observatory Archive at https://www2.keck.hawaii.edu/koa/public/koa.php, and the VLT Archive at http://archive.eso.org/cms.html. Reduced GGG data are available at the CDS via http://cdsarc.u-strasbg.fr/vizbin/qcat?J/MNRAS/445/1745. Other reduced data are available on reasonable request to the corresponding author.   . ESI and X-Shooter spectra used in this work. Panels are labeled with the QSO name and redshift. For each QSO we plot flux per unit wavelength normalized by the continuum flux measured over rest-frame 1270-1380 Å. The spectra have been median filtered using a 3-pixel sliding window. Vertical lines mark the Lyman limit wavelength in the rest frame of the QSO.