Constraints on the Evolution of the Ionizing Background and Ionizing Photon Mean Free Path at the End of Reionization

The variations in Lyα forest opacity observed at z > 5.3 between lines of sight to different background quasars are too strong to be caused by fluctuations in the density field alone. The leading hypothesis for the cause of this excess variance is a late, ongoing reionization process at redshifts below six. Another model proposes strong ionizing background fluctuations coupled to a short, spatially varying mean free path of ionizing photons, without explicitly invoking incomplete reionization. With recent observations suggesting a short mean free path at z ∼ 6, and a dramatic improvement in z > 5 Lyα forest data quality, we revisit this latter possibility. Here, we apply the likelihood-free inference technique of approximate Bayesian computation (ABC) to jointly constrain the hydrogen photoionization rate ΓHI and the mean free path of ionizing photons λ mfp from the effective optical depth distributions at z = 5.0–6.1 from XQR-30. We find that the observations are well-described by fluctuating mean free path models with average mean free paths that are consistent with the steep trend implied by independent measurements at z ∼ 5–6, with a concomitant rapid evolution of the photoionization rate.

The epoch of reionization reflects the cumulative photon output of the first generations of stars and galaxies in the Universe.Determining the precise timing of reionization thus provides crucial clues for understanding structure formation in the first billion years of cosmic time.Observations of the cosmic microwave background suggest a characteristic epoch of z ∼ 7 (Planck Collaboration et al. 2020), and this rough midpoint is supported by constraints derived from the disappearance of Lyα emission from galaxies at z > 6 (e.g.Mason et al. 2018Mason et al. , 2019;;Weinberger et al. 2019;Hoag et al. 2019;Jung et al. 2020, although see Wold et al. 2022) and studies of the Lyα damping wing in the highest redshift quasars known at z ≳ 7 (e.g.Davies et al. 2018c;Wang et al. 2020;Yang et al. 2020a;Greig et al. 2022).
The endpoint of reionization was originally estimated by the disappearance of widespread transmission in the Lyα forest (Becker et al. 2001;Fan et al. 2006), thought to be due to the onset of Gunn-Peterson absorption (Gunn & Peterson 1965) of neutral hydrogen in the intergalactic medium (IGM).That interpretation is complicated, however, by the concomitant decrease in the intensity of the ionizing background (Bolton & Haehnelt 2007;Wyithe & Bolton 2011;Davies et al. 2018b;D'Aloisio et al. 2018) which should lead to a mostlyopaque Lyα forest with strong sightline-to-sightline variations even without significantly neutral gas (e.g.Lidz et al. 2006).With the discovery of the giant Gunn-Peterson (GP) trough at z ∼ 5.5-5.8towards the quasar ULAS J0148+0600 by Becker et al. (2015), however, existing models of the post-reionization IGM were no longer consistent with the distribution of Lyα forest opacity at z ≳ 5.6.This inconsistency was confirmed by subsequent compilations of Lyα forest opacity measurements (Bosman et al. 2018;Eilers et al. 2018) and has recently been carefully quantified to persist to even later times z ∼ 5.3 by Bosman et al. (2022).
Several models were put forward to explain the origin of these excess fluctuations, and in particular the existence of the giant GP trough from Becker et al. (2015).In ionization equilibrium, the Lyα forest opacity is directly connected to the residual neutral hydrogen fraction, which depends on the density of the gas, the strength of the ionizing background, and (through the recombination rate) the gas temperature.D 'Aloisio et al. (2015) suggested that relic temperature fluctuations from a late-ending, extended, and hot reionization process could lead to dramatic contrast in IGM temperatures at z ≲ 6, leading to strong variations in Lyα forest opacity due to the temperature dependence of the hydrogen recombination rate.Such models predict that opaque troughs correspond to large-scale galaxy overdensities (Davies et al. 2018a), however, which is disfavored by galaxy surveys towards the largest GP troughs at z ∼ 5.7 (Becker et al. 2018;Kashino et al. 2020;Christenson et al. 2021).Chardin et al. (2015Chardin et al. ( , 2017) ) showed that the strong fluctuations in the ionizing background would be a natural consequence of a bright and rare source population, i.e. luminous quasars (see also Meiksin 2020), but initial suggestions of an overabundant faint quasar population at z > 5 (Giallongo et al. 2015) are now generally disfavored (e.g.D 'Aloisio et al. 2017;Parsa et al. 2018;Matsuoka et al. 2018).While interest in such quasar-dominated models has begun to resurface in light of the discovery of a substantial population of faint AGNs at z > 5 by JWST surveys (e.g.Kocevski et al. 2023;Harikane et al. 2023;Matthee et al. 2023;Labbe et al. 2023;Maiolino et al. 2023), due to their typically highly reddened nature it is still unknown whether these objects contribute substantially to the ionizing photon budget.Davies & Furlanetto (2016) explored the possibility that galaxy-sourced ionizing background fluctuations could be amplified by a coupling to the mean free path of ionizing photons (see also D'Aloisio et al. 2018) following analytic arguments from McQuinn et al. (2011), but in order to match the broad distribution of Lyα forest opacity the average mean free path would have to be a factor of a few shorter than the extrapolation from lower-redshift measurements (Worseck et al. 2014).All of these models pre-supposed that reionization was more-or-less complete at z ≲ 6, and required adjusting their corresponding parameters to uncomfortable ends of the parameter space to be even qualitatively consistent with Lyα forest observations.
The most successful model thus far was proposed by Kulkarni et al. (2019, see also Nasir & D'Aloisio 2020), wherein reionization is not finished by z ∼ 6, but instead ends at z ≲ 5.5.In this model, the large variations in the Lyα forest opacity at z < 6 are explained by residual "true" (i.e., mostly neutral) GP troughs and shadowing of the ionizing background by large-scale patches of neutral IGM.In particular, late reionization models most readily reproduce observations of rare, extremely large GP troughs (Keating et al. 2020) and so-called "dark gaps" (Zhu et al. 2021(Zhu et al. , 2022) ) even down to redshifts z ≲ 5.5.Such models, however, require careful tuning of the ionizing emissivity evolution to selfconsistently reproduce both the late Lyα forest fluctuations and other reionization-epoch observables (Kulkarni et al. 2019;Keating et al. 2020), and the efficient moment-based radiative transfer method used in these works may suffer from numerical suppression of largescale ionizing background fluctuations inside the ionized regions (Wu et al. 2021).Nevertheless, this explanation has proven to be the most natural one so far, and has begun to deliver powerful constraints on the properties of the reionizing sources (Choudhury et al. 2021;Qin et al. 2021).
An important boundary condition to reionization is the strength of the ionizing background which permeates the Universe after the process is complete, and radiative transfer simulations generally predict a rapid rise in the amplitude of the ionizing background as the ionized bubbles overlap and the last neutral islands disappear (e.g.Gnedin 2000).Constraints on the hydrogen photoionization rate Γ HI derived from Lyα forest observations suggest a rise by a factor of at least a few from z ∼ 6 to z ∼ 5 (Bolton & Haehnelt 2007;Wyithe & Bolton 2011;Calverley et al. 2011;Davies et al. 2018b), followed by a much shallower evolution from z ∼ 5 to z ∼ 2 (Becker & Bolton 2013), which is potentially suggestive of this picture, although we note that the flattening could instead reflect changes in the nature of the gas responsible for the ionizing opacity (Muñoz et al. 2016).Precision measurements of Γ HI at z ∼ 5-6 are complicated by the low mean transmission of the Lyα forest (Eilers et al. 2018;Bosman et al. 2018) and a degeneracy with the relatively unknown thermal state of the IGM (e.g.Becker & Bolton 2013), but recent improvements in data quality (Bosman et al. 2022) and new constraints on the IGM thermal state at z > 5 (e.g.Walther et al. 2019;Boera et al. 2019;Gaikwad et al. 2021) suggest that we can now do much better.
Furthermore, the ionizing background at the very end of reionization provides a census of the ionizing photons being produced throughout the Universe which drove the process to completion, and introduces a strong boundary condition to models of the reionizing sources (e.g.Bouwens et al. 2015a).Connecting the emissivity of ionizing photons ϵ ion to the photoionization rate Γ HI , however, requires an additional ingredient: the mean free path of ionizing photons, λ mfp .While the mean free path can be estimated from the distribution of neutral hydrogen systems at redshifts z ≲ 4 (e.g.Rudie et al. 2013;Prochaska et al. 2014), at higher redshifts the identification of individual absorption lines becomes more challenging.To overcome this difficulty, and to bypass potential uncertainties related to line blending and absorber clustering (Prochaska et al. 2014), the mean free path can instead be measured directly by stacking quasar spectra at their rest-frame Lyman limit (Prochaska et al. 2009).Until recently, such measurements were limited to z ≲ 5 (Prochaska et al. 2009;O'Leary & McQuinn 2012;Fumagalli et al. 2013;Worseck et al. 2014), and so constraints on the ionizing emissivity were limited to a time well after the end of reionization (Becker & Bolton 2013) or required an extrapolation of the mean free path evolution to earlier times (D'Aloisio et al. 2018).
Recently, Becker et al. (2021) measured the mean free path of ionizing photons at z ∼ 6, improving the spectral stacking method by additionally modeling the impact of the quasar ionizing radiation on the line-of-sight Lyman-series and Lyman limit opacity.They measured a mean free path of λ mfp = 0.75 +0.65  −0.45 proper Mpc, much shorter than the extrapolated value implied by the tight power-law evolution measured by Worseck et al. (2014) at z ≲ 5.This mean free path is short enough to strongly affect the progression of reionization and place higher demands on the ionizing output from galaxies (Cain et al. 2021;Davies et al. 2021), and implicitly requires the structure of the gas to trace that of (cold) dark matter down to very small scales (Emberson et al. 2013;Park et al. 2016;D'Aloisio et al. 2020).Notably, an interpolation between the z = 5.1 and z = 6 measurements of Becker et al. (2021) now places the mean free path at the level required by the Davies & Furlanetto (2016) model to explain the Lyα forest fluctuations.
Here we revisit Lyα forest fluctuations in the context of the Davies & Furlanetto (2016) model to constrain parameters of the z > 5 IGM, taking advantage of recent improvements in Lyα forest data quality and new constraints on the IGM thermal state (Gaikwad et al. 2021).We apply a likelihood-free inference methodology based on Davies et al. (2018b) to the extended XQR-30 (D'Odorico et al. 2023) Lyα forest data set (Bosman et al. 2022) to constrain the average photoionization rate Γ HI and the average mean free path of ionizing photons λ mfp from z = 5.0-6.1.We also assess the degree to which the data are consistent with a draw from the model.Finally, we discuss what our constraints imply for the evolution of the ionizing emissivity and the timing of the reionization epoch.
In this work we assume a ΛCDM cosmology with h = 0.685, Ω m = 0.31, Ω b = 0.049, consistent with Planck (Planck Collaboration et al. 2020).Distance units are comoving unless otherwise specified.

SIMULATION METHOD
We simulate post-reionization Lyα forest fluctuations in a manner very similar to Davies et al. (2018b).The general philosophy is to combine skewers through a cosmological hydrodynamical simulation, which can resolve scales small enough for a converged description of the Lyα forest, with skewers from a separate semi-numerical ionizing background simulation, which has a volume sufficiently large to sample the full distribution of largescale background fluctuations.In the process, we unfortunately decouple the density of the gas responsible for Lyα absorption from the intensity of the ionizing radi-ation field.We will revisit this assumption in Section 5 and attempt to quantify the resulting biases.

Nyx hydrodynamical simulation
For the small-scale density, temperature, and velocity structure of the gas, we use a hydrodynamical simulation run with the Nyx code (Almgren et al. 2013).The simulation we employ is 100/h Mpc on a side and has 4096 3 dark matter particles and baryon grid cells, providing a box size and resolution adequate for converged Lyα forest statistics (Lukić et al. 2015;Oñorbe et al. 2019).Our fiducial modeling used three redshift snapshots from z = 5, z = 5.5, and z = 6; for redshifts between these values we use the overdensity field from the snapshot closest in redshift, and scale the physical densities by (1 + z) 3 .We extract 100,000 skewers of gas overdensity, temperature, and peculiar velocity from each snapshot starting from random locations in the simulation box and traversing in one of six random directions along grid axes (±x, y, z).Overdensities are converted into physical densities ∝ (1 + z) 3 according to the evolving redshift across the range of interest.Lyα forest spectra are computed by assuming ionization equilibrium to calculate the corresponding neutral hydrogen density and using the Voigt profile approximation of Tepper-García (2006) to deposit Lyα opacity onto a 2 km/s resolution grid.
One fundamental limitation of our Nyx simulation is that, by virtue of its optically-thin photoionization and photoheating rates from Haardt & Madau (2012), reionization occurs very early (z ∼ 15), and with minimal heat injection (∆T ∼ 10, 000 K).This is in conflict with observational constraints on the timing of reionization (e.g.Davies et al. 2018c) and theoretical models of the expected heating by ionization fronts sweeping through the IGM (Miralda-Escudé & Rees 1994;Abel & Haehnelt 1999;D'Aloisio et al. 2019).As a result, at z = 5-6 the IGM is colder, and the temperaturedensity relation steeper, than our theoretical understanding would otherwise suggest, with further evidence coming from recent constraints on the IGM thermal state at z ∼ 5.5 (Gaikwad et al. 2020).The combination of these two effects will lead to overestimated recombination rates, and thus overestimated values of Γ HI when matched to the observed Lyα forest opacity (e.g.Bolton et al. 2005;Bolton & Haehnelt 2007;Becker & Bolton 2013;D'Aloisio et al. 2018).To overcome this limitation of our hydrodynamical simulations, we adjust the IGM temperatures of our fiducial skewers to match a more realistic average post-reionization thermal state.We model the temperature evolution of the IGM following the analytic method from Davies et al. (2018a)  (see also Upton Sanderbeck et al. 2016) assuming an instantaneous reionization at z re with a heat injection of ∆T = 20, 000 K (e.g.D' Aloisio et al. 2019).
In Figure 1 we show the resulting mean temperature T 0 and slope of the temperature-density relation γ for a range of z re compared to recent measurements of the IGM thermal state (Walther et al. 2019;Boera et al. 2019;Gaikwad et al. 2020).From this comparison, we adopt z re = 7.2 ± 0.5 as our fiducial range of thermal models.However, we will also show the range of inferred Γ HI values corresponding to more extreme thermal models: the cold, native Nyx temperatures corresponding to z re ∼ 15 and a very hot model with z re = 6.2.By manually adjusting the simulated temperatures in postprocessing, we neglect the effect of Jeans smoothing on the Lyα transmission, but this is likely a very small effect (see, e.g., Becker & Bolton 2013;Kulkarni et al. 2015).

Ionizing background fluctuations
The volume of the Nyx simulation described above is not entirely sufficient to capture the large-scale ionizing background fluctuations characteristic of the postreionization IGM, where box sizes of at least a few hundred Mpc are required (e.g.Davies & Furlanetto 2016;D'Aloisio et al. 2018).For our fiducial model, we decided to create a much larger, but independent, cosmological volume within which to simulate the ionizing background -the consequences and biases of this approach will be discussed later.
We computed large-scale simulations of ionizing background fluctuations using the method of Davies & Furlanetto (2016), built off of the semi-numerical framework of the 21cmFAST code from Mesinger et al. (2011).Cosmological initial conditions were generated on a 4096 3 grid in a volume 512 comoving Mpc on a side, then evolved to z = 5.5 onto a coarser 512 3 grid via the Zel'dovich approximation (Zel'dovich 1970).Halos were instantiated from the 4096 3 grid via the excursion set halo-filtering approach of Mesinger & Furlanetto (2007) down to a minimum halo mass M h,min = 2 × 10 9 M ⊙ , and shifted to their evolved positions at z = 5.5 via the Zel'dovich displacement field evaluated on the 512 3 grid at their initial positions from linear theory, a procedure which produces realistic halo clustering down to ∼Mpc scales (Mesinger & Furlanetto 2007).Ultraviolet (UV) luminosities of galaxies corresponding to each halo were assigned by abundance matching to the Bouwens et al. (2015b) UV luminosity functions, resulting in UV magnitudes ranging from −12.6 ≳ M UV ≳ −23.3.For simplicity, we assume that the ionizing luminosity of each galaxy is proportional to its UV luminosity, and leave the constant of proportionality as a free parameter.
We then compute the ionizing background by a brute force radiative transfer algorithm applied to a coarser, 128 3 grid, corresponding to cells 4 Mpc on a side.The ionizing luminosities of the galaxies are deposited onto the coarse grid, and then the ionizing radiation intensity in cell i is computed via where is the intensity due to local sources, d ij is the distance between cells i and j, and τ ij is the ionizing photon optical depth between i and j.The sum over all other cells j is performed assuming locally periodic boundary conditions, i.e. the cell i "sees" the other cells j within a volume of the same size as the full box but centered on cell i.The local source intensity is computed assuming a uniform ionizing emissivity ϵ ν,i inside the cell, where λ ν,i is the mean free path of cell i at frequency ν, l is the side length of the cell, and the constant factor of 0.72 was found to have good resolution convergence properties in Davies & Furlanetto (2016).The optical depth between cells is computed by integrating the opacity κ = dτ /dr on the 128 3 grid, where κ is assumed to vary with the local photoionization rate Γ and overdensity relative to the cosmic mean ∆ = ρ/ρ according to where κ 0 is the opacity of the IGM at the ionizing edge of hydrogen ν HI and at a reference photoionization rate Γ 0 , and we adopt the same power-law dependencies on ∆ and ν as Davies & Furlanetto (2016).Finally, the photoionization rate is calculated by integrating over J ν , where σ HI is the hydrogen photoionization cross-section (from Verner et al. 1996).In practice, for computational efficiency we perform the calculation with a single ionizing photon frequency hν ≈ 17.9 eV, which was found to produce very similar Γ HI fluctuations to a more comprehensive multi-frequency calculation (Davies & Furlanetto 2016).
We computed fluctuating Γ HI fields for 14 different average mean free path values 6, 8, 10, 15, 20, 25, 30, 40, 50, 60, 80, 100, 150 (comoving) Mpc at a fixed redshift z = 5.5.This relatively coarse sampling was necessitated by the computational inefficiency of the Davies & Furlanetto (2016) method; see Gaikwad et al. (2023) for a recent implementation of a more efficient algorithm.We first computed the Γ HI field assuming a uniform κ = κ 0 (ν/ν HI ) −0.9 , and then determine a normalization factor (corresponding to a combination of an ionizing emissivity normalization and a correction for the monochromatic approximation) to initialize the mean Γ HI to the value corresponding to the prescribed ionizing opacity, i.e.Γ 0 in equation ( 4).In the following, we will use only the relative fluctuations in Γ, i.e. the relative fluctuating field ΓHI ≡ Γ HI / ΓHI , so the exact value of Γ 0 is unimportant.Note that due to the dependence of κ on Γ HI , the calculation must be iterated many times to achieve self-consistency.We iterate the calculation until the average change in Γ HI falls below 0.1%, requiring ≳ 20 iterations for λ mfp ≲ 10 Mpc and ∼ 5 iterations for λ mfp ≳ 60 Mpc.
In Figure 2, we show a slice through four of the fluctuating ionizing background fields with λ mfp = 10-80 Mpc.The spatial structure of the fluctuations is very similar, driven by the large-scale clustering of the ionizing sources, but the fluctuations are stronger at shorter average mean free path.We take advantage of this coherence and linearly interpolate in log Γ between adjacent models to produce skewers at arbitrary average mean free path from 5 to 150 Mpc.
As discussed in the following Section, our inference procedure requires many millions of simulated Lyα forest spectra.To reduce the computational requirements for producing the mock spectra, we adopt an approach analogous to the fluctuating Gunn-Peterson approximation (Weinberg et al. 1997) to both adjust the mean Γ HI and introduce the fluctuating Γ HI field.We first compute an initial set of 100,000 Lyα forest spectra from the Nyx simulation skewers with a uniform Γ HI = Γ init , where we set Γ init to the midpoint of the redshiftdependent uniform prior on Γ HI used during the inference procedure (see § 4).We then introduce Γ HI fluctuations by re-scaling the pixel optical depths along the skewer, where Γ new is the new average value of Γ HI and ΓHI is evaluated at z i without taking into account redshiftspace distortions.We find that this procedure is ∼ 50 times faster than a more accurate re-evaluation of the spectra from the neutral hydrogen distribution, without substantially affecting our main results.
While the ionizing background is inherently fluctuating in our simulations, references to "Γ HI " in the rest of the text, particularly in the context of our constraints from observations, will represent the volume-averaged value of Γ HI .

Self-consistent background fluctuations
For comparison and bias assessment purposes, we also produced a similar suite of simulations of ionizing background fluctuations using the halo catalog of the 100 Mpc/h Nyx simulation with the same set of average mean free path values.We abundance-matched dark matter halos down to the same minimum mass, M min = 2 × 10 9 M ⊙ , but we note that due to the coarse force resolution of the uniform grid particle-mesh scheme used to evolve the dark matter particles, the number of halos below M h ∼ 10 10 M ⊙ is significantly underestimated (e.g.Lukić et al. 2007;Almgren et al. 2013).As a result, the UV magnitude range is somewhat restricted compared to our larger volume, −13.9 ≳ M U V ≳ −22.7.
For these complementary simulations, we computed the background fluctuations on a 64 3 grid, corresponding to a spatial resolution of ∼ 2 Mpc.These simulations tend to have weaker large-scale fluctuations than the fiducial ones for relatively short mean free paths, but for long mean free paths they exhibit stronger large-scale fluctuations.The former is due to the suppression of large-scale modes in the density field due to the smaller volume, the latter is due to a limitation of the algorithm used to compute the radiation field.Specifically, the size of the "local volume" seen by any cell sets an upper limit to the effective mean free path of roughly half the box size, or in this case ∼ 73 Mpc.
We will refer to this suite of simulations as the "selfconsistent" model, as it provides ionizing background fluctuations sourced by the same underlying density field as the Lyα forest.As mentioned above, the fiducial model relies on an independent volume for ionizing background fluctuations, which inherently decouples them from the density field.We will compare the constraints obtained in the fiducial model to the selfconsistent model, as well as a de-correlated version of the self-consistent model, where the ionizing background fluctuations are (as in the fiducial model) drawn from a random region in the 100 Mpc/h volume.

LIKELIHOOD-FREE INFERENCE WITH APPROXIMATE BAYESIAN COMPUTATION
Armed with the simulation framework described above, we aim to constrain two parameters: Γ HI and λ mfp .The distribution of Lyα forest opacity at z > 5 is not well-described by Gaussian statistics (e.g.Bosman et al. 2018) that would typically be adopted in Bayesian parameter inference, and at z ≳ 5.6 one must also sensibly account for non-detections that represent the most constraining GP troughs.We adopt a likelihood-free inference technique known as Approximate Bayesian Computation or ABC (Rubin 1984;Tavaré et al. 1997;Pritchard et al. 1999) to overcome these challenges with a principled approach1 .Our method is very similar to that described in Davies et al. (2018b), but with some important differences that we describe below.
The approximation fundamental to the "approximate" nature of ABC is in the definition of the likelihood of the data vector d given the model parameters θ, p(d|θ).In ABC, the likelihood is approximated by (e.g.Marin et al. 2012) where ρ is a (nearly arbitrary) distance metric between the data and a mock data draw x and ϵ is a distance threshold below which the data and mock data are deemed "similar enough".As ϵ approaches zero, the approximate posterior PDF will converge towards the true posterior PDF (e.g.Blum 2010).Typically the raw data are first transformed into a lower dimensional summary statistic s(d) ≡ s d , and the same procedure is applied to make mock observations of the summary statistic s(x) ≡ s x .The ABC procedure involves computing many such mock data sets, with parameter values θ drawn from the prior p(θ), and selecting some number of samples with distances below some threshold ϵthose samples then (approximately) represent samples from the posterior PDF p(θ|d).As discussed below, the threshold ϵ is typically chosen such that one retains a given (small) fraction of the parameter samples.
There is substantial freedom in choosing both the distance metric ρ and the summary statistic s.We choose to summarize the Lyα forest data by their ∆z = 0.1binned effective optical depths s = − ln ⟨exp −τ α ⟩ = τ eff , and compute the Euclidean (L 2 norm) distance between the rank-order set of observed and mock τ eff (see also Davies et al. 2018b).Specifically, we have where τ obs eff,i i and τ mock eff,i are the ith highest τ eff values in the set of N q observed and mock Lyα forest spectra, respectively.While the transformation to effective optical depth provides sensitivity to strongly absorbed regions of the Lyα forest, observational noise can lead to fully opaque GP troughs with mean transmission below zero, leading to undefined τ eff .In both the observed data and mock data, we set τ eff = − ln (2 × σ F ) if the mean transmission is below twice the statistical error σ F ; a common transformation in the z > 5 Lyα forest literature (Becker et al. 2015;Eilers et al. 2018;Bosman et al. 2018Bosman et al. , 2022)).
In Figure 3, we demonstrate an example of the procedure described above applied to a mock Lyα forest data set at z = 5.6.We assume true values of Γ HI = 5 × 10 −13 s −1 and λ mfp = 25 Mpc, and a data set consisting of 50 quasar sightlines with a statistical error σ F = 0.01 on the ∆z = 0.1-binned transmitted flux and a continuum error σ C = 0.1.The cumulative distribution function (CDF) of this "observed" data set is shown by the thick curves in the left panel of Figure 3, where the left (blue) and right (red) curves display undetected transmission with F < 2σ F as τ eff = − ln (2σ F ) or τ eff = ∞, respectively, similar to Bosman et al. (2018Bosman et al. ( , 2022)).We draw 1,000,000 values of θ = {Γ HI , λ mfp } assuming a uniform prior from Γ HI = 0-10 −12 s −1 and λ mfp = 5-150 Mpc, and for each θ we compute a mock data set including the same noise properties as the observed one.The resulting distances ρ(s d , s x ) vs. Γ HI and λ mfp are shown in the (upper) middle and right panels of Figure 3, respectively.We then choose a series of distance thresholds ϵ such that 50%/10%/0.1% of the mock data sets have ρ(s d , s x ) < ϵ, shown by the horizontal lines, and retain the corresponding θ values as samples from the posterior PDFs in the lower panels.The blue and red shaded regions curves in the left panel of Figure 3 show the central 95% of CDFs computed from mock data sets adopting the lowest ϵ posterior mean Γ HI and λ mfp , demonstrating good consistency with the "observed" one.We apply the ABC inference methodology described above to the extended Lyα forest dataset from Bosman et al. (2022).The data consist of 51 optical quasar spectra observed with the X-Shooter spectrograph (Vernet i.e. 67 quasar spectra in total.All quasars are located at z > 5.7 and their spectra possess signal-to-noise ratios (SNRs) larger than 10 per pixel.Intervening damped Lyman-α absorbers are masked based on the detection of associated metal absorbers in the metal catalog of Davies et al. (2023).In a few quasars, portions of the spectra with suspected broad absorption line contamination are also masked.We refer interested readers to (Bosman et al. 2022) for more details of the sample.
In order to measure the effective optical depths, the intrinsic quasar continua are reconstructed based on their un-absorbed emission features at λ rest > 1215.16Å.We use the principal component analysis (PCA) approach of Davies et al. (2018d) with the improvements described in Bosman et al. (2020) and Bosman et al. (2022).The PCA approach reproduces the quasar continua at λ rest < 1180 Å with a wavelength-dependent ∼ 8% accuracy, and sub-percent precision.Finally, the effective optical depths are measured in bins of dz = 0.1 and only bins which are > 50% un-masked are retained.The number of retained quasar sightlines N q at each redshift is shown in Table 1.The masking, observational uncertainties, and wavelength-dependent continuum uncertainties are forward-modeled in all simulated mocks.Bosman et al. (2022) showed that mock data generated in this manner is statistically indistinguishable from the real data in the post-reionization regime (5.0 < z < 5.3), solidifying our confidence in the noise model.
For each redshift bin -the same dz = 0.1 bins defined in Bosman et al. (2022) -we compute 10,000,000 mock data sets, drawing from uniform priors in Γ HI and λ mfp in uniform prior windows2 where the posterior PDF has substantial support.We then define a relative threshold ϵ th by the fraction of retained mocks with distances below the actual ABC threshold, e.g. with ϵ th = 0.1 we would retain the parameter values whose mock data sets resulted in a distance smaller than the 10th percentile of the entire distribution.For our fiducial posterior PDFs we retain the 1,000 samples with the smallest distances in our set of 10,000,000, i.e. ϵ th = 10 −4 .
We show the resulting 2D posterior PDFs of Γ HI and λ mfp from z = 5.0-6.1 in Figure 4.The posteriors show only modest correlations between the two parameters, and both decrease steadily to higher redshift.While Γ HI is well-constrained at all redshifts, along the λ mfp dimension the posterior clearly runs into the edge of the prior (set by our grid of models) at z ≲ 5.3 and z ≳ 5.9.To understand this apparent lack of constraining power, note that our constraints on λ mfp are driven by a difference between the observed Lyα forest fluctuations and the fluctuations one would expect given a uniform ionizing background.At the low redshift end, the weak constraints are due to the fact that the observed Lyα forest fluctuations are fully consistent with a uniform ionizing background (Bosman et al. 2022), while at the high-redshift end the strength of fluctuations is comparable to our shortest λ mfp model.We note that the width of the posterior PDFs in the λ mfp dimension are particularly large at z = 5.5 and z = 5.9 compared to their neighboring redshifts.In general, these variations come about due to the non-linear connection between λ mfp and the width of the τ eff distribution -the strength of radiation field fluctuations varies only modestly for λ mfp ≳ 40 Mpc (e.g. Figure 2), so any upward tail of the posterior will inevitably be elongated.In addition, at z = 5.5 where this apparent difference is most pronounced, it was already noted in Bosman et al. (2022) that the disagreement between the τ eff CDF and simulations assuming a uniform ionizing background was smaller than at any other redshift z > 5.3.We can thus understand the relatively weak constraint, as this disagreement between the observations and uniform background simulations is the source of the constraining power of our analysis.
In Figure 5 we show the posterior medians (black circles) and central 68% credible intervals (black error bars) of Γ HI , marginalized over λ mfp .We have also adjusted the posterior constraints by a few percent to correct for the coarse snapshot sampling employed in the simulated spectra (see Appendix A); the "raw" uncorrected posterior medians are shown as open squares.Note that these statistical error bars apply only to the fiducial IGM thermal model with z re = 7.2.The dark grey shaded regions in Figure 5 show the range of posterior medians obtained for different IGM thermal models with z re = 6.7 and z re = 7.7, while the light shaded region shows the range corresponding to extreme thermal models with z re = 6.2 and z re ∼ 15 (see Section 2).The uncertainty in Γ HI resulting from the IGM thermal state is much larger than the statistical uncertainty except for z ≳ 5.9, where the transmitted flux is much lower and sampled across relatively few sightlines.Tighter statistical constraints could likely be achieved at these redshifts from more informative summary statistics (e.g. the flux PDF on smaller scales and/or Lyβ transmission, cf.Davies et al. 2018b), but we leave a more detailed look at the z ≳ 6 IGM in XQR-30 to future work.The crosses in Figure 5 show the posterior medians in the self-consistent model ( § 2.2.1),where the ionizing background fluctuations are drawn from the same physical location in the Nyx hydrodynamical box as the density field when computing the Lyα forest absorption.We see a gradual positive offset of the Γ HI values in the self-consistent model which increases from a few percent at z = 5 to ∼ 25% at z ≥ 6.This offset comes about due to correlation between ionizing background fluctuations and the density field -regions with high Γ HI tend to have higher density, and regions with low Γ HI tend to have lower density, shifting the mean Lyα forest transmission to lower values.Thus the mean Γ HI must be higher to reproduce the mean Lyα forest opacity.If we randomize the ionizing background fluctuations with respect to the density field in the 100 Mpc/h box, we find that the resulting Γ HI constraints are indistinguishable from the fiducial model.
In the left panel of Figure 6, we compare our constraints on Γ HI (solid points) to literature values (open points) after combining the statistical uncertainties with the thermal state uncertainties in quadrature, along with an additional 0.03 dex of systematic uncertainty to approximate the uncertainty due to Jeans smoothing (e.g.Becker & Bolton 2013) and an additional upward systematic uncertainty given by the bias between the fiducial and self-consistent model constraints.We also show grey points which represent the posterior medians shifted upward by the same amount.Our constraints are consistent with previous measurements by Calverley et al. (2011), Wyithe &Bolton (2011), andD'Aloisio et al. (2018), as well as constraints from the pilot study of Davies et al. (2018b), although we note that we have not carefully accounted for the different IGM thermal states (or ranges of thermal states) assumed by those works.In the right panel of Figure 6, we compare to predictions of Γ HI (z) from empiricallymotivated 1D cosmological radiative transfer models by Haardt & Madau (2012), Khaire & Srianand (2019), and Faucher-Giguère (2020), and the 3D coupled radiationhydrodynamics simulations from Keating et al. (2020), THESAN (Garaldi et al. 2022), and CoDa-III (Lewis et al. 2022).While the 1D radiative transfer models are all in rough agreement with the average Γ HI from z ∼ 5-5.5, they fail to reproduce the steep downturn towards z ∼ 6.In contrast, the Γ HI (z) in the 3D simulations rises more rapidly, and is more consistent with the trend of our constraints.In particular, the evolution from the CoDa-III simulation closely reproduces the rapid rise from z ∼ 6 to z ∼ 5.4, although it overshoots to a more highly ionized IGM at z ∼ 5.
In Figure 7 we show the derived constraints on λ mfp .At z = 5.0-5.2 and z = 6.0-6.1, the posteriors show no clear peak interior to the boundaries of the prior.We define 95% lower and upper limits, respectively, by the λ mfp at which the posterior first falls a factor e −2 from its peak value.At z = 5.3 and z = 5.9, the posterior PDFs are clearly peaked, but the posterior PDF does not quite fall below a factor of e −2 from its peak value at the upper and lower edges of the prior boundary, respectively.We show these strongly prior-influenced constraints as open points with error bars.
The open squares in Figure 7 show the posterior medians for (random) ionizing background fluctuations drawn from the 100 Mpc/h box.At redshifts where the mean free path is constrained to be quite small, λ mfp ≲ 30 Mpc, this model prefers a shorter mean free path.This is due to the general decrease in the amplitude of background fluctuations in the smaller ionizing background simulation at fixed λ mfp .Quantitatively, we find that the width of the central 68% distribution of fluctuations (in log Γ) is roughly a factor of two narrower in the 100 Mpc/h box than the 512 Mpc box for λ mfp = 20 Mpc.For larger mean free paths, the situ-  ation can be reversed -for λ mfp ≳ 50 Mpc, the largerscale background fluctuations in the smaller volume exceed those in the larger one due to the truncation of the periodic boundary conditions in the Davies & Furlanetto (2016) method.
The crosses in Figure 7 show the posterior medians for the self-consistent model.In this case, we see a stark factor of ∼ 2 decrease in the preferred λ mfp at all redshifts.This is due to the ionizing background fluctuations now having to overcome the large-scale fluctuations in the density field in order to increase the scatter in the Lyα forest opacity.That is, regions with a stronger ionizing background will preferentially lie in dense environments that will have a higher baseline Lyα opacity, while regions with weak ionizing background will correspond to voids.This then increases the requirements on the fluctuations in the radiation field relative to the uncorrelated case.However, the difference in λ mfp is likely exaggerated by the relatively small volume of the Nyx simulation, 100 Mpc/h, compared to the fiducial 512 Mpc model, which lacks fluctuations in the source field on the largest coherence scales of the radiation field (≳ 100 Mpc, see Figure 2), thus maintaining a stronger correlation between the background and the gas density on the smaller dz = 0.1 ∼ 50 Mpc scale of the τ eff measurements.As we cannot assess the full strength of this effect, and how a much larger (and thus computationally very expensive) self-consistent model might mitigate this offset, we conservatively opt to extend the lower envelope of the uncertainty on λ mfp by the magnitude of the measured offset between the self-consistent and fiducial method constraints from the 100 Mpc/h box.We also show the posterior medians shifted by this offset as grey points.At z = 5.5 where the fiducial analysis using the 100 Mpc/h box results in a larger value for λ mfp due to the additional effect of the truncated boundary conditions mentioned above, we instead adopt the difference between the self-consistent model and the 512 Mpc model.
In the left panel of Figure 8, we compare our constraints to the measurements from Worseck et al. (2014), Becker et al. (2021), andZhu et al. (2023), as well as the lower limit from Bosman (2021).The single power-law evolution with redshift that Worseck et al. (2014) find to be a good fit at z = 2-5 is disfavored at z ≳ 5.4, with our constraints indicating a more rapid decline to higher redshift.Intriguingly, our constraints are in good agreement with the z = 6-5 trend found by Becker et al. (2021) and Zhu et al. (2023), effectively bridging the gaps between the more direct quasar-stacking measurements.In the right panel of Figure 8, we compare our constraints to theoretical predictions for the evolution of λ mfp .The empirical absorber model from Haardt & Madau (2012), also adopted in a similar form by later ionizing background calculations (e.g.Puchwein et al. 2019;Faucher-Giguère 2020), lies well above our constraints at z > 5.3, suggesting that the distribution of H I absorbers must evolve more rapidly at high red-shift.Meanwhile, the sub-grid opacity model of the Cain et al. (2021) simulations as well as the radiationhydrodynamic simulations from Keating et al. (2020), THESAN (Garaldi et al. 2022), and CoDa-III (Lewis et al. 2022) are in much better agreement, with perhaps a modest over-prediction of λ mfp in the Keating et al. (2020) simulations and THESAN relative to our constraints.
Our constraints on Γ HI and λ mfp , including statisticalonly and total uncertainties, are summarized in Table 1.

Consistency between data and model
As in most other Bayesian parameter inference methods, our ability to place constraints on model parameters does not require that the model actually provides a good description of the data.Here we investigate the degree to which the data are consistent with being a draw from our model.
In Figure 9, we show the CDFs of the observed Lyα forest data from Bosman et al. (2022) as blue and red solid lines, corresponding to the "optimistic" and "pessimistic" definitions from Bosman et al. (2018), similar to Figure 3.In the former case, Lyα mean fluxes below twice the (statistical) noise σ F are set to τ eff = − ln (2 × σ F ), while in the latter case, they are assumed to have τ eff > 10.The value of the red curves at the right-hand edge of the panels in Figure 9 thus corresponds to the fraction of Lyα forest sightlines with detected (i.e.> 2σ F ) flux.The shaded regions corre-  1. Derived constraints on the hydrogen photoionization rate ΓHI in units of 10 −12 s −1 and mean free path λ mfp in units of comoving Mpc.The nominal values represent the median of the posterior PDF, with the statistical error (stat.err.) representing the 16th to 84th percentiles of the posterior PDF.The total error (tot.err.) includes contributions from systematic uncertainties (see text for details), and should be considered highly covariant between redshift bins.The self-consistency correction (s-c corr.) is the ratio of the posterior medians of the self-consistent model and the fiducial model; the grey points in Figure 6 and Figure 8 can be recovered by multiplying the fiducial constraints by the value in this column.Upper and lower limits on λ mfp roughly correspond to 2σ limits, see § 4 for details.spond to the central 95% of the distribution of mock data sets drawn from the posterior mean values of Γ HI and λ mfp , where blue and red similarly correspond to the optimistic and pessimistic CDF definitions.We can see that in the majority of cases, the red and blue curves fit neatly within the model distributions, with the only exceptions being a handful of the most highly opaque sightlines at z = 5.6 and z = 5.8.In the highest redshift bins, z = 6.0 and z = 6.1, the blue curves lie close to the upper end of the model distributions -however, this is the expected behavior in the regime where most observations result in non-detections, as the "optimistic" CDF has an upper bound in τ eff given by the distribution of noise values of the quasar spectra, i.e. the CDF of − ln (2 × σ F ).
To further quantify the goodness-of-fit, we compare the ABC distance between the average simulated CDF at the posterior mean values of Γ HI and λ mfp and the observed data to the distribution of distances to 10,000 mock data sets generated from the same model.We show the distribution of the resulting consistency metric in Figure 10.In general, the distance to the data is well within the range of distances to mock data sets.

CAVEATS AND LIMITATIONS OF OUR ANALYSIS
As mentioned previously, our fiducial model for postreionization ionizing background fluctuations suffers from important limitations.In this section, we reiterate and summarize these limitations, and discuss their potential consequences for the interpretation of our results.

Uncorrelated density and ionization rate fields
In the fiducial model, by virtue of using separate volumes for the ionizing background calculation and the hydrodynamical simulation, the density and ionizing radiation fields of the mock spectra used for inference are decoupled.The two should generally be correlated on large scales (e.g.Mesinger & Furlanetto 2009).Because the Lyα forest opacity roughly scales as τ Lyα ∝ ∆ 2 Γ −1 HI (e.g.Weinberg et al. 1997), the variance of τ Lyα should behave as i.e. the presence of correlated fluctuations in the density and ionization rate (σ ∆ 2 Γ ) should suppress the strength of Lyα forest fluctuations at fixed λ mfp .Because the fluctuations become stronger with decreasing λ mfp , we expect our constraints to be biased high.
We tested this hypothesis by constructing a selfconsistent model of ionizing background fluctuations within the smaller 100 Mpc/h volume of our hydrodynamical simulation ( § 2.2.1).As expected, in the regime of strong Lyα forest fluctuations we recover much shorter mean free paths.A fraction of this difference is due to the suppression of the amplitude of fluctuations by the relatively small volume of the Nyx simulation, 100 Mpc/h on a side, compared to the fiducial ionizing background model, 512 Mpc on a side.Without a much larger self-consistent model we are unable to entirely disentangle these two effects, so we instead opt to conservatively allow for a systematic error that encompasses the constraints derived from the self-consistent model, as seen by the large lower error bars in Figure 8.

Uncertainties in the IGM thermal state
As previously mentioned in § 2, our hydrodynamical simulation was run with an optically thin UV background which reionized and heated the volume at very early times, z re ∼ 15.We have chosen to adjust the temperatures in post-processing to better reflect current constraints on the timing and heat injection of reionization.In particular, we assume a fixed heat injection of ∆T = 20, 000 K, which is a representative temperature of the gas after the passage of the ionization front (Miralda-Escudé & Rees 1994;D'Aloisio et al. 2019), and tuned the range of z re to be consistent with the IGM thermal state measured from the distribution of Lyα transmission spike widths by Gaikwad et al. (2020) (see also Bolton et al. 2012).In principle, one could treat ∆T and z re as additional parameters, introduce priors, and add the Gaikwad et al. (2020) measurements to the computation of the likelihood.Due to the computational demands of our likelihood-free inference method, and because in this work we are not attempting to constrain ∆T or z re , we opted to instead impose a plausible range rather than include them directly in the parameter inference.
We note, however, that the IGM thermal state at z > 5 is still substantially uncertain.Measurements using the 1D Lyα forest flux power spectrum by Walther et al. (2019) and Boera et al. (2019) suggest somewhat lower temperatures with a steeper temperature-density relation at z ∼ 5.If the IGM was reionized earlier, or if the heat injection was much lower, the upper envelope of the light grey region in Figure 5 shows that the Γ HI required to reproduce the observed Lyα forest transmission could be substantially higher, with less evolution required from z ∼ 6 to z ∼ 5.
While our hydrodynamical simulation includes the effect of Jeans smoothing on the gas via its fiducial thermal history, by post hoc altering the IGM thermal state of the simulation we neglect the differences in this smoothing that the different thermal histories would have otherwise imprinted (e.g.Gnedin & Hui Rank-order distance z = 5.9 0.0 0.5 1.0 1.5 2.0 2.5 Rank-order distance z = 6.0 0.0 0.5 1.0 1.5 2.0 2.5 Rank-order distance z = 6.1 1998; Peeples et al. 2010a,b;Kulkarni et al. 2015;Nasir et al. 2016;Oñorbe et al. 2017).The effect of Jeans smoothing on measurements of Γ HI was thoroughly investigated over a wide range of thermal histories by Becker & Bolton (2013), who found that at z ≲ 5 the effect was minor, contributing ≲ 0.03 dex to the error budget.We adopt a somewhat higher 0.03 dex systematic uncertainty to account for the trend of larger error at higher redshift, and add this additional error (in quadrature) to our fiducial constraints.

Lack of post-reionization temperature fluctuations
The thermal state of our hydrodynamical simulation assumes a homogeneous heat injection, but the reionization process is patchy with different reionization timing in different locations (Furlanetto et al. 2004).This should lead to a highly inhomogeneous IGM thermal state immediately after reionization is complete, which persists to much later times (e.g.Trac et al. 2008;Oñorbe et al. 2019;Wu et al. 2019;Keating et al. 2018;D'Aloisio et al. 2019).The thermal state in our simulation is thus too uniform, and lacks any (anti-)correlation with large-scale density.Due to the expected anti-correlation between large-scale density and post-reionization temperature fluctuations (D'Aloisio et al. 2015), our Lyα forest fluctuations are likely overestimated at fixed λ mfp (Davies et al. 2018a;D'Aloisio et al. 2018).This would act to bias our λ mfp mea-surements high; with a more realistic simulation, we would require stronger ionizing background fluctuations to counteract the effect of thermal state fluctuations and reproduce the large variations in the Lyα forest, and thus estimate a lower λ mfp .

Parameter choices in the fluctuating ionizing background simulations
Our model for ionizing background fluctuations has several fixed parameters that we have not explored in detail.For example, we assume that λ mfp ∝ ∆ −1 Γ 2/3 HI , but both power-law indices are uncertain.Chardin et al. (2017) found a shallower density dependence λ mfp ∝ ∆ −0.4 when post-processing the Sherwood simulations (Bolton et al. 2017) with the self-shielding prescription from Rahmati et al. (2013).Incorporating this weaker density dependence would lead to stronger ionizing background fluctuations at fixed λ mfp .Our choice of the Γ HI dependence is motivated by McQuinn et al. (2011), but both higher (up to ∼ 1) and lower values (down to ∼ 1/3) are plausible (see the discussion in Becker et al. 2021), which would increase or decrease the λ mfp required to match the observed Lyα forest fluctuations, respectively.
In addition, our model for the ionizing sources assumes that only halos more massive than 2 × 10 9 M ⊙ produce ionizing photons, and that these halos can be assigned a UV luminosity via abundance matching, and further that the ionizing luminosity is proportional to the UV luminosity.All of these assumptions impact the effective bias of the ionizing emissivity field in the fluctuating ionizing background simulations.Other works have adopted different combinations of these assumptionsfor example, Kulkarni et al. (2019) prescribe ionizing luminosities proportional to halo mass for M h > 10 9 M ⊙ , which results in more ionizing photons coming from lower mass halos and thus a lower bias of the emissivity field relative to our approach.The possible parameter space of source models is quite large, requiring far more efficient approaches to statistical inference than ours to constrain these additional parameters (e.g.Qin et al. 2021).

DISCUSSION
With the caveats above in mind, the constraints presented in this work show a substantial improvement in the statistical precision of Γ HI (z) at z > 5, and provide the first quantitative estimates of λ mfp from the excess fluctuations in the Lyα forest alone.In this section we discuss the implications of these measurements for our understanding of the z > 5 IGM.
Since the successful reproduction of the large-scale z ∼ 5-6 Lyα forest variations by Kulkarni et al. (2019), it has commonly been understood that reionization is incomplete at least down to z ∼ 5.5 (and more recently, z ∼ 5.3, cf.Bosman et al. 2022;Zhu et al. 2021Zhu et al. , 2022)).However, here we have demonstrated good agreement between the data and a model which does not require incomplete reionization, which at first glance goes against this consensus.Crucially, due to the lack of full selfconsistency in our model described in Section 2, we cannot make a strong claim that reionization must be complete at z ≲ 6.That said, previous claims of a complete incompatibility between the Lyα forest opacity distribution at z ∼ 5.5-6 and a fully ionized IGM may not be entirely conclusive.
The IGM models employed by works in the literature that have suggested that incomplete reionization is required at z ∼ 5.5-6 are not without their own limitations.The moment-based radiative transfer method used by Kulkarni et al. (2019) and Keating et al. (2020) has been suggested to exhibit suppressed fluctuations in the radiation field at the end of reionization (Wu et al. 2021, see also Gaikwad et al. 2023), thus they may require more large-scale neutral islands from incomplete reionization in order to achieve strong Lyα forest fluctuations.In addition, the spatial resolution of their radiative transfer models (∼ 100 kpc) is too coarse to resolve self-shielding in Lyman limit systems (e.g.McQuinn et al. 2011), potentially resulting in an artificially elongated mean free path in ionized regions and further reducing the strength of fluctuations.The semi-numerical method of Qin et al. (2021) similarly requires a substantial neutral fraction (x HI ∼ 0.15) to reproduce the Lyα forest transmission statistics at z ∼ 5.8 from Bosman et al. (2018), but they employed an approximate spatial filtering approach to compute ionizing background fluctuations.In contrast, the semi-numerical method of Choudhury et al. (2021) employed a somewhat more sophisticated (but similar) treatment of the radiation field, and required roughly half as much neutral gas to explain the same Lyα forest data set.Similar to our findings, Zhu et al. (2021Zhu et al. ( , 2022) ) found that the early reionization/short mean free path model from Nasir & D'Aloisio (2020) is consistent with the distribution of dark gaps in the Lyα and Lyβ forests.
Finally, we emphasize that our results do not in any way rule out the presence of a significant neutral hydrogen fraction in the z < 6 IGM.Rather, they show that imposing the existence of neutral islands is not explicitly required to match the particular summary statistic we are considering, namely the distribution of τ eff on dz = 0.1 scales.In reality, the λ mfp we measure is only representative of the "actual" mean free path if the universe is fully ionized as assumed in our simulations.Our λ mfp may thus be better interpreted as an effective parameter, reflecting both the existence of neutral islands in the deepest, large-scale voids, but also the strong fluctuations in the ionizing background due to a short, fluctuating mean free path of ionizing photons inside of ionized regions.Indeed, in our shortest mean free path models, large-scale regions with the weakest ionizing background are consistent with the gas remaining mostly neutral.However, compared to Gaikwad et al. (2023) who see a similar effect, we estimate lower volume-averaged neutral fractions of ∼ 0.02-7% at z = 5.8 (∼ 0.02-8% at z = 5.9), with a large systematic uncertainty between our fiducial and self-consistent models.We will explore the nature of these regions more closely in future work.

Constraints on the ionizing emissivity
While the mean free path we measure is highly uncertain, and may be more of an effective parameter as described above, we can nevertheless cautiously explore an interpretation of our Γ HI and λ mfp constraints in terms of the ionizing emissivity at the hydrogen-ionizing edge, ϵ 912 = ϵ ν (λ = 912 Å).Assuming the local source approximation (Meiksin & White 2003), the relationship between Γ HI , λ, and ϵ can be expressed as In practice, at z > 5 the integral over frequency effectively extends only to 4 × ν HI due to the onset of strong He II absorption.Approximating the ionizing emissivity as a power law ϵ ν = ϵ 912 (ν/ν HI ) −α , and recalling that we also treat the mean free path as a power law with frequency (equation 4), we numerically solve this expression for ϵ 912 at each redshift using our measured Γ HI and λ mfp values.We assume a fiducial emissivity spectral index of α = 2, following Becker & Bolton (2013).
Our ionizing emissivity estimates are subject to the statistical and systematic uncertainties in the Γ HI and λ mfp measurements described above, as well as an additional systematic uncertainty in the spectral index α assumed for the ionizing sources.Following Becker & Bolton (2013), we adopt a range of α = 1.0-3.0.For our fiducial uncertainty estimates shown in Figure 11, we propagate the statistical uncertainties assuming ϵ 912 ∝ Γ HI /λ mfp , and combine them in quadrature with the systematic uncertainties on Γ HI and λ mfp as well as the systematic uncertainty from the spectral index variations α = 1.0-3.0.For simplicity, here we treat the self-consistent model corrections for Γ HI and λ mfp by inflating the systematic error terms in the upper and lower directions, respectively.
The resulting ionizing emissivities and their corresponding uncertainties are shown in Figure 11, compared to literature constraints from Becker & Bolton (2013), D'Aloisio et al. (2018), and Becker et al. (2021).Note that ϵ 912 is related to the number of ionizing photons emitted per comoving volume Ṅion by Ṅion = hαϵ 912 (e.g.Becker & Bolton 2013), shown on the right axis assuming α = 2.We find that the ionizing emissivity is consistent with a roughly constant ϵ 912 ∼ 10 25 erg s −1 Hz −1 Mpc −3 from z ∼ 6-5, consistent with previous work.

Comparison to XQR-30 analysis by Gaikwad et al.
In a recent complementary work, Gaikwad et al. (2023) have also constrained the photoionization rate and ionizing photon mean free path from the same XQR-30 Lyα forest data that we employ here.Here we compare our results and discuss the differences in analysis methodology.Gaikwad et al. (2023) use a self-consistent model (in our terminology, cf.§ 2.2.1) consisting of smoothedparticle-hydrodynamics (SPH) simulations with 2048 3 baryon and dark matter particles in a volume 160 Mpc/h on a side.They compute ionizing background fluctuations using an independently developed version of the Davies & Furlanetto (2016) method which is heavily optimized by using a tree decomposition of the emissivity field, allowing them to evaluate the radiation field at a much higher resolution (512 3 , ∼ 0.5 Mpc).The statistical comparison between their dense model grid and the XQR-30 Lyα forest data is performed using the Anderson-Darling test on the τ eff CDF, where they set a p-value threshold to map out a 1 − σ contour around the best-fit model.
The largest difference between the two works is primarily a philosophical one - Gaikwad et al. (2023) perform inference on a different definition of the mean free path λ mfp than adopted here.Specifically, while we treat the mean free path as a fully sub-grid quantity that arises from unresolved gas clumping, Gaikwad et al. (2023) measure λ mfp from the density field of their simulation after applying the fluctuating Γ HI field.That is, for a large number of skewers, they evaluate the density of neutral hydrogen n HI given the local gas density and Γ HI , and then calculate the ionizing opacity at the Lyman limit τ HI by integrating the contribution from each resolution element, τ HI (x) = x 0 n HI σ HI (ν HI ) dl. (10) They then compute the corresponding mean free path by fitting an exponential profile to the average transmis-  2023) (blue crosses), where we have shifted the latter slightly in redshift for clarity.The mean trend agrees very well, but the wider assumed range in IGM thermal state parameters leads to larger error bars in Gaikwad et al. (2023).Right: Comparison to the mean free path constraints from Gaikwad et al. (2023).Despite the substantial philosophical differences in the definition of λ mfp (see text for details), our constraints agree quite well.
sion (F (x) = exp [−τ HI (x)] = exp [−x/λ mfp ]) of the simulated skewers.This procedure imposes a physical prior on the possible combinations of Γ HI and λ mfp via the gas density distribution of their hydrodynamical simulation.In addition, this allows Gaikwad et al. (2023) to obtain tight constraints on λ mfp even without an excess in Lyα forest fluctuations, as it can be obtained from the density field even if the ionizing background is entirely uniform.In contrast, our method infers the mean free path solely from the excess in Lyα forest fluctuations over the uniform case, i.e. from the character of the radiation field fluctuations alone.
In Figure 12, we compare our constraints on Γ HI (left) and λ mfp (right) to Gaikwad et al. (2023).In general the two agree very well, suggesting that despite our rather different statistical and modeling methodologies, the rapid evolution in both Γ HI and λ mfp is robustly indicated by the XQR-30 Lyα forest data.In detail, our method recovers a smaller uncertainty on Γ HI , primarily due to the broader range of IGM thermal parameters marginalized over by Gaikwad et al. (2023), while they achieve tighter constraints on λ mfp more uniformly across the full redshift range of their study, subject to the difference in λ mfp definition described above.

CONCLUSION
In this work, we have constrained the evolution of the photoionization rate Γ HI and the mean free path of ionizing photons λ mfp in the IGM at z = 5.0-6.1 using the extended XQR-30 (D'Odorico et al. 2023) Lyα forest data set from Bosman et al. (2022).We assume that the excess fluctuations in the z > 5.3 Lyα forest are due to a strong coupling between the ionizing background and mean free path (Davies & Furlanetto 2016), and constrain parameters using the likelihood-free inference technique of approximate Bayesian computation, or ABC.
We recover a smooth evolutionary trend in Γ HI , which increases by a factor of ∼ 4 from z = 6 to z = 5; in agreement with past observations but with a finer redshift sampling and smaller statistical uncertainty.The increase is much stronger than predicted by empirically-motivated 1D cosmological radiative transfer models, but is largely consistent with the evolution found by state-of-the-art 3D radiation-hydrodynamic cosmological simulations.We similarly find consistency with recent measurements of the mean free path λ mfp from stacked quasar spectra and recent radiationhydrodynamic models.
We find that the statistical constraining power of the coarsely-binned τ eff measurements is very strong, to the extent that for both Γ HI and λ mfp we are strongly limited by systematic uncertainties.For the former, we are limited by our knowledge of the IGM thermal state, while for the latter (and the former, to a lesser extent) we are limited by the computational resources required to simulate a converged Lyα forest in a large enough volume to capture the full intensity of ionizing background fluctuations.Additionally, the relic thermal fluctuations left by reionization at earlier times (e.g.D' Aloisio et al. 2015) will also confound efforts to precisely constrain the parameters of the z = 5-6 IGM without a complete model of the entire reionization process.Reducing these systematic uncertainties will require extracting information from both the small-scale (e.g.Gaikwad et al. 2020) and large-scale properties of the Lyα forest, and then performing statistical inference on full reionization lightcones with a detailed model for the sources and sinks (e.g.Qin et al. 2021).
At the low redshift end of our data set, z ∼ 5, our ability to constrain λ mfp is limited by the insensitivity of the Lyα forest to long mean free paths, as the resulting fluctuations in the radiation field become weak compared to the intrinsic fluctuations already imprinted by the IGM density field.However at these redshifts the direct stacking of quasar spectra appears to be a far more sensitive probe with fewer underlying assumptions and systematic uncertainties (Worseck et al. 2014;Becker et al. 2021;Zhu et al. 2023).At the high redshift end, z ∼ 6, we are instead limited by the relative sparsity of sightlines and the poor spatial resolution of our ionizing background simulations.The latter can already be rectified by adopting more efficient methods of computing the radiation field (e.g.Gaikwad et al. 2023), while the former will require additional deep spectroscopy of yet higher redshift quasars than the XQR-30 sample (e.g.Yang et al. 2020b

Figure 1 .
Figure 1.Evolution of the temperature at mean density in the IGM (top) and the slope of the temperature-density relation (bottom).The curves show the analytic models described in the text for our fiducial range of reionization heating models, while the data points show measurements from Walther et al. (2019), Boera et al. (2019), and Gaikwad et al. (2020).

Figure 2 .
Figure 2. Slices through a subset of our fluctuating ΓHI fields, 4 Mpc-thick and 512 Mpc on a side, showing a representative range of λ mfp values of 10, 20, 40, and 80 Mpc from left to right.
4. CONSTRAINTS ON Γ HI AND λ MFP FROM THE EXTENDED XQR-30 DATA SET

Figure 3 .
Figure 3. Demonstration of ABC on a mock data set.The blue and red curves in the left panel correspond to the cumulative distribution functions of Lyα forest effective optical depth from a mock XQR-30 data set, with non-detections treated as optimistic (F = 2 × σF ) or pessimistic (F = 0), respectively, following Bosman et al. (2018, 2022).The shaded regions correspond to the central 95% of the scatter of additional mock data sets with ΓHI and λ mfp set to their mean posterior estimates.The grey points in the upper halves of the middle and right panels show the distance metric ρ (equation 8) computed from 1,000,000 mock data sets.The horizontal lines show three different ρ thresholds below which lie 50%, 10%, and 0.1% of the mock data samples from top to bottom.The corresponding posterior PDFs on ΓHI and λ mfp are shown in the lower panels, with the true values indicated by the vertical dashed lines.et al. 2011), primarily sourced from the Extended XQR-30 sample (D'Odorico et al. 2023), as well as 16 spectra taken with the ESI(Sheinis et al. 2002) spectrograph, i.e. 67 quasar spectra in total.All quasars are located at z > 5.7 and their spectra possess signal-to-noise ratios (SNRs) larger than 10 per pixel.Intervening damped Lyman-α absorbers are masked based on the detection of associated metal absorbers in the metal catalog ofDavies et al. (2023).In a few quasars, portions of the spectra with suspected broad absorption line contamination are also masked.We refer interested readers to(Bosman et al. 2022) for more details of the sample.In order to measure the effective optical depths, the intrinsic quasar continua are reconstructed based on their un-absorbed emission features at λ rest > 1215.16Å.We use the principal component analysis (PCA) approach ofDavies et al. (2018d)  with the improvements described inBosman et al. (2020) andBosman et al. (2022).The PCA approach reproduces the quasar continua at λ rest < 1180 Å with a wavelength-dependent ∼ 8% accuracy, and sub-percent precision.Finally, the effective optical depths are measured in bins of dz = 0.1 and only bins which are > 50% un-masked are retained.The number of retained quasar sightlines N q at each redshift is shown in Table1.The masking, observational uncertainties, and wavelength-dependent continuum uncertainties are forward-modeled in all simulated mocks.Bosman et al. (2022) showed that mock data generated in this manner is statistically indistinguishable from the

ΓFigure 4 .
Figure 4. Joint posterior PDF of ΓHI and λ mfp at z = 5.0-6.1.The inner and outer contours contain 68% and 95% of the distribution, respectively.The blue dashed line shows the lower edge of the λ mfp prior.

Figure 5 .
Figure 5. Posterior medians (black circles) and central 68% credible intervals (black thin error bars) on ΓHI from the XQR-30 data set assuming zre = 7.2.The dark grey shaded region shows the deviation of the posterior medians for zre = 6.7 and zre = 7.7, while the light grey shaded region shows the range from more extreme thermal models with zre = 6.2 and zre ∼ 15.The open orange points show the constraints without the correction for the coarse redshift snapshot sampling (see Appendix A).The purple crosses show the posterior medians from the self-consistent model in the L = 100 Mpc/h hydrodynamical simulation volume, see § 2.2.1.

Figure 7 .
Figure 7. Posterior medians (black circles) and central 68% credible intervals (black thin error bars) on λ mfp from the XQR-30 data set, with open circles corresponding to constraints and arrows corresponding to 2σ limits (see text for details).Open orange points show the posterior medians from the 100 Mpc/h simulation with uncorrelated density and ionizing background, while the purple crosses show the posterior medians from the fully self-consistent simulation.

Figure 8 .
Figure 8. Left: Our constraints on λ mfp compared to observations are shown by the black points, while the grey points show the constraints shifted to approximately correct for a bias due to the lack of self-consistency between the density field and ionizing background.Measurements from Worseck et al. (2014) (green squares), Becker et al. (2021) (purple pentagons), Zhu et al. (2023) (red diamonds), and the lower limit from Bosman (2021) (orange arrow) are shown for comparison.The dotted line shows the extrapolated fit to the full redshift range of measurements from Worseck et al. (2014).Right: Comparison to the mean free path in theoretical models, either prescribed by extrapolation of empirical constraints on the column density distribution of hydrogen absorbers (Haardt & Madau 2012, purple) or directly measured in the simulations by Keating et al. (2020) (brown), Cain et al. (2021) (blue), Garaldi et al. (2022) (orange; THESAN), and Lewis et al. (2022) (red; CoDa-III).

Figure 9 .
Figure 9.Comparison between the Lyα forest τ eff CDFs from Bosman et al. (2022) (solid curves) and the central 95% of the distribution of mock CDFs from our simulations when adopting the posterior mean values of ΓHI and λ mfp (shaded regions).The red color corresponds to CDFs where Lyα fluxes detected at less than 2σ significance are assumed to be zero (i.e.τ eff = ∞), while the blue color corresponds to those fluxes being set to twice the statistical noise (i.e., τ eff = − ln(2 × σN )).

Figure 10 .
Figure10.Distribution of distances (equation 8) between the average set of rank-order τ eff and individual sets from 10,000 mock observations, both adopting the posterior mean values of ΓHI and λ mfp at each redshift.The vertical orange lines show the distance between the average simulated rank-order τ eff and the observed Lyα forest data.

Figure 12 .
Figure12.Left: Our constraints on ΓHI (black and grey circles) compared toGaikwad et al. (2023) (blue crosses), where we have shifted the latter slightly in redshift for clarity.The mean trend agrees very well, but the wider assumed range in IGM thermal state parameters leads to larger error bars inGaikwad et al. (2023).Right: Comparison to the mean free path constraints fromGaikwad et al. (2023).Despite the substantial philosophical differences in the definition of λ mfp (see text for details), our constraints agree quite well.

Table
). SEIB is funded by the Deutsche Forschungsgemeinschaft (DFG) under Emmy Noether grant number BO 5771/1-1.GK is partly supported by the Department of Atomic Energy (Government of India) research project with Project Identification Number RTI 4002, and by the Max Planck Society through a Max Planck Partner Group.For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.