Future Constraints on the Reionization History and the Ionizing Sources from Gamma-ray Burst Afterglows

We forecast the reionization history constraints, inferred from Lyman-alpha damping wing absorption features, for a future sample of $\sim 20$ $z \geq 6$ gamma-ray burst (GRB) afterglows. We describe each afterglow spectrum by a three-parameter model. First, L characterizes the size of the ionized region (the"bubble size") around a GRB host halo. Second, $\langle{x_{\rm HI}\rangle}$ is the volume-averaged neutral fraction outside of the ionized bubble around the GRB, which is approximated as spatially uniform. Finally, $N_{\mathrm{HI}}$ denotes the column-density of a local damped Lyman-alpha absorber (DLA) associated with the GRB host galaxy. The size distribution of ionized regions is extracted from a numerical simulation of reionization, and evolves strongly across the Epoch of Reionization (EoR). The model DLA column densities follow the empirical distribution determined from current GRB afterglow spectra. We use a Fisher matrix formalism to forecast the $\langle{x_{\rm HI}(z)\rangle}$ constraints that can be obtained from follow-up spectroscopy of afterglows with SNR = 20 per R=3,000 resolution element at the continuum. We find that the neutral fraction may be determined to better than 10-15\% (1-$\sigma$) accuracy from this data across multiple independent redshift bins at $z \sim 6-10$, spanning much of the EoR, although the precision degrades somewhat near the end of reionization. A more futuristic survey with $80$ GRB afterglows at $z \geq 6$ can improve the precision here by a factor of $2$ and extend measurements out to $z \sim 14$. We further discuss how these constraints may be combined with estimates of the escape fraction of ionizing photons, derived from the DLA column density distribution towards GRBs extracted at slightly lower redshift. This combination will help in testing whether we have an accurate census of the sources that reionized the universe.


INTRODUCTION
Gamma ray bursts (GRBs) are, electromagnetically, the brightest explosions in the universe. The longduration class of GRBs are thought to be powered by the collapse of massive rotating stars (Woosley 1993): as nuclear fuel is exhausted at the end of the star's life, the core loses pressure support and implodes to form a black hole. In this process, jets of energetic particles are produced, which plow through the remaining stellar envelope. This relativistic outflow leads to highly-beamed gamma-ray emission, while subsequent afterglows occur at longer wavelengths as the ejecta runs into surrounding gas in the circumburst and interstellar media (ISM). The incredible luminosity of these events makes longduration GRBs detectable from even some of the earliest phases of cosmic history. Indeed, long-duration GRBs have been observed back into the Epoch of Reionization (EoR) (Tanvir et al. 2018;Melandri et al. 2015;Chornock et al. 2013;Cucchiara et al. 2011;Salvaterra et al. 2009;Tanvir et al. 2009;Totani et al. 2006), when early generations of galaxies and/or accreting black holes emitted ionizing photons and gradually filled the universe with ionized gas (e.g., Loeb & Furlanetto 2013). Therefore, GRB afterglows can be exploited as bright backlights, with follow-up spectroscopy of the afterglows at high spectral resolution and signal-to-noise ratio providing valuable constraints on the neutral fraction of the intergalactic medium (IGM) (e.g., Totani et al. 2006;Hartoog et al. 2015) and the chemical enrichment history of the ISM in the GRB host galaxy (e.g., Lamb & Reichart 2000;Savaglio 2006).
A number of projects are being proposed to study reionization and metal enrichment using GRB afterglows. These include the Gamow Explorer which is expected to find roughly ∼ 20 GRB afterglows in the EoR at z ≥ 6 using an X-ray telescope with a wide field of view, exploiting that the GRB emission at high redshifts is largely shifted into the X-ray band (White 2020;Ghirlanda et al. 2015). The Gamow Explorer will then determine photometric redshifts for these sources, using an on-board infrared telescope, from the drop in flux blueward of Lyman-alpha (Ly-α) at the source redshift. This will enable rapid targeted follow-up spectroscopic campaigns using the James Webb Space Telescope (JWST) in space or existing eight-meter and future thirty-meter class telescopes from the ground. The Gamow Explorer is proposed to launch in 2028, with a nominal three year mission lifetime. The detectable GRB rate during the EoR from the Gamow Explorer should improve on the Swift telescope's (Gehrels et al. 2004) detection rate of such events by roughly an order of magnitude. In addition, HiZ-GUNDAM, which is similar in scope and timescale to the Gamow Explorer, is being planned in Japan (Yonetoku et al. 2014). In the slightly longer term, a larger European-led satellite, Theseus, is being proposed for launch in 2032, and should detect between 30-80 z ≥ 6 GRBs (Amati et al. 2018;Stratta et al. 2018;Amati et al. 2021).
Motivated by these future prospects, the aim of the current paper is to forecast the reionization history constraints expected from these upcoming missions. The key spectral signature of remaining neutral gas in the IGM is the damping wing of the Ly-α line owing to its natural line width (Miralda-Escude 1998). This leads to spectrally extended absorption redward of Ly-α at the source redshift, provided significantly neutral gas -with a neutral hydrogen fraction of order unity -remains in the IGM at the GRB redshift.
In this regard, GRB afterglows offer a few potential advantages over quasars (Miralda-Escude 1998; Barkana & Loeb 2004) which have nevertheless provided our best backlights for studying the IGM thus far. First, the intrinsic unabsorbed spectrum of a GRB afterglow has a simple power-law form around Ly-α at the source redshift. In contrast, in the case of quasars, the damping wing signature needs to be extracted in the presence of a Ly-α emission line (e.g., Davies et al. 2018a). Second, the quasars themselves can ionize significant regions of neutral hydrogen around them. Furthermore, these sources likely reside in highly biased locations where surrounding galaxies form earlier than in typical regions of the universe. Reionization may hence occur earlier in quasar environments than in average parts of the universe (Lidz et al. 2007). Together, the ionizing impact of quasars and surrounding galaxies may lead to weaker damping wing signatures around quasars than around less-biased sources such as long duration GRBs, which are expected to trace more typical regions of massive star formation. This may make the damping wing signature both more difficult to detect around quasar sources, and also complicate the interpretation of quasar-based measurements which require modeling the relationship between the quasar environments and typical parts of the IGM. Third, the quasar luminosity function falls very steeply with increasing redshift and current forecasts suggest that only a single quasar (with UV absolute magnitude M AB ≤ −26) is detectable in the observable universe at z 9 (Fan et al. 2019). GRB afterglows provide the best candidates for observing damping wing signatures in the early phases of the EoR.
The GRB afterglow measurements are not without their challenges, however. First, one needs to detect the GRBs and their afterglows, rapidly determine photometric redshifts, and obtain high quality follow-up spectroscopic observations of the fading afterglow. Second, current samples of GRB afterglows frequently show strong damped Ly-α (DLA) absorption from neutral hydrogen in the ISM of the host galaxy (e.g., Tanvir et al. 2019). In some cases, the host DLA can swamp the signature of damping wing absorption from neutral gas in the IGM (Miralda-Escude 1998). Third, a challenge for interpreting damping wing signatures of the IGM towards both GRB afterglows and quasar backlights is the patchiness of the reionization process (McQuinn et al. 2008;Mesinger & Furlanetto 2008). This leads to a significant sightline-to-sightline scatter in the neutral fraction towards these background sources. More importantly, it implies that the damping wing profile is at least partly sensitive to the precise distribution of neutral gas along each line of sight and it is unclear how to best model the resulting spectra.
Here we consider a simple parameterization (see also, e.g., McQuinn et al. 2008) that partly accounts for the patchiness of the reionzation process, includes the plausible impact of host DLAs, and allows for rapid, fully analytic forecasts of the reionization history constraints obtainable from future GRB afterglow spectral samples. At each redshift, afterglow spectra are described by a bubble size parameter, L, a neutral fraction, x HI (approximated as uniform outside the host bubble), and a host DLA column density, N HI . The bubble size and column densities are drawn from probability distributions, inferred respectively from reionization simulations and observations, to allow for the likely broad sightline-tosightline scatter in these quantities. In this framework, we use the Fisher matrix formalism to determine the neutral fraction constraints across several redshift bins, marginalized over the bubble size and column density parameters.
In addition to providing backlights for studying the IGM damping wings, GRBs also provide valuable information about the escape fraction of ionizing photons (Chen et al. 2007;Tanvir et al. 2019) and the starformation rate density (SFRD) at high redshifts. Specifically, a plausible working assumption is that GRB afterglows provide lines of sight towards massive stars which are representative of the pathways traversed by ionizing photons, as they propagate out into the IGM from the sources of reionization. The fraction of afterglows that have host columns less than roughly 10 18 cm −2 , corresponding to Lyman continuum optical depths at the photoionization threshold of τ HI 6.3, then provides an estimate of the escape fraction. 4 Further, if the conversion between GRB rate and SFRD can be calibrated accurately at lower redshift, the redshift evolution of the GRB rate can be used to determine the SFRD during the EoR, including contributions from lower mass host galaxies that are challenging to observe directly with other star formation indicators (e.g., Blain & Natarajan 2000;Kistler et al. 2008;Robertson & Ellis 2012;Vergani et al. 2015). Measurements of the neutral fraction can be combined with estimates of the escape fraction and internal or external determinations of the SFRD to address whether the neutral fraction constraints are consistent with known source populations. We consider this in the context of our Fisher matrix analyses.
The outline of this paper is as follows. In §2 we review the IGM damping wing signature and the impact of DLA host absorbers. We then describe the key ingredients of our model: the redshift distribution of GRB afterglows ( §2.1), the reionization history model ( §2.2), the bubble size distribution ( §2.3), and the probability distribution of host column densities ( §2.4). §3 describes our Fisher matrix calculations and characterizes parameter degeneracies, while §4 forecasts the resulting neutral fraction constraints. In this section we further discuss the dependence of these constraints on the quality of follow-up spectroscopic observations and the possibilities for obtaining the requisite spectra. §5 discusses the prospects for determining the escape fraction of ionizing photons and the resulting implications. We conclude in §6.

METHOD AND MODEL PARAMETERS
We start by reviewing the physics of the Ly-α damping wing, which is the key feature for probing the IGM neutral fraction with GRB afterglows (Miralda-Escude 1998). The Ly-α optical depth close to line center is large even in the case of a highly ionized universe at the redshifts of interest, as evident in Eq. 2 below. Consequently, the traditional Ly-α forest blueward of Ly-α at the source redshift is expected to be highly opaque at z 6 or so (as also indicated by current observations, e.g., Fan et al. 2006;Yang et al. 2020), regardless of the precise ionization state of the IGM. If significantly neutral gas -with a neutral hydrogen fraction of order unity -remains in the IGM, however, the broad damping wing of the Ly-α line owing to its natural line width, leads to an extended absorption feature redward of Ly-α at the source redshift. In contrast, the probability of absorption in the wing is negligible if the IGM is instead highly ionized.
In what follows, we ignore the impact of peculiar velocities, consider only the wing of the Ly-α line (i.e., we do not treat absorption in the Doppler core), and approximate the gas density by its cosmic mean value. These should be good approximations for our goal of modeling the extended damping wing redward of line center. Further, we adopt the approximation to the line profile in Miralda-Escude (1998), which includes classical Rayleigh scattering and the Lorentzian form for the natural line width (see Lee 2013 for refinements, which we will neglect here as they should be unimportant for our forecasts).
In this case, consider a significantly neutral portion of the IGM with volume-averaged neutral fraction, x HI , between redshift z beg and z end along the line of sight to a GRB afterglow at some higher source redshift, z s . For simplicity, in our modeling we place each GRB at the center of its redshift bin and consider bins of width ∆z = 1, as described further below. The optical depth in the damping wing of the Ly-α line from this neutral patch -at wavelength offset ∆λ from Ly-α at the source redshift -is denoted by τ DW (∆λ). Hence the observed wavelength is λ obs = λ α (1 + z s ) + ∆λ, where λ α = 1215.67Å is the rest-frame wavelength of the Ly-α transition. We further denote the fractional wavelength offset as δ = ∆λ/[λ α (1 + z s )]. The damping wing optical depth is then (see Miralda-Escude 1998; we mostly follow the notation in that work): Here R α is a dimensionless constant related to the natural line-width of the Ly-α line, with R α = Λ α /(4πν α ), where Λ α is the Ly-α decay constant and ν α is the restframe frequency of the Ly-α transition. Numerically, R α = 2.02 × 10 −8 . Here τ GP denotes the Gunn-Peterson optical depth of fully neutral gas at the cosmic mean density (Gunn & Peterson 1965). The Gunn-Peterson optical depth is: where Ω m and Ω b are the usual cosmological matter and baryon density parameters, h in the Hubble parameter in units of 100 km/s/Mpc, and 1 − Y is the hydrogenic mass fraction. Finally, In what follows, we imagine that each GRB host is located within an ionized region and so z beg is set by the distance between the GRB host halo and the first significantly neutral region encountered along the line of sight to the observer. As discussed further below, this distance is drawn from a distribution based on numerical simulations of reionization. Next, we approximate the exterior neutral fraction as uniform outside of the ionized bubble around the GRB host. Our baseline model further adopts the average neutral fraction x HI at the center of the redshift bin in question and includes contributions only from neutral gas within the redshift bin. In other words, z end in Eqs. 1 and 3 is set by the lower redshift end of the redshift bin under consideration. Put differently, we neglect redshift evolution in the neutral fraction across each redshift bin, and ignore contributions from neutral gas at lower redshifts outside of the bin. Although the precise assumptions here are arbitrary and we will test their impact in what follows, note that the damping wing is fairly insensitive to neutral gas at redshifts much smaller than the source redshift (see Appendix B and Figure 15).
Throughout we assume perfect knowledge of the GRB host redshift. We suppose that the GRB redshifts may be determined accurately through measuring the redshifts of host-galaxy metal lines detected in absorption, perhaps associated with a local DLA (e.g., Totani et al. 2006). Alternatively, follow-up observations after the afterglow fades can be used to determine the redshift of the host galaxy using emission lines (e.g., Krühler et al. 2015).

Assumed Redshift Distribution
The first key ingredient in our forecasts is the assumed redshift distribution of GRB afterglows. Provided the GRB rate traces the SFRD, the number of GRBs per unit redshift may be written as (Kistler et al. 2008): whereρ denotes the SFRD, i.e., the star formation rate per co-moving volume per unit time, dV co /dz is the comoving volume per unit redshift, and N 0 is a normalization constant (with units of observing time). The factor of 1 + z in the denominator accounts for time dilation. We adopt the empirical fitting formula for the SFRD from Madau & Dickinson (2014). We ignore any redshift dependence in the ability of the upcoming surveys -Assumed redshift distributions for the rate of longduration GRBs during the EoR. The blue histogram shows our fiducial model in which the GRB rate traces the SFRD, and the total number of GRBs above z = 6 is normalized to 20, as forecast for the Gamow Explorer. The red histogram shows an alternate case, in which the GRB rate is enhanced at high redshift (see text), which is also considered in what follows. Finally, the black histogram shows the upper end of the predicted GRB rate for the planned Thesesus mission.
to detect the GRBs and determine their redshifts in this equation, which should be a good approximation for the bright bursts targeted by the Gamow Explorer and related missions.
Our fiducial model fixes the normalization constant by requiring N GRB (≥ z = 6) = ∞ 6 dz dN GRB /dz = 20. This is based on detailed projections of the number of high redshift GRBs detectable with the Gamow Explorer (White 2020), including models for the GRB luminosity function, prompt emission and afterglow spectral shape, observing time, and the Gamow instrumental sensitivity (Ghirlanda et al. 2015). In order to explore some of the uncertainties involved we also consider a model with an enhanced GRB rate at higher redshifts: in this alternate scenario, we multiply the left-hand side of Eq 4 by (1 + z) 1.5 , normalized to the same dN GRB /dz at z = 6. This case yields 31 GRBs at z ≥ 6. This is intended to roughly represent scenarios where either the SFRD estimates from Madau & Dickinson (2014) are incomplete at z 6, owing for example to star formation in faint galaxies, or possible evolution in the GRB to star-formation rate relationship. Finally, in order to explore the longer term best-case scenario prospects we also consider a case with N GRB (z ≥ 6) = 80 GRBs. This is at the upper end of the range of 30-80 z ≥ 6 GRBs forecast for the Theseus mission (Amati et al. 2018).
The GRB redshift distributions in each of these three models are shown in Figure 1, for redshift bins of width ∆z = 1. In what follows, we adopt the nearest integer to the average afterglow count forecast in each bin. In the blue and red historgram models, 11 afterglows lie in the redshift bin centered around z s = 6, and at least one GRB is detectable out to the z s = 10 bin in the fiducial case and z s = 12 in the more optimistic scenario. In our fiducial reionization history, described below, these redshifts span most of the reionization process (see Figure 2), and should push past the highest redshift that will be probed using quasars, even with futuristic surveys, as mentioned in the Introduction). In the longer term, the optimistic upper-end of the Theseus projections gives 44 GRBs in our z s = 6 redshift bin, 5 in the z s = 10 bin and one GRB, on average, in the z s = 14 bin. 5 2.2. Fiducial reionization history -Fiducial model reionization history compared to some current observational constraints. The solid black curve shows our fiducial theoretical model described in the text below. The red triangles show upper limits on the neutral fraction from the dark pixel fraction in the Ly-α and Ly-β forests (McGreer et al. 2015). The green squares show estimates of the neutral fraction from damping wing measurements towards quasars (Davies et al. 2018b). Finally, the blue pentagon and magenta star, respectively, show a lower limit and a measurement of the neutral fraction using observations of Ly-α emission lines towards Lyman-break galaxies (Mason et al. 2018. In all cases, the error bars are 1 − σ confidence estimates. Our fiducial reionization history model is based on an approximate photon-counting description of reionization, in which the volume-averaged ionization fraction is determined through the following equation (Shapiro & Giroux 1987;Madau et al. 1999): Here the first term describes the rate of ionizing photon production, which is assumed to be proportional to the time derivative of the halo collapse fraction, above a minimum galaxy host halo mass, M min . The second term accounts for recombinations, andt rec gives the average time for ionized gas to recombine, which may be written ast rec = C/[α B (T ) n e (z) ] where C is a clumping factor, α B is the case-B recombination coefficient of hydrogen, and n e (z) is the proper number density of free electrons at redshift z. Our fiducial model adopts an ionizing efficiency parameter, ζ = 20, a minimum galaxy host halo mass of M min = 10 9 M , a gas temperature of T = 2 × 10 4 K, and a clumping factor of C = 3 (see e.g., Lidz 2015 and references therein for a discussion of these parameters). The halo collapse fraction is calculated using the Sheth-Tormen halo mass function (Sheth et al. 2001). The ionization fraction is fixed to unity below the redshift where the solution of Eq. 5 first reaches this limiting value. The volume-averaged neutral fraction is determined through x HI (z) = 1 − x i (z) . The resulting ionization history is shown in Figure 2. We emphasize that the fiducial model here is intended only as a plausible baseline for our ensuing forecasts, and has not been designed or tuned to match current constraints in detail. It is nevertheless in good agreement with current measurements, as illustrated in the figure. Furthermore, the electron scattering optical depth in our fiducial model is τ e = 0.51 (assuming helium is singly ionized along with hydrogen and neglecting the likely small contribution from the second ionization of helium), consistent with recent Planck 2018 measurements: the fiducial model τ e value differs from their TT, TE, EE, lowE, lensing + BAO constraint by less than 1 − σ (Aghanim et al. 2018).

Bubble size distribution
The next important ingredient in our model is to account for the ionized bubbles around GRB host halos. In general, these regions reflect the combined ionizing output of many surrounding galaxies, with the typical size scale growing as ionized bubbles expand and merge over the course of the reionization process (e.g., McQuinn et al. 2007b). In order to appreciate the role of the ionized bubbles here, consider the example of a GRB that occurs in a large ionized region. More precisely, the relevant thing here is the the extent of the ionized skewer traversed along the line of sight from the GRB to the observer before encountering neutral gas in the IGM. 6 In this case, photons initially just redward of Lyα at the source redshift, for instance, will redshift significantly before encountering neutral gas and the source will show only weak damping-wing absorption. In contrast, GRB afterglow photons originating in small bubbles or at the edge of larger ones will redshift less before running into neutral gas and the damping wing feature will be stronger. Therefore the size of the ionized bubbles around GRB host halos plays an important role in determining the strength of the damping wing absorption.
In order to account for this in our model, we turn to the -The mean bubble (or ionized skewer) size and its scatter -as a function of the volume-averaged neutral fraction -extracted from reionization simulations. The black solid line gives the average co-moving distance between a GRB host and the first neutral region, while the blue dashed lines indicate the sightline to sightline scatter around the mean. The mean and scatter grow sharply as reionization proceeds and the average neutral fraction drops. reionization simulations of McQuinn et al. (2008McQuinn et al. ( , 2007a, which treat the radiative transfer of ionizing photons in a post-processing step, on top of an evolved N-body simulation of cosmological structure formation. From a number of simulation snapshots spanning much of the EoR, we extract sightlines in random directions towards dark matter halos with mass between 10 10 M ≤ M ≤ 10 11 M and record the co-moving distance between each simulated dark matter halo and the first significantly neutral cell (defined by when the cell-averaged neutral fraction exceeds 0.9). 7 These halos are well-resolved in the simulations and are plausible hosts for GRBs during the EoR. The reionization history in the simulations differs a little from our fiducial model, and so we record the first-crossing distance as a function of the volumeaveraged neutral fraction rather than redshift. In each GRB redshift bin, we then determine the mean ionized skewer size and the (1-σ) scatter from the simulation output with average neutral fraction closest to the fiducial model value at the desired redshift. Relying on simulation outputs at the correct neutral fraction, but slightly different redshifts, should be a good approximation: previous studies have shown that the spatial structure of reionization is sensitive mainly to the average ionization fraction and is fairly independent of the precise redshift when a given stage of reionization is reached (McQuinn et al. 2007b). This first-crossing distance then sets the upper limit of the integral over redshift, z beg in Eqs 1 and 3. 8 To briefly reiterate, the main approximation in our scheme is to ignore spatial variations in the neutral fraction outside of the first ionized bubble. We expect this to be most accurate in the early stages of reionization when the neutral fraction is large and the size of typical ionized bubbles is small relative to the extent of the damping wing. Especially towards the end of reionization, a more elaborate parameterized model involving the size of the first neutral region, the extent of the second ionized bubble encountered, etc. may be required for more accurate results. We confine our current study, however, to the simpler three-parameter model. Figure 3 shows the resulting distribution of ionized skewer sizes. The mean skewer size (in co-moving units) evolves very strongly during the EoR, ranging from roughly ∼ 1 Mpc/h when the universe is ∼ 95% neutral, to ∼ 10 Mpc/h when it is ∼ 50% neutral, and to ∼ 60 Mpc/h when about ∼ 5% of the IGM volume is neutral. The scatter also grows strongly as reionization proceeds.
The size distribution of ionized skewers may be conveniently approximated by a lognormal distribution. Specifically, where L is the bubble size in units of Mpc/h. The lognormal parameters µ andσ are related to the mean ionized skewer size, L , and sightline to sightline scatter, σ L , via: andσ 2 = ln σ L L 2 + 1 (8) 2.4. Contamination from Host DLAs In addition to absorption from neutral hydrogen in the IGM, we must account for absorption from neutral hydrogen in the host galaxy. Empirically, as discussed further below, many GRB sightlines show strong DLA absorption. We assume that the DLA lies precisely at the redshift of the GRB host galaxy, which is likely an excellent approximation given that long duration GRBs are expected to trace star-forming regions, and most star formation occurs in relatively small mass halos at the redshifts of interest. For instance, the circular velocity of a 10 10 M halo at z s = 8 is v circ = 74 km/s (Barkana & Loeb 2001) and so a DLA moving at the halo circular velocity relative to the host galaxy is offset in wavelength by only ∆λ/λ ∼ v circ /c ∼ 2 × 10 −4 , comparable to the spectral resolution assumed in what follows ( §3). This is hence a negligible effect, although a DLA in a neighboring halo could in principle be more of a concern; an unlikely coincidence would, however, be required for the neighboring absorber to produce a large column density.
Since we are interested in the absorption in the wing of the line, we approximate the optical depth of the DLA at wavelength offset ∆λ using a Lorentzian profile: where σ α = 3λ 2 α Λ α /(8πν α ) = 4.48 × 10 −18 cm 2 is the frequency-integrated cross section. The transmission through the IGM accounting for the IGM damping wing and the DLA is then given from Eqs 1 and 9 as: (10) Fig. 4.-Example damping wing profile. The black solid line shows the IGM damping wing in our fiducial model at zs = 8 for the mean ionized skewer size, while the blue dashed lines show the profiles for skewers of length ±1 − σ above and below the mean. In each case z end = 7.5 (see text). The red dotted line shows the profile of a DLA at the host redshift with a column density of N HI = 10 20 cm −2 . The black dotted line shows the transmission including both the DLA optical depth and the IGM contribution for the mean skewer size. The more gradual IGM profile is distinguishable from the DLA profile for the parameters shown, although the spatial variations in the ionized skewer size lead to fairly strong variations in the transmission close to line center.
We will explore the redshift and parameter dependence of the transmission profiles in more detail in the next section, but Figure 4 gives an illustrative example for starters. Here we consider a GRB source at z s = 8, in which case the IGM is 70% neutral in our fiducial model ( Figure 2). The mean size of ionized skewers at this neutral fraction is L = 4.8 Mpc/h, while the sightline to sightline dispersion is σ L = 4.0 Mpc/h ( Figure 3). In this example, the figure illustrates that the DLA dominates the absorption close to line center for the average ionized skewer size, but ionized skewers with size 1 − σ lower than the mean value produce comparable absorption at line center to the DLA. However, for any skewer size within the ±1 − σ range around the mean in this example, we expect the IGM contribution to the damping wing to dominate over the DLA component far from line center (for N HI = 10 20 cm −2 ). That is, the IGM component, scaling roughly as τ DW ∝ 1/(∆λ) (Eq. 1) 9 , is more extended than the DLA contribution which goes like τ DLA ∝ 1/(∆λ) 2 (Eq. 9). As fleshed out further in §3 below, larger column density DLAs can nevertheless swamp the IGM damping wing contribution, especially at lower redshifts where the neutral fraction is smaller, the universe is less dense, and the ionized skewers from the GRB host halo tend to be larger. A key, yet uncertain, issue then regards the probability distribution of DLA column densities from local host absorption around GRBs during the EoR. The largest empirical sample currently addressing this question is published in Tanvir et al. (2019); these authors determined the column density distribution from a sample of 138 GRB afterglows ranging in redshift from 1.6 < z < 6.7. The probability distribution from the tabulated data in this study is shown in Figure 5. Interestingly, Tanvir et al. (2019) derive an average ionizing photon escape fraction of only f esc = 0.005, with a 98% upper limit of f esc ≤ 0.015. Indeed, the median column density in this sample is N HI = 10 21.59 cm −2 , relatively few sightlines have N HI ≤ 10 20 cm −2 , and only a tiny fraction of sightlines have column densities low enough to allow appreciable ionizing photon escape fractions.
An important question then is how strongly the probability distribution of DLA column densities evolves with redshift. The current sample includes only 7 lines of sight above z > 5, and while there is currently no evidence for redshift evolution in the column density probability distribution function (PDF) (Tanvir et al. 2019), a larger data set will be required to test this more sharply. In fact, provided this distribution is indeed a reliable indicator of the escape fraction of ionizing photons, one can indirectly argue that the distribution must evolve strongly: otherwise it is hard to understand how the universe could be reionized by z 5.5 (as is observed) given the escape fraction estimate in Tanvir et al. (2019). For instance, escape fractions larger than roughly f esc 0.1 are required to reionize the universe with star-forming galaxies, subject somewhat to uncertainties in: faint-end extrapolations down the galaxy UV luminosity function, the ionizing spectral shape, and the clumping factor of the IGM (e.g., Bouwens et al. 2016). This is a factor of ∼ 20 larger than the best fit value from the current GRB samples, which are however concentrated at lower redshifts. The strong evolution implied here is at least broadly consistent with inferences of the escape fraction at z ∼ 2.5 − 5 based on comparing galaxy UV luminosity functions to the ionizing emissivity extracted using the transmission through the Ly-α forest at these redshifts, along with estimates of the mean free path of ionizing photons (Faucher-Giguère 2020; Faucher-Giguere et al. 2008). Nevertheless, our fiducial model adopts the current DLA column density PDF since it is hard to predict how this may evolve towards higher redshift. We should keep in mind that this may be a pessimistic assumption (see also §4.1).

FISHER MATRIX FORMALISM
In order to forecast the constraints on the ionization history, we adopt a Fisher matrix approach. The Fisher matrix calculations assume a quadratic expansion around the logarithm of the likelihood function, which in turn peaks for a set of specified fiducial model parameters. As shown in Appendix A, we can incorporate the expected broad distributions in the ionized skewer size and DLA column density PDFs by first computing a conditional Fisher matrix for sightlines with specified bubble size L and column density N HI . We then find the expected error bar on the neutral fraction for that L and N HI by inverting the conditional Fisher matrix. Finally, the error bar on the neutral fraction for an ensemble of GRB afterglow sightlines can be determined by computing the average of the neutral fraction's inverse-variance, weighted by the bubble size and column density PDFs (see Eqs. 25 and 26). The error bar in each redshift bin is scaled with 1/ N GRB (z s ) where N GRB (z s ) is the number of independent afterglow spectra expected in the bin centered around source redshift z s . For each sightline, we adopt a parameter vector q with three components, (q xHI , q L , q NHI ). Here each component of q describes the fractional variation around the fiducial model value. For example, q xHI = ( x HI − x HI,fid )/ x HI,fid . Note also that the parameters q L and q NHI vary in meaning depending on the value of L and N HI for each sightline.
The conditional Fisher matrix for a sightline with bubble size L and column density N HI may be computed as: Here the sum is over spectral pixels, which are taken to be uniformly spaced in ln(λ) with ∆λ/λ = 1/3, 000, i.e., we consider a spectral resolution of R = λ/∆λ = 3, 000.
We sum over all pixels out to 1, 500Å redward (observed frame) of Ly-α at the source redshift. Appendix B demonstrates that our results are insensitive to this choice (see Figure 16). In general, we assume that the noise in each afterglow spectrum is set by a combination of Poisson fluctuations in the photons from the source itself and a contribution from the sky background. In this case, where N c is the number of photons per spectral resolution element received from the source at the continuum and N sky is the same for the sky background. We adopt a flat continuum spectrum, which should be a very good approximation across the range of wavelengths considered in Eq. 11. 10 The average number of photon counts in spectral pixel k is then given by N c f (k), i.e., the observed photon count is the number of photons at the continuum multiplied by the transmission. The signal-to-noise-ratio squared at the continuum is hence SNR 2 = N 2 c /(N c + N sky ). Our baseline assumption is that the signal-to-noise ratio at the continuum per spectral pixel is SNR = 20, and that N sky N c . That is, we assume that the observations are in the sky background dominated limit, where the source counts make a negligible contribution to the transmission variance. In this case, var [f (k)] = N sky /N 2 c = 1/SNR 2 . We discuss the possibilities for obtaining our baseline SNR value in §4.2 and show that the sky background dominated limit applies for groundbased observations, but that space-based observations with the JWST will be in the source dominated limit. We compare the sky background limited and source dominated cases in Appendix B and Figure 17.
We consider a redshift independent SNR for the purposes of having a single baseline noise value; this helps us understand the overall requirements for spectroscopic follow-up observations. In practice, the sky background is a strongly increasing function of wavelength and so in practice it will be hard to achieve the requisite SNR from the ground for higher redshift bursts (see §4.2). However, in the source-dominated case applicable to JWST spacebased observations, a redshift independent SNR may be a decent approximation. This is the case because given afterglows observed at a fixed time interval after the explosion, higher redshift ones are seen at an earlier time in the source frame -owing to cosmological time dilation -when they are brighter, compensating for the dimming from the increased luminosity distance (e.g., Barkana & Loeb 2004). We will consider variations around the SNR assumptions in §4.2.

Parameter Derivatives
It is instructive to first examine some examples of the transmission derivatives with respect to model parameters, as enter into Eq. 11. Specifically, Figure 6 shows 10 Throughout we assume that the unabsorbed continuum level is determined perfectly, i.e., we neglect continuum fitting errors in our analysis. derivatives with respect to the neutral fraction, bubble size, and DLA column density parameters for the case of our fiducial reionization history with L = L and N HI = 10 20 cm −2 at each of z s = 6, 8, and 10 (top, middle, and bottom panels, respectively). In general, increasing the DLA column density leads to reduced transmission relatively close to the redshift of the source, while increasing the bubble size parameter has the opposite effect, as expected. Increasing the neutral fraction reduces the transmission, and the impact of the neutral fraction parameter is more extended in wavelength than either of the DLA or bubble size parameters, owing to the gradual fall-off in the IGM damping wing profile of Eq. 1. The fairly similar shape, yet opposite sign, of the derivatives with respect to the bubble size and column density parameters suggests that there will be a fairly strong degeneracy direction between these parameters. That is, the enhanced transmission with increasing bubble size can be compensated by raising the DLA column density.
Turning to the redshift dependence in Figure 6, one can see that the neutral fraction parameter becomes more important than the DLA and bubble size parameters as one moves to higher redshift. This results because the strength of the IGM damping wing grows as the universe becomes more neutral and denser towards higher redshift and as the bubble sizes shrink. Note further that the transmission changes less with increasing DLA column density towards high redshift because the derivatives are taken around a fiducial transmission profile, which is more absorbed at high redshift. For the ex-ample DLA column density and bubble sizes, the IGM damping wing signature is swamped at z s ∼ 6 but becomes more prominent at z s 8. On the other hand, we expect to detect more GRB afterglows near z s ∼ 6 than at z s ∼ 8−10 (see Figure 1) and so there is a trade-off between the weaker damping wing signatures towards low redshift and the more abundant data samples available there. Also, we should keep in mind that the IGM damping wing signature is more detectable towards sightlines with smaller ionized bubbles and lower column density DLAs (see §4). Fig. 7.-Parameter degeneracies for the simplified case that all of the sightlines in the zs = 6 bins have bubble size L = L and N HI = 10 20 cm −2 . The ellipses show 1 − σ (red), 2 − σ (blue), and 3 − σ confidence intervals (green). The bottom right-hand corner gives the resulting 1-d likelihood, marginalized over bubble size and DLA column density parameters, for the average neutral fraction parameter in this redshift bin. We refine this error bar estimate in the next section by accounting for the full distribution of ionized skewer sizes and DLA column densities, but this simplified case nevertheless suffices for understanding the degeneracies involved.

Parameter Degeneracies
In order to quantify the ionization history constraints in detail, we must account for the full probability distributions in bubble size and DLA column density. First, however, it is helpful to consider parameter degeneracies for particular realizations of the bubble size and DLA column density. That is, we can invert the conditional Fisher matrix of Eq 11 and inspect the resulting parameter ellipses for particular values of L and N HI . For simplicity of illustration, we suppose here that all of the sightlines in the redshift bin of interest have the specified bubble size and DLA column density. Here we take our fiducial GRB redshift distribution with 20 z ≥ 6 GRBs ( Figure 1). Figure 7 shows an example case at z s = 6 with L = L = 32 Mpc/h and N HI = 10 20 cm −2 . The first conclusion one can draw from this figure is that the frac-tional constraints on the DLA column density in this example are much tighter than those forecast for either the bubble size or the neutral fraction. This results because the derivative with respect to the DLA column density parameter is much stronger than the bubble size and neutral fraction parameter derivatives here, as previously illustrated in the top panel of Figure 6. 11 The strongest degeneracy direction here is between ionized skewer size and neutral fraction, since one can trade-off for an increasing neutral fraction by boosting the bubble size. Similarly, the column density and bubble size parameters are positively correlated. Naively, the positive correlation between DLA column density and neutral fraction seen in the upper left hand panel of the figure is surprising, but this results because the bubble size is also an important parameter in the problem: one can compensate for an increase in both the neutral fraction and DLA column density by increasing the bubble size. Indeed, if we fix the bubble size parameter entirely we find the opposite degeneracy direction between q NHI and q xHI , as expected. In spite of the weak effect from the IGM damping wing in this case, one would still expect an interesting upper bound on the neutral fraction after accumulating the 11 GRB afterglow spectra in this redshift bin, as illustrated in the bottom right hand corner of Figure 7. This panel shows the 1-d likelihood for x HI after marginalizing over the DLA column density and bubble size parameters. This estimate needs to be refined, however, to incorporate the DLA and skewer size distributions. As an aside, we note that metal absorption lines -detected slightly redward of the damping wing feature in the same observations -may help in breaking the degeneracies shown. In particular, low-ionization metal lines should trace the same cold neutral gas in the ISM of the GRB host galaxy that produces the DLA contamination. Specifically, Figure 11 of Mas-Ribas et al. (2017) shows a strong trend between the equivalent width of lowionization metal absorption lines and DLA column densities, as determined from stacking BOSS quasar spectra. Therefore measurements of the equivalent widths of the low-ionization metal lines may provide independent constraints on the DLA column density parameter towards each GRB afterglow. Although we do not consider this further in what follows, it is an interesting direction for future work. Figures 8 and 9 show analogous plots at z s = 8 and z s = 10, respectively. The main change relative to Figure 7 is that the neutral fraction parameter now has a stronger impact than either of the DLA and bubble size parameters. This reflects the increasing neutral fraction and smaller bubble sizes in our model at these redshifts, which make the IGM damping wing more prominent (see the middle and bottom panel of Figure 6). Even though we expect fewer GRB afterglows in these redshift bins than at z s ∼ 6 (for our fiducial GRB redshift distribution there are 3 at z s = 8 and 1 at z s = 10, Figure 1), we expect tighter constraints on the neutral fraction. Indeed, the neutral fraction constraints in this simplified case (bottom right hand panels in Figures 8 and 9) look promising. In principle, the bubble size parameter itself along with its distribution and redshift evolution, provides valuable information about reionization. However, the error bar forecasts on this parameter are fractionally large (Figures 7-9), and so we do not consider the reionization information contained in this parameter further in what follows. Finally, the black points and error bars show the constraints for the Theseus-inspired ∼ 80z ≥ 6 GRBs sample. Each of these results incorporate the ionized bubble size distribution and DLA column density PDFs discussed in the text. Each afterglow spectrum is assumed to have an SNR = 20 per R = 3, 000 resolution element at the continuum. Note that in the red/blue models the error bars are identical in the zs = 6 redshift bin since these cases give an identical GRB count in this bin.

FULL IONIZATION HISTORY FORECASTS
We next compute full ionization history forecasts, assuming the bubble size distribution of Eqs. 6-8 and the DLA column density PDF of Figure 5. The required computations here are described by Eqs. 25 and 26 of Appendix A, and amount to averaging the inverse variance of the neutral fraction over the L, N HI PDFs.
The results of these calculations are shown in Figure 10. These appear quite encouraging: provided the spectroscopic follow-up observations we considered are feasible, a future GRB mission like the Gamow Explorer should help to determine the redshift evolution of the neutral fraction during cosmic reionization in some detail. Fractionally, the error bar on the neutral fraction is largest towards the end of reionization: in the z s = 6 bin we forecast a fractional error of 70% at 1-σ. In the other redshift bins, the error bars are in the 10 − 15% range in our fiducial model. The absolute errors range from σ xHI = 0.08 − 0.14 with the largest values in the first and last redshift bin. The error bars are still smaller in the case of the optimistic GRB redshift distribution model, and two additional redshift bins are populated in this case -out to z s = 12 -at which point the fiducial neutral fraction is around x HI = 0.98.
Finally, the black points and error bars show the more futuristic case in which the Theseus mission obtains high signal-to-noise ratio follow-up spectra for 80 z ≥ 6 GRB afterglows. In this scenario, the statistical precision of the error bars is further improved, with the error bar forecasts approaching 5% level constraints across many redshift bins in the EoR, degrading mainly towards low redshift (where the damping wing feature is typically somewhat feeble), and at very high redshift where there are few GRBs. Nevertheless, a project like this could potentially determine the neutral fraction from z s ∼ 6 − 14 with good statistical precision. The bi-modal nature arises because the N HI PDF is peaked at high column densities (see Figure 5), but the sightlines with low column density systems contain more information about the neutral fraction. The histograms are somewhat noisy because they rely on empirical estimates of the N HI PDF. Bottom panel: The same for the bubble size. The vertical lines indicate the average bubble size at each redshift. The Fisher information here is dominated by sightlines with bubble size smaller than the cosmic mean, with the results turning over as small bubbles bubbles become rare. Here Fx,x in the legend is shorthand for F x HI , x HI .

Contributions to the Fisher information
In order to unpack these results a little further it helps to examine the integrand of the Fisher information for the neutral fraction in Eq. 25. Specifically, Figure 11 plots both dlnF xHI , xHI /dlnL (bottom panel) and dlnF xHI , xHI /dlog 10 N HI (top panel) as a function of L and log 10 N HI , respectively. Put differently, the figure shows how much of the neutral fraction inverse variance is contributed by sightlines with various bubble sizes and column densities.
The bottom panel illustrates that ionized regions smaller than the mean skewer size dominate the Fisher information at each redshift. The information grows towards small bubble size because the damping wing is stronger for small bubbles, and then turns over as sufficiently small ionized regions become rare. This scale grows as reionization proceeds.
Likewise, the top panel shows that the Fisher information peaks for rare sightlines with low column density DLAs, although there is still appreciable information in the bins with 20 log 10 (N HI /cm −2 ) 22. This reflects a trade-off between the fact that sightlines with lower column densities contain more information about the neutral fraction, but are rarer than the high column density sightlines nearer to the peak of the fiducial N HI distribution ( Figure 5), which tend to swamp the IGM damping wing signature. The histograms in Figure 11 are somewhat noisy, because we rely on an empirical model for the N HI PDF. For instance, the dip just below log 10 (N HI /cm −2 ) = 20 reflects a corresponding drop in the observed N HI PDF ( Figure 5).
Another important implication of this figure is that the constraint forecasts may improve appreciably if the DLA column density PDF evolves strongly with redshift, shifting towards lower column densities, as may be expected if this distribution is directly related to the ionizing photon escape fraction (see §2.4). For example, the actual Fisher information in the column density bin around log 10 (N HI /cm −2 ) = 18.6 is a factor of ∼ 7 larger than the information in the bin centered on log 10 (N HI /cm −2 ) = 20.6; they only make comparable contributions to dlnF xHI , xHI /dlog 10 N HI because the lower column density bin here is on the rare tail of the assumed N HI PDF. Note that a factor of ∼ 7 gain in the Fisher information translates into a smaller error bar on the neutral fraction by a factor of ∼ √ 7.8 = 2.8, and so this could be significant.

Dependence on Spectral SNR
Here we consider the prospects for obtaining spectral SNRs of ∼ 20 per R = 3, 000 resolution element at the continuum, as adopted in the calculations of Figure 10. We caution that our goal here is only to obtain rough estimates and to explore the dependence on source redshift, flux density, sky background, and observing time. More detailed calculations will be required to cement the observing strategy and improve on our estimates, which are likely accurate only to within a factor of ∼ 2 in SNR. We also quantify the impact of varying the SNR around the fiducial value of SNR=20. To estimate the expected SNR of an observation, we first calculate the number of photons received by a telescope in each spectral resolution element from the source and/or the night sky background: Here N c denotes the number of photon counts from the source at the continuum, while f c ν is the source flux density at the continuum and at an observed wavelength of λ obs = λ α (1 + z s ). 12 The quantities N sky and f sky ν denote, respectively, the number of photons and the flux density received from the sky background. In this equation, we approximate the telescope aperture as circular, so that the collecting area is πD 2 tel /4 where D tel is the diameter of the telescope, is an efficiency factor to account for photon losses, ∆λ/λ is the size of the spectral pixels, t obs is the observing time, and h p is Planck's constant.We further denote N sky = αN c , noting that α 1 for ground-based observations.
In order to estimate the sky background for groundbased observations, we use the VLT X-Shooter exposure time calculator which returns the specific intensity of the sky background in different frequency bands 13 Specifically, this gives (59, 340, 650) µJy/arcsec 2 in the (I,Y,J) bands which, respectively, include Ly-α from source redshifts of z s = (6, 7 − 8, 9 − 10). These values then determine the flux density of the sky background, further assuming one-arcsecond spatial pixels for ground-based observations. It should be possible to improve on this seeing-limited case using adaptive optics, but we do not consider this further in what follows. The sky background is dominated by atmospheric emission, or "airglow", and is a strong function of wavelength; this makes the highest redshift GRB measurements challenging from the ground. From space, with the JWST, the sky background is expected to be dominated by Zodiacal light (see below) and such observations will be in the sourcedominated regime, i.e., N c N sky , α ∼ 0. Using Eq. 13 we can estimate the spectral SNR (squared), accounting for Poisson fluctuations in the source and sky background photon counts as: We neglect read-out noise and dark current. Plugging in some characteristic numbers, we find: That is, we find that an 8-meter telescope can achieve an SNR of 20 on a 10 µJy GRB afterglow in t obs = 6.1 hours of observing time at a spectral resolution of R=3,000.
Here we use α = 34 appropriate for ground-based observations in the Y-band of a source at z s ∼ 8, where f sky ν = 340 µJy. 14 In other words, an SNR of 20 could be achieved on a 10µJy afterglow at z s = 8 from the ground on an 8-meter telescope in a single night of observations. Although this is a fairly bright burst, note that the Totani et al. (2006) spectrum, for example, was observed 3.4 days after the burst, and so could have been detected as a 10µJy afterglow if it had been followed-up more promptly, a little less than a day after the burst (assuming a typical ∼ t −1 decay). This can also be seen directly in the light curve of this source presented in the discovery paper by Haislip et al. (2006). 13 https://www.eso.org/observing/etc/bin/gen/form?INS. NAME=X-SHOOTER+INS.MODE=spectro. We use the default fixed sky model parameters, a Moon Phase parameter FLI = 0.5, and Airmass = 1.5. We further assume that the sky background is comparable at other good observing sites.
14 We tested the accuracy of Eq. 15 by comparing with the afterglow spectrum in Totani et al. (2006). This was a 2.6µJy burst at z = 6.3, observed for 4 hours. Adopting a sky background estimate of ∼ 80 µJy, in between the I and J band values above, we estimate an SNR of 8 from Eq. 15 which appears consistent with the published spectrum at the continuum near Ly-α.
Moreover, in the future it will be possible to exploit the power of 30-meter class telescopes on the ground and the JWST in space. For example, Eq. 15 shows that using a 30-meter telescope one could observe a burst comparable to the Totani et al. (2006) one, with a flux of f c ν = 3 µJy, at an SNR of 20 in 4.7 hours, assuming the z s = 8 sky background numbers. As mentioned earlier, adaptive optics might help to reduce the sky bacgkround contribution and the required observing time. Using the JWST (with D tel = 6.5m), a f c ν = 2 µJy burst observation would reach an SNR of 20 after 1.3 hours of observing time. 15 Since the JWST observations are sourcecount limited -and not by the wavelength-dependent sky background -the SNR here only depends on the source flux density, with time dilation also counteracting the dimming from increasing luminosity distance towards higher redshifts, as mentioned earlier. JWST will therefore play a crucial role in obtaining high SNR afterglow spectra of the most distant GRBs at z 9, where airglow makes ground-based observation especially challenging, although the prospects here with 30 meter telescopes towards relatively bright afterglows are also fairly promising. Table 1 provides a further summary of the exposure time requirements for obtaining an SNR of 20. Here we show examples at redshifts z = 6, 8 and 10, for afterglows with specific flux densities of f c ν = 1, 5 and 10 µJy. We show results for 8 and 30 meter telescopes from the ground, and for space-based observations with the JWST.
While these estimates suggest there are multiple paths forward for obtaining the requisite SNR ∼ 20 follow-up spectra, further work is needed to refine our SNR estimates, and to assess what fraction of the GRBs detectable with the Gamow Explorer, HiZ-GUNDAM, and the Theseus mission will meet these demands. This depends on the flux distribution of observable GRBs and on how rapidly the follow-up spectroscopic observations can be executed.
We can also illustrate how our Fisher matrix forecasts depend on the quality of the afterglow follow-up observations. This is shown in Figure 12 and demonstrates the improvements that may be possible if the SNR is boosted to 30. It also shows how much the constraints degrade if the typical follow-up spectrum has SNR=10 per R = 3, 000 resolution element.

SOURCES
In this section we briefly illustrate how the combination of GRB afterglow constraints on both the escape fraction and the ionization history can help assess whether we have an accurate census of the ionizing source populations. The redshift distribution of long-duration GRBs also offers a handle on the evolution of the SFRD -provided evolution in the GRB to star-formation rate is understood and/or empirically well calibrated. Essentially, 15 For the case of JWST, we consider NIRSPEC observations with a pixel size of 0.01 arcsec 2 and the highest spectral resolution mode with R = 2, 700.
We use the sky background data from https://jwst-docs.stsci.edu/ jwst-observatory-functionality/jwst-background-model. We assume the same efficiency factor, = 0.25, as for the groundbased observations. -Forecasted constraints on the reionization history of the universe vs spectral signal-to-noise-ratio (SNR). Identical to Figure 10, except here we vary the SNR of the spectral pixels from SNR=10 (magenta error bars) to SNR=20 (blue points, identical to the blue points in Figure 2), and SNR=30 (red points). In all cases, we assume a spectral resolution of R = 3, 000 and the fiducial GRB redshift distribution. the source term in Eq. 5 is determined (at each redshift) by the SFRD, the escape fraction, and the ionizing spectrum of the sources. The recombination term in Eq. 5 is generally less important, since the recombination time is fairly long near the cosmic mean density and current simulation work suggests a fairly small clumping factor, C ∼ 2 − 3 (e.g., McQuinn et al. 2011), during most of the EoR. Hence if the ionizing spectrum may be constrained with e.g., JWST observations, GRB afterglow data may allow one to predict the ionization history from the observed SFRD and escape fraction; this may be compared with the damping wing measurements of the reionization history. 16 The SFRD determinations may of course be compared with that from JWST observations, which probe the bright end of the galaxy populations. Lineintensity mapping offers another potential handle on the cumulative impact of faint source populations (Sun et al. 2020;Kovetz et al. 2017). Taken together, these observa-tions should provide a valuable test of whether we have an accurate census of the sources that reionized the universe.
The escape fraction determinations clearly play a key role here, because this quantity is highly uncertain both theoretically and empirically, and there is relatively little hope (other than through GRB afterglows) for direct constraints during the EoR, especially in the small mass galaxies that may dominate the reionization process. A fairly large sample of afterglows with follow-up spectroscopy is, however, required to determine the escape fraction. The data required motivates follow-up spectroscopy of afterglows at slightly lower redshift, e.g., near z ∼ 5: the higher SFRD rate and slightly less stringent requirements for detecting these more nearby sources will yield larger samples than the 20 z s ≥ 6 afterglows expected for the Gamow Explorer.
This post-reionization sample avoids any damping wing contribution from the IGM; the escape fraction is then determined by fitting for the neutral hydrogen column density associated with each GRB host galaxy, and averaging over the sample of sightlines (e.g., footnote 4, Chen et al. 2007;Tanvir et al. 2019). The likely scenario here is that a fraction 1−f esc of the sightlines have DLAs and hence effectively vanishing ionizing photon leakage, while f esc of the afterglows show negligible neutral hydrogen absorption, i.e., no DLAs or even Lyman-limit systems. In this case, the requirements for follow-up spectroscopy are likely less stringent than in the case of measuring the full shape of the IGM damping wing. In a first pass, one mainly wants to determine whether or not a given sightline has a DLA or not. For estimating the escape fraction, it is not necessary to determine the precise DLA column density since sightlines with DLAs give negligible escape fractions in any case. In the case of sightlines that do not show DLAs, further follow-up may be required to determine whether there is a Lyman-limit system and the associated column density.

Prospects for Measuring the Escape Fraction
For example, suppose that the escape fraction is f esc = 0.1 and that roughly 100% of ionizing photons escape along a fraction f esc of sightlines while no ionizing photons escape from the remaining 1 − f esc of afterglows. In this case, on average one expects N esc = f esc N GRB sightlines to show escaping ionizing photons in a sample of N GRB afterglows. Assuming Poisson fluctuations in N esc around this average, the fractional error on the escape fraction from the sample is σ fesc /f esc ∼ 1/ √ N esc ∼ 1/ √ f esc N GRB . That is, in a sample of N GRB = 100 GRB afterglows, we expect N esc = 10 sightlines to show escaping ionizing photons, and ∼ 30% errors on the escape fraction estimate. Here we assume that the errors are dominated by Poisson fluctuations in the number of sightlines with escaping ionizing photons, and ignore additional errors that may arise from imperfect column density determinations, for example. Note that the escape fraction assumed here implies strong redshift evolution from the Tanvir et al. (2019) sample, but this may be required to reionize the universe with star-forming galaxies (e.g., Bouwens et al. 2016), as discussed in §2.4. It would be interesting to split future samples by redshift and stellar mass (obtained from follow-up host galaxy observations) to explore any correlations or trends in these quantities. Fig. 13.-Illustration of how joint constraints on the ionization history and escape fraction can improve our understanding of ionizing source populations. The red and blue points with error bars are the forecasted constraints on the reionization history from Figure 10. The black solid line shows the fiducial model (with a minimum host halo mass of M min = 10 9 M ). The black dotted line shows the impact of increasing the minimum host halo mass in the model to M min = 5 × 10 9 M , while fixing the escape fraction. The black dashed line also adopts M min = 5×10 9 M , but increases the escape fraction by a factor of two, which largely compensates for the delay in reionization with this M min . A measurement of the escape fraction can help rule out this case, which would otherwise be allowed by the ionization history measurements and determinations of the bright end of the UV luminosity function. Even a modest ∼ 30% level constraint on the escape fraction would be interesting and helpful. In order to illustrate this, Figure 13 considers an alternate model for the ionizing sources, in which our fiducial minimum galaxy hosting halo mass is boosted from M min = 10 9 M to M min = 5 × 10 9 M . This is intended to loosely mimic scenarios in which there are no additional faint galaxy populations below the threshold of UV luminosity function measurements: that is, it is a case where bright galaxy populations produce most of the ionizing photons, termed "reionization by oligarchs" by Naidu et al. (2019). Although different than our fiducial model, this type of scenario is also in agreement with current data, as discussed in that work. Adopting the same value of the escape fraction, or more precisely the same ionizing efficiency factor in Eq. 5, gives the dotted black curve in Figure 13. This delays the reionization history of the universe relative to our fiducial model and would be ruledout by the damping wing measurements given the forecasts of the previous section. In fact, this case is already moderately disfavored by current observations (e.g., Figure 2), but suffices for illustration.
However, one can compensate for this delay by increasing the escape fraction. In particular, the dashed line in Figure 13 boosts the escape fraction by a factor of two and this largely counteracts the effects from increasing the minimum mass. 17 Put differently, there are degeneracies between the minimum host halo mass and the escape fraction (see also, for example, Figure 9 of Sun & Furlanetto 2016). This can be broken by adding the GRB determinations of the escape fraction, which -in our fiducial model -would argue against the factor of two boost in escape fraction. Without the escape fraction measurements, one could be spuriously led into thinking that the brighter sources reionize the universe. Hence the combination of ionization history constraints and escape fraction measurements can play an important role in determining whether we have accurately identified the source populations responsible for reionization. Furthermore, the comparison between the redshift distribution of the GRBs and the SFRD determined from bright populations with e.g., the JWST might provide independent evidence for important contributions from faint source populations, since long-duration GRBs may be produced in such galaxies, while these may remain below the JWST flux limits.

Source Parameter Constraints
We can sharpen this argument by further forecasting the error bars on the parameters of the ionizing sources, which may be obtained from both the GRB afterglow constraints on the ionization history and the escape fraction determinations. To do this, we perform a Fisher matrix forecast for the ionizing source parameters in Eq. 5, given the error bars previously forecast for x HI (z) . Here, for simplicity, we consider a twoparameter model, characterized by q ζ = (ζ − ζ fid )/ζ fid and q Mmin = (M min − M min,fid )/M min,fid . That is, we fix the clumping factor at C = 3 and ignore any redshift or mass dependence in these parameters. Our fiducial parameter values are ζ fid = 20 and M min = 10 9 M . The Fisher matrix for these parameters may be computed as where the prime notation distinguishes this Fisher matrix from the one in Eq. 11. The sum runs over the redshift bins. The parameter derivatives are computed from the ionization equation (Eq. 5), and the variance of the neutral fraction estimates in each bin follow from the calculations in §4. Finally, we can add the independent information provided by an escape fraction estimate. Note that in principle the ionizing efficiency parameter depends on the overall star-formation efficiency -Constraints on the ionizing source properties from the reionization history and escape fraction measurements. The blue, red, and black contours show 2 − σ constraints on the ionizing efficiency parameter, ζ, and the minimum host halo mass of the ionizing sources, M min , from the ionization history constraints for the fiducial, optimistic, and Theseus-inspired GRB samples, respectively. The black dashed lines show 1 and 2 − σ constraints on the escape fraction from a sample of 100 afterglows, assuming that the uncertainty in the ionizing efficiency parameter is dominated by that in the escape fraction (see text). The combination of these measurements is helpful for determining the M min parameter. and the ionizing spectrum, as well as the escape fraction. For simplicity, we suppose that in the JWST-era the escape fraction is the main unknown, so that a ∼ 30% error in the escape fraction translates into a ∼ 30% error in the ionizing efficiency. In this case, assuming a fiducial escape fraction of 10% and a sample of 100 GRBs determines ζ with a fractional accuracy of ∼ 1/ √ 10 as discussed earlier.
The resulting constraint ellipses are shown in Figure 14 for each of the GRB redshift distributions in Figure 1, which result in the corresponding neutral fraction error bars in Figure 10. 18 As anticipated previously, the ionization history measurements alone leave a strong positive degeneracy between q ζ and q Mmin . However, the f esc measurement can help break this degeneracy. Quantitatively, in our baseline N GRB (z ≥ 6) = 20 model, we find that the fractional M min error is 74% for a joint constraint but only 200% without the information from f esc . For the more optimistic GRB redshift distribution with N GRB (z ≥ 6) = 31, the corresponding numbers are 69% and 150%. In the Theseus-inspired scenario, with N GRB (z ≥ 6) = 80, these numbers improve to 59% and 92%, respectively. 19 A good deal of the constraining power is driven by the bursts at high redshift: for instance, if we only include the bursts at z ≥ 9 in the Theseus case, the error bars on the minimum mass only increase by 8%. The high redshift bursts play a valuable role here because the halo mass function becomes exponentially sensitive to M min , as the galaxy hosting halos form only in the rare high-σ peaks of the mass distribution in these highest redshift bins. The numbers quoted are 1 − σ constraints. Hence the escape fraction limits do help noticeably, and they may be still more valuable if f esc is larger than our fiducial value of f esc = 0.1, which would improve the escape fraction measurement error.
Further, we have adopted a very simple model for the ionizing sources in this illustration. In more realistic cases the parameters may depend on halo mass, redshift, or other galaxy host properties. It will be important to combine ionization history constraints with a full suite of measurements to help pin down the properties of the ionizing sources. Determinations of the escape fraction, and any measurements or limits on the redshift dependence of this quantity and/or correlations with host properties should be extremely valuable here. Although our simple illustration shows improvements on the M min constraint from an f esc determination, this may in fact understate the value of such measurements.

CONCLUSIONS
We have forecast the constraints on the reionization history for upcoming GRB afterglow samples during the EoR. Assuming 20 z ≥ 6 GRB afterglow spectra with a typical SNR of 20 per R = 3, 000 spectral resolution element at the continuum, we find that the volume averaged neutral fraction can be determined to 10-15% accuracy from z ∼ 7 − 10, and with ∼ 70% errors near z = 6. Although there appear to be reasonable prospects for obtaining the requisite follow-up observations (see §4.2), more work is required to assess the number of GRBs that will meet these requirements in different redshift bins.
The ionization history forecasts represent a significant improvement on current constraints (e.g., Figure 2), and will provide a useful cross-check on future damping wing measurements towards high redshift quasar samples. Perhaps more importantly, GRBs can help extend the quasar measurements past z 8 − 9 to probe earlier phases of the reionization history before bright quasars formed. Measurements of the redshifted 21 cm power spectrum and its evolution may offer more precise measurements of the reionization history (DeBoer et al. 2017) than our GRB forecasts, but these rely on understanding the mapping between the observable, spatial fluctuations in the 21 cm signal and the average ionization history, and face a range of systematic challenges (Liu & Shaw 2020). While the GRB constraints also rely on understanding the mapping between an observable (the damping wing) and the average ionization history, they provide an independent constraint from the 21 cm efforts, and will be valuable. In the future, it will be interesting to consider joint fits with 21 cm, future cosmic variance limited CMB EE polarization measurements (which provide mostly an integral constraint on the ionization history, e.g., Heinrich et al. 2017), the patchy kinetic Sunyaev-Zel'dovich effect (Ferraro & Smith 2018), the quasar constraints, Lyman-alpha emitter samples (e.g., Mason et al. 2019), and other probes. In the slightly longer term, a sample of as many ∼ 80 z ≥ 6 GRB afterglows from the Theseus mission (Amati et al. 2018), with follow-up spec-troscopy, could return constraints as tight as ∼ 5% across much of the EoR and place some limits out to z ∼ 14. In principle, these measurements could be used to infer the electron scattering optical depth to CMB photons, and thereby sharpen cosmological parameter constraints (e.g. Liu et al. 2016).
We emphasized that an appealing aspect of GRB afterglow probes of reionization is that these constraints may also be combined with escape fraction determinations from the GRB column density distribution and with internal or external measurements of the SFRD. In order, however, to build up a sufficient afterglow sample to determine the escape fraction we recommend a follow-up campaign targeting afterglows at slightly lower redshifts, z ∼ 5 and 50 − 100 afterglows. Taken together, the ionization history, escape fraction, and SFRD constraints will provide a powerful test of our census of the sources that reionized the universe.
It will also be interesting to perform follow-up observations to compare the large-scale environment around GRBs with strong and weak damping wings. Surveys of the galaxy distribution and/or line-intensity maps in various emission lines can potentially identify some of the surrounding galaxies responsible for ionizing the intergalactic gas near GRBs that show weak damping wing absorption. It may also be possible to construct maps of the surrounding neutral hydrogen using 21 cm surveys (Beardsley et al. 2015), which would presumably be brighter in the proximity of GRBs showing strong damping wing absorption. The combination of these observations should help in understanding the interplay between the ionizing sources and the neutral gas in the IGM during reionization.
The biggest caveat with our current work is that we assume a simplified three parameter model for the damping wing signature towards each GRB afterglow. Although patchy reionization is partly incorporated into this model, we assume a uniform neutral fraction outside of the ionized bubble around the GRB host. This description does not fully account for the -generally complicated -distribution of neutral and ionized patches towards each afterglow (McQuinn et al. 2008;Mesinger & Furlanetto 2008). While we expect our simplified description to be mostly adequate for forecasting the constraints that may be obtained with future missions, it is insufficient for interpreting upcoming measurements. The model should be least accurate in the late stages of reionization when the ionized structures become large compared to the effective damping wing length scale. Note that while neglecting scatter in the exterior neutral fraction may cause us to underestimate the error forecasts on this quantity, it also implies that the damping wing feature will be stronger and more detectable towards some lines of sight than in our model; the precise impact of this approximation is unclear and likely sensitive to the reionization model considered.
We encourage further simulation studies to improve on this treatment. One approach may be to develop a generalized version of the three parameter model informed by mock damping wing spectra extracted from reionization simulations. In this case, a description might include the extent of the first neutral region, the size of the second ionized patch, and so forth. In this case, the average neutral fraction may no longer even be a parameter in the model. Instead, one considers only the length of each neutral segment and their distance to the source. These could be fit for independently of any particular model, while the simulations could be used in a second more model-dependent step to calibrate the relationship of these parameters to the volume averaged neutral fraction. It may also be possible to calculate the average transmission profile around GRB hosts of a given halo mass using a halo-model (Cooray & Sheth 2002) type calculation (see Appendix C of Malloy & Lidz 2015 for work in this direction).
The simulations will also enable mock end-to-end tests of parameter extraction pipelines. Such efforts can test the Fisher matrix approximation adopted in this work, which assumes a quadratic expansion in the logarithm of the likelihood around a fiducial input model. A suite of mock spectra may further help in understanding the sample-to-sample variations in the expected error bars; our Fisher forecasts assume typical samples and do not consider the error on the error bar.
Another important issue concerns the redshift evolution of the DLA column density distribution. Our fiducial model adopts the current empirical distribution of DLA host column densities (Tanvir et al. 2019), but this may be a pessimistic assumption. It might be interesting to use numerical simulations of galaxy formation, capable of modeling the escape fraction of ionizing photons (e.g., Ma et al. 2020), to try and predict the redshift evolution in the DLA column density distribution towards GRB hosts. Interesting in its own right, this could also be used in our forecasts in place of the empirical distribution, which mostly consists of lower redshift GRBs.
In summary, the prospects appear promising for future GRB missions to extract valuable constraints on the reionization history. There is a good deal of work to be done to better understand the expected damping wing signatures and to best prepare for the upcoming missions, but we expect these to play an important role in our understanding of cosmic reionization. and this results in Eq. 11 in the main body of the text.
It is further useful to define the likelihood of a neutral fraction estimate, denotedx here for brevity of notation, for a given realization of bubble size and DLA column density, yet marginalized over the estimates of bubble size and column density. We denote this quantity as L(x|L, N HI ). We can compute the width of this distribution, around some fiducial volume-averaged neutral fraction, by determining the inverse of the conditional Fisher matrix. Specifically, where the factor of x HI 2 in the second equality is needed to move from the variance of q xHI (see §3) to that of x HI itself.
Next, we can finally incorporate the bubble size and column density distributions. The likelihood ofx, marginalized over the bubble size and column density PDFs is: The non-conditional Fisher matrix of this one-parameter distribution is given by: We take the second derivative of Eq. 22, using that ∂L(x|L, N HI )/∂x|x = xHI = ∂L(x)/∂x|x = xHI = 0. This gives: where the integration overx, weighted by L(x), takes care of the ensemble averaging. We can rearrange the order of integration and rewrite this as: .
Although the notation is a little cumbersome, we arrive at a simple result. The non-conditional Fisher matrix, F xHI , xHI -i.e., the inverse variance of the neutral fraction estimate -is a weighted-average of the conditional neutral fraction inverse variance, where the weighting is over the bubble size and DLA column density PDFs. That is, our final expression is: In summary, our procedure is to first compute the conditional Fisher matrices of Eq. 11 and Eq. 18 for many different values of L, N HI , spanning the full range of their PDFs. We then invert each conditional Fisher matrix, using Eq. 19 to determine the inverse variance of the neutral fraction estimate for that L, N HI . Finally, we average together the results using Eq. 25 to weight by the PDFs of L, N HI . The final error bar follows simply from this equation as: where the N GRB (z s ) factor accounts for the number of independent GRB sightlines in the redshift bin.

APPENDIX B
Here we show how sensitive our results are to variations around the some of the assumptions in our modeling. First we consider how the results depend on our treatment of the lower limit of integration in our damping wing computations (Eq. 1). Recall that throughout the text we took each GRB to reside at the center of its redshift bin (at redshift z s ), and computed the damping wing using the average neutral fraction at bin center, truncating the integral at the edge of each redshift bin, z s − 0.5. This neglects redshift evolution across the bin and, more importantly, the presence of neutral gas at redshifts below the lower redshift end of the bin in question. While this estimate is imperfect, it simplifies the calculations and is a reasonable approximation given that the damping wing is dominated by relatively nearby gas in the IGM. The blue points and error bars show our fiducial model assumption for the lower redshift limit on the damping wing integral in Eq. 1. In this case, z end = zs − 0.5: that is, we adopt the average neutral fraction at bin center and integrate down to the lower edge of the bin. The magenta points instead adopt the average neutral fraction at bin center but integrate all the way down to z end = 5.5. We adopt the fiducial GRB redshift distribution with 20 z ≥ 6 GRBs. Note that the two cases coincide for the zs = 6 bin. The results are not strongly sensitive to the treatment of z end .
As one simple test of the impact of this approximation, we simply adjust the lower limit of the redshift integration from z end − 0.5 to z end = 5.5, while keeping x HI fixed at its bin center value. This latter redshift value coincides with the edge of the lowest redshift bin in our analysis and is close to the end of reionization in the model. This case deliberately exaggerates the impact of neutral gas in the IGM below the edge of each redshift bin, especially in the highest redshift bins (while it gives identical results for the lowest redshift bin). We then follow the procedure in §4 to compute the error bars on the ionization history for this alternate choice of z end . The results of this test are shown in Figure 15; the error bars tighten slightly in this case especially in the highest redshift bin. However, the difference is relatively small and the case with z end = 5.5 deliberately exaggerates the effect of the neutral gas below the bin edge. Hence our results are not strongly impacted by the choice of z end . A second choice relates to the range of wavelengths included in our Fisher matrix calculations (Eq. 11). In the computations in the text we include the wavelength range between Ly-α at the source redshift, at observed wavelength λ obs = λ α (1 + z s ), and out to 1, 500Å redward of this wavelength (independent of z s ). In Figure 16 we consider variations around this choice. Specifically, for the case of our fiducial GRB redshift distribution we compute the marginalized fractional error on the neutral fraction in each redshift bin, including a range of different choices for the long-wavelength cut-off. In each bin the results converge provided 500Å redward of Ly-α at the source redshift are included (as might also be expected based on Figures 4 and 6), justifying our fiducial choice of 1, 500Å. Fig. 17.-Error bar forecasts for different sky background to source count ratios. In each case, the SNR at the continuum is fixed at SNR=20, but the blue points and error bars take the sky background dominated limit adopted in the main text, while the green error bars assume the sky background is 0.1 times the source flux density, and the magenta errors give the case where the sky background is 100 times smaller than the source. In each case, we assume the fiducial GRB redshift distribution. The results are not sensitive to the assumptions here, although the errors are slightly smaller in the source-dominated case.
Finally, the calculations in the text assume sky background limited measurements, which is appropriate for the case of ground-based follow up spectroscopic observations. As discussed in §4.2, space-based observations with the JWST will instead be in the source-dominated regime. At fixed SNR, the transmission variance in the source-dominated case will be slightly smaller in absorbed regions than in the sky-dominated regime (see Eq. 12). This slightly improves the ability to measure the transmission profile close to Ly-α at the source redshift, which helps in determining the bubble size and column density parameters. This then leads to less strong parameter degeneracies than in the skydominated case of Figures 7-9, and slight improvements on the marginalized neutral fraction error forecasts. Figure 17 compares the sky dominated case with alternate scenarios where the sky flux is 0.1 and 0.01 times the source flux. The improvements are relatively small, and so we expect rather similar constraints at fixed SNR in each of the sky and source dominated cases.