Robustness of direct measurements of the mean free path of ionizing photons in the epoch of reionization

Measurements of the mean free path of Lyman-continuum photons in the intergalactic medium during the epoch of reionization can help constrain the nature of the sources as well as sinks of hydrogen-ionizing radiation. A recent approach to this measurement has been to utilize composite spectra of multiple quasars at $z\sim 6$, and infer the mean free path after correcting the spectra for the presence of quasar proximity zones. This has revealed not only a steep drop in the mean free path from $z=5$ to $z=6$, but also potentially a mild tension with reionization simulations. We critically examine such direct measurements of the mean free path for biases due to quasar environment, incomplete reionization, and quasar proximity zones. Using cosmological radiative transfer simulations of reionization combined with one-dimensional radiative transfer calculations of quasar proximity zones, we find that the bias in the mean free path due to overdensities around quasars is minimal at $z\sim 6$. Patchiness of reionization at this redshift also does not affect the measurements significantly. Fitting our model to the data results in a mean free path of $\lambda_{\mathrm{mfp}}=1.49^{+0.47}_{-0.52}$~pMpc at $z=6$, which is consistent with the recent measurements in the literature, indicating robustness with respect to the modelling of quasar proximity zones. We also compare various ways in which the mean free path has been defined in simulations before the end of reionization. Overall, our finding is that recent measurements of the mean free path appear to be robust relative to several sources of potential bias.


INTRODUCTION
The epoch of reionization marks the era during which the neutral hydrogen in the Universe became largely ionized, most likely due to UV photons from stars in galaxies across the cosmic volume, with potentially an additional minor contribution from quasars.To understand how reionization happened, we need to know how the photons that escaped from the various sources of radiation interacted with the intergalactic medium (IGM) to ionize the universe.An important characteristic of the propagation of these ionizing photons is determined by their mean free path (MFP), defined as the average distance a photon travels before getting absorbed (Rybicki & Lightman 1985).In a homogeneous IGM, the MFP increases with the background UV photoionization rate, as this lowers the IGM opacity.More generally, the MFP  MFP varies with the background photoionization rate Γ bg as a power law (Miralda-Escudé et al. 2000), and, furthermore, both of these quantities can be expressed as functions of the cosmological emissivity  of ionizing radiation.Knowing the MFP therefore puts constraints on the ionizing sources as well as the sinks of ionizing photons.At higher redshifts, the rate of increase of the MFP with redshift indicates how rapid the progress of reionization is.Measure-★ E-mail: sindhu@theory.tifr.res.inment of the MFP of hydrogen-ionizing photons is thus important for characterizing reionization.
During reionization, the sources of radiation carve out regions of ionized hydrogen in the IGM around themselves.These ionized regions later coalesce, at which point a background of ionizing radiation gets established throughout the Universe (Gnedin & Kaurov 2014;Bauer et al. 2015).This is considered to be the end of reionization and is thought to occur at  ∼ 5.3 (Kulkarni et al. 2019;Bosman et al. 2022;Zhu et al. 2022).At lower redshifts, in the post-reionization Universe, the MFP is dominated by the residual neutral hydrogen systems in the otherwise fully ionized IGM.These systems retain neutral hydrogen due to self-shielding thanks to their high density.These appear as high-column-density absorbers of the ionizing radiation background.The corresponding MFP is then set by the average spacing between such absorbers, obtained by counting the number of optically-thick absorbers or Lyman-limit systems (LLS) along the lines of sight in quasar spectra (Storrie-Lombardi et al. 1994;O'Meara et al. 2007O'Meara et al. , 2013)).At higher redshifts, before reionization has ended, however, the boundaries of the ionized regions also play an important role in setting the MFP.At these times, the MFP is the average distance between LLSs or the typical size of ionized regions, depending on which is smaller (Shukla et al. 2016;Madau 2017).
Observationally, measurements of the MFP have been inferred using direct and indirect means.One way to obtain a direct measurement of the MFP is to use stacked quasar spectra bluewards of restframe 912Å.The MFP is then computed as the distance at which the effective Lyman-continuum optical depth becomes unity (Prochaska et al. 2009;Fumagalli et al. 2013;O'Meara et al. 2013;Worseck et al. 2014).This has resulted in a measurement of  mfp = 10.3±1.6 pMpc at  = 5.16, with the MFP increasing as (1 + ) −5.4 down to  = 4.56.At higher redshifts, the MFP can however become comparable or smaller than the size of the proximity zone of the quasars whose spectra are used in the stacks.This biases the MFP towards lower values (Worseck et al. 2014;D'Aloisio et al. 2018).Becker et al. (2021, B21) reported a measurement of the MFP in which they sought to correct the bias due to the quasar proximity zones by modifying the direct measurement method to account for the excess flux due to the quasar.The resultant value of the MFP inferred by B21 is 0.75 +0.65  −0.45 pMpc at  = 6.The MFP thus shows a sharp drop at  > 5 relative to the (1 + ) −5.4 behaviour seen at redshifts lower than 5.Moreover, the value of the MFP at  = 6 is smaller than the value in several of the latest reionization simulations by a factor of 2 or more (Keating et al. 2020a).Davies et al. (2021) discussed the implications of this tension for the ionizing budget of galaxies to argue that a shorter MFP requires an ionizing emissivity that is up to six times larger than the typically assumed values.Cain et al. (2021) found that the short MFP reported by B21 is consistent with a late reionization scenario powered by fainter galaxies with a high ionizing photon production efficiency (Lewis et al. 2022;Garaldi et al. 2022).More recently, Zhu et al. (2023b) have measured the MFP for quasars at redshifts between 5.1 <  < 6 using the B21 method and found the value at redshift 6 to agree with that measured by B21.They also find the decrease in the MFP to be steeper with increasing redshift, with a sharp drop in MFP by nearly 75% between  ∼ 5.6 and 5.9.Such a sharp decline in the MFP points towards a rapid end to reionization.Another direct measurement technique is to define the MFP by averaging free paths measured from the distribution of absorbers along individual quasar sightlines (Romano et al. 2019;Songaila & Cowie 2010).This approach suggested a milder evolution of the MFP with redshift, between 3 <  < 6, when compared to the direct stacking method.Bosman (2021) have extended this approach to detect absorption due to LLS in the six lowest-order Lyman-series transitions and put a lower limit based on the average of individual free paths defined this way to be  MFP > 0.31 pMpc at redshift 6, consistent with the measurement of B21.
Indirect measurements of the MFP have been obtained by comparing the observed Ly opacity of the IGM with that in numerical simulations.The mean and the scatter of the cumulative distribution function (CDF) of the effective Ly opacity has been used to simultaneously constrain both the MFP and the photoionization rate for a given emissivity (e.g., Gaikwad et al. 2023;Wolfson et al. 2023, andDavies et al. 2023, in preparation).Overall, the indirect MFP measurements based on the Ly opacity (Gaikwad et al. 2023) are consistent with the direct MFP measured from the quasar stack beyond the Lyman limit, within the error bars (Becker et al. 2021;Zhu et al. 2023b).However, the slope of evolution of the MFP with redshift differs beyond redshift  ∼ 5.5, with the direct measurements yielding a steeper slope than the indirect measurements.A gradual slope agrees better with the reionization simulations of D' Aloisio et al. (2018), Keating et al. (2020a), andCain et al. (2021), as compared to the steeper slope seen in the models of Lewis et al. (2022) and Garaldi et al. (2022).Possible reasons for the mismatch between the indirect and direct measurements of the MFP could be the assumptions made by the two methods.B21 assume an analytic model to measure the effective optical depth, where the opacity  is proportional to the photoionization rate as Γ −  .In the QSO proximity zone, the total Γ is a sum of both Γ bg and Γ qso , the latter predominantly varying with distance from the quasar due to geometric dilution.B21 solve for Γ qso numerically, assigning average parameters for their QSO stack.They thereby keep the background value of the photoionization rate, Γ bg fixed while keeping the MFP as a free parameter.In reality, the photoionization rate and MFP will both depend on the ionizing emissivity and will co-evolve.B21 also fix the opacity due to higher-order Lyman series absorption to the values obtained using an optically thin simulation while fitting their model.Roth et al. in prep discuss the effect of this on the inference of B21.The slope of the variation of the MFP with redshift is also different between the counting LLS method, indirect inference method and the direct stacking method, the latter showing the steepest evolution.These differences reflect the dependence of the MFP on the nature of absorbing sources, each sensitive to a different measurement technique (Inoue et al. 2014).
In this work, we critically examine the direct measurement method of B21 for possible biases due to (a) higher cosmological densities around high-redshift quasars, (b) incomplete reionization at  > 5.3, and (c) differences between the structure of quasar proximity zones as computed using the analytical model of B21 and that obtained via radiative transfer calculations.We also investigate the challenge of defining the MFP during the epoch of reionization, when a cosmological UV radiation background is not yet established uniformly.This paper is structured as follows.We discuss our data and models in Section 2. We discuss the definition of the MFP before the end of reionization, and the effect of overdense quasar environments on the MFP in Section 3. Section 4 presents a discussion on the effect of incomplete reionization on the B21 measurements by discussing mock measurements of the MFP in our simulations using the B21 method.Finally, we present a new measurement of the MFP using the B21 data in Section 5, incorporating a radiative transfer model of quasar proximity zones.We end with a discussion and a summary of our conclusions in Section 6.Our ΛCDM cosmological model has Ω b = 0.0482, Ω m = 0.308, Ω Λ = 0.692, ℎ = 0.678,  s = 0.961,  8 = 0.829, and  He = 0.24 (Planck Collaboration XVI 2014).All distances are in proper (physical) units unless specified otherwise.

DATA AND MODELS
Our strategy in this paper is to use a cosmological radiative transfer (RT) simulation of reionization.Such a simulation provides us with a realistic model for a partially ionized IGM at  > 5.3, on top of which we model quasar proximity zones using one-dimensional radiative transfer.We then use the resultant one-dimensional spectra to examine the MFP.The cosmological RT simulation used here is described by Kulkarni et al. (2019).It consists of a cosmological hydrodynamical simulation performed using P-GADGET-3 (which is a modified version of GADGET-2, described by Springel 2005), that is postprocessed with three-dimensional radiative transfer using the ATON code (Aubert & Teyssier 2008, 2010).The box size of the simulation is 160 cMpc/ℎ with 2048 3 gas and dark matter particles.This model is calibrated to reproduce the observed mean Ly transmission at  > 5, and agrees with several high-redshift observations such as the spatial fluctuations in the Ly effective opacity (Becker et al. 2015;Bosman et al. 2022), the Thomson scattering optical depth to the last scattering surface (Planck Collaboration et al. 2020), the large-scale radial distribution of galaxies around opaque Ly troughs (Keating et al. 2020a;Becker et al. 2015), quasar damping wings (Greig et  of the IGM temperature (Gaikwad et al. 2020), and the luminosity function and clustering of Ly emitters (Weinberger et al. 2018(Weinberger et al. , 2019)).The mid-point of reionization in this simulation is at  = 7, and reionization ends at  = 5.3.At  = 5.95, the range of halo masses resolved in the simulation is in between 2.32 × 10 8 M ⊙ and 4.59 × 10 12 M ⊙ .We direct the reader to Kulkarni et al. (2019) and Keating et al. (2020a) for further details.
To construct mock spectra including the proximity zone of the QSO, we post-process our simulations with a one-dimensional radiative transfer code that is described by Satyavolu et al. (2023).All quasars are assigned a magnitude of  1450 = −27.0,corresponding to the mean magnitude of the quasars in the B21 stack.In our fiducial model, the quasars have a 'lightbulb' lightcurve with a lifetime of 10 6 yr and are placed in halos with masses between 10 11 M ⊙ ≲  halo ≲ 10 12 M ⊙ .We assume a broken power-law SED, given by Lusso et al. (2015), such that the spectral slopes above and below 912Å are −0.61 and −1.7 respectively.We then solve the thermochemistry equations for the hydrogen and helium ionized fractions along a sightline for a given quasar magnitude and lifetime at each redshift.We also simultaneously solve for the gas temperature, which is regulated by photoionization heating, adiabatic cooling, radiative cooling due to recombinations, collisional excitations, Compton cooling due to CMB photons, and Bremsstrahlung.
For our measurement of the MFP, we use the observed data from B21.This is in the form of a composite of quasar spectra with mean redshift  qso = 5.97, constructed with 13 quasar spectra obtained using the Keck/ESI and VLT/XSHOOTER spectrographs.The redshifts used for the quasars in the stack were obtained using the [C II] line, from the Ly halo emission, or from the apparent start of the Ly  1).This curve is also shown in Figure 1.The solid pink curve shows the the value of the MFP computed as the mean of the higher-value subpopulation of the bimodal distribution of free-path lengths in the simulation box.The dashed pink curve shows the MFP when computed as the mean of the full bimodal distribution of the free-path lengths.
forest absorption.While selecting quasars in the stack, a minimum cut-off of / ≳ 20 per 30 km s −1 was applied for the signal-tonoise ratio near the rest-frame 1285Å in the continuum.The mean brightness of the quasars in the sample is  1450 = −27, with a range of −27.8 <  1450 < −25.7.The typical resolution for the quasars in the stack is FWHM ∼ 45 km s −1 for spectra obtained using ESI and FWHM ∼25 km s −1 for spectra obtained using XSHOOTER.While stacking, the quasar flux is normalised by the median flux between 1270-1380 Å, although the normalisation with respect to the intrinsic continuum is kept a free parameter for our analysis, similar to B21.In order to study the Lyman-continuum MFP, we focus on the wavelength range of 912-825Å.

MFP OF HYDROGEN-IONIZING PHOTONS BEFORE THE END OF REIONIZATION
The post-reionization ( ≲ 5.5) MFP measurements are in good agreement with several reionization simulations (Keating et al. 2020b;D'Aloisio et al. 2020;Cain et al. 2021;Lewis et al. 2022).At higher redshifts  ≳ 5.5, when reionization is still believed to be ongoing, both observations and simulations differ in the measured MFP values.While some simulations agree better with the MFP measured by B21 (Lewis et al. 2022;Garaldi et al. 2022), some models (e.g., Cain et al. 2021), including ours (Kulkarni et al. 2019), agree better with the MFP measured by Gaikwad et al. (2023).Keeping this in mind, we look at possible biases in the computation of the MFP from the simulations.We use the standard definition for the MFP (Kulkarni et al. 2016) where the average Lyman-continuum transmission across random sightlines in the comoving frame is fit by an exponential with an -folding length scale of  MFP , ⟨⟩ =  0 exp −   MFP . (1) Figure 1 shows the evolution in our model of the MFP defined in this manner.The orange curve shows the MFP obtained by fitting Equation (1) to composite spectra of 13 randomly drawn sightlines from our simulation, to match the number of sightlines used in B21.
In order to obtain the sample variance, we repeat this computation for 10,000 randomly drawn samples.This method for measuring the MFP is analogous to the approach of Prochaska et al. (2009) 1 .The solid curves in Figure 1 show the median of the resultant distribution of MFP values, with the shaded regions showing the one-and two-sigma scatter.We see that the disagreement between the MFP measured in this manner from the measured value of B21 and Zhu et al. ( 2023b) is more than 2 at  = 6.Also missing in the simulations is the steep decline in the MFP observed by B21 and Zhu et al. (2023b) between  ∼ 5.5 and  ∼ 6.We also see that the evolution of the MFP in our model becomes even more gradual once reionization is complete at  ∼ 5.3.At these post-reionization redshifts, the optical depth is controlled by the overdense regions responsible for LLSs.
Direct measurements infer MFP from sightlines towards QSOs, which tend to reside in overdensities (e.g.Wang et al. 2023).To check for the bias due to large scale structure, we draw random sightlines originating only from halos with the highest masses sampled by our simulation volume, between 10 12 and 10 13 M ⊙ , otherwise applying Equation (1) as above.The resulting median MFP with the one-and two-sigma scatter is shown in Figure 1 in purple.The MFP along biased sightlines shows qualitatively similar behavior with redshift when compared to the MFP along random sightlines, but with a slightly shallower slope.The difference between the two MFPs is insignificant at  = 6.The random sightlines include sightlines that start from any location, including halos.Any difference seen between the two cases in Figure 1 is because sightlines originating from massive halos are ionized earlier than other regions in the IGM.Once reionization is complete, the over-density along sightlines originating from halos will also lead to a higher optical depth or lower MFP when compared to the measurement of random sightlines.We see in Figure 1 that this biased MFP is in better agreement with the measurement by Zhu et al. (2023b) and B21 than that by Worseck et al. (2014) at  ≲ 5.5.We discuss this further in Section 4 below.
Another potential source of bias in the MFP as formally defined by Equation ( 1) is that this exponential attenuation assumes a spatially constant opacity.In reality, reionization is not yet complete at  = 6.Consequently, the Lyman-continuum opacity of the IGM has large spatial variations.An alternative definition of the MFP that addresses this complexity is one that obtains a distribution of the free-path lengths of photons in a simulation box and computes its mean (Rahmati & Schaye 2018).We implement this in our simulation by following Rahmati & Schaye (2018) and measuring the free path as the distance at which the Lyman-continuum optical depth along the sightline becomes unity.The distribution of such free-path lengths is bimodal (see Appendix A).One part of the distribution describes the mean free path in neutral regions that are on average not overdense.The sightlines with larger free paths on the other hand also include those which encounter neutral islands along the path.The fraction of sightlines with small free paths becomes smaller than those with larger free paths at lower redshifts, until the free path distribution becomes unimodal Gaussian.We measure the MFP as the average of the dominant part of the bimodel distribution, which represents the ionized IGM.The MFP defined in this manner is shown by the solid pink curve in Figure 2. We see that this results in a value of the MFP that is very close to the value obtained by using Equation (1), shown by the orange curve in Figure 2.
In their CoDa III simulation, Lewis et al. (2022) compute the MFP using a similar free-path length method.The free-path lengths are defined in this case by doing an exponential fit to the flux along individual sightlines, which they find to be similar to the free paths measured following Rahmati & Schaye (2018).The dashed pink curve in Figure 2 shows the result from our simulation when we average over all free paths, similar to one of the approaches used in Lewis et al. (2022) to measure the MFP.As expected from the shape of the distribution of the free-path lengths, the MFP is now biased towards lower values.We find that this bias is small at  ≲ 6.This could be a reflection of the inadequate spatial resolution of our simulation, due to which a large number of free paths become smaller than our cell size.

EFFECT OF INCOMPLETE REIONIZATION ON THE B21 MEASUREMENT
We have seen above that neither the large-scale-structure bias nor the variations in the definition of the MFP cause a significant change in the value of the MFP in our simulation at  = 6.We now investigate the effect of the residual neutral hydrogen 'islands' in the IGM at this redshift on the B21 measurement of the MFP.Given the good agreement of our simulations with the measurements of the MFP by The black curve shows the average Γ qso for QSOs with lifetimes ranging between 10 5 <  q < 10 7 yr with magnitude  1450 = −27 at redshift  = 5.95, computed using our 1D radiative transfer simulation.The green, orange and blue curves show the average Γ qso computed analytically for the same QSO magnitude at redshift 6, for  = 0.33, 0.67, and 1.0 respectively.The value of  eq is computed using Γ bg = 3 × 10 −13 s −1 .Worseck et al. (2014), Gaikwad et al. (2023) and Bosman (2021), it is important to consider if the absence of patchy reionization in the models of B21 and Zhu et al. (2023a) could potentially bias their measurements.In order to do this, we construct mock data from our simulation and apply the B21 method to it.We then compare the resultant measurement of the MFP with the 'true' value of the MFP in the simulation given by Equation 1.A similar test was also performed by B21 themselves.However, the mock spectra used in their test were created from numerical simulations where reionization is assumed as instantaneous from a uniform ionizing background and post-processed to include fluctuations following the approach of Davies & Furlanetto (2016).

The B21 method
We begin by briefly reviewing the B21 method for measuring the MFP.The effective Lyman-continuum optical depth for a photon that is emitted at  qso and redshifts to the Lyman limit at  912 is given by where  912 is the opacity to 912 Å photons.The dependence on redshift is as follows: the exponent 2.75 is a result of the dependence of cross-section on frequency as  LL  ∝  −2.75 , while the exponent −5.25 comes from the conversion of comoving distance to redshift as  ∝  ∝ (1 + ) −1 × (1 + ) −3/2 in the matter-dominated era.The opacity  912 is assumed to scale with the photoionization rate as a power law, so that at a distance  from the quasar, where  bg 912 is the opacity to the ionized background and Γ tot () = Γ qso () + Γ bg () is the sum of the photoionization rate due to the QSO and the local background photoionization rate.The opacity and photoionization rate are related to each other through their mutual dependence on the shape and number density of neutral gas absorbers.This information is parameterised using the power-law index , which has been studied using analytic arguments as well as numerical simulations (?McQuinn et al. 2011).We discuss this parameter in detail in Section 4. Equation 4 is applied to a stack of QSO spectra, so the terms in the equation are all averaged quantities, and the average of the local background photoionization is assumed to be uniform and equal to Γ bg .The average photoionization rate due to all quasars at a location  is computed by iteratively solving for with the initial condition . (5) Here,  eq is the distance at which the photoionization rate due to the quasar's radiation is equal to the background photoionization rate in the absence of any absorption. eq depends on the QSO magnitude and spectrum in addition to the background photoionization rate.
The free parameters in this model are therefore,  bg 912 ,  eq ,  and Γ bg .The mean flux of the stack is fit by the above model of opacity and the MFP  MFP is inferred as the distance from the quasar at which  LL eff ( 912 ) becomes equal to 1.

Mock data
To construct stacks of mock spectra at redshift 6, we consider QSOs with magnitude of  1450 = −27 with a broken power law spectral profile for the specific luminosity as in Lusso et al. (2015), Although it is not a deciding factor, we choose the magnitude of our mock stack to be of the value same as the mean magnitude of quasars in the B21 stack at redshift 6.The power law index values are slightly different from those assumed in B21 (−0.5 and −1.5 for wavelengths greater than and less than 912 Å respectively).After placing QSOs in halos, we post-process the sightlines with our 1D radiative transfer code as described in Section 2 to obtain the distribution of the neutral hydrogen density.Using this, we compute the ionizing optical depth along a sightline as, where  LL HI = 6.3×10 −22 m 2 is the hydrogen ionization cross-section.Since Equation 4 does not explicitly assume the QSO lifetime or the host halo mass, we randomize over them in our stack, with halo masses between 10 10 M ⊙ ≲  halo ≲ 10 12 M ⊙ and lifetimes between 10 5 <  q < 10 7 Myr.Each of our composite stacks comprises 13 sightlines.
Figure 3 shows composite stacks obtained in this manner with and without quasars for the above range of quasar lifetimes and host halo masses.The top panel shows the variation of the mean continuum flux along quasar sightlines originating from ∼ 10 12 M ⊙ halos and having lifetimes between 10 5 <  q < 10 7 Myr.The mean flux computed along sightlines drawn at random directions is shown In each case, the MFP is determined using the B21 method, where  eq and Γ bg are chosen from our simulation and  = 0.67.As a comparison, the hatched region shows the 68% scatter in the MFP in the simulation, obtained using Equation (1). in grey.Once the quasar is on, the opacity near the quasar decreases as the quasar emits ionizing photons that reduce the neutral hydrogen density.The mean flux thus increases with increase in the quasar lifetime.While the decrease in flux is gradual when the quasar turns on in a uniformly ionized medium, we notice that the flux fall off is steeper in our patchy model.The bottom panel shows the variation of the mean continuum flux when the quasar lifetime is fixed to be  q = 10 6 yr and the quasar host halo masses are varied to be in three ranges between 10 10 M ⊙ ≲  halo ≲ 10 12 M ⊙ .In our simulation, more massive halos are ionized earlier and therefore the mean flux is slightly higher along these sightlines.Figure 3 shows that the mean flux changes almost doubles in the presence of the quasar, while changing the lifetime or host halo masses of the quasars by an order of magnitude changes the mean flux by less than around 10%.For the purpose of recovery of the MFP using the B21 method, we do not include the absorption due to higher order Lyman series transitions in the mock spectra.We correspondingly fit the mock spectra without including the contribution from the Lyman series opacity, as shown in Equation 8. We also do not add a zero-point correction to our mock stack and hence do not include this parameter while fitting using Equation 8. .Open circles show the MFP inferred using the B21 method from a mock composite spectrum drawn from our simulation.The curves shows the MFP in the simulation, obtained using Equation (1), using sightlines starting from randomly chosen halos (orange curve), and sightlines starting from only the most massive halos (purple curve).When we use the B21 method for measuring the MFP, we use  = 0.67 and use values of  eq and Γ bg derived from the ionized regions in the simulation.

Choice of 𝑅
with  bg 912 as free parameter.The other parameters are  eq , Γ bg and , which we discuss below.We chose to keep these parameters fixed, as B21 did for their nominal measurement.We have inspected the results while keeping  as a free parameter, and found them to be consistent with what was reported in Zhu et al. (2023b).We leave self-consistent parametrization of Γ bg with  bg 912 for future work.The distance  eq is given by Calverley et al. ( 2011) The analytic expression for  eq is computed under the approximation that the quasar photoionization rate falls as 1/ 2 , where  is the distance from the source.This equation can be used to further compute Γ qso in the presence of absorption by iterating over Equation 4. For the mock stacks, we chose a value of  eq = 13 pMpc, for  912 corresponding to magnitude of  1450 = −27.0 and Γ bg of 3 × 10 −13 s −1 as discussed below.B21 chose a value of 3 × 10 −13 s −1 for Γ bg in their analysis.The background photoionization rate in our simulations is inhomogeneous.In ionized regions, we find the average Γ bg to be around 3 × 10 −13 s −1 at redshift  ∼ 6.However, the average value of Γ bg is 10 −13 s −1 , almost 3 times smaller than the optically thin value assumed by B21.Given that at redshift  ∼ 6, the volume averaged ionized hydrogen fraction is close to 80%, it would suffice to use this value of background photoionization rate to compute Γ qso .We note that Zhu et al. (2023a) use a different Γ bg of 1.5 × 10 −13 s −1 at  ∼ 6 from Gaikwad et al. (2023).They find that at  ∼ 5.93, assuming a value of 3 × 10 −13 s −1 for Γ bg increases their mean free path from the measured value of 0.81 pMpc to ∼1 pMpc.The parameter  encodes the nature of the density of the absorbers that set the local MFP (Furlanetto & Oh 2005).For an isothermal absorber, the theoretical prediction is  = 2/3, while the scaling inferred from simulations is dependent on the self-shielding systems and can range between 0.33 and 1.0 (McQuinn et al. 2011;D'Aloisio et al. 2020).For , we use the nominal parameter from B21 where it is assumed to be 0.67 based on the arguments presented in Miralda-Escudé et al. (2000).In Figure 4, we compute the photoionization rate using Equation 4 for a quasar at redshift 6 and magnitude  1450 = −27. eq was computed using Equation 9, assuming Γ bg of ∼ 3 × 10 −13 s −1 .The resulting Γ qso is compared to the Γ qso from our simulations averaged over 1000 lines of sight to the QSO.The analytic Γ qso is in agreement with our simulated value within 25%, for a  value of 0.67.The analytic method however, under-estimates Γ qso by almost two orders of magnitude close to the quasar for all values of .Keeping  as a free parameter can change the MFP measurement by up to a factor of 2, as has been discussed in Zhu et al. (2023a).Roth et al. in prep discuss that keeping  as a free parameter yields a better fit to their mock spectra in the presence of neutral islands.They however find that such a fit might result in a MFP that is biased high with respect to the true value in their simulation.

Recovery of the 'true' MFP using the B21 method
We use the B21 method to measure the MFP in our simulation using the mock data generated using the methods discussed in Section 4. We find that the resulting MFP agrees with the MFP of our simulation measured in Kulkarni et al. (2019).This shows that the MFP measured using the B21 method is not biased due to the residual neutral hydrogen islands, quasar lifetimes or due to the overdensities around quasars.
In Figure 5, we show the distribution of MFP measured in our simulations along 10000 stacks of 13 QSOs each.In the fiducial model, the MFP is measured along random sightlines without a QSO by simply fitting an exponential to the average flux as discussed in Section 3. To compute MFP from mock stacks including a QSO, we use the B21 method with parameters discussed in the previous subsection.We consider mock stacks with different QSO properties, varying QSO lifetimes and varying host halo masses.The resulting distribution of MFP for 10000 stacks each of 13 QSOs for different QSO lifetimes and host halo masses in shown in the top and bottom panels of Figure 5, respectively.We find that irrespective of the quasar lifetime, we are able to recover the MFP from the mock stacks using the B21 method reasonably well, up to within ∼ 10% of our fiducial value.Similarly, the MFP computed from mock stacks with lower mass halos is in good agreement with the fiducial MFP.If we use a quasar stack constructed along sightlines originating in the heaviest mass halos with 10 11 M ⊙ ≲  halo ≲ 10 12 M ⊙ , the measured MFP distribution is offset from the fiducial value by ∼ 25%.
We show the MFP measured from the 10,000 mock stacks for nominal values of quasar lifetime (10 6 yr) and host halo mass (∼ 10 12 M ⊙ ) at redshifts  = 5.1 and 6 in Figure 6.Also shown are the MFPs along random sightlines and biased sightlines originating in halos but without a QSO, as discussed in Section 3. The MFP recovered from our mock stacks using the B21 method matches better with the MFP measured along biased sightlines, which is particularly evident at post-reionization redshifts.This would suggest that the analytic computation of Γ qso is robust in accounting for the QSO ionization flux, but is biased to measure the MFP along sightlines originating in overdensities.Note however that at higher redshifts, this bias becomes negligible.Another concern with the analytic model is the uncertainties introduced by the assumptions about the free parameters in the model,  eq ,  and Γ bg .Zhu et al. (2023a) have shown that increasing or decreasing the value of  can half or double the MFP estimates relative to the estimates with nominal values.This raises a need for a better constrained value of the MFP independent of the choice of .For this purpose, we use simulated models as discussed in the following section.

A DIRECT MEASUREMENT OF THE MFP USING RADIATIVE TRANSFER MODELS
We now proceed to analyse the composite spectrum of B21 using our radiative transfer models of patchy reionization and quasar proximity zones.Figure 7 shows the composite spectra from 13 randomly drawn sightlines along a QSO with magnitude  1450 = −27.0 and lifetime of 1 Myr in our simulation for a range of values of the mean free path.As our reionization simulation is calibrated to the observed Ly mean transmission, there is no freedom to change the mean free path at a given redshift.Instead, to construct Figure 7, we use the values of the ionized hydrogen fraction, the gas density, and temperature from different redshifts, in each case scaling the densities by (1 + ) 3 to  = 6.In order to fit the model stacks in Figure 7 to the data, two modifications need to be made.First, we need to account for the absorption due to higher-order Lyman-series transitions that will cause absorption bluewards of the rest-frame quasar Lyman-limit.Second, to match the observed flux, we need to consider the absorption due to the Lyman-series and continuum photons of the intrinsic QSO spectrum.Following B21 we assume the intrinsic quasar spectrum is a power law with  SED  =  912 (/912 Å) −  ion , and keep  912 as a free parameter.
We include the foreground absorption due to higher-order Lymanseries transitions as follows.At a given  obs , the absorption due to transition of the jth Lyman series line will happen due to an absorber at a redshift   =  obs /  − 1.The effective Lyman series optical depth due to all transitions is then, To obtain the optical depth due to the jth transition, we rescale the corresponding Ly optical depths at   by the product of their oscillator strength and rest-frame wavelength such that   =       /     .This scaling is valid in the optically thin regime corresponding to most of the redshifts of our interest, while in the damping wings the optical depths scale roughly as square of oscillator strengths (Malloy & Lidz 2015).We use the post-processed sightlines to compute the Ly optical depths in the presence of a quasar.The optical depths outside the proximity zone were assumed to be from the 40-2048 Sherwood simulation (Bolton et al. 2017), rescaled to the match the observed mean flux at low  as discussed in B21.From the mean flux at   , computed as negative exponential of   , the corresponding effective optical depth is obtained as   eff = −log⟨⟩.Including the absorption from Lyman-limit and 38 higher order Lyman-series terms, the observed flux is where  0 is the zero-point correction factor, that accounts for the uncertainty in the estimated zero-point of the data in B21.We fit our models to the data by sampling the posterior distribution in a Bayesian manner.We use the likelihood where the column matrices m and d denote the model and data vectors, respectively, both with  elements corresponding to the number of wavelength bins used.The covariance matrix C is obtained as follows.The data covariance matrix C data over the wavelength range of our interest (912-850 Å) is computed from the data by B21, using bootstrap realisations of the mean flux.Due to the limited data set of 13 quasars, this matrix is noisy and may underestimate the sample variance.We therefore compute a separate bootstrapped model covariance matrix for each of the model parameters using 30,000 simulated stacks and obtain a model correlation matrix .We then regularise the data covariance matrix C data using the model correlation matrix  to obtain the covariance matrix C used in Equation ( 12) at each point in the parameter space as where See applications by ??? for a discussion of this approach.We assumed uniform priors in the ranges 0.3 <  MFP < 2, 0 <  912 < 1 and −0.001 <  0 < 0.01.The model stacks and the likelihood function were linearly interpolated over a grid of models wherever necessary.We used the emcee package (?) to perform the inference by sampling the posterior distribution using MCMC via Bayes' theorem, (|d) ∝ L (d|)p(). (15) The resultant best-fit curve, given by the mean of the posterior, is shown in Figure 8, along with the 1 uncertainty.The marginalised posterior distributions of our model parameters is shown in Figure 9.We expect our assumptions for quasar lifetimes and host halo masses to introduce additional uncertainty on these parameters.We approximate these uncertainties by adding them in quadrature to the uncertainty derived from the posterior distributions.We saw in Figure 5 that varying quasar lifetimes and host halo masses leads to a ∼10% and ∼25% change, respectively, in the mean free path.With this additional error included, our inferred value of MFP is  MFP = 1.49+0.47 −0.52 pMpc.Figure 10 shows the comparison of our new measurement with previous measurements and our simulation.While our new measurement is slightly higher than the MFP measured by B21, it is consistent with their result.Our measurement is also consistent with the other estimates by Gaikwad et al. (2023) and Zhu et al. (2023a).Our new measurement remains lower than the MFP of our simulations by more than 2.Nonetheless, the consistency between our inference and that  of B21 suggests that the B21 method is robust with respective to the quasar proximity zone modelling.

CONCLUSIONS
In this paper, we have taken a closer look at the mean free path of hydrogen-ionizing photons at  = 6 in order to examine potential sources of bias in recent measurements.Our findings are as follows: • At least for  ≲ 6.5, there is no significant difference in the value of the mean free path obtained by using any of the multiple definitions of this quantity that have been used in the recent literature for redshifts at which reionization is still incomplete.
• The bias in the mean free path due to the overdensities around quasars is also minimal at  ∼ 6.At lower redshifts, with  ≲ 5, the overdensities can bias the mean free path towards lower values by about 50%.
• Due to the short mean free path at  ∼ 6, the Lyman-continuum composite spectra used in the literature for a direct MFP measurement are affected by the quasar proximity zones.The dependence of the inferred mean free path on the variations in quasar lifetime and host halo masses appears, however, only to be small and less than 25%.
• The patchiness of reionization also has minimal effect on the direct measurement of the MFP reported by B21.
• Using radiative transfer modelling of patchy reionization and quasar proximity zones, we reanalyze the data obtained by B21 to find a MFP of  MFP = 1.49+0.47 −0.52 pMpc at  = 6.This is nearly two times larger than the B21 measurement and smaller than the MFP in our reionization simulation by less than a factor of ∼ 2.
We fit to the data while self-consistently modelling quasar proximity effect on the MFP in our simulations.We include the quasars through post-processing and hence do not consider the response of gas densities to the photoheating caused by the quasar ionization.However, the dynamical timescale for the gas to respond (∼ 100 Myr) is usually much larger than the average episodic lifetime of the quasar, which is around ∼ 1 Myr (D'Aloisio et al. 2020;Morey et al. 2021;Satyavolu et al. 2023).We also do not include feedback from AGN which could result in gas heating as well as AGN driven clumpy outflows, although the former effect is shown to be not major (Trebitsch et al. 2018).The response of gas densities to heating from reionization might nonetheless be not captured in our main cosmological RT simulation, which was also run in post-processing.This would require running either hybrid or fully coupled cosmological hydrodynamical simulations (e.g.Puchwein et al. 2023;Garaldi et al. 2022).We did not check the dependence of our result on the quasar spectral profile due to computational costs.The error on MFP due to this however is expected to be small (B21).Cain et al. (2021) have modeled sinks in a subgrid fashion and found that the small-scale clumping reduces the MFP.There has also been evidence for presence of excess LLSs along the line-of-sight to quasars (Hennawi et al. 2006).We furthermore do not resolve the dense gas within mini-halos in our simulation.Recent works show that the Ly flux is suppressed by 10% on average due to mini-halos, with the suppression being enhanced along lines of sight in the vicinity of large halos, above redshifts  ≳ 5.5 (Park et al. 2023).Inclusion of dense absorbing gas such as LLSs or mini-halos in our simulations therefore will not only reduce the MFP along the random sightlines shown in orange in Figure 1, but also further reduce the continuum flux along quasar sightlines.This would push our measured MFP in Figure 10 to higher than its current value.
In summary, the shortness of the mean free path relative to radiative transfer simulation could be due to unresolved excess Lyman-limit systems (LLSs) along quasar sightlines.Consequently, the overall picture that emerges is consistent with one in which reionization ends much later than redshift 6, with the photon budget dominated by faint galaxies with high ionizing efficiency (Davies et al. 2021;Cain et al. 2021). .Distribution of free paths across redshifts.At redshifts  < 6, the free path length distribution becomes bimodal with a large number of free paths < 0.01 pMpc, smaller than the resolution of our base simulation at these redshifts.

APPENDIX A: DISTRIBUTION OF FREE PATHS
Figure A1 shows the distribution of free paths in our simulation between redshifts 8 and 5.The divergence seen at redshifts  < 6 are numerical artefacts arising from the short free paths limited by the resolution of our simulation.

APPENDIX B: EFFECT OF OFF-DIAGONAL TERMS IN THE COVARIANCE MATRIX
To understand the effect of the off-diagonal terms in the covariance matrix, we repeat our analysis of Section 5 by fitting our model to the stacked data of B21 using only the diagonal elements of the data covariance matrix provided by B21.The full, unregularised, data covariance matrix is nearly diagonal, but has some non-zero off-diagonal elements, particularly around ∼ 912Å.Setting all off-diagonal terms to zero results in an inference shown by the purple diamond in Figure B1.The magenta diamond in this figure, on the other hand, shows the result of the analysis from Section 5. We see that the off-diagonal elements lead to a larger mean value of the posterior distribution and a narrower posterior distribution, resulting in lower 1 uncertainty.The purple diamond shows the same for the case in which only the diagonal elements of the covariance matrix are used while computing the likelihood.
In this case the off-diagonal terms are set to zero.The black points show the B21 measurements.We see that the off-diagonal elements lead to a larger mean value of the posterior distribution and a narrower posterior distribution, resulting in lower 1 uncertainty.The three MFP measurements shown at  ∼ 6 are at  = 6 but we have shown them with small displacements in the redshift directions for better legibility.

Figure 2 .
Figure 2. A comparison of three definitions of the MFP in an incompletely reionized IGM.The solid orange curve shows the MFP obtained from our simulation using Equation (1).This curve is also shown in Figure1.The solid pink curve shows the the value of the MFP computed as the mean of the higher-value subpopulation of the bimodal distribution of free-path lengths in the simulation box.The dashed pink curve shows the MFP when computed as the mean of the full bimodal distribution of the free-path lengths.

1Figure 3 .
Figure 3.Effect of quasar lifetime (top panel) and host halo mass (bottom panel) on the Lyman-continuum composite spectrum of 1000 sightlines in our simulation.In both cases, the grey curve shows the transmission without the effect of the quasar proximity zone.The volume averaged neutral hydrogen fraction in our simulation at  = 5.95 is ⟨  HI ⟩ = 0.13.

Figure 4 .
Figure 4. HI photoionization rate due to QSO: The black curve shows the average Γ qso for QSOs with lifetimes ranging between 10 5 <  q < 10 7 yr with magnitude  1450 = −27 at redshift  = 5.95, computed using our 1D radiative transfer simulation.The green, orange and blue curves show the average Γ qso computed analytically for the same QSO magnitude at redshift 6, for  = 0.33, 0.67, and 1.0 respectively.The value of  eq is computed using Γ bg = 3 × 10 −13 s −1 .

Figure 5 .
Figure5.Distribution of the MFP using 10000 draws of 13-sightline stacks from our simulation for different quasar ages (top panel), and different quasar host halo mass (bottom panel).In each case, the MFP is determined using the B21 method, where  eq and Γ bg are chosen from our simulation and  = 0.67.As a comparison, the hatched region shows the 68% scatter in the MFP in the simulation, obtained using Equation (1).
Figure 6.Open circles show the MFP inferred using the B21 method from a mock composite spectrum drawn from our simulation.The curves shows the MFP in the simulation, obtained using Equation (1), using sightlines starting from randomly chosen halos (orange curve), and sightlines starting from only the most massive halos (purple curve).When we use the B21 method for measuring the MFP, we use  = 0.67 and use values of  eq and Γ bg derived from the ionized regions in the simulation.

Figure 7 .
Figure 7. Stacked quasar spectra bluewards of rest-frame 912Å in our model for a range of values of the MFP at  = 6.In each case, the solid curve shows the composite of 13 randomly chosen sightlines from the simulation box.The shaded region shows the 68% scatter across 10,000 such samples.

Figure 8 .
Figure 8.Our best-fit model for the B21 composite spectrum, with the 1 uncertainty indicated by the shaded region.

Figure 9 .
Figure 9. Corner showing the posterior distributions of the normalization parameters (  912 ,  0 ) and the mean free path ( MFP ).The rightmost panels in each row indicate the marginalized posteriors.The black, green, and yellow dotted lines represent the median, 1, and 2 levels of the distribution, respectively.
Figure A1.Distribution of free paths across redshifts.At redshifts  < 6, the free path length distribution becomes bimodal with a large number of free paths < 0.01 pMpc, smaller than the resolution of our base simulation at these redshifts.

Figure B1 .
Figure B1.Effect of the off-diagonal terms in the covariance matrix C on our inferred value of the mean free path.Magenta diamond shows the mean of the posterior distribution computed using the full, regularised covariance mastrix as described in the text.Also shown is the one-sigma uncertainty.The purple diamond shows the same for the case in which only the diagonal elements of the covariance matrix are used while computing the likelihood.In this case the off-diagonal terms are set to zero.The black points show the B21 measurements.We see that the off-diagonal elements lead to a larger mean value of the posterior distribution and a narrower posterior distribution, resulting in lower 1 uncertainty.The three MFP measurements shown at  ∼ 6 are at  = 6 but we have shown them with small displacements in the redshift directions for better legibility.
Bosman (2021).(2023), andZhu et al. (2023b)e MFP using the B21 data with our model.The error bar on this point shows the 1 uncertainty.We see that this measurement is consistent with that of B21 which is shown by the black point.This suggests that the B21 measurement is robust against the sources of potential bias that we have considered in this paper.The orange curve with shaded regions shows the mean free path in the reionization simulations ofKulkarni et al. (2019).We see that our new measurement prefers a factor-of-two smaller MFP than this simulation.We also show other measurements of the MFP, byWorseck et al. (2014),Bosman (2021),Gaikwad et al. (2023), andZhu et al. (2023b).Measurements byBecker et al. (2021);Gaikwad et al. (2023)andBosman (2021)have been displaced in redshift by  = ± 0.03 and 0.05 respectively for legibility.