# Revisiting the Intergalactic Medium Around GRB 130606A and Constraints on the Epoch of Reionization

H.M Fausey, A.J van der Horst, N.R Tanvir, K Wiersema, J.P.U Fynbo, D

Hartmann, A de Ugarte Postigo

► To cite this version:

H.M Fausey, A.J van der Horst, N.R Tanvir, K Wiersema, J.P.U Fynbo, et al. Revisiting the Intergalactic Medium Around GRB 130606A and Constraints on the Epoch of Reionization. The Astrophysical Journal, 2025, 985 (1), pp.28. 10.3847/1538-4357/adc5fc hal-04891919

HAL Id: hal-04891919

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

Submitted on 10 Aug 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
OPEN ACCESS

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

CrossMark

# Revisiting the Intergalactic Medium around GRB 130606A and Constraints on the Epoch of Reionization

H. M. Fausey<fnref n="1" />, A. J. van der Horst<fnref n="1" />, N. R. Tanvir<fnref n="2" />, K. Wiersema<fnref n="3" />, J. P. U. Fynbo<fnref n="4" />,<fnref n="5" />, D. Hartmann<fnref n="6" />, and A. de Ugarte Postigo<fnref n="7" />,<fnref n="8" />

Received 2024 December 6; revised 2025 March 3; accepted 2025 March 24; published 2025 May 13

## Abstract

Gamma-ray bursts (GRBs) are excellent probes of the high-redshift Universe due to their high luminosities and the relatively simple intrinsic spectra of their afterglows. They can be used to estimate the fraction of neutral hydrogen (i.e., the neutral fraction) in the intergalactic medium (IGM) at different redshifts through the examination of their Ly$\alpha$ damping wing with high-quality optical-to-near-infrared spectra. Neutral fraction estimates can help trace the evolution of the Epoch of Reionization (EoR), a key era of cosmological history in which the IGM underwent a phase change from neutral to ionized. We revisit GRB 130606A, a $z \sim 5.9$ GRB for which multiple analyses, using the same damping-wing model and data from different telescopes, found conflicting neutral fraction results. We identify the source of the discrepant results to be differences in assumptions for key damping-wing model parameters and data range selections. We perform a new analysis implementing multiple GRB damping-wing models and find a $3\sigma$ neutral fraction upper limit ranging from $x_{H\alpha} \lesssim 0.20$ to $x_{H\alpha} \lesssim 0.23$. We present this result in the context of other neutral fraction estimates and EoR models, discuss the impact of relying on individual GRB lines of sight, and highlight the need for more high-redshift GRBs to effectively constrain the progression of the EoR.

Unified Astronomy Thesaurus concepts: Reionization (1383); Intergalactic medium (813); Gamma-ray bursts (629); Infrared spectroscopy (2285)

## 1. Introduction

The Epoch of Reionization (EoR) is a key era of cosmological history in which the hydrogen in the intergalactic medium (IGM) underwent a phase change from neutral to ionized. It was likely driven by ionizing radiation from the first galaxies, with possible contributions from faint active galactic nuclei (AGN; J. Arons & D. W. Wingert 1972; C.-A. Faucher-Giguère et al. 2008; R. J. Bouwens et al. 2012; G. D. Becker & J. S. Bolton 2013; M. McQuinn 2016; J. Matthee et al. 2024).

To date, a lot of work has been done to model the progression of the EoR (e.g., B. E. Robertson et al. 2015; M. Ishigaki et al. 2018; S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; A. Lidz et al. 2021; Y. Qin et al. 2021; S. Bruton et al. 2023). However, there is still much uncertainty surrounding many key aspects of EoR models, such as the fraction of ionizing photons that escape from a galaxy (i.e., the escape fraction), the UV luminosity distribution of galaxies at different redshifts (i.e., the UV luminosity function), and the primary sources of ionizing radiation.

Recent models agree that the EoR likely ended around a redshift of $z \sim 5.5-6$ (B. E. Robertson et al. 2015; M. Ishigaki et al. 2018; S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; A. Lidz et al. 2021; Y. Qin et al. 2021; S. Bruton et al. 2023),

but due to uncertainty around some key model components, the overall progression is still debated (J. B. Muñoz et al. 2024; J. Witstok et al. 2025). Measuring the fraction of neutral hydrogen (i.e., the neutral fraction) in the IGM at different redshifts can help to constrain which models best represent the EoR progression, and better understand the underlying aspects of cosmological history. This can be done using a variety of sources and methods.

Ly$\alpha$ emitters (LAEs), star-forming galaxies with a high amount of Ly$\alpha$ emission, can be used to infer the neutral fraction at different redshift by examining the evolution of their luminosity functions and the Ly$\alpha$ equivalent width (A. K. Inoue et al. 2018; A. Konno et al. 2018; C. A. Mason et al. 2018, 2019; J. Witstok et al. 2025). LAEs are also increasingly clustered at high redshifts due to the clumpiness of reionization, since Ly$\alpha$ emission from LAEs residing in large ionized bubbles is more likely to be transmitted than emission from LAEs in neutral regions (M. Ouchi et al. 2018). The Planck survey uses cosmic microwave background anisotropies to estimate the total integrated optical depth of reionization, and infer midpoint of the EoR (Planck Collaboration et al. 2020). The fraction of dark pixels and length of dark gaps in the Ly$\alpha$ and Ly$\beta$ forests in high-redshift spectra can be used to estimate the fraction of neutral hydrogen in the IGM. This method does not have any dependency on the modeling of the source (Y. Zhu et al. 2022; X. Jin et al. 2023). For other high-redshift probes such as quasars, galaxies, and gamma-ray bursts (GRBs), examining the Ly$\alpha$ damping wing, a spectral feature created by the presence of neutral hydrogen, can provide

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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

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

<footnote n="3">³Centre for Astrophysics Research, University of Hertfordshire, Hatfield, AL10 9AB, UK</footnote>

<footnote n="4">⁴Cosmic Dawn Center (DAWN), Copenhagen, Denmark</footnote>

<footnote n="5">⁵Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

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

<footnote n="7">⁷Observatoire de la Côte d'Azur, Université Côte d'Azur, Artemis Boulevard de l'Observatoire, 06304 Nice, France</footnote>

<footnote n="8">⁸Aix Marseille Univ, CNRS, CNES, LAM Marseille, France</footnote>
estimates for the fraction of neutral hydrogen in the IGM along the line of sight (J. Miralda-Escude 1998; T. Totani et al. 2006; E. Bañados et al. 2018; F. B. Davies et al. 2018; B. Greig et al. 2019; T. Y.-Y. Hsiao et al. 2024; H. Umeda et al. 2024).

GRBs are valuable probes of the high-redshift Universe and the EoR (J. Miralda-Escude 1998; T. Totani et al. 2006; M. McQuinn et al. 2008; O. E. Hartoog et al. 2015; A. Lidz et al. 2021). They have beamed outflows (i.e., jets; R. Sari et al. 1999), and two distinct phases: the prompt emission, where shocks due to internal variability in the outflow (M. J. Rees & P. Meszaros 1992) or magnetic reconnection (C. Thompson 1994; H. C. Spruit et al. 2001; D. Giannios & H. C. Spruit 2007; Y. Lyubarsky 2010; P. Beniamini & J. Granot 2016) produce a burst of gamma rays; and the afterglow, multi-wavelength emission that arises when the front of the jet collides with the surrounding medium (R. Sari et al. 1998). GRBs fall into two categories: short GRBs, which have prompt emission that is usually shorter than 2 s (E. P. Mazets et al. 1981; C. Kouveliotou et al. 1993), tend to have harder spectra (C. Kouveliotou et al. 1993), and are thought to arise from compact object mergers (D. Eichler et al. 1989; R. Narayan et al. 1992; B. P. Abbott et al. 2017); and long GRBs, which have prompt emissions that are usually longer than 2 s (E. P. Mazets et al. 1981; C. Kouveliotou et al. 1993), tend to have softer spectra (C. Kouveliotou et al. 1993), and originate from the core collapse of Wolf-Rayet stars (S. E. Woosley 1993; T. J. Galama et al. 1998; R. A. Chevalier & Z.-Y. Li 1999; J. Hjorth et al. 2003). However, recently some long GRBs associated with compact object mergers (H. Gao et al. 2022; J. C. Rastinejad et al. 2022; A. J. Levan et al. 2024a) and a short GRB associated with a collapsar (A. Rossi et al. 2022) have been observed.

GRB afterglows have a relatively simple power-law spectrum (R. Sari et al. 1998) as compared to other high-redshift probes like quasars, LAEs, and Lyman-break galaxies (LBGs), making GRBs ideal for studying early star formation and the initial mass function (N. M. Lloyd-Ronning et al. 2002; C. L. Fryer et al. 2022), Population III stars (N. M. Lloyd-Ronning et al. 2002; M. A. Campisi et al. 2011), the chemical evolution of the Universe (S. Savaglio 2006; C. C. Thöne et al. 2013; M. Sparre et al. 2014; A. Saccardi et al. 2023), and the EoR (J. Miralda-Escude 1998; R. Barkana & A. Loeb 2004; T. Totani et al. 2006; A. Lidz et al. 2021). Long GRBs are particularly useful as probes of the EoR due to their extreme luminosities, with some long GRBs' isotropic equivalent luminosities exceeding ~$10^{54}$ erg s⁻¹ (D. D. Frederiks et al. 2013; E. Burns et al. 2023), allowing them to be seen out to high redshifts (B. Ciardi & A. Loeb 2000; D. Q. Lamb & D. E. Reichart 2000; L. J. Gou et al. 2004; R. Salvaterra 2015).

GRBs have been detected up to redshifts of $z \sim 8-9$ (R. Salvaterra et al. 2009; N. R. Tanvir et al. 2009; A. Cucchiara et al. 2011; N. R. Tanvir et al. 2018). However, to date only a handful of high-redshift GRBs have optical-to-near-IR spectra of a high enough quality to enable neutral fraction measurements using the Lyα damping wing (see T. Totani et al. 2006; M. Patel et al. 2010; R. Chornock et al. 2013; T. Totani et al. 2014; O. E. Hartoog et al. 2015; A. Melandri et al. 2015; T. Totani et al. 2016; H. M. Fausey et al. 2025). Of this already limited sample, the neutral fraction around one key GRB, GRB 130606A, sparked controversy. Four analyses, across three separate groups and using data from three different telescopes, found a wide range of neutral fraction

results, from upper limits as low as $x_{H_i} < 0.05$ (O. E. Hartoog et al. 2015) to detections as high as $x_{H_i} = 0.47_{-0.07}^{+0.08}$ (T. Totani et al. 2014). For GRB damping-wing analyses to effectively contribute to the study of the EoR, we need to understand how different assumptions impact the neutral fraction estimate. We revisit the GRB 130606A spectrum taken with the X-shooter instrument (J. Vernet et al. 2011) on the Very Large Telescope (VLT) in Chile to reproduce each previous result using the assumptions from the corresponding analysis, which will provide a better understanding of the cause of the controversy around the neutral fraction measurement for this GRB. We also analyze the X-shooter spectrum using more recent models to obtain an improved estimate for the neutral fraction along the line of sight of GRB 130606A.

In Section 2, we describe the data, modeling, and fitting methodology. In Section 3, we reproduce the neutral fraction results from previous analyses. In Section 4, we obtain new estimates for the neutral fraction around GRB 130606A using assumptions based on new insights, and a range of models. In Sections 5 and 6, we discuss the implications of the result, and summarize our findings. We use cosmological parameters $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.315$, $\Omega_b h^2 = 0.0224$, and $Y_P = 0.2454$ throughout this analysis (Planck Collaboration et al. 2020).

## 2. Data, Modeling, and Methodology

### 2.1. Data

We use the VLT X-shooter (J. Vernet et al. 2011) data for GRB 130606A, which provides a higher resolution than the Gemini Multi-Object Spectrograph (GMOS) and Subaru Faint Object Camera and Spectrograph (FOCAS) data (~0.2 Å for the X-shooter VIS spectra as compared to ~1.38 Å for Gemini GMOS and ~0.74 Å for Subaru FOCAS). The VLT X-shooter observation of GRB 130606A started at 03:57:41 UT on 2013 June 7, about 7 hr after the Neil Gehrels Swift Observatory detection (T. N. Ukwatta et al. 2013). We use the same data reduction as the X-shooter damping-wing analysis performed by O. E. Hartoog et al. (2015), which was reduced using the X-shooter pipeline version 2.2.0 (P. Goldoni 2011). The X-shooter VIS data were binned at 0.2 Å px⁻¹, and the spectrum was corrected for telluric absorption using the spectra of the telluric standard star Hip095400 and the Spextool software (W. D. Vacca et al. 2003). For an in-depth description of the X-shooter data reduction, see O. E. Hartoog et al. (2015).

To obtain an accurate fit to the Lyα damping wing it is important to remove any absorption lines and bad regions of data, which would affect estimates of the continuum or damping-wing profile shape. For the reproduction of each result, we use roughly the same line removal of each corresponding analysis. For the R. Chornock et al. (2013) reproduction, we use data between 8404.41 Å (i.e., $\lambda_{\alpha}(1+z)$, where $z = 5.9314$ is the redshift assumed in the R. Chornock et al. 2013 analysis, and $\lambda_{\alpha_e} = 1215.67$ Å is the Lyα wavelength in a vacuum) to 8902 Å. We remove all absorption lines listed in Table 1 of R. Chornock et al. (2013) as well as any regions with a transmission of <90% in the ESO SKYCALC sky model to account for any telluric absorption (S. Noll et al. 2012; A. Jones et al. 2013). For the T. Totani et al. (2014, 2016) reproductions, we use data between 8426 and 8902 Å, and remove the metal lines listed in T. Totani et al. (2014) and roughly any additional regions omitted in Figure 1.
of their paper. For the O. E. Hartoog et al. (2015) reproduction, we use data between 8403.74 Å (i.e., $\lambda_{\alpha}(1+z)$, where $z=5.91285$) and 8462 Å. We remove any metal lines discussed in O. E. Hartoog et al. (2015), as well as any regions with transmission of <90% in the ESO SKYCALC sky model (S. Noll et al. 2012; A. Jones et al. 2013).

For our new analysis, we assume the same redshift and line-removal method as was used for the O. E. Hartoog et al. (2015) reproduction, but with the wavelength data range extended out to 8902 Å. We also remove a spectral feature between 8469 and 8480 Å, which remains unexplained by absorption lines but is observed in the VLT, Gemini, and Subaru spectra, as well as a Gran Telescopio CANARIAS (GTC) spectrum from A. J. Castro-Tirado et al. (2013).

## 2.2. Modeling

We first attempt to reproduce each result using the J. Miralda-Escude (1998) model, which was used for all previous analyses. The J. Miralda-Escude (1998) model assumes a uniform neutral fraction between two fixed redshifts, $z_{IGM,u}$ and $z_{IGM,l}$, and no neutral hydrogen below $z_{IGM,l}$. The presence of neutral hydrogen increases the optical depth and alters the shape of the Ly$\alpha$ damping wing, allowing the neutral fraction to be estimated.

Since GRBs are relatively short-lived, they do not contribute to the ionization of the medium around their host galaxies, so they do not have a large-scale “proximity effect” like some other high-redshift probes like quasars (L. Christensen et al. 2011; S. D. Vergani 2013; O. E. Hartoog et al. 2015). However, there will still be an ionized bubble around the GRB host galaxy due to ionizing radiation from the galaxy itself or the combined ionizing output of multiple nearby galaxies (A. Lidz et al. 2021). To account for an ionized bubble, we also obtain neutral fraction estimates using the M. McQuinn et al. (2008) model, which is an approximation of the J. Miralda-Escude (1998) model but with the addition of an ionized bubble around the GRB host galaxy with radius $R_b$. This allows for the size of the bubble to be fit as a free parameter rather than assumed with a fixed choice of $z_{IGM,u}$. We assume no neutral hydrogen within the bounds of the ionized bubble, aside from the neutral hydrogen present in the host galaxy itself.

Finally, we use a shell implementation of the J. Miralda-Escude (1998) model to better account for patchiness and evolution in the IGM. We use shells of width $\Delta z = 0.1$ each with a different neutral fraction. We apply this model in two ways. For one implementation, we use completely independent neutral fractions in each shell, as was done in H. M. Fausey et al. (2025). For the other implementation, the neutral fraction in the highest-redshift shell acts as a free parameter, and the neutral fraction in the subsequent shells are tied together by an equation that describes the evolution of the neutral fraction as a function of redshift. Given the short range of redshifts being examined, we assume a linear evolution with slope $dx_{H1}/dz$, which is also treated as a free parameter. Since the coupled version of the shell implementation requires the neutral fraction to decrease with decreasing redshift, we perform the fits at a range of different $z_{IGM,u}$, and in some cases allow $z_{IGM,u}$ to vary as a free parameter, so that $x_{H1}$ is not driven lower by the presence of an ionized bubble around the GRB host galaxy. For an in-depth discussion of each model implementation, see H. M. Fausey et al. (2025).

## 2.3. Methodology

For the J. Miralda-Escude (1998) model, we use four free parameters: the normalization at 8730 Å, A; the spectral index of the underlying power-law continuum, $\beta$; the column density, $N_{H1}$; and the neutral fraction, $x_{H1}$. The M. McQuinn et al. (2008) model also uses these parameters, but includes an additional parameter $R_b$ for the radius of the ionized bubble around the GRB host galaxy. The shell implementation of the J. Miralda-Escude (1998) model with independent shells does not include $R_b$ as a free parameter, but includes separate $x_{H1}$ parameters for the neutral fraction in each shell. The dependent-shell implementation only treats the neutral fraction of the highest-redshift shell as a free parameter, but also includes a parameter for the slope of $x_{H1}$ as a function of redshift which relates the neutral fraction in each shell. In some fits, we also allow $z_{IGM,u}$ of the highest shell to act as a free parameter to better account for an ionized bubble around the GRB redshift.

We use the emcee implementation of a Markov Chain Monte Carlo method (D. Foreman-Mackey et al. 2013). For all fits performed using the J. Miralda-Escude (1998) and M. McQuinn et al. (2008) models, we use 50 walkers with a 2000 step burn-in phase and 4000 step production phase, which upon visual inspection is sufficiently long for the walkers to converge to a preferred region of parameter space. For the shell implementations of the J. Miralda-Escude (1998) model, we increase the number of walkers to 100, the burn-in to 10,000 steps, and the production chain to 20,000 steps to adjust for the additional complexity of the model. We use a Bayesian likelihood function $log(L) = -\chi^2/2$, where $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 \ge 0$, the spectral index to $0 \le \beta \le 3$ (where $F_\nu \propto \nu^{-\beta}$), the host neutral hydrogen column density to $18 < \log(N_{H1}/\text{cm}^{-2}) < 23$ (N. R. Tanvir et al. 2019), and the neutral fraction to $0 \le x_{H1} \le 1$. For fits using the M. 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 $\sim 90 \text{ Mpc}$), since the latter corresponds to the A. Lidz et al. (2021) prediction for ionized bubble size for a largely ionized ($x_{H1} \sim 0.05$) IGM. For the independent-shell implementation of the J. Miralda-Escude (1998) model, the neutral fraction of each shell is restricted to $0 \le x_{H1} \le 1$. For the dependent-shell implementation, the neutral fraction of the highest-redshift shell is restricted to $0 \le x_{H1} \le 1$ and the slope is restricted to $0 \le dx_{H1}/dz \le 2$ (i.e., decreasing neutral fraction with decreasing redshift, or increasing neutral fraction with increasing redshift).

For model comparison, we use marginal likelihoods and the Bayes factor using harmonic (J. D. McEwen et al. 2021) in addition to the using $\chi^2$ and reduced $\chi^2$. For additional discussion about marginal likelihoods and the Bayes factor, see R. E. Kass & A. E. Raftery (1995). For a description of its use for model comparison in the context of GRB damping-wing analyses, see H. M. Fausey et al. (2025).

# 3. Reconstructing Previous Results

All previous results (R. Chornock et al. 2013; T. Totani et al. 2014; O. E. Hartoog et al. 2015; T. Totani et al. 2016) used the J. Miralda-Escude (1998) model for approximation of the neutral fraction in the IGM around GRB 130606A (see Section 2.2). However, each result was found using data from
<table><caption>**Table 1**
Comparison of All Assumptions and Results Reported in Each Previous Analysis</caption><thead><tr><th></th><th>R. Chornock et al. (2013)</th><th>T. Totani et al. (2014)</th><th>O. E. Hartoog et al. (2015)</th><th colspan="2">T. Totani et al. (2016)</th></tr></thead><tbody><tr><td>Instrument</td><td>GMOS</td><td>FOCAS</td><td>X-shooter</td><td>FOCAS</td><td>X-shooter</td></tr><tr><td>Range</td><td>...</td><td>8426–8902 Å</td><td>8406–8462 Å</td><td colspan="2">8426–8902 Å</td></tr><tr><td>z<sub>GRB</sub></td><td>5.9134</td><td>5.9131</td><td>5.91285</td><td colspan="2">5.9131</td></tr><tr><td>z<sub>IGM,u</sub></td><td>z<sub>GRB</sub></td><td>z<sub>GRB</sub></td><td>5.83</td><td>z<sub>GRB</sub></td><td>z<sub>GRB</sub></td></tr><tr><td>z<sub>IGM,i</sub></td><td>...</td><td>5.67</td><td>5.8</td><td colspan="2">5.8</td></tr><tr><td>β value</td><td>...</td><td>0.94 ± 0.04</td><td>0.74<sup>+0.09</sup><sub>-0.07</sub></td><td>1.02 (fixed)</td><td>1.02 (fixed)</td></tr><tr><td>log (N<sup>host</sup><sub>H<sub>I</sub></sub>)</td><td>19.93 ± 0.07</td><td>19.719 ± 0.04</td><td>19.801 ± 0.023</td><td>19.91 ± 0.02</td><td>19.768<sup>+0.032</sup><sub>-0.032</sub></td></tr><tr><td>x<sub>H<sub>I</sub></sub></td><td>&lt;0.11(2σ)</td><td>0.086<sup>+0.012</sup><sub>-0.011</sub></td><td>0.47<sup>+0.08</sup><sub>-0.07</sub></td><td>&lt;0.05(3σ)</td><td>0.061<sup>+0.007</sup><sub>-0.007</sub></td></tr></tbody></table>

different telescopes and using various underlying assumptions about the GRB redshift, the values of $z_{IGM,u}$ and $z_{IGM,i}$, and what ranges of the spectrum should be used. We discuss the differences between analyses below, and attempt to reproduce the results. A summary of all assumptions and results is presented in Table 1.

## 3.1. Previous Analyses and Results

R. Chornock et al. (2013) used spectra from the GMOS (I. M. Hook et al. 2004) on Gemini North, and found a GRB redshift of $z_{host} = 5.9134$. They performed a joint fit of the host column density $log(N_{H}/\text{cm}^{-2})$ and the neutral fraction $x_{H_{I}}$, and found that the highest likely neutral fraction values were slightly above zero, and the highest-likelihood log-column density values were slightly below 19.9 (see R. Chornock et al. 2013, their Figure 8). They reported a $2\sigma$ upper limit of $x_{H_{I}} < 0.11$. While not included in their fit to the Ly$\alpha$ damping wing, they found evidence of a potential DLA at $z = 5.806 \pm 0.001$ based on the presence of Si II, O I and C II lines (see R. Chornock et al. 2013, their Figure 2). They also noted a dark trough in Ly$\alpha$ transmission around this redshift (see R. Chornock et al. 2013, their Figure 5).

T. Totani et al. (2014) used the FOCAS (N. Kashikawa et al. 2002) on the Subaru Telescope and found a slightly different GRB redshift of $z = 5.9131$. They fitted the spectral data from 8426 to 8902 Å. This data range provided a lever arm to constrain the spectral index, but only included the top portion of the Ly$\alpha$ damping wing. They chose to omit data blueward of 8426 Å because they are dominated by host-galaxy absorption (T. Totani et al. 2014). The spectral index was treated as a free parameter for all fits. They performed multiple fits to the damping wing, including two fits assuming a host H I component and an IGM contribution to the Ly$\alpha$ damping wing. For both of these fits $z_{IGM,i}$ was set to 5.67, the lower edge of the dark trough in Ly$\alpha$ transmission identified by R. Chornock et al. (2013) and associated with the potential DLA at $z = 5.806$. For $z_{IGM,u}$, one fit used the host redshift, and the other used $z_{IGM,u} = 5.83$, which T. Totani et al. (2014) identified as the most likely value through a chi-squared analysis. They noted that $z = 5.83$ also corresponded to the upper bound of the dark-pixel region identified in R. Chornock et al. (2013). They reported a neutral fraction of $x_{H_{I}} = 0.086^{+0.012}_{-0.011}$ and $x_{H_{I}} = 0.47^{+0.08}_{-0.07}$ for models using $z_{IGM,u} = z_{host}$ and $z_{IGM,u} = 5.83$, respectively. T. Totani et al. (2014) also included fits that assume no IGM contribution to the damping wing. For the purposes of our

reanalysis, we focus on the fits that include the neutral fraction as a free parameter, but additional details on other fits can be found in T. Totani et al. (2014).

O. E. Hartoog et al. (2015) used data from the X-shooter spectrograph (J. Vernet et al. 2011) on the VLT, and found a GRB redshift of $z = 5.91285$. They performed their fit of the damping wing assuming $z_{IGM,u} = z_{host}$ and $z_{IGM,i} = 5.8$. They found that the neutral fraction result was not sensitive to their choice of $z_{IGM,i}$. The authors performed the fit from 8406 to 8462 Å. Given the limited lever arm, they used a fixed spectral index of $\beta = 1.02$, based on a joint fit of data from X-shooter, the Gamma-Ray Burst Optical/Near-Infrared Detector (GROND; P. Afonso et al. 2013), and the Swift X-Ray Telescope (XRT) repository (P. A. Evans et al. 2007, 2009). They also included fits with the spectral index fixed to $\beta = 0.96$ and $\beta = 1.08$ (corresponding to $\pm 3\sigma$ spectral index range) to check the impact of their choice of spectral index. Using these assumptions, O. E. Hartoog et al. (2015) reported a $3\sigma$ upper limit of $x_{H_{I}} < 0.05$.

Finally, T. Totani et al. (2016) revisited both the VLT X-shooter and Subaru FOCAS data sets for GRB 130606A using similar assumptions to O. E. Hartoog et al. (2015). They took $z_{IGM,u} = z_{host} = 5.9131$ and fixed the spectral index to $\beta = 1.02$. They also used the same data range as in T. Totani et al. (2014). Using these assumptions, T. Totani et al. (2016) found a neutral fraction of $x_{H_{I}} = 0.087^{+0.017}_{-0.029}$ and $x_{H_{I}} = 0.061 \pm 0.007$ for the VLT X-shooter and Subaru FOCAS data sets, respectively.

## 3.2. Results Reconstruction

We attempt to reproduce each result using only the VLT X-shooter data, but the same corresponding data ranges and underlying assumptions for each result. This allows us to test whether the different assumptions were the main cause of the conflicting neutral fraction values.

### 3.2.1. R. Chornock et al. (2013) Reconstruction

The values of $z_{IGM,u}$ and $z_{IGM,i}$ used in the R. Chornock et al. (2013) analysis are unspecified, so we attempt the fit with $z_{IGM,u}$ set to the GRB redshift, and $z_{IGM,i}$ set to 5.8 and 5.7. The exact ranges of spectra used for the neutral fraction fit are also unspecified, so we choose to use the X-shooter spectrum out to ~8900 Å, but remove any absorption lines identified in Table 1 of R. Chornock et al. (2013) as well as any atmospheric absorption and emission lines. Using these assumptions, we
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Example fit to the X-shooter spectrum of GRB 130606A using the assumptions from R. Chornock et al. (2013), and $z_{IGM,1} = 5.8$ (left) along with a zoomed-in examination of the damping-wing fit and residuals (right). Regions with metal or telluric lines are shaded gray and excluded from the fit. Top: spectral data (black) with the 100 final positions of each walker (blue). Middle: residual plot in erg s⁻¹ cm⁻² Å⁻¹; the dashed black line represents zero flux. Bottom: residual plot in σ; the dashed and dotted lines represent 0 and 3σ, respectively.</figcaption></figure>

find a column density of $\log(N_H/\text{cm}^{-2}) = 19.84 \pm 0.02$ and a neutral fraction of $x_{H\i} = 0.03 \pm 0.02$ when assuming $z_{IGM,1} = 5.7$, and $\log(N_H/\text{cm}^{-2}) = 19.85 \pm 0.02$ and $x_{H\i} = 0.03 \pm 0.02$ when assuming $z_{IGM,1} = 5.8$ (see Figure 1). Both of these results are consistent with the $x_{H\i} < 0.11~2\sigma$ upper limit found by R. Chornock et al. (2013). The distributions of our posteriors from both of our fits to the X-shooter data are also similar to the distribution found in the R. Chornock et al. (2013) analysis (see Figure 2 and R. Chornock et al. 2013, their Figure 8, for comparison).

### 3.2.2. *T. Totani et al. (2014) Reconstruction*

While T. Totani et al. (2014) perform four different fits on the Subaru FOCAS spectrum, we choose to focus on reproducing the two fits that include an IGM contribution. We first attempt to reproduce the T. Totani et al. (2014) fit with $z_{IGM,1} = z_{\text{Hsost}}$ and $z_{IGM,1} = 5.67$. When we allow the spectral index to vary freely, we find a column density of

$\log(N_H/\text{cm}^{-2}) = 19.47_{-0.21}^{+0.13}$ and $x_{H\i} = 0.23 \pm 0.05$ (as compared to $x_{H\i} = 0.08_{-0.011}^{+0.012}$). While this value is higher than the one found by T. Totani et al. (2014), we agree that these assumptions lead to a positive neutral fraction detection, and the two values are within $3\sigma$ of each other. However, the spectral index is unusually low at $\beta = 0.23 \pm 0.12$ (see Figure 3, left), which is also inconsistent with the T. Totani et al. (2014) fit. We note that the spectral index and host column density appear to have a slight positive correlation. We would normally expect a slight negative correlation between the two parameters since a higher (shallower) spectral slope requires a lower $\log(N_H)$ to create the same damping-wing shape. However, the column density has a stronger negative correlation with the neutral fraction, which also has a strong negative correlation with the spectral index, making the column density appear to have a slight positive correlation with the spectral index in some cases. This is especially true for the T. Totani et al. (2014) fits, where large portions of the damping wing are omitted from the analysis. The bottom of the Ly$\alpha$ damping
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Posterior distribution corresponding to a fit using the assumptions from R. Chornock et al. (2013), with $z_{\text{IGM},1} = 5.8$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Posterkraft for the T. Totani et al. (2014) results reconstruction. Left: results reconstruction for $z_{\text{IGM,u}} = z_{\text{host}}$, and a free spectral index. Right: results reconstruction for $z_{\text{IGM,u}} = 5.83$, and a Gaussian spectral index prior.</figcaption></figure>

wing is a key region for constraining the host column density, and without it the column density posteriors are more affected by correlations with other parameters.

If we implement a Gaussian prior following the T. Totani et al. (2014) spectral index for this fit ($\beta = 0.94 \pm 0.04$), with the likelihood fixed to zero outside of the $3\sigma$ range, the neutral
fraction no longer significantly deviates from zero, with a peak at $x_{H_i} \sim 0.01$ and a $3\sigma$ upper limit of $x_{H_i} \sim 0.08$. This value is also not consistent with the original T. Totani et al. (2014) value of $0.086^{+0.012}_{-0.011}$.

When we attempt to reproduce the fit using $z_{\text{IGM},u} = 5.83$, we find a higher neutral fraction than the previous fit, with $x_{H_i} = 0.79^{+0.11}_{-0.14}$. It is also significantly higher than the T. Totani et al. (2014) result ($0.47^{+0.08}_{-0.07}$). This fit also results in a most likely spectral index of zero with a $3\sigma$ upper limit of $\beta < 0.47$, which is unusually low (L. Li et al. 2015). When we implement a Gaussian spectral index prior according to the T. Totani et al. (2014) spectral index ($\beta = 0.074 \pm 0.09$) with the likelihood set to zero outside of the $3\sigma$ range, we instead find a neutral fraction of $x_{H_i} = 0.28 \pm 0.13$, which is lower than the T. Totani et al. (2014) result but still within $3\sigma$ (see Figure 3, right). We note that the posterior distribution for the spectral index displays a bimodal distribution, with a small peak around $\beta \sim 0.7$ and a large peak at the boundary of the spectral index prior ($\beta \sim 0.47$). The spectral index in this fit also appears to have an anticorrelation with the neutral fraction, with spectral indices around $\beta \sim 0.7$ resulting in a neutral fraction of $x_{H_i} \sim 0.15$, and spectral indices around $\sim 0.5$ resulting in a spectral index closer to $\sim 0.4$. However, this anticorrelation is expected as a smaller spectral index requires a stronger absorption to create the same-shaped damping wing.

It is important to note that the T. Totani et al. (2014) data range often results in a poor fit of the damping wing (see Figure 4). This can be explained by the fact that the majority of the damping-wing data are not included in the fit. The omission of the damping wing could also be the cause of the volatility of the spectral index and neutral fraction result. Such a limited range of data affected by the neutral fraction can lead to a wide spread of damping-wing profiles (see Figure 4, right), and a small change in spectral index could have a large impact on the neutral fraction result.

### 3.2.3. O. E. Hartoog et al. (2015) Reconstruction

When using the assumptions from O. E. Hartoog et al. (2015), we find a column density of $\log(N_H/\text{cm}^{-2}) = 19.91 \pm 0.01$ and a neutral fraction $3\sigma$ upper limit of $x_{H_i} \lesssim 0.07$ (see Figure 5). Both the column density and neutral fraction are consistent with the original O. E. Hartoog et al. (2015) result.

### 3.2.4. T. Totani et al. (2016) Reconstruction

Following the neutral fraction result from O. E. Hartoog et al. (2015), T. Totani et al. (2016) performed a reanalysis of both the Subaru FOCAS and VLT X-shooter spectra using the assumptions from the O. E. Hartoog et al. (2015) analysis ($z_{\text{IGM},u} = 5.8$, fixed spectral index of $\beta = 1.02$), but with the same data ranges as the T. Totani et al. (2014) analysis (omission of data below 8426 Å). Using these assumptions and the X-shooter data, we find a column density of $\log(N_H/\text{cm}^{-2}) = 19.73 \pm 0.06$ and $x_{H_i} = 0.07 \pm 0.03$ (see Figure 6). This result is similar to the T. Totani et al. (2016) FOCAS and X-shooter neutral fraction results ($x_{H_i} = 0.061 \pm 0.007$ and $x_{H_i} = 0.087^{+0.017}_{-0.029}$, respectively). The column density is also consistent with those from the T. Totani et al. (2016) fits.

## 4. New Analysis with Updated Models and Assumptions

Most results from previous analyses can be reproduced using only the X-shooter data and each paper's assumptions, which suggests that the main source of the discrepancies in results from the different papers stems from the assumptions each paper made in their respective analyses. It is therefore important to carefully examine what assumptions we choose and how they impact the neutral fraction result. Here we present a new analysis using the J. Miralda-Escude (1998) and T. Totani et al. (2006) models using motivated assumptions from the previous analyses. We also explore the neutral fraction result when using more realistic models that better account for the patchiness of the EoR.

### 4.1. J. Miralda-Escude (1998) Methodology

We first perform an analysis using the original J. Miralda-Escude (1998) model. We allow the spectral index to vary freely, and use the GRB redshift estimate from the X-shooter analysis ($z_{\text{GRB}} = 5.91285$; O. E. Hartoog et al. 2015) since it provides the highest-resolution spectrum of GRB 130606A ($\sim 0.2$ Å for the X-shooter VIS spectra as compared to $\sim 1.38$ Å for Gemini GMOS and $\sim 0.74$ Å for Subaru FOCAS). We fit the VLT X-shooter spectrum from 8403.74 to 8902 Å to include both the Ly$\alpha$ damping wing and a long-wavelength lever arm to help constrain the continuum spectral index. When assuming $z_{\text{IGM},u} = z_{\text{GRB}}$ and $z_{\text{IGM},I} = 5.75$, we find a $3\sigma$ upper limit of $x_{H_i} \lesssim 0.04$ with a column density of $N_{H_i} \sim 19.91 \pm 0.01$, which is consistent with the column densities found in the R. Chornock et al. (2013) and O. E. Hartoog et al. (2015) analyses. We find a spectral index of $\beta = 0.63 \pm 0.06$. This is consistent with estimates from an optical-to-near-IR spectral energy distribution using GROND data, which suggest a spectral index of $\beta \sim 0.7$ (P. Afonso et al. 2013). This spectral index is also consistent with results from the Swift-XRT spectrum repository, which reports an X-ray photon index of $\Gamma = 1.71^{+0.11}_{-0.10}$ (90% uncertainties; P. A. Evans et al. 2007, 2009), where $n(E)dE \propto E^{-\Gamma}$. This photon index corresponds to a spectral index with a $1\sigma$ uncertainty of $\beta = 0.71 \pm 0.07$, which means that there is no spectral break between the X-ray and optical regimes, as suggested in O. E. Hartoog et al. (2015).

We also perform a fit with $z_{\text{IGM},u} = 5.8$ and $z_{\text{IGM},I} = 5.65$ to examine the dark trough in Ly$\alpha$ forest emission identified in R. Chornock et al. (2013), as was done in the T. Totani et al. (2014) analysis. In this case, we find a neutral fraction of $x_{H_i} \lesssim 0.53$ with a column density of $N_{H_i} = 19.91 \pm 0.01$. The upper limit on the neutral fraction increases significantly, which may indicate some presence of neutral hydrogen in the system around $z \sim 5.8$. However, this increase is likely partially due to the increased distance between $z_{\text{IGM},u}$ and $z_{\text{host}}$, since neutral hydrogen at a redshift further from the source has a less discernible impact on the Ly$\alpha$ damping wing. We also note that the spectral index estimate for this case is a bit low, with $\beta = 0.57^{+0.08}_{-0.10}$.

If we use the Swift-XRT photon index estimate to implement a Gaussian spectral index prior of $\beta = 0.71 \pm 0.07$ with likelihood set to zero outside of the $3\sigma$ range, we instead find a $3\sigma$ upper limit of $x_{H_i} \lesssim 0.23$, with a column density of $N_{H_i} = 19.91 \pm 0.01$ and spectral index of $\beta = 0.69^{+0.03}_{-0.04}$ (see Figure 7). This neutral fraction upper limit is still higher than the one found using $z_{\text{IGM},u} = z_{\text{host}}$, but provides a tighter constraint than the one found using a uniform spectral index.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left: example fit to the X-shooter spectrum of GRB 130606A using the assumptions from the T. Totani et al. (2014) fit using $z_{\text{IGM,u}} = z_{\text{host}}$. Right: a zoomed-in examination of the damping-wing fit and residuals. See Figure 1 for panel descriptions, and Section 3.2.2 for a discussion of the poor fit to the damping wing in this case.</figcaption></figure>

prior and $z_{\text{IGM,u}} = 5.8$. We also note that the log marginal likelihood for this case is slightly better than the others. The log Bayes factor for comparing model 1 to model 0 is defined as $\ln(B_{10}) = \ln(\text{ML}_1) - \ln(\text{ML}_0)$. According to R. E. Kass & A. E. Raftery (1995), a Bayes factor of $6 < 2\ln(B_{10}) < 10$ is strong evidence in favor of model 1, and $2\ln(B_{10}) > 10$ is very strong evidence in favor of model 1. From the marginal-likelihood values in Table 2, we find that the model using $z_{\text{IGM,u}} = 5.8$ and a Gaussian spectral index prior has strong evidence ($2\ln(B_{10}) = 8.8$) when compared to the model with $z_{\text{IGM,u}} = z_{\text{host}}$ and a uniform $\beta$ prior, and very strong evidence ($2\ln(B_{10}) = 15.0$) when compared to the model with $z_{\text{IGM,u}} = 5.8$ and a uniform $\beta$ prior. This finding is consistent with results from T. Totani et al. (2014), who found a best-fit value of $z_{\text{IGM,u}} = 5.83$ through a comparison of $\chi^2$ values.

## 4.2. M. McQuinn et al. (2008) Methodology

We also analyze the spectrum using the M. McQuinn et al. (2008) model, which is an approximation of the

J. Miralda-Escude (1998) model but includes a parameter for the size of an ionized bubble around the host galaxy, $R_b$, rather than using the assumed $z_{\text{IGM,u}}$ and $z_{\text{IGM,l}}$ values. The M. McQuinn et al. (2008) model assumes that the IGM is fully ionized within the bounds of the ionized bubble. We first perform a fit with $R_b$ as a free parameter with a uniform prior between $0 \le R_b \le 60$ Mpc $h^{-1}$ (or $\sim 90$ Mpc), as is predicted for the ionized bubble size for a largely ionized ($x_{H_i} \sim 0.05$) IGM (A. Lidz et al. 2021). In this case, we find a neutral fraction $3\sigma$ upper limit of $x_{H_i} \lesssim 0.76$, with an unconstrained bubble radius that tends toward $\sim 60$ Mpc $h^{-1}$, indicating a large ionized bubble around the host galaxy. The column density ($\log(N_{H_i}/\text{cm}^{-2}) \sim 19.91$) is consistent with other results, but the spectral index is lower than expected ($\beta = 0.52_{-0.16}^{+0.11}$).

If we again use a Gaussian prior for the spectral index according to the Swift-XRT photon index estimate ($\beta = 0.71 \pm 0.07$, with likelihood fixed to zero outside of the $3\sigma$ range), we instead find a neutral fraction $3\sigma$ upper limit of $x_{H_i} \lesssim 0.20$. The bubble radius is still unconstrained but tends
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: zoomed-in fit of the X-shooter spectrum of GRB 130606A using the assumptions from O. E. Hartoog et al. (2015). See Figure 1 for panel descriptions. Right: associated posteriors for the O. E. Hartoog et al. (2015) reconstruction. See Figure 1 for panel descriptions.</figcaption></figure>

toward ~60 Mpc h⁻¹, and the column density and spectral index results are consistent with those found in previous analyses (see Figure 8). The M. McQuinn et al. (2008) fit with a Gaussian spectral index prior has a slightly lower χ² value, and has strong evidence in its favor when comparing the marginal likelihoods of the Gaussian and uniform spectral index prior results (2 ln(B₁₀) ~ 6.0; see Table 3). The Gaussian spectral index prior was also strongly preferred for the J. Miralda-Escude (1998) model, which had a similar neutral fraction upper limit of xHᵢ ≲ 0.23.

### 4.3. Shell Implementation of the J. Miralda-Escude (1998) Model

To better account for the patchiness of the EoR, we also use a shell implementation of the J. Miralda-Escude (1998) model. We first attempt fits with independent neutral fraction parameters for each shell. We use four shells with widths of Δz ~ 0.1 (or ~7 proper Mpc) starting at the GRB redshift and ending at z ~ 5.5. For each shell we use a uniform neutral fraction prior of 0 ≲ xHᵢ ≲ 1. We find that for all shells the neutral fraction does not deviate significantly from zero, but the upper limit on xHᵢ increases for shells further from the GRB host galaxy (see Figure 10). This behavior was also observed in the analysis of GRB 210905A (H. M. Fausey et al. 2025). This effect is attributed to neutral hydrogen in the IGM having a diminishing impact on the shape of the Lyα damping wing the farther it is from the GRB. We find xHᵢ ≲ 0.03 between z = 5.91285–5.8, xHᵢ ≲ 0.51 for z = 5.8–5.7, and an unconstrained neutral fraction for all other shells.

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

We also perform fits for which the neutral fraction in each shell is coupled with a slope for the neutral fraction as a function of redshift, dxHᵢ/dz. For these fits, the neutral fraction in the closest shell, xHᵢ,₀, and the slope are both treated as free parameters with uniform priors of 0 ≲ xHᵢ,₀ ≲ 1 and 0 ≲ dxHᵢ/dz ≲ 2. The neutral fraction values in the other shells are determined by these two parameters. We assume a total of four shells of width Δz = 0.1.

To account for an ionized bubble, we first allow zIGM,u, the upper redshift boundary of the nearest shell to the GRB redshift, to vary as a free parameter. zIGM,u is given a uniform prior between 5 < zIGM,u < zhost. We find that when zIGM,u is treated as a free parameter, it tends toward lower redshift, with a flat distribution between z ~ 5.0–5.6 that drops off at higher redshifts. This distribution indicates a large ionized bubble around the GRB host galaxy. The flat distribution for z < 5.6 is likely because beyond this redshift neutral hydrogen no longer has any impact on the damping-wing shape, so there is no way to distinguish between the effects of the choice of these redshifts (see Figure 11). For this fit, both the neutral fraction and slope of the neutral fraction as a function of redshift are also unconstrained, with a flat posterior distribution across their allowed ranges. This is likely because zIGM,u tends toward redshifts where the IGM no longer impacts the shape of the damping wing.

We also perform fits with zIGM,u fixed to a range of values between z = 5.85–5.70. All fits still use four shells of width Δz = 0.1, but with the start of the first shell at different redshifts. In all cases, the neutral fraction does not deviate significantly from zero, but the upper limit on xHᵢ in the nearest neutral shell
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Left: zoomed-in fit of the X-shooter spectrum of GRB 130606A using the assumptions from T. Totani et al. (2016). See Figure 1 for panel descriptions. Right: associated posteriors for the T. Totani et al. (2016) reconstruction. See Figure 1 for panel descriptions.</figcaption></figure>

to the GRB increases as $z_{IGM,u}$ decreases. This behavior was also seen for the independent-shell implementation. However, now that the neutral fraction of each shell is coupled according to some slope $dx_{HI}/dz$, within each individual fit the upper limit in farther shells decreases with redshift. For example, for $z_{IGM,u} = 5.8$, we find $x_{HI} \lesssim 0.48$ for $z = 5.8-5.7$, $x_{HI} \lesssim 0.38$ for $z = 5.7-5.6$, $x_{HI} \lesssim 0.34$ for $z = 5.6-5.5$, and $x_{HI} \lesssim 0.31$ for $z = 5.5-5.4$. The neutral fraction estimate in the highest-redshift shell ($z = 5.8-5.7$) is also consistent with the results when using the original J. Miralda-Escude (1998) model with $z_{IGM,u}$. When implementing the same Gaussian spectral index prior from Sections 4.1 and 4.2, the neutral fraction upper limits in each shell also decrease, as they did for the original J. Miralda-Escude (1998) model (see Figure 12), with $x_{HI} \lesssim 0.22$ for $z = 5.8-5.7$, $x_{HI} \lesssim 0.13$ for $z = 5.7-5.6$, $x_{HI} \lesssim 0.10$ for $z = 5.6-5.5$, and $x_{HI} \lesssim 0.09$ for $z = 5.5-5.4$. For the dependent-shell model with $z_{IGM,u} = 5.8$, we find strong evidence in favor of the fit using a Gaussian prior (see Table 4), which gives a $3\sigma$ neutral fraction upper limit of $x_{HI} \lesssim 0.22$, and is consistent with findings of the statistically preferred fits from the J. Miralda-Escude (1998) and M. McQuinn et al. (2008) models (see Sections 4.1 and 4.2).

For all values of $z_{IGM,u}$ the slope $dx_{HI}/dz$ is unconstrained, which can be explained by the fact that the neutral fraction is already nearly zero, so the slope would not have an impact on lower-redshift shells. These results also point to a neutral fraction that does not significantly deviate from zero. For posteriors associated with the independent and dependent shell fits, see Figures 10, 11, and 12 in the Appendix.

# 5. Discussion

In Section 3, we reproduced the previous results from other analyses using only the X-shooter spectrum, which points to the assumptions and data ranges being the source of the discrepant results in each paper, in agreement with the findings from the T. Totani et al. (2016) reanalysis. In Section 4, we performed a new analysis with assumptions based in new information and a range of models. We found that the preferred results for each model all point to a neutral fraction $3\sigma$ upper limit of $x_{HI} \lesssim 0.20 - 0.23$. In this section, we discuss the potential for a system at $z \sim 5.8$, compare the analysis of GRB 130606A to other GRB damping-wing analyses, and explore the implications of the new results in the broader context of EoR measurements and models.

## 5.1. Potential System at $z \sim 5.8$

The R. Chornock et al. (2013) analysis of GRB 130606A identified a potential DLA at $z \sim 5.8$ using metal lines, and noted that it seemed to correspond to a dark trough in $\text{Ly}\alpha$ transmission from $z \sim 5.72$ to 5.79. T. Totani et al. (2014) noted that their best-fit $z_{IGM,u}$ value corresponded to the same redshift as the dark trough in $\text{Ly}\alpha$ transmission. We do not find sufficient evidence for a DLA or a significant neutral fraction at $z \sim 5.7-5.8$ within the damping-wing analysis. However, $z \sim 5.8$ is already $\sim 50$ Mpc $h^{-1}$ from the GRB redshift. The farther away neutral hydrogen is from the GRB redshift, the higher the neutral fraction must be to have a discernible impact.
<figure><img src="image_11.png" /><figcaption>Figure 7. Posteriors for the J. Miralda-Escude (1998) model with $z_{IGM,u}$ = 5.8, and $z_{IGM,I}$ = 5.65 and a Gaussian spectral index prior.</figcaption></figure>

<table><caption>**Table 2**
Comparison of the $\chi^2$, Reduced $\chi^2$, and Log Marginal Likelihoods of Each New Result Using the J. Miralda-Escude (1998) Model</caption><thead><tr><th>z<sub>IGM,u</sub></th><th>β Prior</th><th>β Result</th><th>x<sub>H</sub> Result</th><th>χ<sup>2</sup></th><th>Red. χ<sup>2</sup></th><th>ln(ML)</th></tr></thead><tbody><tr><td>z<sub>host</sub></td><td>Uniform</td><td>0.63 ± 0.06</td><td>&lt;0.04</td><td>2714.4</td><td>1.52</td><td>-1362.1</td></tr><tr><td>5.8</td><td>Uniform</td><td>0.57<sup>+0.08</sup><sub>-0.10</sub></td><td>&lt;0.53</td><td>2723.2</td><td>1.53</td><td>-1365.2</td></tr><tr><td>5.8</td><td>Gaussian</td><td>0.69<sup>+0.03</sup><sub>-0.04</sub></td><td>&lt;0.23</td><td>2714.7</td><td>1.52</td><td>-1357.7</td></tr></tbody></table>

**Note.** The $\chi^2$ and reduced-$\chi^2$ values for $z_{IGM,u} = z_{host}$ with a uniform $\beta$ prior and $z_{IGM,u} = 5.8$ with a Gaussian $\beta$ prior are the same. However, the marginal likelihood for $z_{IGM,u} = 5.8$ and a Gaussian $\beta$ prior is slightly higher than the others.

It is possible that there is some neutral hydrogen around $z \sim 5.8-5.7$, but not in a high enough quantity to be effectively measured with the Ly$\alpha$ damping wing.

## 5.2. Comparison with Other GRB Results

For GRB 130606A, we find a $3\sigma$ neutral fraction upper limit of $x_{H} \lesssim 0.20 - 0.23$. This result is roughly in agreement with current EoR models and neutral fraction measurements (e.g., M. Ishigaki et al. 2018; S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; Y. Zhu et al. 2022; S. Bruton et al. 2023; X. Jin et al. 2023). However, the neutral fraction upper limit for GRB 130606A is higher than or equal to that of GRB 210905A, a $z \sim 6.3$ GRB with a $3\sigma$ neutral fraction upper limit of $x_{H} \lesssim 0.15 - 0.23$ (H. M. Fausey et al. 2025). There are a number of reasons why the neutral fraction upper limit for GRB 130606A may be larger than that of GRB 210905A.

Some previous analyses of GRB 130606A identified a potential DLA and/or a potentially neutral system at $z \sim 5.8$

(R. Chornock et al. 2013; T. Totani et al. 2014; see Section 5.1). While we did not find clear evidence for this potential system in our damping-wing analysis, it is possible that there is neutral hydrogen between $z \sim 5.8-5.7$ in high enough quantities to drive up the neutral fraction upper limit, but not in high enough quantities to allow for a clear detection. GRB 210905A may have also had an overionized sightline for its redshift (H. M. Fausey et al. 2025), which could explain why it has a lower neutral fraction upper limit than GRB 130606A. This discrepancy between GRB damping-wing results highlights the potential impact that the line of sight can have on a GRB neutral fraction estimate.

Finally, assuming a homogeneous neutral fraction can introduce significant scatter in neutral fraction measurements, particularly for low global neutral fractions. For a global neutral fraction $\lesssim 0.25$, a “picket fence” model with thin walls of neutral hydrogen between large ionized bubbles is more realistic than a uniform distribution or a combination of neutral shells (L. C. Keating et al. 2024). An analysis by
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Posterior distributions associated with the M. McQuinn et al. (2008) model, a Gaussian spectral index prior, and an $R_b$ upper limit of 60 Mpc h⁻¹ or ~90 Mpc.</figcaption></figure>

<table><caption>**Table 3**
Comparison of the $\chi^2$, Reduced $\chi^2$, and Log Marginal Likelihoods of Each New Result Using the M. McQuinn et al. (2008) Model</caption><thead><tr><th>β Prior</th><th>β Result</th><th>x<sub>H</sub> Result</th><th>χ<sup>2</sup></th><th>Red. χ<sup>2</sup></th><th>ln(ML)</th></tr></thead><tbody><tr><td>Uniform</td><td>0.52<sup>+0.11</sup><sub>-0.16</sub></td><td>&lt;0.76</td><td>2730.8</td><td>1.53</td><td>-1358.8</td></tr><tr><td>Gaussian</td><td>0.69<sup>+0.03</sup><sub>-0.04</sub></td><td>&lt;0.20</td><td>2716.4</td><td>1.52</td><td>-1355.8</td></tr></tbody></table>

**Note.** The Gaussian prior model is strongly preferred when comparing the marginal likelihoods of the two models.

A. Mesinger & S. R. Furlanetto (2008) found that assuming a uniform distribution can induce a large scatter in neutral fraction measurements for a low global neutral fraction, and can introduce a bias in which the neutral fraction can be overestimated by up to ~0.3 for a neutral fraction of $x_{H_i} \sim 0.5$. More recent analyses found that the average damping-wing profile for a homogeneous neutral fraction is very similar to a patchy reionization model, but that there is also significant overlap in the scatter of damping-wing profiles for $\Delta x_{H_i} \sim 0.2$ (H. Chen 2024; L. C. Keating et al. 2024). This scatter is caused by differences in the locations of the host galaxies in their ionized bubbles, and the distributions of ionized bubbles along the line of sight. The additional scatter from assuming a homogeneous distribution in a patchy IGM could also explain the discrepancy between the

<table><caption>**Table 4**
Comparison of the $\chi^2$, Reduced $\chi^2$, and Log Marginal Likelihoods of Each New Result Using the Dependent-shell Implementation of the J. Miralda-Escande (1998) Model</caption><thead><tr><th rowspan="2">β Prior</th><th rowspan="2">β Result</th><th colspan="4">x<sub>H<sub>i</sub>,z</sub>=5.8-5.7</th></tr><tr><th>Result</th><th>χ<sup>2</sup></th><th>Red. χ<sup>2</sup></th><th>ln(ML)</th></tr></thead><tbody><tr><td>Uniform</td><td>0.58 ± 0.1</td><td>&lt;0.48</td><td>2720.3</td><td>1.52</td><td>-1362.1</td></tr><tr><td>Gaussian</td><td>0.69 ± 0.03</td><td>&lt;0.22</td><td>2715.2</td><td>1.52</td><td>-1357.5</td></tr></tbody></table>

**Note.** The Gaussian prior model is strongly preferred when comparing the marginal likelihoods.

GRB 210905A and GRB 130606A results. It will be vital to increase the number of high-redshift GRBs with high-quality spectroscopic observations so that we can minimize the impact of scatter in $x_{H_i}$ and avoid relying on the sight lines of just a few GRBs to obtain a neutral fraction estimate at different redshifts.

## 5.3. Comparison with Different Methods for Estimating $x_{H_i}$ Evolution

There are still multiple sources of uncertainty in EoR modeling. The escape fraction, $f_{esc}$, denotes the average fraction of ionizing photons that escape from the galaxies in which they are produced, and is important for understanding the evolution of
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Recent EoR models (M. Ishigaki et al. 2018; S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; S. Bruton et al. 2023) and neutral fraction estimates as a function of redshift using a variety of methods. The cyan stars show results from the analysis of the GRB 130606A damping wing (this paper) and the GRB 210905A damping wing (H. M. Fausey et al. 2025). Other neutral fraction results obtained from GRB damping wings are marked with a black star (GRB 050904; T. Totani et al. 2006). Red circles represent neutral fraction results from dark-pixel fractions/troughs (Y. Zhu et al. 2022; X. Jin et al. 2023); orange triangles from Ly$\alpha$ emitter clustering (M. Ouchi et al. 2018); yellow plus signs from LAE luminosity functions (A. K. Inoue et al. 2018; A. Konno et al. 2018; A. M. Morales et al. 2021); green squares from Ly$\alpha$ equivalent widths (C. A. Mason et al. 2018, 2019; A. Hoag et al. 2019; I. Jung et al. 2020; L. R. Whitler et al. 2020; P. Bolan et al. 2022; S. Bruton et al. 2023; T. Morishita et al. 2023); light blue crosses from the Planck survey (Planck Collaboration et al. 2020); dark blue diamonds from galaxy damping wings (T. Y.-Y. Hsiao et al. 2024; H. Umeda et al. 2024); and purple pentagons from quasar damping wings (E. Bañados et al. 2018; F. B. Davies et al. 2018; B. Greig et al. 2019; F. Wang et al. 2020; J. Yang et al. 2020; B. Greig et al. 2022).</figcaption></figure>

the EoR. It has been measured using a variety of sources with redshifts $z \lesssim 4$ (R. E. Mostardi et al. 2015; M. J. Rutkowski et al. 2016; E. Vanzella et al. 2016; C. C. Steidel et al. 2018; N. R. Tanvir et al. 2019; J. B. Vielfaure et al. 2020; Y. I. Izotov et al. 2021; A. J. Pahl et al. 2021), and even up to $z \sim 5$ with a GRB afterglow (A. J. Levan et al. 2024b). However, determining $f_{esc}$ at higher redshifts is increasingly difficult due to an increase in intergalactic attenuation at higher redshifts (P. Madau 1995; A. K. Inoue et al. 2014; B. E. Robertson 2022). While studies have been done to indirectly estimate the escape fraction at higher redshifts (e.g., K. Kakiichi et al. 2018; N. R. Tanvir et al. 2019; R. A. Meyer et al. 2020), more work is required for a complete understanding of the escape fraction at different redshifts (B. E. Robertson 2022). The UV luminosity function is another key component of reionization models, which describes the distribution of galaxy UV luminosities as a function of redshift (N. R. Tanvir et al. 2012). It has changed significantly with the launch of JWST, which detected more high-luminosity galaxies at high redshifts than previously expected (S. L. Finkelstein et al. 2023; Y. Harikane et al. 2023; J. B. Muñoz et al. 2024), and which could have impacted the early progression of the EoR (B. E. Robertson 2022; S. Bruton et al. 2023). There is also still debate as to whether bright or faint galaxies are the primary sources of ionizing radiation (R. P. Naidu et al. 2020; S. Bruton et al. 2023; Z. Wu & A. Kravtsov 2024), and whether or not AGN also played a role (S. L. Finkelstein et al. 2019).

There are a wide range of methods and probes for estimating the neutral fraction at different redshifts. Ly$\alpha$ damping wings of quasars and LBGs can be used to obtain neutral fraction

estimates, with some additional considerations for their more complex continua, and the impact of continuous ionizing radiation from quasars (E. Bañados et al. 2018; F. B. Davies et al. 2018; B. Greig et al. 2019; F. Wang et al. 2020; J. Yang et al. 2020; B. Greig et al. 2022; T. Y.-Y. Hsiao et al. 2024; H. Umeda et al. 2024). LAEs are also useful for neutral fraction estimation. LAEs are clustered in the sky rather than isotropically distributed, and the amount of clustering is expected to increase at higher redshift due to the patchiness of the IGM, since LAEs in large ionized bubbles are less impacted by Ly$\alpha$ absorption (M. Ouchi et al. 2018). Examining the clustering of LAEs as a function of redshift can provide insight into the neutral fraction at different redshifts (M. Ouchi et al. 2018). The evolution of the LAE luminosity function in comparison with the UV luminosity function can help estimate the change in Ly$\alpha$ transmission, which can be related to the neutral fraction and ionized bubble sizes (A. K. Inoue et al. 2018; A. Konno et al. 2018; A. M. Morales et al. 2021). Finally, the evolution of equivalent widths of LAE Ly$\alpha$ emission lines can provide insight into the evolution of the EoR (C. A. Mason et al. 2018). Dark pixels and troughs in Ly$\alpha$ and Ly$\beta$ transmission can also indicate the presence of neutral hydrogen in the IGM at different redshifts (Y. Zhu et al. 2022; X. Jin et al. 2023). Recently, Y. Zhu et al. (2024) stacked Ly$\alpha$ transmission profiles according to gaps in Ly$\beta$ transmission in search of Ly$\alpha$ damping-wing features. The Planck survey also estimated the midpoint of reionization using electron scattering optical depth estimates (Planck Collaboration et al. 2020). A compilation of results from these methods of neutral fraction estimation are presented in Figure 9 along with theoretical
curves for four different models (M. Ishigaki et al. 2018; S. L. Finkelstein et al. 2019; R. P. Naidu et al. 2020; S. Bruton et al. 2023). For EoR theoretical curves, we show only one line rather than the entire model ranges.

There is still a large amount of uncertainty in the progression of the EoR. There are a wide range of neutral fraction estimates at each redshift, making it difficult to resolve the EoR progression. Increasing the number of neutral fraction measurements from a large range of probes will be vital to understanding the EoR and its evolution. Since GRBs fade rapidly, quick spectral follow-up can greatly improve the data quality. However, it can be difficult to quickly determine which GRBs are high redshift, as they often require near-IR imaging for their identification. Proposed missions like the Gamow Explorer (N. E. White et al. 2021) and Transient High-Energy Sky and Early Universe Surveyor (L. Amati et al. 2021) are designed to quickly identify high-redshift GRBs and alert the community, so they can aid in decreasing the time between GRB detection and observation. New missions such as Einstein Probe (W. Yuan et al. 2022), and the Space Variable Objects Monitor (J. L. Atteia et al. 2022) will also likely increase the sample of high-redshift GRBs. Additionally, JWST (M. A. Greenhouse 2016) and a new generation of 30 m telescopes (B. Neichel et al. 2018), along with new instruments such as SCORPIO (M. Robberto et al. 2020) on the Gemini Telescope, will provide more high-quality optical-to-near-IR spectra for GRB damping-wing analyses, enabling better constraints on the progression of the EoR. In particular, the simultaneous channels of SCORPIO will be easier to calibrate than an instrument like X-shooter, which has curved orders and three separate arms for UV, optical, and near-IR observations, so it may cut down on correlated noise and uncertainties in the spectrum and allow for more precise estimates of the neutral fraction.

## 6. Conclusions

GRBs are excellent probes of the high-redshift Universe. The Lyα damping wing of high-redshift GRBs can provide insight into the neutral fraction at different redshifts and track the progression of the EoR. GRB 130606A is a high-redshift

GRB for which multiple analyses using data sets from different telescopes and varying assumptions found different neutral fraction results. We reproduce all results using the VLT X-shooter spectrum and the corresponding assumptions of each analysis, highlighting the notable impact that assumptions can have on neutral fraction results. We present new analyses using assumptions motivated by new insights and multiple models, to ensure the robustness of the results. For the original J. Miralda-Escude (1998) model, the M. McQuinn et al. (2008) model, and a shell implementation of the J. Miralda-Escude (1998) model, the statistically preferred results give a $3\sigma$ neutral fraction upper limit of $x_{H_i} \lesssim 0.28$, $x_{H_i} \lesssim 0.24$, and $x_{H_i} \lesssim 0.26$, respectively. We compare these results to the neutral fraction analysis of GRB 210905A, which resides at a slightly higher redshift, and present both GRB damping-wing neutral fraction estimates in the context of neutral fraction measurements from other probes and EoR models. More high-redshift GRBs will be vital for probing the EoR at different redshifts, and reducing the reliance on the lines of sight of individual GRBs.

## Acknowledgments

We thank the anonymous referee for their constructive feedback. We also thank Ryan Chornock for providing us with the Gemini spectrum of GRB 130606A. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

*Software:* astropy (Astropy Collaboration et al. 2013, 2018, 2022), emcee (D. Foreman-Mackey et al. 2013), corner (D. Foreman-Mackey 2016), harmonic (J. D. McEwen et al. 2021).

# Appendix
## Shell Implementation Fits Posteriors

In this appendix, we show the posterior distribution associated with the independent- and dependent-shell implementations of the J. Miralda-Escude (1998) fit (see Figures 10, 11 and 12).
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Posteriors for the independent-shell implementation of the J. Miralda-Escude (1998) model for shells with widths of $\Delta z \sim 0.1$ and $z_{IGM,\text{u}}$ fixed to $z_{\text{host}}$. The neutral fraction posterior is most densely populated around zero, with increasing $3\sigma$ upper limits for redshifts further from the GRB.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Posterior distributions for the dependent-shell implementation of the J. Miralda-Escude (1998) model for shells of width $\Delta z = 0.1$ and $z_{\text{IGM,u}}$ treated as a free parameter. $z_{\text{IGM,u}}$ tends toward farther redshifts, indicating a large ionized bubble around the GRB host galaxy. The neutral fraction, and the slope of the neutral fraction with redshift, are unconstrained.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Posteriors for the dependent-shell implementation of the J. Miralda-Escude (1998) model for shells of width $\Delta 0.1$, $z_{IGM,u} = 5.8$, and a Gaussian spectral index prior ($\mu = 0.71$, $\sigma = 0.07$). $x_{H_i}$ does not significantly deviate from zero for all shells, and the slope of the neutral fraction as a function of redshift is unconstrained. The correlation between the posteriors of the neutral fractions in different shells is due to coupling between them according to the slope $dx_{H_i}/dz$.</figcaption></figure>

## ORCID iDs

H. M. Fausey https://orcid.org/0000-0002-2927-2398

A. J. van der Horst https://orcid.org/0000-0001-9149-6707

N. R. Tanvir https://orcid.org/0000-0003-3274-6336

K. Wiersema https://orcid.org/0000-0002-9133-7957

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJL, 848, L13

Afonso, P., Kann, D. A., Nicuesa Guelbenzu, A., et al. 2013, GCN, 14807, 1

Amati, L., O'Brien, P. T., Götz, D., et al. 2021, ExA, 52, 183

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

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Atteia, J. L., Cordier, B., & Wei, J. 2022, IJMPD, 31, 2230008

Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Natur, 553, 473

Barkana, R., & Loeb, A. 2004, ApJ, 601, 64

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

Beniamini, P., & Granot, J. 2016, MNRAS, 459, 3635
Bolan, P., Lemaux, B. C., Mason, C., et al. 2022, MNRAS, 517, 3263

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJL, 752, L5

Bruton, S., Lin, Y.-H., Scarlata, C., & Hayes, M. J. 2023, ApJL, 949, L40

Burns, E., Svinkin, D., Fenimore, E., et al. 2023, ApJL, 946, L31

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

Castro-Tirado, A. J., Sánchez-Ramírez, R., Ellison, S. L., et al. 2013, arXiv:1312.5631

Chen, H. 2024, MNRAS, 528, L33

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

Chornock, R., Berger, E., Fox, D. B., et al. 2013, ApJ, 774, 26

Christensen, L., Fynbo, J. P. U., Prochaska, J. X., et al. 2011, ApJ, 727, 73

Ciardi, B., & Loeb, A. 2000, ApJ, 540, 687

Cucchiara, A., Levan, A. J., Fox, D. B., et al. 2011, ApJ, 736, 7

Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018, ApJ, 864, 142

Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Natur, 340, 126

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177

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

Fausey, H. M., Vejlgaard, S., van der Horst, A. J., et al. 2025, MNRAS, 536, 2839

Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJL, 946, L13

Finkelstein, S. L., D'Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36

Foreman-Mackey, D. 2016, JOSS, 1, 24

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

Frederiks, D. D., Hurley, K., Svinkin, D. S., et al. 2013, ApJ, 779, 151

Fryer, C. L., Lien, A. Y., Fruchter, A., et al. 2022, ApJ, 929, 111

Galama, T. J., Vreeswijk, P. M., van Paradisj, J., et al. 1998, Natur, 395, 670

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

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

Goldoni, P. 2011, AN, 332, 227

Gou, L. J., Mészáros, P., Abel, T., & Zhang, B. 2004, ApJ, 604, 508

Greenhouse, M. A. 2016, Proc. SPIE, 9904, 990406

Greig, B., Mesinger, A., & Bañados, E. 2019, MNRAS, 484, 5094

Greig, B., Mesinger, A., Davies, F. B., et al. 2022, MNRAS, 512, 5390

Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5

Hartog, O. E., Malesani, D., Fynbo, J. P. U., et al. 2015, A&A, 580, A139

Hjorth, J., Sollerman, J., Møller, P., et al. 2003, Natur, 423, 847

Hoag, A., Bradač, M., Huang, K., et al. 2019, ApJ, 878, 12

Hook, I. M., Jørgensen, I., Allington-Smith, J. R., et al. 2004, PASP, 116, 425

Hsiao, T. Y.-Y., Abdurro'uf, Coe, D., et al. 2024, ApJ, 973, 8

Inoue, A. K., Hasegawa, K., Ishiyama, T., et al. 2018, PASJ, 70, 55

Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73

Izotov, Y. I., Worseck, G., Schäerer, D., et al. 2021, MNRAS, 503, 1734

Jin, X., Yang, J., Fan, X., et al. 2023, ApJ, 942, 59

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

Jung, I., Finkelstein, S. L., Dickinson, M., et al. 2020, ApJ, 904, 144

Kakiichi, K., Ellis, R. S., Laporte, N., et al. 2018, MNRAS, 479, 43

Kashikawa, N., Aoki, K., Asai, R., et al. 2002, PASJ, 54, 819

Kass, R. E., & Raftery, A. E. 1995, JASA, 90, 773

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

Konno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, 70, S16

Kouveliotou, C., Meegan, C. A., Fishman, G. J., et al. 1993, ApJL, 413, L101

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

Levan, A. J., Gompertz, B. P., Salafia, O. S., et al. 2024a, Natur, 626, 737

Levan, A. J., Jonker, P. G., Saccardi, A., et al. 2024b, arXiv:2404.16350

Li, L., Wu, X.-F., Huang, Y.-F., et al. 2015, ApJ, 805, 13

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, ApJL, 725, L234

Madau, P. 1995, ApJ, 441, 18

Mason, C. A., Fontana, A., Treu, T., et al. 2019, MNRAS, 485, 3947

Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, ApJ, 856, 2

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Mazets, E. P., Golenetskii, S. V., Ilinskii, V. N., et al. 1981, Ap&SS, 80, 3

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

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

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

Melandri, A., Bernardini, M. G., D'Avanzo, P., et al. 2015, A&A, 581, A86

Mesinger, A., & Furlanetto, S. R. 2008, MNRAS, 385, 1348

Meyer, R. A., Kakiichi, K., Bosman, S. E. I., et al. 2020, MNRAS, 494, 1560

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

Morales, A. M., Mason, C. A., Bruton, S., et al. 2021, ApJ, 919, 120

Morishita, T., Roberts-Borsani, G., Treu, T., et al. 2023, ApJL, 947, L24

Mostardi, R. E., Shapley, A. E., Steidel, C. C., et al. 2015, ApJ, 810, 107

Muñoz, J. B., Mirocha, J., Chisholm, J., Furlanetto, S. R., & Mason, C. 2024, MNRAS, 535, L37

Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109

Narayan, R., Paczynski, B., & Piran, T. 1992, ApJL, 395, L83

Neichel, B., Mouillet, D., Gendron, E., et al. 2018, in SF2A-2018: Proc. Annual Meeting of the French Society of Astronomy and Astrophysics, ed. P. Matteo, F. Billebaud, F. Herpin et al. (Paris: SF2A), 3

Noll, S., Kausch, W., Barden, M., et al. 2012, A&A, 543, A92

Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13

Pahl, A. J., Shapley, A., Steidel, C. C., Chen, Y., & Reddy, N. A. 2021, MNRAS, 505, 2447

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

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

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

Rastinejad, J. C., Gompertz, B. P., Levan, A. J., et al. 2022, Natur, 612, 223

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

Robberto, M., Roming, P. W. A., van der Horst, A. J., et al. 2020, Proc. SPIE, 11447, 1144774

Robertson, B. E. 2022, ARA&A, 60, 121

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

Rossi, A., Rothberg, B., Palazzi, E., et al. 2022, ApJ, 932, 1

Rutkowski, M. J., Scarlata, C., Haardt, F., et al. 2016, ApJ, 819, 81

Saccardi, A., Vergani, S. D., De Cia, A., et al. 2023, A&A, 671, A84

Salvaterra, R. 2015, JHEAp, 7, 35

Salvaterra, R., Della Valle, M., Campana, S., et al. 2009, Natur, 461, 1258

Sari, R., Piran, T., & Halpern, J. P. 1999, ApJL, 519, L17

Sari, R., Piran, T., & Narayan, R. 1998, ApJL, 497, L17

Savaglio, S. 2006, NJPh, 8, 195

Sparre, M., Toogood, O. E., Krühler, T., et al. 2014, ApJ, 785, 150

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

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123

Tanvir, N. R., Fox, D. B., Levan, A. J., et al. 2009, Natur, 461, 1254

Tanvir, N. R., Fynbo, J. P. U., de Ugarte Postigo, A., et al. 2019, MNRAS, 483, 5380

Tanvir, N. R., Laskar, T., Levan, A. J., et al. 2018, ApJ, 865, 107

Tanvir, N. R., Levan, A. J., Fruchter, A. S., et al. 2012, ApJ, 754, 46

Thompson, C. 1994, MNRAS, 270, 480

Thöne, C. C., Fynbo, J. P. U., Goldoni, P., et al. 2013, MNRAS, 428, 3590

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

Totani, T., Aoki, K., Hattori, T., et al. 2014, PASJ, 66, 63

Totani, T., Kawai, N., Kosugi, G., et al. 2006, PASJ, 58, 485

Ukwatta, T. N., Stamatikos, M., Maselli, A., et al. 2013, GCNR, 444, 1

Umeda, H., Ouchi, M., Nakajima, K., et al. 2024, ApJ, 971, 124

Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389

Vanzella, E., de Barros, S., Vasei, K., et al. 2016, ApJ, 825, 41

Vergani, S. D. 2013, in SF2A-2013: Proc. Annual Meeting of the French Society of Astronomy and Astrophysics, ed. L. Cambresy et al. (Paris: SF2A), 377

Vernet, J., Dekker, H., D'Odorico, S., et al. 2011, A&A, 536, A105

Vielfaure, J. B., Vergani, S. D., Jabelj, J., et al. 2020, A&A, 641, A30

Wang, F., Davies, F. B., Yang, J., et al. 2020, ApJ, 896, 23

White, N. E., Bauer, F. E., Baumgartner, W., et al. 2021, Proc. SPIE, 11821, 1182109

Whitler, L. R., Mason, C. A., Ren, K., et al. 2020, MNRAS, 495, 3602

Witstok, J., Jakobsen, P., Maiolino, R., et al. 2025, Natur, 639, 897

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

Wu, Z., & Kravtsov, A. 2024, OJAp, 7, 56

Yang, J., Wang, F., Fan, X., et al. 2020, ApJL, 897, L14

Yuan, W., Zhang, C., Chen, Y., & Ling, Z. 2022, in Handbook of X-ray and Gamma-ray Astrophysics, ed. C. Bambi & A. Sangangelo (Singapore: Springer), 86

Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2022, ApJ, 932, 76

Zhu, Y., Becker, G. D., Bosman, S. E. I., et al. 2024, MNRAS, 533, L49