First constraints on the local ionization topology in front of two quasars at z ~ 7.5

TL;DR

This paper introduces a Bayesian inference framework to analyze JWST spectra, providing new constraints on the local ionization topology and quasar properties at z ~ 7.5, revealing short quasar lifetimes and partial IGM neutrality.

Contribution

It presents a novel Bayesian HMC method to extract local HI distribution and quasar parameters from damping wing signatures in high-redshift quasar spectra.

Findings

J1007+2115 likely in a 32% neutral IGM

J1342+0928 likely in a 58% neutral IGM

Quasars have short to intermediate lifetimes

Abstract

Thus far, Lyman-$\alpha$ damping wings towards quasars have been used to probe the \textit{global} ionization state of the foreground intergalactic medium (IGM). A new parameterization has demonstrated that the damping wing signature also carries \textit{local} information about the distribution of neutral hydrogen (HI) in front of the quasar before it started shining. Leveraging a recently introduced Bayesian \texttt{JAX}-based Hamiltonian Monte Carlo (HMC) inference framework, we derive constraints on the Lorentzian-weighted HI column density $N_\mathrm{HI}^\mathrm{DW}$, the quasar's distance $r_\mathrm{patch}$ to the first neutral patch and its lifetime $t_\mathrm{Q}$ based on JWST/NIRSpec spectra of the two $z \sim 7.5$ quasars J1007+2115 and J1342+0928. After folding in model-dependent topology information, we find that J1007+2115 (and J1342+0928) is most likely to reside in a $\langle x_\mathrm{HI} \rangle = 0.32_{-0.20}^{+0.22}$ ($0.58_{-0.23}^{+0.23}$) neutral IGM while shining for a remarkably short lifetime of $\log_{10} t_\mathrm{Q} /\mathrm{yr} = 4.14_{-0.18}^{+0.74}$ (an intermediate lifetime of $5.64_{-0.43}^{+0.25}$) along a sightline with $\log_{10} N_\mathrm{HI}^\mathrm{DW} /\mathrm{cm}^{-2} = 19.70_{-0.86}^{+0.35}$ ($20.24_{-0.22}^{+0.25}$) and $r_\mathrm{patch} = 28.9_{-14.4}^{+54.0} \,\mathrm{cMpc}$ ($10.9_{-5.9}^{+5.6} \,\mathrm{cMpc}$). In light of the potential presence of local absorbers in the foreground of J1342+0928 as has been recently suggested, we also demonstrate how the Lorentzian-weighted column density $N_\mathrm{HI}^\mathrm{DW}$ provides a natural means for quantifying their contribution to the observed damping wing signal.

Equations

$$N_{\mathrm{HI}}^{\mathrm{DW}}\equiv$$

$$\times {\left(\frac{1+z_{\mathrm{QSO}}}{1+7.54}\right)}^2{\int }_{r_{\min }/r_{\mathrm{T}}}^{r_{\max }/r_{\mathrm{T}}}\frac{x_{\mathrm{HI}}(r)\cdot \Delta (r)}{(r/r_{\mathrm{T}}+1{)}^2}\mathrm{d}(r/r_{\mathrm{T}}),$$

$$L(\mathbf{f}|\mathbf{\sigma },\mathbf{\theta },\mathbf{\eta })=\mathcal{N}\left(\mathbf{f};⟨\mathbf{t}⟩\circ ⟨\mathbf{s}⟩,\mathbf{\Sigma }+⟨\mathbf{S}⟩{\mathbf{C}}_{\mathbf{t}}⟨\mathbf{S}⟩+⟨\mathbf{T}⟩{\mathbf{C}}_{\mathbf{s}}⟨\mathbf{T}⟩\right),$$

$$P_{\mathrm{top}}$$

$$\times P_{\mathrm{top}}({\log }_{10}{N}_{\mathrm{HI}}^{\mathrm{DW}}/{\mathrm{cm}}^{-2},r_{\mathrm{patch}}|⟨x_{\mathrm{HI}}⟩)/P(\mathbf{f}|\mathbf{\sigma }).$$

$$N_{\mathrm{HI}}^{\mathrm{DW}}=N_{\mathrm{HI}}^{\mathrm{DW},\mathrm{IGM}}+N_{\mathrm{HI}}^{\mathrm{DW},\mathrm{abs}}$$

$$\log N_{\mathrm{HI}}^{\mathrm{DW},\mathrm{abs}}=\log N_{\mathrm{HI}}^{\mathrm{abs}}-2\log (r_{\mathrm{abs}}/r_{\mathrm{T}}+1),$$

Full text

††pubyear: 2025††pagerange: First constraints on the local ionization topology in front of two quasars at $z\sim 7.5$–References

1 Introduction

The highly sensitive Lyman-$\alpha$ transition observed in the spectra of high-redshift sources carries a wealth of information about the epoch of reionization, most particularly when the absorption imprint from the foreground neutral intergalactic medium (IGM) saturates, manifesting in a characteristic damping wing signature redward of Lyman-$\alpha$ line center (Miralda-Escudé, 1998). Quasars (e.g. Mortlock et al., 2011; Greig et al., 2017, 2019; Greig et al., 2022; Davies et al., 2018; Wang et al., 2020; Yang et al., 2020; Ďurovčíková et al., 2024) as well as galaxies (e.g. Curtis-Lake et al., 2023; Mason et al., 2025) have been used as background sources to infer constraints on the global volume-averaged fraction $\langle x_{\mathrm{HI}}\rangle$ of neutral hydrogen (HI) at the redshift of the source. Stacking spectra at a similar redshift, or statistically combining the inferred $\langle x_{\mathrm{HI}}\rangle$ constraints makes IGM damping wings a natural probe of the timing of reionization. But recent studies have constructed parameterizations that capture the shape of the IGM damping wing significantly more tightly than the global IGM neutral fraction $\langle x_{\mathrm{HI}}\rangle$, and even carry information about the local ionization topology in front of a given source (Chen, 2024; Keating et al., 2024; Mason et al., 2025; Kist et al., 2025b).

Specifically, Kist et al. (2025b) put forward a three-parameter model applicable in the context of quasars, consisting of two summary statistics of the $x_{\mathrm{HI}}$ field before the quasar started shining (which we dub the pre-quasar $x_{\mathrm{HI}}$ field), and the quasar’s lifetime $t_{\mathrm{Q}}$ encoding the effects of its ionizing radiation. The former two statistics of the pre-quasar ionization topology are 1) the column density $N_{\mathrm{HI}}^{\mathrm{DW}}$ of neutral material in front of the quasar, weighted by a Lorentzian profile mimicking the frequency dependence of the Lyman-$\alpha$ cross section $\sigma_{\alpha}$ appearing in the optical depth integral; and 2) the distance $r_{\mathrm{patch}}$ from the quasar to the first neutral patch. In a companion paper, Kist et al. (2025a) introduced a fully Bayesian pipeline to infer these two parameters along with the quasar lifetime $t_{\mathrm{Q}}$ in a reionization model-independent fashion based on observed high-redshift quasar spectra, accounting for all relevant sources of uncertainty such as the unknown intrinsic continuum of the quasar, IGM transmission fluctuations and spectral noise. Statistical tests on hundreds of mock spectra demonstrated that the pipeline allows them to constrain $t_{\mathrm{Q}}$ to $0.58_{-0.13}^{+0.13}\,\mathrm{dex}$, $N_{\mathrm{HI}}^{\mathrm{DW}}$ to $0.72_{-0.25}^{+0.08}\,\mathrm{dex}$ and $r_{\mathrm{patch}}$ to $35.7_{-31.9}^{+7.8}\,\mathrm{cMpc}$ in case a noticeable damping wing is present in the spectrum. Here we take advantage of this framework to infer the first local IGM damping wing constraints from JWST/NIRSpec spectra of two of the highest-redshift quasars known to date, that is, J1007+2115 at $z=7.51$ (Yang et al., 2020) and J1342+0928 at $z=7.54$ (Bañados et al., 2018).

We briefly summarize in Section 2 the underlying theory and our local damping wing analysis framework, and proceed in Section 3 by introducing the data and presenting our analysis results for the two aforementioned objects. We put these measurements into context with existing literature constraints and conclude in Section 4.

2 Methods

We conduct our analysis in the context of the local damping wing parameterization put forward by Kist et al. (2025b), harnessing the inference pipeline originally introduced in Hennawi et al. (2025b) and adapted by Kist et al. (2025a) for the use with this parameterization. Kist et al. (2025a) also conceptualized a theoretical framework for tying the inferred local parameter constraints which are agnostic to the underlying reionization topology to a specific reionization model to obtain a global constraint on the timing of reionization. We restrict this section to a brief but self-contained summary of each of these analysis components and refer the reader to the works above for more comprehensive descriptions of the respective parts.

The first of the two local summary statistics of the pre-quasar ionization topology defined by Kist et al. (2025b) that we aim to constrain in this work is the Lorentzian-weighted HI column density

[TABLE]

where $r$ denotes the comoving distance, and the denominator in the integrand is the Lorentzian weighting kernel which accounts for the frequency dependence of the Lyman-$\alpha$ cross section that the neutral hydrogen field $x_{\mathrm{HI}}(r)$ and the overdensity field $\Delta(r)$ are integrated against in the optical depth integral. The three distances $r_{\mathrm{min}}$, $r_{\mathrm{max}}$ and $r_{\mathrm{T}}$ are hyperparameters which we fix to $r_{\mathrm{min}}=4\,\mathrm{cMpc}$, $r_{\mathrm{max}}=r_{\mathrm{min}}+100\,\mathrm{cMpc}$ and $r_{\mathrm{T}}=18\,\mathrm{cMpc}$ throughout (for details see Kist et al., 2025b). Our second statistic $r_{\mathrm{patch}}$ can be straightforwardly defined as the (comoving) distance between the source and the first neutral patch. Recall that both $N_{\mathrm{HI}}^{\mathrm{DW}}$ and $r_{\mathrm{patch}}$ are summaries of the pre-quasar ionization topology $x_{\mathrm{HI}}$ as this is the cosmological field we are aiming to constrain. Our third parameter, the lifetime $t_{\mathrm{Q}}$ of the quasar, summarizes the effects of the quasar’s ionizing radiation.

We simulate IGM transmission profiles $\boldsymbol{t}$ based on a hybrid approach following Davies et al. (2018) where we combine density, velocity and temperature skewers from the Nyx hydrodynamical simulations (Almgren et al., 2013; Lukić et al., 2015) and synthetic $x_{\mathrm{HI}}$ skewers (Kist et al., 2025b) generated at $21$ column density values between $17.48\leq\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2}\leq 21.08$ and $18$ different neutral patch distances between $0.3\,\mathrm{cMpc}\leq r_{\mathrm{patch}}\leq 143.0\,\mathrm{cMpc}$. We then perform one-dimensional radiative transfer (Davies et al., 2016) assuming $51$ quasar lifetimes between $10^{3}\,\mathrm{yr}\leq t_{\mathrm{Q}}\leq 10^{8}\,\mathrm{yr}$, and forward-model realistic instrumental effects (Hennawi et al., 2025b). All our models are generated assuming $z_{\mathrm{QSO}}=7.54$ with an ionizing photon emission rate of $Q=10^{57.14}\,\mathrm{s}^{-1}$, explicitly resembling J1342+0928 but also appropriate for J1007+2115. Further, to tie our local constraints to the global IGM neutral fraction $\langle x_{\mathrm{HI}}\rangle$ in the way described in Kist et al. (2025a), we determine the stochastic mapping $P_{\mathrm{top}}(\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2},r_{\mathrm{patch}}\,|\,\langle x_{\mathrm{HI}}\rangle)$,111Throughout, we use a ’top’ subscript to denote probability distributions that are defined in the context of a specific reionization model. by combining the same hydrodynamical Nyx sightlines with $x_{\mathrm{HI}}$ skewers extracted from semi-numerical reionization topologies simulated with a modified version of 21cmFAST (Mesinger et al., 2011; Davies & Furlanetto, 2022) at $21$ different global IGM neutral fractions between $0\leq\langle x_{\mathrm{HI}}\rangle\leq 1$.

In combination with a principal component analysis (PCA) model for the intrinsic quasar continuum $\boldsymbol{s}$ constructed based on a dataset of $44,587$ quasar spectra at $1.878<z<3.427$ from the SDSS-III Baryon Oscillation Spectroscopic Survey (BOSS) and SDSS-IV Extended BOSS (eBOSS) (Hennawi et al., 2025a; Kist et al., 2025a), this allows us to write down the likelihood of an observed spectrum $\boldsymbol{f}$ given the astrophysical parameters $\boldsymbol{\theta}$, the PCA coefficients $\boldsymbol{\eta}$ and its observational noise vector $\boldsymbol{\sigma}$ as the following multivariate Gaussian distribution (Hennawi et al., 2025b):

[TABLE]

where $\boldsymbol{t}\circ\boldsymbol{s}$ is the element-wise (Hadamard) product of the two mean vectors $\langle\boldsymbol{t}\rangle$ and $\langle\boldsymbol{s}\rangle$, and $\boldsymbol{C}_{\boldsymbol{t}}$ and $\boldsymbol{C}_{\boldsymbol{s}}$ are the covariance matrices of $\boldsymbol{t}$ and $\boldsymbol{s}$, and we write $\boldsymbol{\Sigma}\equiv\mathrm{diag}(\boldsymbol{\sigma})$, $\boldsymbol{T}\equiv\mathrm{diag}(\boldsymbol{t})$ and $\boldsymbol{S}\equiv\mathrm{diag}(\boldsymbol{s})$.

We follow Kist et al. (2025a) and initially infer the local parameters $\boldsymbol{\theta}\equiv(\log_{10}t_{\mathrm{Q}}/\mathrm{yr},\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2},r_{\mathrm{patch}})$ under the assumption of a logarithmically flat prior on the lifetime between $3\leq\log_{10}t_{\mathrm{Q}}/\mathrm{yr}\leq 8$, and a constant two-dimensional topology-agnostic prior on $(\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2},r_{\mathrm{patch}})$ with a non-trivial boundary enclosing all physically permitted regions of parameter space (for details see Kist et al., 2025a). After marginalizing out the nuisance parameters $\boldsymbol{\eta}$, we can probabilistically tie these local constraints on $\boldsymbol{\theta}$ to the aforementioned semi-numerical reionization topology via the conditional distribution $P_{\mathrm{top}}(\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2},r_{\mathrm{patch}}\,|\,\langle x_{\mathrm{HI}}\rangle)$ and obtain a joint constraint on the global IGM neutral fraction:

[TABLE]

For comparison, we also perform the inference following the conventional approach without local summary statistics where $\boldsymbol{\theta}\equiv(\langle x_{\mathrm{HI}}\rangle,\log_{10}t_{\mathrm{Q}}/\mathrm{yr})$ such that no conversion according to Eq. (2) is necessary. We henceforth refer to this as the global parameterization (as opposed to our local one). Here we assume a flat neutral fraction prior between $0\leq\langle x_{\mathrm{HI}}\rangle\leq 1$.

Practically, we sample from the respective posterior distribution via the NumPyro Hamiltonian Monte-Carlo (HMC) implementation with No U-Turn Sampler (NUTS). Each inference run consists of $8$ HMC chains with $1000$ warm-up and $2000$ sampling steps per chain (for details, see Hennawi et al., 2025b; Kist et al., 2025a). We reweight these samples based on the coverage tests performed in Kist et al. (2025a) to ensure that we are quoting statistically faithful constraints.

3 Local IGM damping wing constraints at ${z}\sim 7.5$

We now proceed by leveraging the framework summarized in the previous section to infer the first local IGM damping wing constraints for two of the highest-redshift quasars known to date. That is, J1007+2115 at $z=7.51$ (Yang et al., 2020) and J1342+0928 at $z=7.54$ (Bañados et al., 2018) with absolute magnitudes of $M_{1450}=-26.82$ and $M_{1450}=-26.34$ at $1450\,\text{\AA }$ in the rest frame, respectively, determined based on Euclid photometry (Euclid Collaboration: Yang et al., 2025).

This work presents the first damping wing analysis of JWST/NIRSpec spectra of these two objects (JWST GO 1764, PI: Fan, and JWST GTO 1219, PI: Luetzgendorf; see also Christensen et al., 2023; Hennawi et al., 2025a). The spectra have signal-to-noise ratios of $\mathrm{S}/\mathrm{N}\simeq 5-20$ at a resolution of $R\simeq 2700$ and were taken with the NIRSpec Fixed Slits together with the G140H/F070LP and G235H/F170LP grating and filter combinations. The data were reduced via a combination of the JWST Science calibration pipeline CALWEBB (version 1.13.4) and the python-based semi-automated reduction pipeline PypeIt (Prochaska et al., 2020). The full details of the reduction procedure are discussed in Hennawi et al. (2025a).

The spectral coverage of these observations allows us to operate on a rest-frame wavelength range of $1175\,\text{\AA }\leq\lambda\leq 3000\,\text{\AA }$. After removing narrow absorbers, we rebin both spectra to a velocity pixel scale of $\mathrm{d}v=500\,\mathrm{km}/\mathrm{s}$ (and likewise are our forward-modeled IGM transmission profiles) as our current likelihood prescription does not allow us to extract all information from the full resolution data at its native pixel scale in a statistically faithful manner (Kist et al., 2025c, a) which is pending improvements in our likelihood prescription, for example through simulation-based inference (Chen, 2024).

3.1 J1007+2115

We start by analyzing the quasar J1007+2115 at $z=7.51$. Previous studies have inferred IGM neutral fractions of $\langle x_{\mathrm{HI}}\rangle=0.39_{-0.13}^{+0.22}$ (Yang et al., 2020) and $\langle x_{\mathrm{HI}}\rangle=0.27_{-0.17}^{+0.21}$ (Greig et al., 2022) from ground-based spectra of this object. Note that the former analysis also obtained a constraint of $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.89_{-0.76}^{+1.15}$ on its lifetime (Eilers et al., 2021), and used, upon some minor differences, the same simulation models as we do in this work, however with a significantly different analysis pipeline. Our re-analysis is based on the new JWST/NIRSpec spectrum of this object, and we perform the analysis both in the context of the conventional global as well as our local IGM damping wing parameterization. We depict the final rebinned spectrum that forms the input to our pipeline as the black line in Figure 1 with associated noise vector in yellow. The upper panel depicts our reconstruction in the context of the local parameterization, and the lower panel in the context of the global one. The median reconstructed continua are shown in blue, and the full models including IGM absorption in red. The shading around these lines highlights the associated $16\,\%$ and $84\,\%$ percentile variations, reflecting parameter uncertainty, continuum reconstruction errors, as well as spectral noise (for further details see Hennawi et al., 2025b).

Overall, the models show a remarkable agreement, both providing a good fit to the spectrum across the entire spectral range, not only matching the shape of the proximity zone and the Lyman-$\alpha$ damping wing, but also the smooth emission lines redward of Lyman-$\alpha$. The median models slightly undershoot the observed spectrum around $\lambda\approx 1225\,\text{\AA }$ (somewhat more pronounced in the global parameterization), and slightly overshoot at $\lambda\approx 1245\,\text{\AA }$. However, in both regions the observed value is still within the $68\,\%$ region. Somewhat more redward, at $\lambda\approx 1265\,\text{\AA }$, three pixels got rejected by the sigma-clipping procedure we apply prior to the inference (see Kist et al., 2025a). These pixels are possibly affected by a mild broad absorption line (BAL) system but this does not appear to impact our conclusions as we obtain the same results when no clipping is applied.

The only main difference between the global and the local model curves is a reduced degree of scatter in the local parameterization as this parameterization excludes the scatter due to the stochastic distribution of neutral patches along the line-of-sight from the inference task. Both models only show a very mild damping wing. The black contours in Figure 2 depict the associated local parameter constraints. The lifetime posterior peaks relatively sharply at $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.09_{-0.15}^{+0.36}$, implying that we are concerned with a comparably young object, as already hinted at by its small proximity zone (see e.g. also the objects in Eilers et al., 2020, 2021). The data are not overly constraining with respect to the two local parameters $N_{\mathrm{HI}}^{\mathrm{DW}}$ and $r_{\mathrm{patch}}$. In essence, only the highest HI column densities (and correspondingly short neutral bubble distances that are still located in the physically permitted region of parameter space enclosed by the dashed line in the $(N_{\mathrm{HI}}^{\mathrm{DW}},r_{\mathrm{patch}})$ panel) are excluded as these would lead to a noticeable damping wing imprint.

We obtain a global $\langle x_{\mathrm{HI}}\rangle$ constraint on the timing of reionization based on these local constraints by following the procedure introduced in Kist et al. (2025a) and summarized in Section 2, folding in the distribution of our local parameters within the realistic semi-numerical reionization topology as a prior. The green contours in Figure 2 show the resulting posterior distribution. We see that the prior on $N_{\mathrm{HI}}^{\mathrm{DW}}$ and $r_{\mathrm{patch}}$ induced by this topology (shown explicitly in purple in the extra panel of the plot) disfavors the (low-$N_{\mathrm{HI}}^{\mathrm{DW}}$, low-$r_{\mathrm{patch}}$) region of the posterior and therefore gives more weight to the (low-$N_{\mathrm{HI}}^{\mathrm{DW}}$, high-$r_{\mathrm{patch}}$) and (high-$N_{\mathrm{HI}}^{\mathrm{DW}}$, low-$r_{\mathrm{patch}}$) peaks.222Note that the green converted posterior seemingly has support in regions where the black unconverted one does not. In fact, both posteriors have support in these regions, but before performing the conversion, these are not included in the $95\,\%$ contour shown in the plot. Further, due to smoothing effects and the statistical reweighting procedure we apply, the contours marginally extend beyond the physical $(N_{\mathrm{HI}}^{\mathrm{DW}},r_{\mathrm{patch}})$ boundary enclosed by the dashed line. The resulting $\langle x_{\mathrm{HI}}\rangle$ posterior suggests a comparably low to intermediate neutral fraction of $\langle x_{\mathrm{HI}}\rangle=0.32_{-0.20}^{+0.22}$ at this redshift. The lifetime constraint largely remains unaffected by folding in the topology information, still peaking at $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.14_{-0.18}^{+0.74}$, on the lower end of the literature value of $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.89_{-0.76}^{+1.15}$ (Yang et al., 2020; Eilers et al., 2021). Our local parameter constraints, on the other hand, are fully prior-dominated, as a comparison of the green and purple contours in the extra $(N_{\mathrm{HI}}^{\mathrm{DW}},r_{\mathrm{patch}})$ panel in the same figure suggests.

We test the robustness of the $\langle x_{\mathrm{HI}}\rangle$ constraint by also comparing it to the one we obtain when inferring $\langle x_{\mathrm{HI}}\rangle$ and $t_{\mathrm{Q}}$ directly in the context of the global parameterization. The corresponding model curves are shown in the lower panel of Figure 1. The resulting posterior, depicted by the yellow contours in Figure 2, is in good agreement with the one obtained by converting our local constraints (green). The only differences between the two are that the converted posterior exhibits a somewhat more extended tail towards longer lifetimes. Since these longer lifetimes go hand in hand with a higher IGM neutral fraction, the marginal $\langle x_{\mathrm{HI}}\rangle$ posterior obtained via the local parameterization (green) shows a very mild peak at around $\langle x_{\mathrm{HI}}\rangle\simeq 0.3$, in contrast to the directly inferred one which peaks very close to $\langle x_{\mathrm{HI}}\rangle\simeq 0$. Overall, however, both posteriors are in clear statistical agreement, underlining the robustness of our constraints despite the highly different ways in which they were obtained.

3.2 J1342+0928

The second object of our study, the quasar J1342+0928 at $z=7.54$, has already been extensively analyzed in previous literature (Bañados et al., 2018; Davies et al., 2018; Greig et al., 2019; Ďurovčíková et al., 2020; Reiman et al., 2020; Greig et al., 2022), resulting in varying constraints on the IGM neutral fraction, with median values ranging from $\langle x_{\mathrm{HI}}\rangle\approx 0.2-0.6$. We again apply both versions of our inference framework to the JWST/NIRSpec spectrum of this object and obtain the fits and parameter constraints shown in Figures 3 and 4, again in good statistical agreement with a reduced degree of scatter in the local parameterization. Unlike in the case of J1007+2115, we infer the clear presence of a damping wing with a correspondingly high HI column density of $\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2}=19.99_{-0.40}^{+0.27}$ and a rather short distance of $r_{\mathrm{patch}}=12.0_{-7.4}^{+5.0}\,\mathrm{cMpc}$ to the first neutral patch, increasing/reducing further to $\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2}=20.24_{-0.22}^{+0.25}$ and $r_{\mathrm{patch}}=10.9_{-5.9}^{+5.6}$ after folding in the topology dependence (green contours). The lifetime posterior shows a clear degeneracy with $N_{\mathrm{HI}}^{\mathrm{DW}}$ in the local parameterization, peaking at $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=5.64_{-0.43}^{+0.25}$, and with $\langle x_{\mathrm{HI}}\rangle$ in the global one, preferring a somewhat lower value of $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.69_{-0.33}^{+1.40}$. However, the global posterior shows a long axis of degeneracy with a second peak closer to the higher value preferred in the local model. Previous analyses inferred $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=5.38_{-1.30}^{+0.72}$ with a similarly wide degeneracy (Davies et al., 2018, 2019; Eilers et al., 2021). The corresponding $\langle x_{\mathrm{HI}}\rangle$ constraints suggest an intermediate to high IGM neutral fraction of $\langle x_{\mathrm{HI}}\rangle=0.58_{-0.23}^{+0.23}$ in the context of the local parameterization, and a very high one of $\langle x_{\mathrm{HI}}\rangle=0.78_{-0.30}^{+0.12}$ in the context of the global model.

These results ought to be treated with caution since Davies et al. (2025) recently pointed out the possibility of contamination of the IGM in front of J1342+0928 by proximate DLA absorbers located at $r_{\mathrm{abs}}\simeq 21\,\mathrm{cMpc}$, $31\,\mathrm{cMpc}$ and/or $56\,\mathrm{cMpc}$ based on weak Mg II absorption lines they identified in the JWST/NIRSpec spectrum that is also the subject of this analysis. While the authors concluded that such proximate absorbers would have to be unusually metal-poor, we use this possibility to highlight an additional virtue of our local parameterization: including the (Lorentzian-weighted) HI column density $N_{\mathrm{HI}}^{\mathrm{DW}}$ as an explicit model parameter allows for a straightforward comparison of our constraints to the (unweighted) column density $N_{\mathrm{HI}}^{\mathrm{abs}}$ of a putative absorber. The Lorentzian-weighted version $N_{\mathrm{HI}}^{\mathrm{DW,abs}}$ of the absorber’s column density then combines with the column density $N_{\mathrm{HI}}^{\mathrm{DW,IGM}}$ sourced by the neutral IGM to the total column density of

[TABLE]

giving rise to the practically observed damping wing signal. By approximating the localized absorber as a Dirac delta peak at position $r_{\mathrm{abs}}$, we can straightforwardly use Eq. (1) to determine $N_{\mathrm{HI}}^{\mathrm{DW,abs}}$ based on its classical column density $N_{\mathrm{HI}}^{\mathrm{abs}}$:

[TABLE]

provided the absorber is located within the integration range for $N_{\mathrm{HI}}^{\mathrm{DW}}$ enclosed by $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$. Any inferred value of $N_{\mathrm{HI}}^{\mathrm{DW}}$ therefore only constitutes an upper limit on the column density $N_{\mathrm{HI}}^{\mathrm{DW,IGM}}$ sourced by the IGM (and thus relevant to reionization), from which the contributions from any potential absorber would have to be subtracted. Vice versa, the column density $N_{\mathrm{HI}}^{\mathrm{DW,abs}}$ of a given absorber places a lower limit on the full $N_{\mathrm{HI}}^{\mathrm{DW}}$. We plot these limits based on the values determined by Davies et al. (2025) assuming oxygen abundances relative to solar of $[\mathrm{O}/\mathrm{H}]=-2.5$ ($-3.5$) as light (dark) blue arrows in the $(N_{\mathrm{HI}}^{\mathrm{DW}},r_{\mathrm{patch}})$ panel of Figure 4. Note that the length of the arrows carries no information and we place these constraints at $r_{\mathrm{patch}}=r_{\mathrm{abs}}$ for each absorber even though the two distances are not necessarily related given that we define $r_{\mathrm{patch}}$ as the first extended neutral patch in front of the quasar. The fact that the constraint from the closest absorber in the $[\mathrm{O}/\mathrm{H}]=-3.5$ case aligns well with the peak of the posterior is in excellent agreement with the fact that a such metal-poor absorber at this distance would solely account for the observed damping wing imprint without any additional IGM contribution. The Lorentzian-weighted HI column density of an $[\mathrm{O}/\mathrm{H}]=-2.5$ absorber, on the other hand, would largely be negligible compared to our inferred $N_{\mathrm{HI}}^{\mathrm{DW}}$, and hence not bias our conclusions about the ionization state of the IGM. Our local parameterization thus also paves the way for a principled, joint inference of the damping wing signal due to the IGM and potential local absorbers which we leave to future work.

4 Summary and Conclusions

We summarize in Figure 5 our constraints on the reionization history obtained by analyzing the JWST/NIRSpec spectra of the two $z\sim 7.5$ quasars J1007+2115 and J1342+0928 in the context of the local damping wing parameterization (green stars) put forward in Kist et al. (2025b, a) and the conventional global one (yellow pentagons). We compare these constraints to the CMB constraint on the reionization optical depth inferred by the Planck Collaboration et al. (2020) whose $68\,\%$- and $95\,\%$-contours are marked as gray swathes, as well as previous IGM damping wing measurements of the two objects studied in this work (Bañados et al., 2018; Davies et al., 2018; Yang et al., 2020; Reiman et al., 2020; Ďurovčíková et al., 2020; Greig et al., 2022). Our constraints are in clear statistical agreement with all these literature constraints, also owing to the comparably large statistical uncertainties that are unavoidable for individual objects due to the stochastic nature of reionization (Kist et al., 2025c, b, a). Note that despite the methodological improvements in our pipeline as compared to previous approaches (Hennawi et al., 2025b; Kist et al., 2025a), our uncertainties are not necessarily smaller than those quoted in the literature. This is explained by the fact that we rigorously folded in all relevant uncertainties due to stochastic reionization, the unknown quasar lifetime, continuum reconstruction, and spectral noise.

We conclude that J1007+2115 prefers a somewhat lower neutral fraction of $\langle x_{\mathrm{HI}}\rangle=0.32_{-0.20}^{+0.22}$, compared to $\langle x_{\mathrm{HI}}\rangle=0.58_{-0.23}^{+0.23}$ based on J1342+0928 (both inferred in the framework of our local parameterization). The directly inferred global $\langle x_{\mathrm{HI}}\rangle$ posterior of the latter object ($\langle x_{\mathrm{HI}}\rangle=0.78_{-0.30}^{+0.12}$) peaks at the high end compared to our local and most literature damping wing constraints (though closer to Planck Collaboration et al. (2020)), but note here the extended axis of degeneracy of the full posterior of this object (c.f. Figure 4). In addition, all damping wing-based $\langle x_{\mathrm{HI}}\rangle$ constraints would be biased high in case an unusually metal-poor absorber was indeed present in the foreground of this object (Davies et al., 2025) which would pull the actual corrected IGM neutral fraction closer to the value inferred from J1007+2115. Our local parameterization provides a natural framework for studying this possibility more carefully in future work, particularly relevant also to JWST observations of IGM damping wings towards galaxies (Mason et al., 2025; Huberty et al., 2025). More conclusive statements about the ionization state of the IGM at $z\sim 7.5$ will be enabled by applying our robust inference approach to the spectra of additional objects identified at these redshifts by the Euclid wide field survey (Euclid Collaboration et al., 2019; Bañados et al., 2025; Euclid Collaboration: Yang et al., 2025).

In addition, we constrained the lifetimes of the two objects of this study and found J1007+2115 to be remarkably young with $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=4.14_{-0.18}^{+0.74}$, whereas J1342+0928 (with $\log_{10}t_{\mathrm{Q}}/\mathrm{yr}=5.64_{-0.43}^{+0.25}$) is closer to the average expected lifetimes of $t_{\mathrm{Q}}\sim 10^{6}\,\mathrm{yr}$ (Khrykin et al., 2021; Morey et al., 2021). Most importantly, we also obtained the first constraints on the local ionization topology in front of these two sources. With a Lorentzian-weighted HI column density of $\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2}=19.70_{-0.86}^{+0.35}$ and a first neutral patch at $r_{\mathrm{patch}}=28.9_{-14.4}^{+54.0}\,\mathrm{cMpc}$, we remained in the prior-dominated regime for the pre-quasar sightline originating from J1007+2115, but we measured $\log_{10}N_{\mathrm{HI}}^{\mathrm{DW}}/\mathrm{cm}^{-2}=20.24_{-0.22}^{+0.25}$ and $r_{\mathrm{patch}}=10.9_{-5.9}^{+5.6}\,\mathrm{cMpc}$ for J1342+0928. While this work was focused on the conversion of these measurements to constraints on the global timing on reionization, constraining $N_{\mathrm{HI}}^{\mathrm{DW}}$ and $r_{\mathrm{patch}}$ for larger statistical samples of objects will also open the doors to direct constraints on the topology of reionization with quasar IGM damping wings (Sharma et al., 2025).

Acknowledgements

We acknowledge helpful conversations with the ENIGMA group at UC Santa Barbara and Leiden University. This work made use of NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), JAX (Bradbury et al., 2018), NumPyro (Bingham et al., 2018; Phan et al., 2019), sklearn (Pedregosa et al., 2011), Astropy (Astropy Collaboration et al., 2013, 2018, 2022), PypeIt (Prochaska et al., 2020), SkyCalc_ipy (Leschinski, 2021), h5py (Collette, 2013), Matplotlib (Hunter, 2007), corner.py (Foreman-Mackey, 2016), and IPython (Pérez & Granger, 2007). TK and JFH acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 885301). JFH acknowledges support from NSF grant No. 2307180. SEIB is supported by the Deutsche Forschungsgemeinschaft (DFG) under Emmy Noether grant number BO 5771/1-1. FW acknowledges support from NSF award AST-2513040.

Data Availability

The derived data generated in this research will be shared on reasonable requests to the corresponding author.