# Neutral Fraction of Hydrogen in the Intergalactic Medium Surrounding High-Redshift Gamma-Ray Burst 210905A

H.M Fausey, S Vejlgaard, A.J van der Horst, K.E Heintz, L Izzo, D.B

Malesani, K Wiersema, J.P.U Fynbo, N.R Tanvir, S.D Vergani, et al.

► To cite this version:

H.M Fausey, S Vejlgaard, A.J van der Horst, K.E Heintz, L Izzo, et al., Neutral Fraction of Hydrogen in the Intergalactic Medium Surrounding High-Redshift Gamma-Ray Burst 210905A. Monthly Notices of the Royal Astronomical Society, 2025, 536 (3), pp.2839-2856. 10.1093/mnras/stae2757. hal-04537829

HAL Id: hal-04537829

https://hal.science/hal-04537829v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Neutral fraction of hydrogen in the intergalactic medium surrounding high-redshift gamma-ray burst 210905A

H. M. Fausey<sup>1,★</sup> S. Vejlgaard,<sup>2,3</sup> A. J. van der Horst,<sup>1</sup> K. E. Heintz<sup>1,2,3</sup> L. Izzo,<sup>4</sup> D. B. Malesani<sup>1,2,3,5</sup>

K. Wiersema,<sup>6,7</sup> J. P. U. Fynbo,<sup>2,3</sup> N. R. Tanvir,<sup>8</sup> S. D. Vergani,<sup>9,10,11</sup> A. Saccardi,<sup>9</sup> A. Rossi<sup>1,2</sup>

S. Campana,<sup>10</sup> S. Covino,<sup>10</sup> V. D'Elia,<sup>13,14</sup> M. De Pasquale,<sup>15</sup> D. Hartmann,<sup>16</sup> P. Jakobsson,<sup>17</sup>

C. Kouveliotou,<sup>1</sup> A. Levan,<sup>5,18</sup> A. Martin-Carrillo,<sup>19</sup> A. Melandri,<sup>14</sup> J. Palmerio,<sup>9,11</sup> G. Pugliese<sup>20</sup>
and R. Salvaterra<sup>21</sup>

Accepted 2024 December 10. Received 2024 December 6; in original form 2024 March 18

## ABSTRACT

The Epoch of Reionization (EoR) is a key period of cosmological history in which the intergalactic medium (IGM) underwent a major phase change from being neutral to almost completely ionized. Gamma-ray bursts (GRBs) are luminous and unique probes of their environments that can be used to study the timeline for the progression of the EoR. Here, we present a detailed analysis of the European Southern Observatory Very Large Telescope X-shooter spectrum of GRB 210905A, which resides at a redshift of $z \sim 6.3$. We focus on estimating the fraction of neutral hydrogen, $x_{\text{H}\text{I}}$, on the line of sight to the host galaxy of GRB 210905A by fitting the shape of the Lyman-$\alpha$ damping wing of the afterglow spectrum. The X-shooter spectrum has a high signal-to-noise ratio, but the complex velocity structure of the host galaxy limits the precision of our conclusions. The statistically preferred model suggests a low neutral fraction with a $3\sigma$ upper limit of $x_{\text{H}\text{I}} \lesssim 0.15$ or $x_{\text{H}\text{I}} \lesssim 0.23$, depending on the absence or presence of an ionized bubble around the GRB host galaxy, indicating that the IGM around the GRB host galaxy is mostly ionized. We discuss complications in current analyses and potential avenues for future studies of the progression of the EoR and its evolution with redshift.

**Key words:** methods: statistical – techniques: spectroscopic – gamma-ray burst: individual: GRB 210905A – early Universe.

## 1 INTRODUCTION

Gamma-ray bursts (GRBs) are extremely bright transient events that can be observed from across the Universe. GRBs have two phases:

the prompt emission, in which internal shocks due to variability in the outflow (Rees & Meszaros 1992) or magnetic reconnection (Thompson 1994; Spruit, Daigne & Drenkhahn 2001; Giannios & Spruit 2007; Lyubarsky 2010; Beniamini & Granot 2016) produce a flash of gamma-rays; and the afterglow, multiwavelength emission due to the relativistic jet interacting with the surrounding medium (Sari, Piran & Narayan 1998).

<footnote n="1"><sup>1</sup>Department of Physics, George Washington University, 725 21st Street NW, Washington, DC 20052, USA</footnote>

<footnote n="2"><sup>2</sup>Cosmic Dawn Center (DAWN), Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="3"><sup>3</sup>Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark</footnote>

<footnote n="4"><sup>4</sup>DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen, Denmark</footnote>

<footnote n="5"><sup>5</sup>Department of Astrophysics/IMAPP, Radboud University, NL-6525 AJ Nijmegen, the Netherlands</footnote>

<footnote n="6"><sup>6</sup>Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hertfordshire AL10 9AB, UK</footnote>

<footnote n="7"><sup>7</sup>Physics Department, Lancaster University, Lancaster LA1 4YB, UK</footnote>

<footnote n="8"><sup>8</sup>School of Physics and Astronomy, University of Leicester, University Road, Leicester LE1 7RH, UK</footnote>

<footnote n="9"><sup>9</sup>GEPI, Observatoire de Paris, Université PSL, CNRS, 5 pace Jules Janssen, F-92190 Meudon, France</footnote>

<footnote n="10"><sup>10</sup>INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807 Merate, LC, Italy</footnote>

<footnote n="11"><sup>11</sup>Institut d'Astrophysique de Paris, UMR 7095, CNRS-SU, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="12"><sup>12</sup>INAF – Osservatorio di Astrofisica e Scienza dello Spazio, Via Piero Gobetti 93/3, I-40129 Bologna, Italy</footnote>

<footnote n="13"><sup>13</sup>Space Science Data Center (SSDC) – Agenzia Spaziale Italiana (ASI), I-00133 Roma, Italy</footnote>

<footnote n="14"><sup>14</sup>INAF – Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monte Porzio Catone, Italy</footnote>

<footnote n="15"><sup>15</sup>MIFT Department, Polo Papardo, University of Messina, Via Ferdinando Stagno d’Alcontres 31, I-98166 Messina, Italy</footnote>

<footnote n="16"><sup>16</sup>Department of Physics & Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634, USA</footnote>

<footnote n="17"><sup>17</sup>Center for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107 Reykjavik, Iceland</footnote>

<footnote n="18"><sup>18</sup>Department of Physics, University of Warwick, Coventry CV4 7AL, UK</footnote>

<footnote n="19"><sup>19</sup>School of Physics and Centre for Space Research, University College Dublin, Belfield, Dublin 4, Ireland</footnote>

<footnote n="20"><sup>20</sup>Astronomical Institute Anton Pannekoek, University of Amsterdam, NL-1090 GE Amsterdam, the Netherlands</footnote>

<footnote n="21"><sup>21</sup>INAF – Istituto di Astrofisica Sapiziale e Fisica Cosmica, Via Alfonso Corti 12, I-20133 Milano, Italy</footnote>

<footnote>*E-mail: hfausey@gwu.edu</footnote>
GRBs can be divided into two classes: short-hard and long-soft GRBs (Mazets et al. 1981; Kouveliotou et al. 1993). Short-hard GRBs are thought to arise from compact object mergers (Eichler et al. 1989; Narayan, Paczynski & Piran 1992). They have a prompt emission that usually lasts less than 2 s, and tend to have harder spectra (Kouveliotou et al. 1993). Long-soft GRBs are usually associated with the core collapse of Wolf–Rayet stars (Woosley 1993; Galama et al. 1998; Chevalier & Li 1999; Hjorth et al. 2003). Their prompt emission generally lasts longer than 2 s and they tend to have softer spectra (Kouveliotou et al. 1993). However, recently some long-duration GRBs with kilonova counterparts have been associated with compact object mergers (Gao, Lei & Zhu 2022; Rastinejad et al. 2022; Levan et al. 2024), and a short-duration GRB with a soft spectrum and associated supernova (Rossi et al. 2022b) have been observed.

While both classes of GRBs are extremely luminous, long GRBs can be particularly bright, with some having isotropic equivalent luminosities larger than $10^{54}$ erg s⁻¹ (Frederiks et al. 2013; Burns et al. 2023). This allows long GRBs to potentially be detected out to z ~ 20 (Lamb & Reichart 2000), making them cosmological probes of the high-redshift Universe (Campana et al. 2022). The afterglows of GRBs have a simple power-law spectrum (e.g. Sari et al. 1998) as compared to other high-redshift probes such as quasars, Lyman-α (Lyα) emitters, and Lyman break galaxies, which have relatively complex continuum spectra. GRBs are also short-lived compared to quasars, which emit a steady stream of ionizing radiation into the surrounding medium, potentially biasing them towards higher ionization states (Totani et al. 2006). Because GRBs do not emit ionizing radiation for an extended period of time, there is no need to correct for additional intergalactic medium (IGM) ionization from the GRB. For these reasons, GRBs are ideal for studying the chemical evolution of galaxies (Savaglio 2006; Thöne et al. 2013; Sparre et al. 2014; Heintz et al. 2023; Saccardi et al. 2023) and the Epoch of Reionization (EoR; Miralda-Escude 1998; Totani et al. 2006; McQuinn et al. 2008; Hartoog et al. 2015; Lidz et al. 2021), and could potentially contribute to understanding early star formation and the initial mass function (Lloyd-Ronning, Fryer & Ramirez-Ruiz 2002; Fryer et al. 2022) as well as population III stars (Lloyd-Ronning et al. 2002; Campisi et al. 2011).

The EoR is a key era in cosmological history, yet is still poorly understood. It was likely driven by ionizing radiation from the first galaxies, with possible contributions from faint active galactic nuclei (Arons & Wingert 1972; Tegmark, Silk & Blanchard 1994; Haiman, Rees & Loeb 1996; Miralda-Escude 1998; Madau, Haardt & Rees 1999; Faucher-Giguère et al. 2008; Bouwens et al. 2012; Becker & Bolton 2013; Chardin et al. 2015; Giallongo et al. 2015; Madau & Haardt 2015; McQuinn 2016; Matthee et al. 2024). Models and observations suggest that the EoR ended around z ~ 5.5–6 (Totani et al. 2006; Robertson et al. 2015; Ishigaki et al. 2018; Finkelstein et al. 2019; Naidu et al. 2020; Qin et al. 2021), but its progression at higher redshifts is poorly understood. Quasars have been used to track the progression of the EoR by estimating the fraction of neutral hydrogen, i.e. the neutral fraction $x_{H1}$, in the IGM along the line of sight (Mortlock et al. 2011; Simcoe et al. 2012; Bosman & Becker 2015; Greig et al. 2017, 2022; Bañados et al. 2018; Davies et al. 2018; Ďurovčíková et al. 2020; Wang et al. 2020; Yang et al. 2020; Fan, Bañados & Simcoe 2023). There have been fewer than 10 GRBs with spectroscopic or photometric redshifts greater than 6 to date (Tanvir et al. 2009, 2018; Cucchiara et al. 2011; Salvaterra 2015), and a few of these have been used as cosmological probes of the EoR. The latter is done by examining the shape of the red Lyα damping wing, which provides insight into the ionization state

of the IGM surrounding the host galaxy by estimating the fraction of neutral hydrogen (Miralda-Escude 1998; Totani et al. 2006, 2014, 2016; Hartoog et al. 2015).

GRB 210905A was a z ~ 6.3 GRB with one of the most luminous late-time optical afterglows ever observed, which allowed for detailed multiwavelength analyses (Rossi et al. 2022a). Its redshift was determined using metal and fine structure lines associated with the GRB host galaxy (Tanvir et al. 2021; Rossi et al. 2022a; Saccardi et al. 2023). Spectroscopy with the X-shooter instrument on the European Southern Observatory (ESO) Very Large Telescope (VLT) showed two absorption systems around z = 6.3118 and z = 6.3186, as well as additional intervening absorbers at z = 5.739 and z = 2.830 (Saccardi et al. 2023). The two velocity systems at z ≈ 6.3, with a velocity difference of 278 km s⁻¹ and both associated with the host galaxy, were well identified in the metal lines, allowing for an analysis of their metallicity and chemical composition (Saccardi et al. 2023). The metallicity of the z ~ 6.3 system was found to be [M/H]$_{tot}$ = −1.72 ± 0.13 (Saccardi et al. 2023), which is consistent with observations for other GRB damped Ly systems (DLAs) at z ~ 6, and is on the higher end of metallicities for quasar DLAs when extrapolated to high redshift (De Cia et al. 2018; Saccardi et al. 2023).

Determining the neutral fraction of the IGM at different redshifts is important to estimate the progression of the EoR. In this paper, we use the X-shooter spectrum of GRB 210905A to decouple the interstellar medium (ISM) and IGM contributions to the shape of the red Lyα damping wing, and obtain an estimate of the neutral fraction in the IGM surrounding its z ~ 6.3 host galaxy, using the model developed by Miralda-Escude (1998). We perform detailed modelling using the multiple velocity systems found by Saccardi et al. (2023), taking them to be DLAs (systems of gas and dust with a large amount of Lyα absorption due to neutral hydrogen) for this analysis. We also include considerations of the McQuinn et al. (2008) damping wing model to investigate the robustness of the result. We use $H_0 = 67.4$ km s⁻¹ Mpc⁻¹, $\Omega_m = 0.315$, $\Omega_b h^2 = 0.0224$, and $Y_P = 0.2454$ for the cosmological parameters throughout this paper (Planck Collaboration VI 2020). In Section 2, we discuss the observations and data reduction. In Section 3, we present the models and fitting methodology used in this study. In Section 4, we present the results of our damping wing analysis for different assumptions regarding the DLA systems and spectral index. In Section 5, we use robust statistical methods to reach conclusions regarding the model parameters, in particular the neutral fraction. In Section 6, we compare the results for an ionized bubble model. In Sections 7 and 8, we discuss the implications of the results for future high-redshift GRB analyses and summarize our findings.

## 2 OBSERVATIONS AND MODELLING

### 2.1 Data reduction and selection

GRB 210905A was observed with the X-shooter spectrograph (Ver- net et al. 2011), mounted on Unit Telescope 3 (UT3; Melipal), at the VLT observatory in Cerro Paranal (Chile). Data acquisition started 2.53 h after the GRB detection. It consisted of four nod- on-slit exposures of 20 min in the three arms of the instrument: ultraviolet-blue (UVB), visual (VIS), and near-infrared (NIR) with K- band blocking filter, which cover the wavelength range from 3000 to 21 000 Å (Tanvir et al. 2021). The data set was reduced using standard esorex prescriptions (Modigliani et al. 2010), although for the UVB and VIS arms each single exposure was first reduced in STARE mode,
<figure><img src="image_2.png" /><figcaption>**Figure 1.** GRB 210905A spectrum used for the Lyα damping wing analysis (solid black) with error spectrum (dotted red). Regions with metal or telluric lines are shaded in grey and are excluded from the analysis.</figcaption></figure>

after which all exposures were aligned and subsequently stacked to obtain the final 2D spectrum. Post-processing scripts described in Selsing et al. (2019) were applied to correct the final spectrum for slit losses, improve the residual skyline correction, and extract the final 1D spectrum. The VIS arm was further corrected for telluric lines using a PYTHON wrapper for the line-by-line radiative transfer model (Clough, Iacono & Moncet 2014), which provides a best-fitting telluric spectrum obtained from a model of the atmospheric conditions at Paranal and a fit of the telluric lines in the observed data set. This spectrum is the same as those used in the Rossi et al. (2022a) and Saccardi et al. (2023) analyses.

For the damping wing analysis, we chose to use only the VIS data, which range from 5595 to 10 240 Å (Vernet et al. 2011) with 0.2 Å bins, to avoid any issues related to calibration offsets between the X-Shooter VIS and NIR arms. Based on the GRB redshift, the bottom of the red wing of the Lyα forest ends at 8896 Å [i.e. $\lambda_{\alpha}(1+z)$, where $\lambda_{\alpha} = 1215.67$ Å is the Lyα wavelength in a vacuum], and we focus on analysing the spectrum redward of this wavelength. We removed any wavelength bins with telluric absorption greater than 0.9 in either the esorex prescriptions or the ESO Skycalc Model Calculator (Noll et al. 2012; Jones et al. 2013), as well as the metal absorption lines identified in Saccardi et al. (2023). We also excluded data between 9300 and 9530 Å due to the high density of absorption and emission lines in those regions. Finally, for our modelling we only used the data up to 9700 Å, which is the transition between the two largest wavelength orders in the VIS arm, to avoid systematic errors that arise from transitions between these highly curved orders. The final spectrum used for the damping wing analysis is shown in Fig. 1. We note that there is still a deviation in the flux of this GRB compared to the error spectrum, which results in a relatively high chi-squared per degree of freedom in our fitting. For the data reduction, we have followed the ESO-recommended strategies, which are informed by processing many years of data in quality control. We are confident that the data were processed in a way that minimizes any systematic behaviour in error estimates, but

any fit bias due to inaccuracy in the noise model cannot be completely excluded.

# 3 MODEL AND FITTING METHODOLOGY

In this section, we describe the different models and fitting methodologies used to fit the Lyα damping wing.

## 3.1 Model

### 3.1.1 ISM contribution

For the ISM contribution to the damping wing, we consider the two velocity systems identified at $z_{DLA} = 6.3118$ and $z_{DLA} = 6.3186$. Saccardi et al. (2023) take the GRB redshift to be $z = 6.3118$, and conclude that the two velocity systems are likely both from the ISM of the GRB host galaxy. To ensure that we account for any potential absorption from each velocity system, for this analysis we will consider the two main absorption systems identified by Saccardi et al. (2023) as two DLAs. We adopt the model used in Totani et al. (2006) for DLA contributions to the damping wing:

$$ \tau_{DLA}(v_{obs}) = N_{H1}\sigma_{\alpha}[v_{obs}(1+z_{DLA})], \quad (1) $$

where $v_{obs} = c/\lambda_{obs}$, and $\sigma_{\alpha}(v)$ is defined as

$$ \sigma_{\alpha}(v) = \frac{3\lambda_{\alpha}^2 f_{\alpha} \Lambda_{\text{cl},\alpha}}{2\pi} * \frac{\Lambda_{\alpha} \left(\frac{v}{v_{\alpha}}\right)^4}{16\pi^2(v - v_{\alpha})^2 + \Lambda_{\alpha}^2 \left(\frac{v}{v_{\alpha}}\right)^6}. \quad (2) $$

The rate of spontaneous radiative decay between levels in the Lyα transition is given by $\Lambda_{\alpha} = 3 (g_{u}/g_{l})^{-1} f_{\alpha} \Lambda_{\text{cl},\alpha}$, where $g_{u}/g_{l} = 3$ and $f_{\alpha} = 0.4162$ are the ratio of the weights of the upper and lower atomic levels and the absorption oscillator strength for the first Lyα transition, respectively (Weisskopf 1933; Peebles 1993; Madau & Rees 2000; Morton 2003; Bach & Lee 2014), and
$\Lambda_{cl,\alpha} = 8\pi^2 e^2 / 3m_e c \lambda_{\alpha}^2$ is the classical damping constant (Weisskopf 1933; Totani et al. 2006; Bach & Lee 2014).

### 3.1.2 Miralda-Escude (1998) IGM model

One model for the red Ly$\alpha$ damping wing is based on the one developed by Miralda-Escude (1998) and used in previous analyses of other high-redshift GRB damping wings (Totani et al. 2006, 2014, 2016; Chornock et al. 2013; Hartoog et al. 2015). It has also been used to estimate the neutral fraction in other fields (e.g. Bañados et al. 2018; Andika et al. 2023; Umeda et al. 2024). The model assumes a uniform distribution of neutral hydrogen between a fixed upper and lower redshift, $z_{IGM,u}$ and $z_{IGM,l}$. There is assumed to be no neutral hydrogen below $z_{IGM,l}$. The presence of neutral hydrogen increases the optical depth of the medium surrounding the host galaxy, and alters the shape of the red damping wing, which allows the neutral fraction to be estimated. The Miralda-Escude (1998) model does not account for patchiness in the IGM. However, Chen (2024) and Keating et al. (2024) find that the Miralda-Escude (1998) model does sufficiently well at estimating the neutral fraction from the damping wing profile when compared to more complex simulations that account for a patchy Universe using the Cosmic Reionization on Computers (CROC; Gnedin & Kaurov 2014) and Sherwood-Relics (Bolton et al. 2017; Puchwein et al. 2023) simulations, respectively. They argue that, while patchiness can introduce some scatter, the shape of the damping wing is most strongly affected by the average neutral hydrogen density along the line of sight. We also implement the Miralda-Escude (1998) model in multiple shells with independent neutral fractions to examine the neutral fraction results in different redshift ranges (see Section 7.2).

In the model developed by Miralda-Escude (1998) and further refined by Totani et al. (2006), the IGM optical depth as a function of observed wavelength is described as

$$ \tau_{\text{IGM}}(\lambda_{\text{obs}}) = \frac{x_{\text{H1}} \Lambda_{\alpha} \lambda_{\alpha} \tau_{\text{GP}}(z_{\text{host}})}{4\pi^2 c} \left[ \frac{1+z_{\text{obs}}}{1+z_{\text{host}}} \right]^{3/2} \times \left[ I\left(\frac{1+z_{\text{IGM,u}}}{1+z_{\text{obs}}}\right) - I\left(\frac{1+z_{\text{IGM,l}}}{1+z_{\text{obs}}}\right) \right], \quad (3) $$

where $x_{H1}$ is the neutral fraction, $z_{host}$ is the host galaxy redshift, $z_{obs}$ is defined by $1+z_{obs} \equiv \lambda_{obs}/\lambda_{\alpha}$, and $z_{IGM,u}$ and $z_{IGM,l}$ are the upper and lower redshifts between which the distribution of neutral hydrogen is assumed to be constant. Outside of this redshift range the neutral fraction is assumed to be negligible. Totani et al. (2006) model the Gunn-Peterson optical depth using

$$ \tau_{\text{GP}}(z) = \frac{3 f_{\alpha} \Lambda_{\text{cl},\alpha} \lambda_{\alpha}^3 \rho_{\text{crit}} \Omega_{\text{B}} (1 - Y_{\text{p}})}{8 \pi m_{\text{p}} H_0 \Omega_{\text{M}}^{1/2}} (1 + z)^{3/2}, \quad (4) $$

where $\rho_{crit} = 3H_0^2/8\pi G$. The function $I(x)$ is defined as (Miralda-Escude 1998; Totani et al. 2006)

$$ I(x) = \frac{x^{9/2}}{(1-x)} + \frac{9}{7}x^{7/2} + \frac{9}{5}x^{5/2} + 3x^{3/2} + 9x^{1/2} - \frac{9}{2}\ln\left(\frac{1+x^{1/2}}{1-x^{1/2}}\right). \quad (5) $$

We set $z_{IGM,u}$ to the lower DLA redshift $z_{IGM,u} = 6.3118$, since it is the lowest redshift system that is likely to be part of the GRB host galaxy, and $z_{IGM,l} = 6$ since this is the redshift at which the EoR is theorized to have ended (Robertson et al. 2015; Ishigaki et al. 2018; Naidu et al. 2020). Some models suggest that the EoR could have ended at $z \sim 5.5$ (e.g. Finkelstein et al. 2019; Qin et al. 2021; Bosman et al. 2022). However, the choice of $z_l$ does not have an effect on the neutral fraction result so long as it is assumed to be

sufficiently far from the GRB redshift (see Section 7.1 for further discussion on the effect of the $z_l$ parameter).

### 3.1.3 McQuinn et al. (2008) IGM model

To assess the robustness of our results, we also performed fits using the McQuinn et al. (2008) model, which accounts for patchiness in the IGM by assuming a comoving bubble of ionized hydrogen around the host galaxy with radius $R_b$. This model is an approximate version of the Miralda-Escude (1998) model, but allows the radius of an ionized bubble to be fit as a free parameter rather than assumed from the value of $z_{IGM,u}$. The McQuinn et al. (2008) model optical depth due to neutral hydrogen in the IGM is

$$ \tau_{\text{IGM}} \approx 900 \text{ km s}^{-1} \times x_{\text{H1}} \left(\frac{1+z_{\text{DLA}}}{8}\right)^{\frac{3}{2}} \times \left[ \frac{H(z_{\text{DLA}})R_b}{1+z_{\text{DLA}}} - c \left( \frac{v_{\text{obs}}(1+z)-v_{\alpha}}{v_{\alpha}} \right) \right]^{-1}, \quad (6) $$

where $H(z) = H_0 \sqrt{\Omega_M(1+z)^3 + \Omega_{\Lambda} + \Omega_{\text{k}}(1+z)^2}$ (Farooq & Ratra 2013; Verde, Protopapas & Jimenez 2014; Chen, Kumar & Ratra 2017; Wei & Wu 2017). When fitting with the McQuinn et al. (2008) IGM model, we use the same parameters as for the Miralda-Escude (1998) model, but without $z_{IGM,u}$ and $z_{IGM,l}$, and with the additional parameter $R_b$. This model does not account for additional patchiness outside of the ionized region directly around the host galaxy. However, McQuinn et al. (2008) test their model against three radiative transfer simulations and find that the resulting $R_b$ estimates always roughly match the actual $R_b$ value, and the measured $x_{H1}$ along the line of sight is always within 0.3 of the global neutral fraction value.

## 3.2 Fitting

For the fitting, we use four parameters: the normalization $A$ at 10000 Å; the continuum spectrum power-law index $\beta$; the column density $N_{H1}$ for the DLAs; and our main parameter of interest, the neutral fraction $x_{H1}$. We do not fit for extinction since previous analyses found dust extinction to be negligible for GRB 210905A (Saccardi et al. 2023). The column density and neutral fraction are not independent of one another, so we perform a joint fit of the DLA column densities and neutral fraction in an attempt to disentangle their contributions to the damping wing. Because there are two DLAs, we use two column density parameters, one for each DLA (but we also show one-DLA fits for comparison). The column densities of different elements in the two systems are very similar, so we perform some fits of the data with a fixed difference between their H I column densities according to $\log(N_{H1,6.3118}) = \log(N_{H1,6.3186}) + 0.38$, based on the difference of their average Si II column densities, 14.36 and 14.74 at $z = 6.3186$ and $z = 6.3118$, respectively, as found by Saccardi et al. (2023). We also perform modelling with the column densities of the two DLAs uncoupled to see the impact on the results when attempting an independent measurement of the two DLA column densities.

We are only fitting the VIS data to avoid any offsets between the VIS and NIR arms, and we model the data up to 9700 Å, as discussed in the previous section. While this approach results in a more reliable estimate of our key parameters of interest, namely $N_{H1}$ and $x_{H1}$, it also means that there is not enough of a lever arm to confidently constrain the spectral index. We fix the value of the spectral index to results reported in the literature from spectral energy distribution
(SED) fits to multiwavelength data, but include fits with the spectral index as a free parameter as well for completeness.

There are conflicting results for the optical spectral index based on X-ray to NIR analysis. The *Swift* X-Ray Telescope (XRT) spectrum repository reports a late-time X-ray photon index Γ of 1.90 ± 0.15 (with 90 per cent uncertainties; Evans et al. 2009), where $n(E)dE \propto E^{-\Gamma}$. This photon index corresponds to an X-ray spectral index $\beta_X$ of 0.90 ± 0.09 (with $1\sigma$ uncertainties). Multiwavelength analysis reveals that the SED is consistent with slow cooling, with the synchrotron cooling break falling between the X-ray and optical regimes of the spectrum (Rossi et al. 2022a). Assuming that this break is the cooling break, $\beta_{opt} = \beta_X - 1/2$ leads to an optical spectral index of 0.40 ± 0.09 (Sari et al. 1998). In their X-ray to NIR analysis, Rossi et al. (2022a) find a spectral index (with $1\sigma$ uncertainty) of 0.60 ± 0.04. We perform fits fixing the spectral index to both $\beta = 0.40$ and $\beta = 0.60$ to account for both possible spectral indices in our modelling (see Table 1). A similar approach for the spectral index has been taken for other GRB damping wing analyses (e.g. Hartoog et al. 2015). We note that a spectral fit of the X-shooter NIR arm returns a spectral index of ~0.45, but we note that the NIR arm is much noisier than the VIS arm in X-shooter.

The spectral fitting is performed using the EMCEEPYTHON package implementation of the Markov chain Monte Carlo (MCMC) method (Foreman-Mackey et al. 2013). MCMC methods estimate parameter values by constructing posterior distributions by combining a likelihood function with a set of priors, and probing the resulting parameter space using a certain number of walkers (Spade 2020). To fit the $\text{Ly}\alpha$ damping wing using the Miralda-Escude (1998) model, we use 50 walkers with a 1500-step burn-in and 3000-step production chain, which is sufficiently long for the walkers to settle into a stable region of parameter space. We increase the length of the burn-in and production chains (5000 and 10 000 steps, respectively) for the McQuinn et al. (2008) model to adjust for the additional parameters. We use a Bayesian likelihood function $log(\mathcal{L}) = -\chi^2/2$, where $\mathcal{L}$ is the likelihood and $\chi^2$ is the standard definition of chi-squared. All parameters have linearly uniform priors, but we restrict the normalization to $A > 0$, the DLA column density to $18 < \log(N_{H1}/\text{cm}^{-2}) < 23$ (Tanvir et al. 2019), and the neutral fraction to $0 \le x_{H1} \le 1$. For fits in which $\beta$ is a free parameter, we use a linearly uniform prior restricting the spectral index to $0 \le \beta \le 2$. For fits using the McQuinn et al. (2008) model, we use a linearly uniform prior for $R_b$, but restrict the bubble size to $0 \text{ Mpc h}^{-1} \le R_b \le 60 \text{ Mpc h}^{-1}$ (or ~90 Mpc), since the latter corresponds to the Lidz et al. (2021) prediction for ionized bubble size for a largely ionized ($x_{H1} \sim 0.05$) IGM. We focus on fits using this prior, but also include results with $R_b \le 130 \text{ Mpc h}^{-1}$ (or ~193 Mpc, corresponding to $z \sim 6.0$) and $R_b \le 355 \text{ Mpc h}^{-1}$ (or 527 Mpc, corresponding to $z \sim 5.5$) to ensure that our choice of prior is not biasing the results (see Appendix A).

Fitting results using the Miralda-Escude (1998) model are reported in Sections 4 and 5; fitting results using the McQuinn et al. (2008) are reported in Section 6; and fitting results from the Miralda-Escude (1998) shell implementation are discussed in Section 7.2.

# 4 RESULTS

To measure the fraction of neutral hydrogen in the IGM surrounding the GRB host galaxy, we perform a fit of the $\text{Ly}\alpha$ damping wing using spectroscopic data from VLT/X-shooter instrument. As described in the previous section, we perform multiple fits with $\beta$ fixed to two different values, as well as an additional free $\beta$ fit. We consider both coupled and uncoupled column densities for the two velocity systems

in the GRB host galaxy, and explore how the results change when we fit the damping wing with and we include fits for one DLA at the GRB redshift for completeness and comparison. The results for the two-DLA fits are presented in Table 1 and discussed in the following sections. The fits for one DLA are discussed in Section 7.3.

## 4.1 Parameter results

In our analysis, we perform a joint fit of the host ISM column density and neutral fraction parameters. If the posterior distribution for a parameter is most densely populated around a non-zero value, we estimate the best-fitting parameter to be the 50 per cent centile, and use the 16 per cent and 84 per cent centiles for the $1\sigma$ uncertainties. If instead the posterior distribution is most densely populated around 0, we report the $3\sigma$ upper limit using the 99.7 per cent centile.

When fitting the $\text{Ly}\alpha$ damping wing of GRB 210905A with two coupled DLAs, we find that the posterior distribution for $x_{H1}$ is most densely populated around 0, and obtain upper limits for the neutral fraction for all spectral indices. When fixing the spectral index to $\beta = 0.60$, as obtained by the SED analysis, we find $\log(N_{H1,6.3186}/\text{cm}^{-2}) \sim 20.560 \pm 0.002$, corresponding to $\log(N_{H1,6.3118}/\text{cm}^{-2}) \sim 20.940 \pm 0.002$, and a $3\sigma$ upper limit of $x_{H1} < 0.03$ for the neutral fraction. We find similar results when setting the spectral index to $\beta = 0.40$, with $\log(N_{H1,6.3186}/\text{cm}^{-2}) \sim 20.565 \pm 0.003$, corresponding to $\log(N_{H1,6.3118}/\text{cm}^{-2}) \sim 20.945 \pm 0.003$, and a neutral fraction of $x_{H1} < 0.04$. Including the spectral index as a free parameter gives similar results for the two column densities, and a slightly higher upper limit of $x_{H1} < 0.08$ for the neutral fraction as compared to the fixed spectral index fits. Fig. 2 shows the best fit for the coupled DLAs and the spectral index fixed to $\beta = 0.60$, along with the residuals in both flux and $\sigma$ units. Fig. 3 shows the posterior distributions and correlations for all parameters of this particular fit, showing that the neutral fraction does not significantly deviate from 0.

When uncoupling the two DLAs, the column densities of the two DLAs are similar, with $\Delta \log(N_{H1}/\text{cm}^{-2}) \sim 0.1-0.4$. For example, for a fixed spectral index of $\beta = 0.6$, the column densities for the DLAs at $z = 6.3186$ and $z = 6.3118$ are $\log(N_{H1,6.3186}/\text{cm}^{-2}) \sim 20.88 \pm 0.04$ and $\log(N_{H1,6.3118}/\text{cm}^{-2}) \sim 20.48_{-0.18}^{+0.12}$. While the $H1$ column densities are similar, we consistently obtain a higher column density for the DLA at $z = 6.3186$, while the Saccardi et al. (2023) analysis finds that the metal column density for the DLA at $z = 6.3118$ is higher. However, this does not have a significant effect on the results for the neutral fraction.

When uncoupling the DLAs, the posterior distribution for the neutral fraction does still not significantly deviate from 0, but the upper limits increased, ranging from a $3\sigma$ upper limit of $x_{H1} < 0.13$ for a fixed spectral index of 0.60 (see Figs 4 and 5) up to as high as $x_{H1} < 0.28$ when including the spectral index as a free parameter (see Figs 6 and 7). All corner plots were generated using the CORNER package (Foreman-Mackey 2016). All results are summarized in Table 1.

## 4.2 Parameter correlations

For both DLA models, the neutral fraction $3\sigma$ upper limit tends to decrease as the spectral index increases. In the case of the coupled DLAs, the column densities do not seem to depend on spectral index. In the uncoupled DLA case, the column density of the highest redshift DLA increases with spectral index, while the lowest redshift DLA shows the opposite trend. However, the sum of the column densities is stable. The anticorrelation between the column densities of the two uncoupled DLAs can also be seen in Fig. 3. This anticorrelation is
<table><caption>**Table 1.** Modelling results for two DLAs and different spectral indices, with uncertainties and upper limits on the parameters corresponding to $1\sigma$ and $3\sigma$ levels, respectively. The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ at 10 000 Å. The $\chi^2$, reduced $\chi^2$, and degrees of freedom associated with each set of best-fitting parameters are given as well.</caption><thead><tr><th></th><th></th><th>β = 0.40</th><th>β = 0.60</th><th>β free</th></tr></thead><tbody><tr><td rowspan="6">Coupled DLAs</td><td>A</td><td>2.643 ± 0.003</td><td>2.684 ± 0.003</td><td>2.608 ± 0.010</td></tr><tr><td>β</td><td>0.40</td><td>0.60</td><td>0.23 ± 0.05</td></tr><tr><td>log (N<sub>H1,6.3186</sub>/cm<sup>-2</sup>)</td><td>20.565 ± 0.003</td><td>20.560 ± 0.002</td><td>20.568 ± 0.003</td></tr><tr><td>log (N<sub>H1,6.3118</sub>/cm<sup>-2</sup>)</td><td>20.945 ± 0.003</td><td>20.940 ± 0.002</td><td>20.948 ± 0.003</td></tr><tr><td>x<sub>H1</sub></td><td>&lt;0.04</td><td>&lt;0.03</td><td>&lt;0.08</td></tr><tr><td>χ<sup>2</sup></td><td>5596.87</td><td>5640.9</td><td>5587.0</td></tr><tr><td rowspan="3">Uncoupled DLAs</td><td>Red. χ<sup>2</sup></td><td>3.55</td><td>3.58</td><td>3.55</td></tr><tr><td>d.o.f.</td><td>1569</td><td>1569</td><td>1568</td></tr><tr><td>A</td><td>2.634 ± 0.004</td><td>2.670 ± 0.004</td><td>2.621 ± 0.011</td></tr><tr><td rowspan="6">Uncoupled DLAs</td><td>β</td><td>0.40</td><td>0.60</td><td>0.32 ± 0.06</td></tr><tr><td>log (N<sub>H1,6.3186</sub>/cm<sup>-2</sup>)</td><td>20.81 ± 0.04</td><td>20.88 ± 0.04</td><td>20.79 ± 0.06</td></tr><tr><td>log (N<sub>H1,6.3118</sub>/cm<sup>-2</sup>)</td><td>20.67<sup>+0.08</sup><sub>-0.11</sub></td><td>20.48<sup>+0.12</sup><sub>-0.18</sub></td><td>20.71<sup>+0.08</sup><sub>-0.12</sub></td></tr><tr><td>x<sub>H1</sub></td><td>&lt;0.15</td><td>&lt;0.13</td><td>&lt;0.28</td></tr><tr><td>χ<sup>2</sup></td><td>5583.5</td><td>5597.23</td><td>5595.6</td></tr><tr><td>Red. χ<sup>2</sup></td><td>3.54</td><td>3.55</td><td>3.55</td></tr><tr><td></td><td>d.o.f.</td><td>1568</td><td>1568</td><td>1567</td></tr></tbody></table>

expected as the two DLAs have the same effect on the shape of the damping wing.

## 5 MODEL COMPARISON

We performed fits of the X-shooter spectrum using models that differ in spectral index and DLA coupling. As shown in the previous section, we obtain qualitatively similar results for the different model assumptions, but there are quantitative differences for the neutral fraction depending on our model input parameters. Here, we use the PYTHON software package HARMONIC (McEwen et al. 2021) to compare the evidence for the different models using the Bayes factor.

The Bayes factor is used for identifying a preferred model by comparing posterior probabilities (Jeffreys 1935, 1939; Kass & Raftery 1995), and is essentially a comparison of evidence for two models. The PYTHON software package HARMONIC is an implementation of the learnt harmonic mean estimator (Newton & Raftery 1994; McEwen et al. 2021). Newton & Raftery (1994) found that the harmonic mean of the likelihoods of posterior chains could be used to estimate the marginal likelihood, or evidence, for a model. While the original harmonic mean estimator could be unreliable, HARMONIC uses the alternative Gelfand & Dey (1994) harmonic mean estimator, which improved upon the original one by introducing a target distribution (McEwen et al. 2021). HARMONIC uses the posteriors of a fit to estimate the optimal target distribution using machine learning, and to compute the Bayesian evidence and posterior probability for a model (McEwen et al. 2021).

### 5.1 Model selection and preferred model results

The marginal likelihood is an average likelihood for a model based on the probability of generating the data given the prior (McEwen et al. 2021). Using HARMONIC, we find the estimated natural-log of the marginal likelihood, ln(ML), for each model using the posterior chains. The natural-log of the Bayes factor comparing model 1 to

model 0 is defined as $ln(B_{10}) = ln(ML_1) - ln(ML_0)$. The estimated marginal likelihoods for each of the two DLA models can be found in Table 2.

Based on the estimated Bayes factors, we find that the model with the coupled DLA column densities is generally not preferred. When taking the model with coupled column densities to be model 0 and the model with uncoupled column densities to be model 1, $ln(B_{10})$ ranges from ~11.6 to ~28.6. According to Kass & Raftery (1995), $2ln(B_{10}) > 10$ is very strong evidence against model 0, so the model with coupled column densities is strongly disfavoured. A lower spectral index is also generally preferred. When comparing the different uncoupled DLA fits, we find $ln(B_{10}) \sim 6.3$ when comparing the models with their spectral indices fixed to 0.60 and 0.40, and $ln(B_{10}) \sim 2.4$ when comparing the model with the spectral index fixed to 0.40 to the model with a free spectral index, which found a spectral index of $\beta = 0.32 \pm 0.06$. The fit with the spectral index fixed to 0.60 is strongly not preferred, but the models with a fixed spectral index of $\beta = 0.40$ and a free spectral index are harder to distinguish as $2 < 2ln(B_{10}) < 6$ is considered positive, but not strong, evidence against model 0 (Kass & Raftery 1995). Since we are not confident in our ability to constrain the spectral index given the short lever arm, and since there is no strong evidence that the free spectral index fit is preferred over a simpler version of the model where the spectral index is fixed to 0.40, we take the model with $\beta = 0.40$ and uncoupled DLAs to be the preferred model. For this model, the neutral fraction does not significantly deviate from 0 and has $3\sigma$ upper limits of $x_{H1} \lesssim 0.15$. The fit and posteriors for this fit are shown in Figs 8 and 9, respectively.

## 6 COMPARISON TO IONIZED BUBBLE MODEL RESULTS

The McQuinn et al. (2008) model is an approximate version of the Miralda-Escude (1998) model, but it differs from the Miralda-Escude
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example fit for the model with coupled DLA column densities. This fit assumes two DLAs at $z = 6.3118$ and $z = 6.3186$, and a spectral index of $\beta = 0.60$. Regions with metal or telluric lines are shaded in grey and were excluded from analysis. **Top:** VIS spectral data (black) with the 100 final positions of the walkers (blue). **Middle:** Residual plot in erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$. The dashed black line represents 0 flux. **Bottom:** Residual plot in $\sigma$; the dashed and dotted black lines represent 0 and $3\sigma$, respectively.</figcaption></figure>

(1998) and Totani et al. (2006) models in the sense that it allows for an ionized bubble of radius $R_b$ around the host galaxy, rather than relying on the assumed values for the $z_{IGM,u}$ and $z_{IGM,1}$. We performed fits with the McQuinn et al. (2008) model to check the robustness of our result.

In all cases, the McQuinn et al. (2008) model also results in a neutral fraction posterior that does not deviate significantly from 0. An example fit and posterior distribution using the statistically preferred set-up (uncoupled DLAs, a fixed spectral index of $\beta = 0.40$) for the McQuinn et al. (2008) model give a $3\sigma$ upper limit of $x_{H1} \le 0.23$ (see Figs 10 and 11). This is slightly higher than the result obtained with the Miralda-Escude (1998) model with the same set-up ($3\sigma$ upper limit of $x_{H1} \le 0.15$), but consistent. We note that the DLA column density results are also very similar to those found for uncoupled DLAs when using the Miralda-Escude (1998) model. A comparison of the neutral fraction results for the Miralda-Escude

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Corner plot for the fit shown in Fig. 2. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the $1\sigma$ uncertainty ranges (16 per cent and 84 per cent centiles), or the $3\sigma$ upper limit (99.7 per cent centile). The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.</figcaption></figure>

(1998) and McQuinn et al. (2008) models with uncoupled DLAs can be seen in Table 3.

The radius of the comoving ionized bubble is not well constrained, but tends towards the largest allowed radius of $\sim 60$ Mpc $h^{-1}$ or $\sim 90$ Mpc. The behaviour of this parameter could indicate that there is no clear boundary between the ionized bubble around the host galaxy and the IGM, which would indicate that the IGM is highly ionized along the line of sight.

When performing these fits for larger $R_b$ upper bounds, the posteriors display the same behaviour, with an unconstrained $R_b$ tending towards the largest allowed radius. This indicates that the behaviour of the $R_b$ parameter is not caused by our chosen $R_b$ range. When increasing $R_b$, the neutral fraction still does not significantly deviate from 0. The upper limits on $x_{H1}$ do increase with $R_b$, but this is likely because the farther away neutral hydrogen is from the redshift associated with the red Ly$\alpha$ damping wing, the less noticeable its impact on the shape of the damping wing and the less precisely it can be constrained. This same effect can be seen when using a shell implementation of the Miralda-Escude (1998) model (see Section 7.2). Fits using $0 \text{ Mpc } h^{-1} \le R_b \le 130 \text{ Mpc } h^{-1}$ and $0 \text{ Mpc } h^{-1} \le R_b \le 355 \text{ Mpc } h^{-1}$ priors can be found in Appendix A.

For completeness, we also performed a fit with the Miralda-Escude (1998) model with $z_{IGM,u}$ fixed to a lower redshift to create an ionized region around the host galaxy similar to the McQuinn et al. (2008) model. We use $z_{IGM,u} = 6.23$ corresponding to a bubble radius of $\sim 32$ Mpc $h^{-1}$. This is the bubble radius that is expected for a neutral fraction of $\sim 0.25$, which is the estimated value at $z \sim 6.3$ based on recent EoR models (e.g. Ishigaki et al. 2018; Naidu et al. 2020; Lidz et al. 2021). For uncoupled DLAs and a fixed spectral index of $\beta = 0.4$, we find a $3\sigma$ upper limit of $x_{H1} \lesssim 0.26$, which is consistent with the result obtained with the McQuinn et al. (2008) model. Given
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Example fit for a fit with uncoupled column densities for the two DLAs, and a fixed spectral index of $\beta = 0.60$. See Fig. 2 for description of the three panels.</figcaption></figure>

that this result requires a strong assumption about the size of the ionized bubble around the host galaxy, we adopt the 3$\sigma$ upper limit of $x_{H1} \lesssim 0.23$ obtained with the McQuinn et al. (2008) model as the preferred result.

# 7 DISCUSSION

In the previous sections, we have shown that the preferred models point to a negligible neutral fraction of hydrogen with an upper limit of $x_{H1} \lesssim 0.15$ and $x_{H1} \lesssim 0.23$ when using the Miralda-Escude (1998) and McQuinn et al. (2008) models, respectively. Here, we compare our results with those using a shell implementation of the Miralda-Escude (1998) model, explore how the result would be different if we did not obtain a high-quality spectrum (i.e. we would have only one DLA), and discuss the effects of potential spectral curvature and absorption lines on the neutral fraction estimate. We also explore the implications of these results in the broader context of the EoR, and for future high-redshift GRB analyses.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Corner plot for the fit shown in Fig. 4. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the 1$\sigma$ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ Å}^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.</figcaption></figure>

## 7.1 Effect of $z_{IGM_1}$ for the Miralda-Escude (1998) model

The model parameter $z_{IGM_1}$ represents the redshift at which the fraction of neutral hydrogen in the IGM becomes negligible. For redshifts between $z_{IGM_1}$ and $z_{IGM_2}$, the fraction of neutral hydrogen is assumed to be constant. The choice of $z_{IGM_1}$ can have an impact on the modelling results, as a $z_{IGM_1}$ value closer to the GRB redshift gives a higher neutral fraction measurement. However, beyond a certain redshift, $z_{IGM_1}$ becomes sufficiently large that the neutral fraction result is no longer impacted. This is likely because beyond a certain redshift, the presence of neutral hydrogen no longer has a discernible effect on the damping wing.

To test the impact of our choice of $z_{IGM_1}$ we simulated a Ly$\alpha$ damping wing using $z_{IGM_1} = 6.0$, $z_{DLA} = 6.3118$, $z_{Host} = 6.3186$, $\beta = 0.5$, $\log(N_{H1}/\text{cm}^{-2}) = 21$, and $x_{H1} = 0.10$. We fit the simulated data with fixed $z_{IGM_1}$ values ranging from $z = 5.25$ to $z = 6.31$, to assess how the result for the neutral fraction changes with $z_{IGM_1}$ (see Fig. 12). It is evident that our choice of $z_{IGM_1}$ does not impact the retrieved neutral fraction value, and choosing $z_{IGM_1} = 5.5$ to correspond to a late ending reionization would still result in a very similar neutral fraction value. Therefore, our initial choice of $z_{IGM_1}$ did not impact the analysis presented in this paper, and neutral fraction results using the Miralda-Escude (1998) method are independent of the chosen reionization model as long as $z_{IGM_1}$ is sufficiently far from the GRB redshift ($z_{IGM_1} \lesssim 6.2$).

## 7.2 Miralda-Escude (1998) model with shells

To better account for the patchiness of reionization, we implemented the Miralda-Escude (1998) model in a series of shells around the GRB host galaxy. The shell implementation has independent contributions to the shape of the damping wing from shells with widths of $\Delta z = 0.1$ from the GRB redshift out to $z = 5.5$. This allows more flexibility
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Example fit for a fit with uncoupled column densities for the two DLAs, and a free spectral index with a prior of $0 \le \beta \le 2$. See Fig. 2 for description of the three panels.</figcaption></figure>

to explore any evolution in the ionization state of the IGM, but may also be pushing the limits of what can be done with the data.

We perform a fit using the statistically preferred assumptions (uncoupled DLAs and a fixed spectral index of $\beta = 0.4$), but with a 10 000-step burn-in and 20 000-step production chain. We find that the neutral fraction always tends towards 0, but the posterior gradually becomes flatter for shells further away from the GRB redshift. This is expected as neutral hydrogen farther away from the GRB has a less pronounced effect on the $Ly\alpha$ damping wing and becomes harder to constrain. For shells with $z = 6.3186-6.2$, $z = 6.2-6.1$, and $z = 6.1-6.0$, we find $3\sigma$ upper limits of $x_{H1} \lesssim 0.23$, $x_{H1} \lesssim 0.39$, and $x_{H1} \lesssim 0.79$, respectively. For $z = 6.0-5.7$, the neutral fraction is most densely populated around 0, but no meaningful upper limit can be given. For lower redshifts, the posterior still tends slightly towards 0, but is mostly flat. The fit and corresponding corner plot can be seen in Figs B1 and B2, respectively. These results point towards a largely ionized medium around GRB 210905A, as was also found using the original Miralda-Escude (1998) model and the McQuinn et al. (2008) model.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Corner plot for the fit shown in Fig. 6. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the $1\sigma$ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ Å}^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.</figcaption></figure>

<table><caption>**Table 2.** The estimated marginal likelihoods of each two-DLA model using the HARMONICPYTHON software package.</caption><thead><tr><th></th><th>β = 0.4</th><th>β = 0.6</th><th>β free</th></tr></thead><tbody><tr><td>Coupled DLAs</td><td>-2809.8</td><td>-2832.6</td><td>-2804.9</td></tr><tr><td>Uncoupled DLAs</td><td>-2797.7</td><td>-2804.0</td><td>-2795.3</td></tr></tbody></table>

## 7.3 Modelling with one DLA

The modelling with two DLAs shows that the neutral fraction does not significantly deviate from 0. However, if the spectrum had a significantly lower spectral resolution and/or signal-to-noise ratio, the absorption lines from the two DLAs would have been indistinguishable, and we would have fitted the damping wing assuming one DLA. When fitting for only one DLA, the neutral fraction deviates from 0 for all spectral indices, with results ranging from $x_{H1} = 0.08 \pm 0.04$ for a fixed spectral index of 0.60 to $x_{H1} = 0.29 \pm 0.06$ when the spectral index is treated as a free parameter; see Table 4 for the full results. However, it makes sense that the neutral fraction increases when one of the DLAs is removed, since there is now one less component contributing to the damping wing, so the neutral fraction must increase to create the same shape through the wing (see Fig. 13). The posterior distribution for one DLA at $z = 6.3186$ and a spectral index fixed to $\beta = 0.40$ is shown in Fig. 14.

This highlights the importance of using instruments capable of obtaining a decent spectral resolution at wavelengths covering the $Ly\alpha$ damping wing for GRBs at $z \sim 6-10$, such as X-shooter, for GRB follow-up observations. It is clear in the analysis by Saccardi et al. (2023) that there are two velocity systems belonging to the same GRB host galaxy complex, and both are close enough to impact the shape of the $Ly\alpha$ damping wing. Had we not been aware of the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Fit for the statistically preferred model with uncoupled column densities for the two DLAs and a fixed spectral index of $\beta = 0.40$. See Fig. 2 for description of the three panels.</figcaption></figure>

second DLA, we would have derived an inflated neutral fraction in the IGM around the GRB host galaxy. Future instruments such as the Spectrograph and Camera for Observations of Rapid Phenomena in the Infrared and Optical (SCORPIO; Robberto et al. 2020), which includes a long-slit spectrograph, will be even better for GRB follow-up. Since long-slit spectrographs are easier to flux-calibrate, it could eliminate some of the issues that may arise when fitting for $\beta$.

## 7.4 Spectral indices and curvature

In our analysis, we fit the data with multiple fixed spectral index values, because the lack of a sufficient lever arm makes it difficult to fit for the spectral index. In Table 1, we show that in general both the fit quality and the neutral fraction increase with a lower spectral index. Based on the latter, one may prefer the $\beta \sim 0.4$ fits, which is based on extrapolation from the X-ray spectrum with an assumed cooling break. However, Rossi et al. (2022a) find an optical spectral index of $\beta_{opt} = 0.60 \pm 0.04$, and argue that the cooling break occurs below but close to the X-ray energy range. We would like to note that the cooling break does not result in a sharp transition between two

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Corner plot for the fit shown in Fig. 8. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the 1σ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ Å}^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.</figcaption></figure>

spectral indices, but instead is a more gradual change potentially over multiple orders of magnitude in wavelength (Granot & Sari 2002; Van Eerten & Wijers 2009; Uhm & Zhang 2014). Therefore, even if the cooling break is relatively far away in wavelength from either the optical or X-ray regime, it may still have an effect on the derived spectral indices, in which case the neutral fraction estimate would be less certain.

## 7.5 Impact of telluric spectral lines

For GRB 210905A, there is a large number of absorption lines throughout the damping wing, particularly in the upper part of the wing where a higher neutral fraction has the biggest impact. This makes it more difficult to constrain the neutral fraction confidently. For GRBs with a similar redshift, a significant portion of the upper damping wing will be contaminated by telluric emission and absorption lines present around ~9000 Å. Having a larger sample of high-redshift GRBs at a range of different redshifts will alleviate this issue of telluric spectral lines in a crucial part of the damping wing. This is obviously not an issue for space-based observatories such as JWST (Gardner et al. 2006).

## 7.6 Broader context

The neutral fraction around GRB 210905A is generally in agreement with other neutral fraction analyses (e.g. Fan et al. 2006; Totani et al. 2006; Ouchi et al. 2010; Chornock et al. 2013; Hartoog et al. 2015; Ouchi et al. 2018). However, some more recent results are in moderate disagreement. For example, Zhu et al. (2024) stack Ly$\alpha$ transmission profiles based on dark gaps in Ly$\beta$ transmission to find evidence of damping wing-like features. They find a neutral fraction of $x_{H1} \ge 0.061 \pm 0.039$ at $z \sim 5.8$, suggesting that some dark gaps at
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Fit using the McQuinn et al. (2008) model with uncoupled column densities for the two DLAs, a fixed spectral index of $\beta = 0.4$, and an $R_b$ upper limit of 60 Mpc $h^{-1}$ or ~90 Mpc. See Fig. 2 for description of the three panels.</figcaption></figure>

this redshift are caused by neutral islands. A similar analysis by Spina et al. (2024) finds $x_{H1} < 0.44$ at $z \sim 5.9$. The estimate for the fraction of neutral hydrogen around the host galaxy of GRB 210905A was complicated by a range of possible spectral indices and absorption lines in a key portion of the Ly$\alpha$ damping wing. There were also multiple DLAs near the GRB redshift, which can make it more difficult to disentangle the contributions of each component (Heintz et al. 2024). Measuring the neutral fraction around the end of the EoR is also difficult. Lidz et al. (2021) discuss how estimating the neutral fraction is easier at higher redshifts where the ionized bubble around the host galaxy is smaller and changes in $x_{H1}$ have a greater impact on the shape of the damping wing. Finally, it is always possible that the host galaxy of GRB 210905A falls in a region of space with an unusually large ionized region along the line of sight, and does not fully represent the global neutral fraction at this redshift. This is commonly a concern for quasars since they are usually found in massive dark matter haloes, which tend to fall in regions of space that are already highly ionized (Alvarez & Abel 2007; Lidz et al. 2007; Davies et al. 2018). The result for the neutral fraction around

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Corner plot for the fit shown in Fig. 10. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the $1\sigma$ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ and the ionized bubble radius is given in Mpc $h^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.</figcaption></figure>

<table><caption>**Table 3.** Table comparing results from the Miralda-Escude (1998) and McQuinn et al. (2008) models. The upper limit when fitting for $\beta$ is much larger for the McQuinn et al. (2008) model, but for a fixed spectral index the upper limits from the two models are consistent.</caption><thead><tr><th></th><th>β = 0.4</th><th>β = 0.6</th><th>β free</th></tr></thead><tbody><tr><td>Miralda-Escude (1998)</td><td>&lt;0.15</td><td>&lt;0.13</td><td>&lt;0.28</td></tr><tr><td>McQuinn et al. (2008)</td><td>&lt;0.23</td><td>&lt;0.12</td><td>&lt;0.57</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Figure showing impact of $z_{IGM_1}$ on neutral fraction. In this case, $z_{IGM_1}$ does not have a large impact on the neutral fraction result as long as $z_{IGM_1} \lesssim 6.2$.</figcaption></figure>
<table><caption>**Table 4.** Modelling results for different spectral indices assuming one DLA at $z = 6.3186$, with parameter uncertainties given at $1\sigma$ levels. The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ at 10 000 Å. The $\chi^2$, reduced $\chi^2$, and degrees of freedom associated with each set of best-fitting parameters are given as well.</caption><thead><tr><th></th><th>β = 0.40</th><th>β = 0.60</th><th>β free</th></tr></thead><tbody><tr><td>A</td><td>2.639 ± 0.005</td><td>2.668 ± 0.005</td><td>2.619 ± 0.012</td></tr><tr><td>β</td><td>0.40</td><td>0.60</td><td>0.26 ± 0.08</td></tr><tr><td>log (N<sub>H1</sub>6.3186 / cm<sup>-2</sup>)</td><td>20.968 ± 0.007</td><td>20.982 ± 0.006</td><td>20.957 ± 0.009</td></tr><tr><td>x<sub>H1</sub></td><td>0.20 ± 0.04</td><td>0.08 ± 0.04</td><td>0.29 ± 0.06</td></tr><tr><td>χ<sup>2</sup></td><td>5594.1</td><td>5609.6</td><td>5590.9</td></tr><tr><td>Red. χ<sup>2</sup></td><td>3.55</td><td>3.56</td><td>3.55</td></tr><tr><td>d.o.f</td><td>1569</td><td>1569</td><td>1568</td></tr></tbody></table>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Fit comparing the contributions of different absorbers for the one DLA and two uncoupled DLA cases. **Top:** Fit for one DLA at $z = 6.3186$ and a fixed spectral index of $\beta = 0.40$ with the original power law (dotted blue), the individual contributions from the IGM (dash-dotted red) and the DLA (dashed purple), and the total absorption (solid green). **Bottom:** Fit for two uncoupled DLAs and a fixed spectral index of $\beta = 0.40$ with the original power law (dotted blue), the individual contributions from the DLA at $z = 6.3118$ (dash-dotted red) and the DLA at $z = 6.3186$ (dashed purple), and the total absorption (solid green).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Corner plot for the fit shown in Fig. 13. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the $1\sigma$ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$.</figcaption></figure>

GRB 210905A is one data point among many that can help improve our understanding of the EoR, but it cannot in itself be used to support or rule out any reionization models. Increasing the sample of high-redshift GRBs and other high-redshift probes will be vital to increasing our understanding of the EoR (see Section 7.8).

## 7.7 Sightline analysis at $z \sim 6.3$

To determine the likelihood of encountering an overionized sightline for GRB 210905A, we perform a basic simulation to determine the number of fully ionized sightlines for an IGM with $x_{H1} \sim 0.25$, as would be expected for a redshift of $\sim 6.3$ based on recent EoR models (e.g. Ishigaki et al. 2018; Naidu et al. 2020; Lidz et al. 2021). When using a shell implementation of the Miralda-Escude (1998) model, we found that regions below $z \lesssim 6$ were too far away to have a discernible effect on the damping wing (see Section 7.2). Since the sightline of GRB 210905A would only have to be clear out to a redshift of $z \sim 6$ to get a non-detection, we run the simulation for a spherical region of radius $R = 125$ Mpc $h^{-1}$, the comoving distance between $z = 6.3$ and $z = 6$.

We assume spherical neutral islands with $x_{H1} = 1$, and $x_{H1} = 0$ everywhere else, as is sometimes assumed for toy models of neutral gas distributions (e.g. Spina et al. 2024). We determine the radius of each neutral island by pulling from a distribution based on the corresponding lengths of Ly$\beta$ gaps from Zhu et al. (2024). These neutral islands are then given a random position that does not overlap with any other neutral regions, or the GRB host galaxy at the centre of the simulation. The neutral islands continue to populate the simulation until their combined volume occupies at least 25 per cent of the simulation space. 10 000 random sightlines from the simulation are examined and used to estimate the percentage of overionized sightlines. This process is re-run 10 times, and the percentage of
overionized sightlines is averaged, to ensure that the result is not
dictated by an unusual distribution of neutral islands.

We first run this simulation with the constraint that the centre
of all neutral islands must be more than 30 Mpc h⁻¹ to account for
the expected ionized zone around the GRB host galaxy at xH1 ~ 0.25
(Lidz et al. 2021). This configuration results in ~11.75 per cent of all
sightlines being overionized. While the 30 Mpc h⁻¹ ionized region
may be more realistic, we also perform the simulation without it to
provide a more conservative estimate of ~8.85 per cent of sightlines
being overionized. If 8.85 per cent (roughly 1 in 12) to 11.75 per cent
(roughly 1 in 9) of sightlines at z ~ 6.3 are overionized, then there is a
reasonable chance that the sightline of GRB 210905A is overionized,
and not representative of global neutral fraction at z ~ 6.3.

Given the simplicity of the simulation, there are multiple assump-
tions that may lead to overestimates or underestimates of the fraction
of ionized sightlines. For one, the assumptions to make the neutral
islands fully spherical, and to make them either fully ionized or fully
neutral, reduce how broadly the neutral hydrogen is distributed in
the simulations, which could potentially lead to an overestimate.
There is also a limited amount of information about neutral island
size distributions at redshifts greater than 6. The neutral island size
distribution used in this simulation is largely based on dark troughs
in Lyα transmission between z ~ 5.5 and 6 (Zhu et al. 2024). Since
lower redshifts are more ionized, we may be underestimating the size
of the neutral islands at z ~ 6.3. However, even a small amount of
neutral hydrogen (xH1 ~ 10⁻⁴ to 10⁻³) can create or contribute to
dark gaps (Zhu et al. 2022; Gaikwad et al. 2023), so the full length
of the dark troughs may not be exclusively caused by highly neutral
regions. If residual neutral hydrogen in highly ionized regions is
also contributing to the lengths of the dark troughs from Zhu et al.
(2024), then we may instead be overestimating the radii of the neutral
islands. The simulation also takes sightlines with any interaction
with a neutral island, regardless of how short the distance, to be non-
ionized sightlines. Short segments of neutral hydrogen, especially
ones farther away from the GRB redshift, may not have a strong
enough impact on the damping wing to be detectable, and could still
appear to be an overly ionized sightline in the context of a damping
wing analysis, even though it is not fully ionized. In this way, the
simulation may also be underestimating the fraction of overionized
sightlines.

In comparison, the mean free path λ_mfp of photons through the
IGM has been estimated using Lyα and Lyβ transmission profiles,
with estimates of λ_mfp = 0.75<sup>+0.65</sup><sub>-0.45</sub> proper Mpc (pMpc) at z = 6 (or
~3.54 Mpc h<sup>-1</sup>; Becker et al. 2021) and λ_mfp = 0.81<sup>+0.73</sup><sub>-0.48</sub> pMpc
at z = 5.93 (or 3.78 Mpc h<sup>-1</sup>; Zhu et al. 2023). While no direct
measurements of the mean free path have been made beyond z ~ 6,
hydrodynamical simulations suggest a mean free path of ~0.2 pMpc
(~1 Mpc h<sup>-1</sup>; Lewis et al. 2022). For any of these mean free paths,
the probability of a photon reaching 125 Mpc h<sup>-1</sup> is virtually 0,
which is not consistent with the findings of our simulation. One
cause of this discrepancy is that the mean free path is affected by
small amounts of neutral hydrogen in highly ionized regions, so
it is not exclusively dependent on interactions with mainly neutral
regions. If the mean free path at z ~ 6 is mainly limited by residual
neutral hydrogen in ionized regions, then many photons travelling
along overionized sightlines would still likely be absorbed, and
the fraction of photons reaching 125 Mpc h<sup>-1</sup> would be much
lower than the fraction of overionized sightlines. Another potential
contributor to the inconsistency are large-scale fluctuations in the
ionizing ultraviolet (UV) background. Both neutral islands and large-
scale fluctuations in the UV background can cause dark gaps in
transmission that would impact the optical depth, but the latter would

have no impact on a damping wing (Zhu et al. 2023, 2024). This could
explain why the probability for an overionized sightline is so much
higher than the probability for a photon to be unattenuated between
z ~ 6 and 6.3. However, this large discrepancy still highlights the
fact that our simulation is idealized, and may be overestimating the
number of clear sightlines at z ~ 6.3.

7.8 Increasing the sample of high-redshift GRBs

There have been a few other GRB damping wing analyses to measure
the fraction of neutral hydrogen in the IGM. In an analysis of
GRB 050904, another z ~ 6.3 GRB, Totani et al. (2006) found a
best-fitting neutral fraction of 0 with upper limits of <0.17 and <0.60
at the 1σ and 2σ confidence levels, respectively. A similar analysis
was performed by Patel et al. (2010) on GRB 080913 at z ~ 6.7,
but they were unable to constrain xH1. Multiple analyses were
performed on the red damping wing of GRB 130606A at z ~ 5.9,
resulting in varying results. An analysis using data from the Multi-
Mirror Telescope (MMT) Blue Channel Spectrograph (Schmidt,
Weymann & Foltz 1989) and the Gemini Multi-Object Spectrograph
(Hook et al. 2004) on Gemini-North resulted in a 2σ upper limit of
xH1 < 0.11 (Chornock et al. 2013), while damping wing analysis of
data from the Faint Object Camera and Spectrograph (Kashikawa
et al. 2002) on Subaru resulted in xH1 ~ 0.1–0.5 depending on
model assumptions (Totani et al. 2014). Two teams using data from
X-Shooter found a 3σ upper limit of xH1 < 0.05 (Hartoog et al.
2015) and a detection of xH1 = 0.061 ± 0.007 (Totani et al. 2016).
Finally, Melandri et al. (2015) examined the red damping wing of
GRB 140515A at a redshift of z ~ 6.33, and found an upper limit
of xH1 ≲ 0.002. Future analyses of GRB Lyα damping wings will
require a careful treatment of the data. It will also be important to
increase the number of high-redshift GRBs with high-quality spectra
to perform more robust measurements of xH1.

Increasing the sample of high-redshift GRBs will impact our
understanding of the EoR, as the more neutral fraction estimates
we can obtain at a given redshift, the more confidence we can
have in our understanding of the ionization state as a function of
redshift. Future missions such as the *Space Variable Objects Monitor*
(Wei et al. 2016; Atteia, Cordier & Wei 2022), as well as proposed
missions such as the *Gamow Explorer* (Gamow; White et al. 2021)
and *Transient High-Energy Sky and Early Universe Surveyor* (Amati
et al. 2021), will be crucial to find a large sample of high-z GRBs.
New instruments and facilities, such as SCORPIO (Robberto et al.
2020) on *Gemini* and a new generation of 30-m telescopes (Neichel
et al. 2018), with spectroscopy capabilities over a wide wavelength
range, will be important for detailed analyses along the lines of what
we present in this paper.

8 CONCLUSIONS

GRBs are excellent probes of the high-redshift Universe. They
are useful for studying the EoR, as they can be used to track its
progression by measuring the fraction of neutral hydrogen in the
IGM at different redshifts. GRB 210905A was a z ~ 6.3 burst that
had one of the most luminous late-time optical afterglows ever
observed (Rossi et al. 2022a), and was used to probe the chemical
composition of its host galaxy and multiple intervening absorbing
systems (Saccardi et al. 2023). We have presented a detailed analysis
of the afterglow spectrum of GRB 210905A to estimate the fraction of
neutral hydrogen in the IGM surrounding its host galaxy. Based on the
preferred model with two uncoupled DLAs, we find a 3σ upper limit
on the neutral fraction of $x_{HI} \lesssim 0.15$ and $x_{HI} \lesssim 0.23$ for the Miralda-
Escude (1998) and McQuinn et al. (2008) models, respectively, which
indicates that the IGM around the GRB host galaxy is highly ionized.
We note that performing fits with a coupling between the DLAs,
based on information from the high-quality spectra, leads to tighter
upper limits of $x_{HI} \lesssim 0.04$. This result is generally in agreement
with other neutral fraction measurements performed with other high-
redshift probes of the IGM. We discussed the complications of this
particular analysis brought about by possible curvature in GRB
afterglow spectra and telluric lines. We also consider avenues for
pushing forward this field of study. More high-redshift GRBs will be
crucial to the continued study of the EoR and its progress at different
points in cosmological history.

ACKNOWLEDGEMENTS

We thank the referee for their timely and constructive feedback.
This study is based on observations collected at the European
Southern Observatory under ESO programme 106.21T6.014. This
work made use of data supplied by the UK Swift Science Data
Centre at the University of Leicester. AR acknowledges sup-
port from Progetti di Ricerca di Interesse Nazionale - Ministero
dell'Istruzione, dell'Università e della Ricerca (PRIN - MIUR) 2017
(grant 20179ZF5KS).

DATA AVAILABILITY

The data presented in this paper are available upon request to the first
author.

REFERENCES

Alvarez M. A., Abel T., 2007, MNRAS, 380, L30

Amati L. et al., 2021, *Exp. Astron.*, **52**, 183

Andika I. T. et al., 2023, ApJ, 943, 150

Arons J., Wingert D. W., 1972, ApJ, 177, 1

Atteia J. L., Cordier B., Wei J., 2022, *Int. J. Mod. Phys. D*, **31**, 2230008

Bach K., Lee H.-W., 2014, MNRAS, 446, 264

Bañados E. et al., 2018, ApJ, 856, L25

Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023

Becker G. D., D'Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853

Beniamini P., Granot J., 2016, MNRAS, 459, 3635

Bolton J. S., Puchwein E., Sijacki D., Haehnelt M. G., Kim T.-S., Meiksin A., Regan J. A., Viel M., 2017, MNRAS, 464, 897

Bosman S. E. I., Becker G. D., 2015, MNRAS, 452, 1105

Bosman S. E. I. et al., 2022, MNRAS, 514, 55

Bouwens R. J. et al., 2012, ApJ, 752, L5

Burns E. et al., 2023, ApJ, 946, L31

Campana S. et al., 2022, *Nat. Astron.*, **6**, 1101

Campisi M. A., Maio U., Salvaterra R., Ciardi B., 2011, MNRAS, 416, 2760

Chardin J., Haehnelt M. G., Aubert D., Puchwein E., 2015, MNRAS, 453, 2943

Chen H., 2024, MNRAS, 528, L33

Chen Y., Kumar S., Ratra B., 2017, ApJ, 835, 86

Chevalier R. A., Li Z.-Y., 1999, ApJ, 520, L29

Chornock R., Berger E., Fox D. B., Lunnan R., Drout M. R., Fong W.-f., Laskar T., Roth K. C., 2013, ApJ, 774, 26

Clough S. E., Iacono M. J., Moncel J.-L., 2014, Astrophysics Source Code Library, record ascl:1405.001

Cucchiara A. et al., 2011, ApJ, 736, 7

Davies F. B. et al., 2018, ApJ, 864, 142

De Cia A., Ledoux C., Petitjean P., Savaglio S., 2018, A&A, 611, A76

Ďurovčíková D., Katz H., Bosman S. E. I., Davies F. B., Devriendt J., Slyz A., 2020, MNRAS, 493, 4256

Eichler D., Livio M., Piran T., Schramm D. N., 1989, *Nature*, **340**, 126

Evans P. A. et al., 2009, MNRAS, 397, 1177

Fan X. et al., 2006, AJ, 132, 117

Fan X., Bañados E., Simcoe R. A., 2023, ARA&A, **61**, 373

Faroog O., Ratra B., 2013, ApJ, 766, L7

Faucher-Giguère C.-A., Lidz A., Hernquist L., Zaldarriaga M., 2008, ApJ, 688, 85

Finkelstein S. L. et al., 2019, ApJ, 879, 36

Foreman-Mackey D., 2016, *J. Open Source Softw.*, **1**, 24

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Frederiks D. D. et al., 2013, ApJ, 779, 151

Fryer C. L., Lien A. Y., Fruchter A., Ghirlanda G., Hartmann D., Salvaterra R., Upton Sanderbeck P. R., Johnson J. L., 2022, ApJ, 929, 111

Gaikwad P. et al., 2023, MNRAS, 525, 4093

Galama T. J. et al., 1998, *Nature*, **395**, 670

Gao H., Lei W.-H., Zhu Z.-P., 2022, ApJ, 934, L12

Gardner J. P. et al., 2006, *Space Sci. Rev.*, **123**, 485

Gelfand A. E., Dey D. K., 1994, *J. R. Stat. Soc. B*, **56**, 501

Giallongo E. et al., 2015, A&A, 578, A83

Giannios D., Spruit H. C., 2007, A&A, 469, 1

Gnedin N. Y., Kaurov A. A., 2014, ApJ, 793, 30

Granot J., Sari R., 2002, ApJ, 568, 820

Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239

Greig B., Mesinger A., Davies F. B., Wang F., Yang J., Hennawi J. F., 2022, MNRAS, 512, 5390

Haiman Z., Rees M. J., Loeb A., 1996, ApJ, 467, 522

Hartoog O. E. et al., 2015, A&A, 580, A139

Heintz K. E. et al., 2023, A&A, 679, A91

Heintz K. E. et al., 2024, *Science*, **384**, 890

Hjorth J. et al., 2003, *Nature*, **423**, 847

Hook I. M., Jørgensen I., Allington-Smith J. R., Davies R. L., Metcalfe N., Murowinski R. G., Crampton D., 2004, PASP, 116, 425

Ishigaki M., Kawamata R., Ouchi M., Oguri M., Shimasaku K., Ono Y., 2018, ApJ, 854, 73

Jeffreys H., 1935, *Proc. Camb. Phil. Soc.*, **31**, 203

Jeffreys H., 1939, *Theory of Probability*. Clarendon Press, Oxford

Jones A., Noll S., Kausch W., Szyszka C., Kimeswenger S., 2013, A&A, 560, A91

Kashikawa N. et al., 2002, PASJ, 54, 819

Kass R. E., Raftery A. E., 1995, *J. Am. Stat. Assoc.*, **90**, 773

Keating L. C., Puchwein E., Bolton J. S., Haehnelt M. G., Kulkarni G., 2024, MNRAS, 531, L34

Kouveliotou C., Meegan C. A., Fishman G. J., Bhat N. P., Briggs M. S., Koshut T. M., Paciesas W. S., Pendleton G. N., 1993, ApJ, 413, L101

Lamb D. Q., Reichart D. E., 2000, ApJ, 536, 1

Levan A. J. et al., 2024, *Nature*, **626**, 737

Lewis J. S. W. et al., 2022, MNRAS, 516, 3389

Lidz A., McQuinn M., Zaldarriaga M., Hernquist L., Dutta S., 2007, ApJ, 670, 39

Lidz A., Chang T.-C., Mas-Ribas L., Sun G., 2021, ApJ, 917, 58

Lloyd-Ronning N. M., Fryer C. L., Ramirez-Ruiz E., 2002, ApJ, 574, 554

Lyubarsky Y., 2010, ApJ, 725, L234

McEwen J. D., Wallis C. G. R., Price M. A., Docherty M. M., 2021, preprint (arXiv:2111.12720)

McQuinn M., 2016, ARA&A, 54, 313

McQuinn M., Lidz A., Zaldarriaga M., Hernquist L., Dutta S., 2008, MNRAS, 388, 1101

Madau P., Haardt F., 2015, ApJ, 813, L8

Madau P., Rees M. J., 2000, ApJ, 542, L69

Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648

Matthee J. et al., 2024, ApJ, 963, 129

Mazets E. P. et al., 1981, Ap&SS, 80, 3

Melandri A. et al., 2015, A&A, 581, A86

Miralda-Escude J., 1998, ApJ, 501, 15

Modigliani A. et al., 2010, in Silva D. R., Peck A. B., Soifer B. T., eds, Proc. SPIE Conf. Ser. Vol. 7737, Observatory Operations: Strategies, Processes, and Systems III. SPIE, Bellingham, p. 773728
Mortlock D. J. et al., 2011, Nature, 474, 616

Morton D. C., 2003, ApJS, 149, 205

Naidu R. P., Tacchella S., Mason C. A., Bose S., Oesch P. A., Conroy C., 2020, ApJ, 892, 109

Narayan R., Paczynski B., Piran T., 1992, ApJ, 395, L83

Neichel B., Mouillet D., Gendron E., Correia C., Sauvage J. F., Fusco T., 2018, in Di Matteo P., Billebaud F., Herpin F., Lagarde N., Marquette J. B., Robin A., Venot O., eds, Overview of the European Extremely Large Telescope and its instrument suite, SF2A-2018: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. preprint (arXiv:1812.06639)

Newton M. A., Raftery A. E., 1994, J. R. Stat. Soc. B, 56, 3

Noll S., Kausch W., Barden M., Jones A. M., Szyszka C., Kimeswenger S., Vinther J., 2012, A&A, 543, A92

Ouchi M. et al., 2010, ApJ, 723, 869

Ouchi M. et al., 2018, PASJ, 70, S13

Patel M., Warren S. J., Mortlock D. J., Fynbo J. P. U., 2010, A&A, 512, L3

Peebles P. J. E., 1993, Principles of Physical Cosmology. Princeton Univ. Press, Princeton, NJ

Planck Collaboration VI, 2020, A&A, 641, A6

Puchwein E. et al., 2023, MNRAS, 519, 6162

Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390

Rastinejad J. C. et al., 2022, Nature, 612, 223

Rees M. J., Meszaros P., 1992, MNRAS, 258, 41

Robberto M. et al., 2020, in Evans C.J. Bryant J.J. Motohara K. Proc. SPIE Conf. Ser. Vol. 11447, Ground-based and Airborne Instrumentation for Astronomy VIII. SPIE, Bellingham, p. 1144774

Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, ApJ, 802, L19

Rossi A. et al., 2022a, A&A, 665, A125

Rossi A. et al., 2022b, ApJ, 932, 1

Saccardi A. et al., 2023, A&A, 671, A84

Salvaterra R., 2015, J. High Energy Astrophys., 7, 35

Sari R., Piran T., Narayan R., 1998, ApJ, 497, L17

Savaglio S., 2006, New J. Phys., 8, 195

Schmidt G. D., Weymann R. J., Foltz C. B., 1989, PASP, 101, 713

Selsing J. et al., 2019, A&A, 623, A92

Simcoe R. A., Sullivan P. W., Cooksey K. L., Kao M. M., Matejek M. S., Burgasser A. J., 2012, Nature, 492, 79

Spade D., 2020, Markov Chain Monte Carlo Methods: Theory and Practice , Elsevier, Amsterdam, Netherlands

Sparre M. et al., 2014, ApJ, 785, 150

Spina B., Bosman S. E. I., Davies F. B., Gaikwad P., Zhu Y., 2024, A&A, 688, L26

Spruit H. C., Daigne F., Drenkhahn G., 2001, A&A, 369, 694

Tanvir N. R. et al., 2009, Nature, 461, 1254

Tanvir N. R. et al., 2018, ApJ, 865, 107

Tanvir N. R. et al., 2019, MNRAS, 483, 5380

Tanvir N., Rossi A., Xu D., Zhu Z. P., Izzo L., Kann D. A., Levan A. J., Stargate Collaboration, 2021, GCN Circ., 30771, 1

Tegmark M., Silk J., Blanchard A., 1994, ApJ, 434, 395

Thompson C., 1994, MNRAS, 270, 480

Thöne C. C. et al., 2013, MNRAS, 428, 3590

Totani T., Kawai N., Kosugi G., Aoki K., Yamada T., Iye M., Ohta K., Hattori T., 2006, PASJ, 58, 485

Totani T. et al., 2014, PASJ, 66, 63

Totani T., Aoki K., Hattori T., Kawai N., 2016, PASJ, 68, 15

Uhm Z. L., Zhang B., 2014, ApJ, 780, 82

Umeda H., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Isobe Y., Zhang Y., 2024, ApJ, 971, 124

Van Eerten H. J., Wijers R. A. M. J., 2009, MNRAS, 394, 2164

Verde L., Protopapas P., Jimenez R., 2014, Phys. Dark Universe, 5, 307

Vernet J. et al., 2011, A&A, 536, A105

Wang F. et al., 2020, ApJ, 896, 23

Wei J.-J., Wu X.-F., 2017, ApJ, 838, 160

Wei J. et al., 2016, preprint (arXiv:1610.06892)

Weisskopf V. F. I., 1933, The Observatory, 56, 291

White N. E. et al., 2021, in Siegmund O. H., ed., Proc. SPIE Conf. Ser. Vol. 11821, UV, X-Ray, and Gamma-Ray Space Instrumentation for Astronomy XXII. SPIE, Bellingham, p. 1182109

Woosley S. E., 1993, ApJ, 405, 273

Yang J. et al., 2020, ApJ, 897, L14

Zhu Y. et al., 2022, ApJ, 932, 76

Zhu Y. et al., 2023, ApJ, 955, 115

Zhu Y. et al., 2024, MNRAS, 533, L49

APPENDIX A: ADDITIONAL IONIZED BUBBLE MODEL FITS

In Section 6, we show that the size of the ionized bubble around the GRB, $R_b$, is not well constrained and tends towards the largest possible bound allowed by the prior. We also performed fits with $R_b$ upper bounds of $R_b < 130$ Mpc (corresponding to $z \sim 6.0$; see Figs A1 and A2) and $R_b < 355$ Mpc (corresponding to $z \sim 5.5$; see Figs A3 and A4), to ensure that this behaviour is not specific to the

<figure><img src="image_17.png" /></figure>

**Figure A1.** Fit using the McQuinn et al. (2008) model with uncoupled column densities for the two DLAs, a fixed spectral index of $\beta = 0.40$, and an $R_b$ upper bound of 130 Mpc $h^{-1}$ or 193 Mpc. See Fig. 2 for description of the three panels.
<figure><img src="image_18.png" /></figure>

**Figure A2.** Corner plot for the fit shown in Fig. A1. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the 1σ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ and the ionized bubble radius is given in Mpc. $N_{H,1}$ corresponds to the DLA at z = 6.3186, and $N_{H,2}$ corresponds to the DLA at z = 6.3118.

chosen prior. We find that the $R_b$ posterior is always unconstrained but tends towards the largest possible $R_b$ regardless of our choice of upper limit. This indicates that the shape of the $R_b$ posterior is not a product of our choice of bounds for the $R_b$ prior.

As $R_b$ increases, the upper limit on $x_{H1}$ also increases, but this is expected because the farther neutral hydrogen becomes from the GRB, the less it impacts the shape of the damping wing and the harder it becomes to constrain its value. This behaviour can also be seen in the $x_{H1}$ versus $R_b$ corner plot panels (see Figs A2 and A4), as there is a wider spread of neutral fraction values at higher values of $R_b$. In all cases, the posterior is still most densely populated around 0.

<figure><img src="image_19.png" /></figure>

**Figure A3.** Fit using the McQuinn et al. (2008) model with uncoupled column densities for the two DLAs, a fixed spectral index of $\beta = 0.40$, and an $R_b$ upper bound of 355 Mpc $h^{-1}$ or 527 Mpc. See Fig. 2 for description of the three panels.
<figure><img src="image_20.png" /></figure>

**Figure A4.** Corner plot for the fit shown in Fig. A3. The solid grey lines represent the 50 per cent centile, while the dashed grey lines show the $1\sigma$ error ranges (16 per cent and 84 per cent centiles). The normalization is given in $10^{-17}$ erg s$^{-1}$ cm$^{-2}$ Å$^{-1}$ and the ionized bubble radius is given in Mpc. $N_{H,1}$ corresponds to the DLA at $z = 6.3186$, and $N_{H,2}$ corresponds to the DLA at $z = 6.3118$.

<figure><img src="image_21.png" /></figure>

**Figure B1.** Fit using the shells implementation of the Miralda-Escude (1998) model with uncoupled column densities for the two DLAs, a fixed spectral index of $\beta = 0.40$, and shells of width $\Delta z = 0.1$ out to $z = 5.5$. See Fig. 2 for description of the three panels.

## APPENDIX B: SHELL IMPLEMENTATION

To better account for clumpiness in the IGM, we also implement the Miralda-Escude (1998) model in shells of width $\Delta z = 0.1$ from the GRB redshift out to $z = 5.5$ (see Section 7.2). The fit and corner plot for this implementation are shown in Figs B1 and B2, respectively.
<figure><img src="image_22.png" /><figcaption>**Figure B2.** Corner plot associated with Fig. B1. The upper limit on $x_{H1}$ gradually increases with redshift. At redshifts sufficiently far from the GRB redshift ($z \sim 5.7–5.5$), the posterior for $x_{H1}$ flattens out, and any impact of these regions on the shape of the Ly$\alpha$ damping wing can no longer be distinguished.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.