1 Introduction

A major goal of contemporary astrophysics and cosmology is to achieve a detailed understanding of the formation of the first collapsed objects (Pop-III and early Pop-II stars, black holes and the primordial galaxies in which they were born) during the first billion years in the life of the Universe. The growth of these early luminous sources is intimately connected to the reionization of the intergalactic medium (IGM) and build-up of heavy chemical elements [1]. Cosmic chemical evolution is very poorly constrained at high redshift, and even in the JWST era, metallicity estimates at z > 6 will rely on crude, or model dependent, emission line diagnostics for only a limited number of the brightest galaxies (typically, MUV < − 19 mag). Regarding reionization, measurement of the Thomson-scattering optical depth to the microwave background by the Planck satellite now suggests that it substantially occurred in the redshift range \(z \sim 7 \)– 8.5 [2], but the timeline and the nature of the driving sources remain contentious [3,4,5].

The timeline of reionization has been constrained chiefly from two directions. The first entails observations of the Gunn-Peterson trough [6] and estimates of the redshift evolution of the Lyman α (henceforth Ly α) optical depths in the spectra of distant quasars, which indicate that reionization was nearly complete at \(z\sim 6\). Due to its sensitivity, the Gunn-Peterson test, which shows total absorption of Ly α photons, therefore saturates for z > 6.1 quasars (see, e.g., [7,8,9]; also Section 8 for supporting evidence from GRB afterglows). A second approach has proven successful at probing the neutral fraction of the IGM at higher redshifts, namely searches for Ly α emission in spectroscopic observations of Lyman-break selected galaxies. These studies yield a drop in the fraction of Ly α emitting galaxies at z > 6 (e.g. [10,11,12,13,14,15,16,17,18,19,20,21,22,23]), although there exist exceptions such as in observations of UV-faint galaxies behind galaxy cluster lenses [24], massive galaxies [25], and possibly in cases of ionised bubbles around bright galaxies [26]. There is also some evidence for a late reionization scenario (e.g. [27]), with the process possibly still ongoing as late as \(z\sim 5.5\) [28,29,30]. Statistical measurements of the fluctuations in the redshifted 21 cm line of neutral hydrogen by future experiments, notably SKA-Low, are expected ultimately to provide further constraints on the time history of reionization [31].

The central question, however, remains whether it was predominantly radiation from massive stars that both brought about and sustained this phase change, as is conventionally thought [32], or whether more exotic mechanisms must be sought. Even in the JWST/ELT era, the numerous galaxies populating the faint-end of the luminosity function (\(M_{\text {UV}} \gtrsim -16\) mag) in the re-ionization epoch (z > 6) will be very hard to access. Gauging star formation occurring in these faint galaxies, the average fraction of ionizing radiation that escapes from galaxies, together with the parallel build-up of metals and their spreading throughout the galactic ISM and the IGM by stellar and galactic winds, will still remain highly challenging problems.

GRBs offer a unique route to detecting the death-throes of individual massive stars to very high redshifts (Fig. 1; [33,34,35,36,37]), which in turn provides multiple powerful probes of early star formation, metal enrichment and galaxy evolution, potentially even before the main phase of reionization (e.g. [38]). Indeed, they are detectable independently of the luminosity of their underlying hosts and so can pinpoint the presence of massive star formation in distant galaxies below the sensitivity limit of even the most powerful facilities foreseen in the long-term future. Additionally, their afterglow counterparts can be used as bright background lighthouses probing in absorption both the IGM and the interstellar medium (ISM) of faint galaxies otherwise not accessible with other observing techniques. The key to exploiting this new window on the early universe will be the detection and follow-up of samples of tens of GRBs in the era of reionization. The proposed THESEUS mission [39] is expected to identify and locate between 40 and 50 GRBs at z > 6 in 3.45 years of scientific operations [40], and allow on-board determination of photometric redshifts more accurate than 10% thanks to the identification of the Lyman break feature being shifted to the imaging sensitivity range of its Infrared Telescope (IRT, see [41]). In this white paper we outline the range of high redshift science that will be enabled by this mission.

Fig. 1
figure 1

A realisation of the expected redshift distribution of long-duration GRBs detected by THESEUS in 3.45 yr of scientific operations (orange and purple histograms) based on GRB population modeling [40], compared with the redshift distribution obtained to date (blue and black histograms)

Full size image

2 Global star formation rate from GRB rate as a function of redshift

Long-duration GRBs are produced by massive stars, and so track star formation, and in particular the populations of UV-bright stars responsible for the bulk of ionizing radiation production. This makes them important probes of global star formation, and significant changes in stellar populations, such as variations in the initial mass function (IMF).

However, there is growing observational evidence that the situation is more complicated in the low-z universe. Population studies based on well selected, complete sub-samples of long GRBs, have reported that their rate increases more rapidly with z and peaks at higher redshift compared to cosmic star formation (e.g. [42,43,44]). Moreover, analyses of the properties of galaxies hosting low-z GRB events have shown that they are typically less luminous, less massive and more metal poor than the normal field galaxy populations [45,46,47,48,49,50]. These findings are consistent with a scenario in which GRBs form preferentially from moderately metal-poor progenitors as indeed required by some theoretical models [51, 52], and inefficiently in super-solar metallicity environments. While the exact metallicity threshold at which GRB form efficiently is still debated, ranging from 0.3Z⊙ to nearly solar [53], there is general consensus that at high redshifts, where most of the galaxies are still metal-poor, GRBs can be used as fair tracers of the cosmic star formation (e.g. [54,55,56,57]). Recent hydrodynamical simulations accounting for complex chemical enrichment (based on [58]) including self-consistent dust production support this view at z ≥ 4 [59], and are consistent with the observed evolution of the dust-to-metals ratio in high-redshift galaxies probed by GRB-DLAs [60]. Under this assumption, their detected rate at z > 6 can be used to estimate the cosmic star formation in an independent way, provided that the proportionality factor between stars and GRB events is calibrated. Even more directly, the evolution of the GRB rate with z should be closely related to the decline of the cosmic star formation rate with increasing redshift, after the cosmic noon.

Intriguingly, analyses of this sort have consistently indicated a higher star formation rate (SFR) density at redshifts z > 6 than traditionally inferred from UV rest-frame galaxy studies ([49, 61]; see Fig. 2), which rely on counting star-forming galaxies and attempting to account for galaxies below the detection threshold as well as the fraction of UV light missed because of dust obscuration. Although this discrepancy has been alleviated by the growing realisation of the extremely steep faint-end slope of the galaxy Luminosity Function (LF) at z > 6 (e.g. [62], it still appears that this steep slope must continue to very faint magnitudes (e.g. MAB > − 11), in order to provide consistency with GRB counts and indeed to achieve reionization (something that can only be quantified via a full census of the GRB population). Interestingly, some numerical models of galaxy formation performed at different cosmic scales, and capable of resolving smaller galaxies clustering around Milky Way-like galaxies, support a similar view indicating that significant SFR is still hidden in undetected faint star forming systems clustering around over-dense ones at high redshift (see [65] Fig. 1). Alternative, and equally interesting, possibilities are that the stellar IMF becomes more top-heavy at early times (GRB progenitors being drawn from populations with birth masses \(\sim 25-40\) M⊙), suggestive of very low metallicity populations [66], or that there are multiple binary pathways to producing long-GRB progenitors with different metallicity dependence [67].

Fig. 2
figure 2

Comoving Star Formation Rate density as a function of redshift as derived from rest-frame UV surveys (shaded curve based on the observed galaxy population [63], blue hatched region accounts for galaxies below the detection limit) and from GRBs based on different assumptions (GRB rate to SFR ratio, GRB progenitor metallicity [64]). Red points show current constraints from the available GRB sample, which is very small at z > 6 [49], but shows some tension compared to estimates based on galaxy counts. Cyan points show corresponding estimates expected from a representative THESEUS sample where the GRB rate has been enhanced at high redshift to illustrate a situation in which the current tension is maintained

Full size image

However, up to now, these studies have been severely limited by the very small size of the high-z GRB samples [55, 61]. THESEUS will establish the GRB redshift distribution, N(z), much more reliably at z > 5.5 than previous missions, thanks to significantly larger numbers detected by the Soft X-ray Imager (SXI) and also more uniform selection via rapid on-board redshift estimates with the IRT, combined in many cases with a refinement of the redshift through follow-up with other facilities. Thus, accounting for the observational selection function, we will obtain the evolution of the global star formation rate density at redshifts where uncertainties from traditional galaxy observations begin to rise. This will firmly establish whether the tension between GRB counts and extrapolated galaxy counts remains, demanding a more radical reappraisal of high-z star formation (Fig. 2).

As another way to demonstrate the power of the THESEUS sample, in Fig. 3 we show the accuracy at which will measure the decline slope of the star formation rate density with redshift (specifically, \(-d\text {log(SFRD)}/d\log (1+z)\)) for the cosmic star formation at z > 6 after 1 and 3.45 years of operations. Already after the first year, the measured number of high-z GRBs will allow us to discard a slope as shallow as 2.4 at 3σ confidence level, while discarding slopes as steep as 3.3 will require the full 3.45-year sample. A slope as steep as 2.7 [68] will be measured with an error of 0.15 (1σ) while slightly tighter (looser) constraints are expected for shallower (steeper) declines.

Fig. 3
figure 3

Accuracy in determining the decline slope of the SFR density as measured by the number of GRBs detected by Theseus/SXI at z > 6. Green (red) line and points represent 1 year (3.45 years) constraints. The grey shaded area reports the 3σ confidence level: so in this example the actual slope is set to 2.9, and the first year results give a 3σ range of about 2.4 to in excess of 3.4, whereas the nominal mission would constrain the result to between about 2.6 and 3.3. (Note, the SFR density slope measured by [68] is \(\sim 2.7\))

Full size image

Thus THESEUS provides an independent way, complementary to deep-field galaxy surveys, to measure this important parameter at very high-z. Even more interestingly, the identification of any discrepancy between GRB-based measurements and those obtained by different methods could provide important insights about role of faint galaxies, missed even in future deep surveys, in shaping the cosmic star formation history, about the possible evolution of the stellar IMF in the very high-z Universe, and about the existence of an important population of bright X-ray transients, such as those expected from Pop-III stars (Section 9).

3 The galaxy luminosity function: detecting undetectable galaxies

The intrinsically very small galaxies, which appear to increasingly dominate star formation at z > 6, are very hard to detect directly. GRBs circumvent this difficulty, sign-posting the existence and redshifts of their hosts, no matter how faint.

The faint end of the galaxy luminosity function is a key issue for our understanding of reionization since, to the depth achieved in the Hubble Ultra-deep Field (HUDF), it appears that the faint-end of the LF steepens significantly with redshift approaching a power-law of slope \(\alpha \sim 2\) at z > 6 (i.e. where the number of galaxies per unit luminosity interval scales as ϕ(L) ∝ L−α for faint galaxies; e.g. [62]). Studies, which take advantage of gravitational lensing by galaxy clusters in the Hubble Frontier Fields, can reach even fainter magnitudes, albeit probing increasingly smaller volumes in the higher magnification regime, and are also consistent with a steep faint-end slope [69,70,71]. The value of the total luminosity integral depends sensitively on the choice of low-luminosity cut-off (and indeed the assumption of continued power-law form for the LF). By conducting deep searches for the hosts of GRBs at high-z we can directly quantify the ratio of star-formation occurring in detectable and undetectable galaxies, with the sole assumption that GRB-rate is proportional to star-formation rate (Fig. 4). Although currently limited by small-number statistics, early application of this technique has confirmed that the majority of star formation at z > 6 occurred in galaxies below the effective detection limit of HST [72, 73] with expected magnitudes fainter than \(m_{\text {AB}} \sim 30\), at the limit of what is reachable with JWST and the ELTs. Since the exact position and redshift of the galaxy is known from the GRB afterglow, follow-up observations to measure the host UV continuum are much more efficient than equivalent deep field searches for Lyman-break galaxies. If the luminosity function is modelled as a Schechter function with a sharp faint-end cut-off, then this analysis allows us to constrain that cut-off magnitude, even though the galaxies are too faint to be observed. This technique applied to a sample of ∼ 40 GRBs detected by THESEUS at z > 6 will yield much tighter constraints on that cut-off magnitude than obtained with GRBs so far (Fig. 8, left).

Fig. 4
figure 4

Mosaic of deep HST imaging of the locations of known GRBs at z > 6, obtained when the afterglows had faded (red boxes are centered on the GRB locations, and are 2 arcsec on a side). Only in 2–3 cases is the host galaxy detected, confirming that the bulk of high-z star formation is occurring in galaxies below current limits. This approach allows us to quantify the contribution of the faint end of the galaxy luminosity function to the star formation budget, even in the absence of direct detections

Full size image

4 The build-up of metals, molecules and dust

Bright GRB afterglows, with their intrinsic power-law spectra, provide ideal backlights for measuring not only the hydrogen column, but also obtaining exquisite abundances and gas kinematics, probing into the hearts of their host galaxies [74]. Thus, they can be used to map cosmic metal enrichment and chemical evolution to early times, and search for evidence of the nucleosynthetic products of even earlier generations of stars.

Based on the specifications of the IRT spectrograph with its spectral resolution \(R\sim 400\), we simulate a set of possible afterglows with varying host galaxy hydrogen column densities and metallicities, and analyse the precision that we will be able to achieve as a function of the signal-to-noise ratio of the spectra (Fig. 5, right hand panel). First, follow-up of the brightest afterglows (\(H_{\text {AB}} \lesssim 17.5\) mag) with the IRT spectroscopic mode will provide constraints within \(\sim 0.2\) dex on the hydrogen column density along the GRB line of sight in the host galaxy. The identification of associated metal absorption lines will also enable spectroscopic redshift determinations to < 1% precision, which will help refining the IRT photo-z estimates and distinguishing cleanly between the GRBs at z > 6 and contaminants from dusty afterglows at lower redshifts. Even though IRT spectra will be able to detect the usual rest-frame UV metal absorption lines [75] for afterglows with metallicities above \(\sim \)1% solar for an afterglow detected at \(H_{\text {AB}}\lesssim 16.5\) mag, accurate abundance, metallicity and abundance pattern determination necessitates follow-up deeper, higher-spectral resolution data from space- or ground-based facilities.

Fig. 5
figure 5

Left: Simulated ELT 30-minute spectrum of a z = 8.0 GRB afterglow with J(AB) = 20 (typical after \(\sim 0.5\) day). The S/N provides exquisite abundance determinations from metal absorption lines (in this example, 1% solar metallicity), while fitting the Ly α damping wing simultaneously fixes the IGM neutral fraction and the host Hi column density, as illustrated by the two overlaid models, a pure 100% neutral IGM (green) and a \(\log (N_{\text {HI}}/\text {cm}^{-2})=21.2\) host absorption with a fully ionized IGM (orange). A well-fitting combined model is shown in red. Right: The same afterglow with a magnitude of 16 in a simulated IRT spectrum with a realistic spectral observing sequence at a total integration time of 1800 s, illustrating that the most prominent metal lines are also clearly detected

Full size image

Taking advantage of the availability of 30 m class ground-based telescopes in the 2030s superb abundance determinations will be possible through simultaneous measurement of metal absorption lines and modelling the red-wing of Ly-α to determine host HI column density, potentially even many days post-burst (Fig. 5, left panel), thanks to the effects of cosmological time dilation for the highest redshift events. This may be supported by ATHENA observations to quantify the high ionization gas content (e.g. [76]) and potentially radio for other atomic fine structure and molecular lines (e.g. [77]).

Using the sample of GRBs discovered by THESEUS to trace the ISM in galaxies at z > 6 will be the only way to map in detail exact metallicities and abundance patterns across the whole range of star forming galaxies in the early Universe, including those at the very faint end of the LF (Fig. 6). Thus, for example, we will be able to search directly for evidence of alpha-enhancement, hinted at by modelling of emission line spectra at intermediate redshifts [78], and evidence of pop III chemical enrichment (see Section 9).

Fig. 6
figure 6

Absorption-line based metallicities relative to solar, corrected for dust depletion as a function of redshift for Quasar Damped Ly α absorbers (DLAs, grey symbols) and GRB-DLAs (blue symbols) (adapted from [97], [98]). Open square symbols show representative expectations for THESEUS, assuming continued evolution of the mass-metallicity relationship, and a dominant population of low mass galaxies at z > 6 (green triangles and red squares assume faint-end slopes of − 1.8 and − 1.4 for the galaxy luminosity function, respectively). GRBs represent the unique way for probing evolution of ISM absorption–based metallicities in the first billion years of cosmic history

Full size image

From the observed LF for galaxies at z > 6 [62], we can estimate the distributions of metallicities for GRBs at the highest redshifts by assuming that GRBs trace galaxy luminosity and that a relation between the stellar mass and gas-phase metallicity is valid. While this mass-metallicity relation has been established at z < 3 for galaxies in emission [79] and based on absorption line metallicities [80], a considerable uncertainty remains on the validity at higher redshifts, and we have to rely on an extrapolation that will be tested with future analyses of galaxies with the JWST. Since the low-mass end of the LF at the highest redshifts is debated, we vary the slope of the LF faint end, which impacts the resulting metallicity distributions. With a significant sample of GRB afterglow metallicities at z > 6 THESEUS will provide an independent test of the slope of the LF at z > 6. Most importantly, this test is entirely independent of being able to detect the host galaxy in emission.

There is a rising theoretical interest in both the cosmic early chemical enrichment and the properties of the ISM of high redshift galaxies, also supported by recent observations of dusty galaxies in the epoch of reionization [81, 82], and complemented by recent results from the ALPINE survey [83], down to \(z\sim 4\). ISM constraints provided by THESEUS will be critical to further improve both cosmological models investigating the early assembly of dusty galaxies [59, 84, 85] and the detailed evolution of the simulated ISM only accessible in zoom-in simulations of isolated galaxies [86, 87]. Many of the processes regulating dust evolution in the galactic environment are in fact strictly depending on the multi-phase ISM in which both metal atoms and dust grains co-evolve.

The imprint of the dust in the host can be seen in the (rest-frame UV/optical) broad-band spectral energy distributions of GRB afterglows. Such studies have found a variety of dust laws, many reasonably approximated by the three canonical local laws (SMC, LMC, Milky Way; e.g. [88, 89]), but along some sight-lines the extinction is unusual and harder to explain [90]. Thus, GRBs offer a remarkable route to assessing the dust content of even low mass galaxies in the early Universe [91]. At moderate redshifts, molecular hydrogen can be detected through electronic transitions (Lyman and Werner bands). Due to its specific formation and destruction processes, H2 provides an excellent probe of the cold neutral medium (CNM) in galaxies [92]. Being inclined to Jeans instability, this phase plays a crucial role in the star-formation process, but is not detectable in emission. In addition, the detection of various rotational levels of H2, together with fine-structure levels of companion species (as Ci) provides very sensitive tools to determine the prevailing physical conditions (in particular density, temperature, and UV field) [93, 94].

Although this becomes unfeasible at high redshift due to the Gunn-Peterson trough obscuring the Lyman-Werner bands, other absorption tracers, notably vibrationally excited H\(_{2}^{*}\), and neutral carbon, Ci, provide reliable indicators of high molecular hydrogen columns [95, 96].

We emphasize that using the THESEUS on-board NIR spectroscopy capabilities will provide, in addition to arcsec accurate location, the redshift estimates and luminosity measurements that are essential to optimising the time-critical follow-up observations using highly expensive next-generation facilities, allowing us to select the highest priority targets and deploy the most appropriate telescope and instrument.

Finally, we note that also in the X-ray band, SXI could allow for minimal studies of the medium surrounding the GRB site. As observed in GRB050904, X-ray data holds the potential to reveal a decrease in the X-ray absorbing column density and to constrain the metallicity [99]. SXI could follow GRB050904-like events up to \(\sim 3000\) s during which the full photoionisation of the medium will occur.

5 The Lyman-continuum escape fraction

GRB afterglow spectroscopy allows us to measure the column density of neutral hydrogen in the host galaxy, thus providing a powerful probe of the opacity of the interstellar medium to EUV photons.

A key issue for the reionization budget is quantifying the fraction of ionizing radiation that escapes the galaxies in which it is produced. Even in the JWST/ELT era, direct observations of escaping Lyman continuum radiation at z > 5 will remain unfeasible for the low mass galaxies responsible for the bulk of star formation in the era of reionization. However, GRB afterglow spectroscopy allows us to measure the column density of neutral hydrogen in the host galaxy, thus providing a strong lower limit to the opacity of the interstellar medium to FUV photons (further EUV attenuation due to dust can also be estimated by modelling the afterglow SEDs). A statistical sample of afterglows can be used to infer the average escape fraction over many lines of sight, specifically to the locations of massive stars dominating global ionizing radiation production. Useful constraints have so far only been possible at z = 2–5, indicating an upper limit of fesc < 2% (Fig. 7; [100, 101]), although note this is lower than found in the recent study of stacked LBG galaxy spectra at \(z\sim 3\) by [102]. However, future observations of the population of z > 5 GRBs detected by THESEUS, with both on-board spectroscopy and with 30 m class ground-based telescope follow-up observations, will provide much more precise constraints on the fraction of ionizing radiation that escaped galaxies during the epoch of reionization.

Fig. 7
figure 7

The host neutral hydrogen column densities measured from the Ly α absorption line in the afterglows of 143 GRBs spanning a wide range of redshift. Typically, the columns imply a high optical depth to ionizing Ly-continuum radiation, and thus a low overall escape fraction (fesc < 2%). Strikingly, there is little evolution in the running median value (red line) or inter-quartile range (pink region) over the redshift for which the distribution is well sampled, but there is a hint of a possible downturn at the highest redshifts. THESEUS will enable us to extend this plot with good statistics at z > 5 providing a direct estimate of the escape fraction on the lines of sight to massive stars in this era

Full size image

6 Did stars reionize the Universe?

The evolution of the IGM from a completely neutral to a fully ionized state is intimately linked to early structure formation, and thus a central issue for cosmology (e.g. [3]). Answering the key question of whether this phase change was primarily brought about stars hinges on two subsidiary issues: how much massive star formation was occurring as a function of redshift, and, on average, what proportion of the ionizing radiation produced by these massive stars escaped from the immediate environs of their host galaxies? Both will be addressed through THESEUS GRB observations.

The former problem, quantifying massive star formation as a function of redshift, can be extrapolated based on observed candidate z > 7 galaxies found in HST deep fields, but two very significant uncertainties are, firstly, the completeness and cleanness of the photometric redshift samples at z > 7, and, secondly, the poorly constrained faint-end behaviour of the galaxy luminosity function (at stellar masses \(\lesssim 10^{8}\) M⊙), especially since galaxies below the HST (and potentially even the JWST) detection limit very likely dominate the star-formation budget. Even though some constraints on fainter galaxies can be obtained through observations of lensing clusters [69], which will be improved further by JWST, simulations suggest that considerable star formation was likely occurring in fainter systems still [103].

As discussed in Section 5, the second problem, that of the Lyman-continuum escape fraction, is even more difficult since it cannot be determined directly at these redshifts, and studies at lower redshifts have found conflicting results. Recent stacked spectroscopic analyses have suggested escape fractions as high as 10% [102], which could be sufficient to drive reionization [64], but it is unclear whether the samples of galaxies studied, at \(z\sim 3\), are representative of all star- forming galaxies, and in particular of the typical, intrinsically fainter galaxies at z > 7. As seen in Fig. 7, GRB studies find that sight-lines to massive stars are generally highly opaque to ionizing radiation, at least up to \(z\sim 5\).

The improvements in both the census of star formation and the escape fraction from THESEUS GRB studies will provide a strong test of the hypothesis that reionization was brought about by star light. Our detailed simulations indicate that THESEUS is expected to detect between 40 and 80 GRBs at z > 6 over a 3.45-year mission, with between 10 and 25 of these at z > 8 (and several at z > 10) [40]. The on-board follow-up capability will mean that redshifts are estimated for almost all of these, and powerful next generation ground- and space-based telescopes available in this era will lead to extremely deep host searches and high-S/N afterglow spectroscopy for many (e.g. using ELT, ATHENA etc.). To illustrate the potential of such a sample, we simulate in Fig. 8 (right) the precision in constraining the product of the UV luminosity density and average escape fraction, ρUVfesc, that would be obtained with samples of 20, 30 and 50 GRBs at 7 < z < 9 having high-S/N afterglow spectroscopy and (∼3hr) ELT depth host searches (for definiteness the ρUV axis corresponds to z = 8). This will allow us to confidently distinguish between conventional models in which starlight brings about reionization, and models that fail (e.g. if the escape fraction remains as low as we find at lower redshifts).

Fig. 8
figure 8

Left: constraints on the faint end cut-off of the galaxy luminosity function (assuming a Schechter function with abrupt cut-off and bright end constrained by galaxy observations). Black curve shows current constraints based on nine GRBs at z ≥ 6; red curve shows a simulation after THESEUS, assuming deep searches for a sample of \(\sim 45\) hosts from 30 m class telescopes. Right: the UV luminosity density from stars at \(z\sim 8\) and average escape fraction 〈fesc〉 are insufficient to sustain reionization [104] unless the galaxy luminosity function steepens to magnitudes fainter than MUV = − 13 (hatched region), and/or 〈fesc〉 is much higher than that found from GRB studies at \(z\sim 2\)–5 (shaded region). Even in the late 2020s, 〈fesc〉 at z > 6 will be largely unconstrained by direct observations. The contours show the 2-σ expectations for samples of 20, 30 and 50 GRBs (thin to thick contours) at \(z\sim 7\)–9 for which deep spectroscopy provides the host neutral column and deep imaging constrains the fraction of star formation occurring in hosts below ELT limits. The cyan contours illustrate a model with high escape fraction (fesc ≈ 20%) and SFRD, whereas the blue contours are for a model with lower values of both (fesc ≈ 4%)

Full size image

7 Topology and timeline of reionization

In practice, it is expected that reionization should proceed in a patchy way, for example, ionized bubbles may grow initially around the highest density peaks where the first galaxies form, expanding and ultimately filling the whole IGM. The topology of the growing network of ionized regions reflects the character of the early structure formation and the ionizing radiation field. The Ly α scattering cross section is composed of a principal, exponential-profile component arising from thermal motion, and a tail component due to natural broadening. This latter damping wing absorption is weaker than the principal absorption component by several orders of magnitude, and since it follows a power-law profile, it ventures into longer wavelengths than the principal component. Effectively, this damping wing absorption is far less sensitive than Gunn-Peterson absorption, and can therefore be used to probe higher neutral fractions where the Gunn-Peterson test saturates. Via this damping wing absorption, therefore, we may constrain the neutral fraction of the IGM for values upward of 10%.

To date, this has been challenging both because of the sparsity of events, and the S/N ratio that is required, such that the experiment has been limited to a very small number of GRBs (e.g. GRB 050904 [105] and GRB 130606A [74, 106]). With high-S/N afterglow spectroscopy, the red damping wing of the hydrogen Ly α line can be decomposed into contributions due to the host galaxy and the IGM, with the latter providing a measure of the progress of reionization local to the burst.

GRBs themselves do not ionise the surrounding IGM, but previously occurring star formation in the host galaxy as well as contributions from nearby galaxies will form Hii regions. A patchy reionization scenario, which sees galaxies carving out Hii bubbles which later merge, would imply that measurements will vary across a sample of GRBs. By carrying out follow-up high-resolution spectroscopy of a sample of several tens of GRBs at z > 6–9, we can begin to investigate statistically the redshift-dependence of the average and variance of the reionization process [107, 108]. To illustrate the capability of follow-up spectroscopy to quantify both host Hi column, and IGM neutral fraction, Fig. 9 shows the results of simulating an afterglow with the characteristics of GRB 090423 (cf. Fig. 5) with a range of different input parameters, for which the measured parameters over a large number of realisations are shown (green points). This demonstrates that the IGM neutral fraction can be recovered from such spectra, even in the presence of a fairly high host column, and similarly the host column can generally be well characterised except when very low. While the situation may be complicated further due to the additional effect of any local ionized bubble in which the host resides, but again this can be modelled with sufficient S/N spectra [107], as is likely to be obtained for many afterglows with 30 m telescope spectroscopy.

Fig. 9
figure 9

Results of simulations showing how well the IGM neutral fraction χHI and host Hi column density can be retrieved from typical ELT-class afterglow spectra (cf. Fig. 5). The input values are shown as green points, and the resultant fits by the red contours. In nearly all cases the IGM neutral fraction is well recovered (within the symbol size), as are host columns for all but the lowest cases, although even here results should be sufficient to indicate a partially transparent sight line. The variance between sight-lines at the same redshift of IGM neutral fraction is expected to be up to \(\delta \chi _{\text {HI}}\sim 0.3\) at mid reionization [107, 108], depending on the reionization model, which can be determined with the THESEUS sample

Full size image

Using GRBs for this endeavour offers advantages over reionization studies that employ Ly α emitters (LAEs; e.g. [24]) for two main reasons: firstly, GRBs suffer less from selection bias compared to LAEs; secondly, one can often get an accurate redshift measurement via the presence of absorption lines in the afterglow spectrum. This may be contrasted to LAEs, and also Lyman Break Galaxies (LBGs), where redshift determination can sometimes be ambiguous when it depends solely on Ly α emission itself (in the former) and the presence of the break (in the latter), with potential confusion with low-z interlopers.

Observations of the damping wing in GRB afterglow spectra also offer advantages over the same type of exercise in quasars (cf. [109]). In the first instance, one expects that GRBs do not represent such biased sources as quasars; in this respect, with GRBs one would not be probing ”special” regions, e.g. as in higher-density environments where quasars are regularly located. Moreover, on account of GRBs being short-lived events, one does not expect their emitted radiation to have a large impact on the ionisation state of the surrounding IGM, which means that one would obtain a more representative picture than with quasars. Finally, quasars often exhibit complicated continuum emission and bright emission lines such as Ly α, making the damping wing modelling and the analysis of the foreground Lyman forest particularly challenging. By contrast, the featureless continuum of GRB afterglows makes this exercise much more simple and leads to Hi fraction and column density estimates with much smaller uncertainties (see also Fig. 5).

8 Population III stars and primordial galaxies

High redshift GRBs provide several routes to exploring the earliest populations of metal-free and ultra-low metallicity stars, from direct spectroscopic determinations of ISM abundances, to changes in the numbers and properties of the bursts themselves, reflective of the increase in average and maximum masses expected for Population III stars [66].

In the current ΛCDM paradigm, the very first stars (the so-called Population III stars) are expected to form from pristine gas, primarily at very high redshifts, \(z\sim 10\)–30 [110]. The absence of heavy elements and consequent inefficiency of cooling at these early cosmic times, is expected to lead to masses that largely exceed those of Pop-I and Pop-II stars (possibly reaching several hundreds of solar masses). When these first stars reach their final stage of evolution, their low-opacity envelope combined with limited mass loss from stellar winds may thus keep large amounts of gas bound until the final collapse, favouring the conditions for jet breakout and for the launch of a very energetic long GRB [111]. Such models predict that the total equivalent isotropic energy released by such Pop-III star explosions could exceed by several orders of magnitude that of GRBs from Pop-I/II progenitors, possibly reaching \(\sim 10^{56} -10^{57}\) erg and making them detectable up to the highest possible redshifts [112]. Given the peculiar energetics and chemistry associated with Pop-III star GRBs, their observed properties should differ from GRBs at lower redshift. In particular, their prompt emission could extend over much longer timescales, potentially lasting up to a month [113].

Similarly, the energy released by the jetted explosion could imply significantly longer times to dissipate. This could give rise to a luminous afterglow emission peaking much later and at higher fluxes than observed for Pop-II GRBs, which in particular would be extremely bright at radio wavelengths if happening in a high-density ISM [114]. A luminous thermal component could also be produced if the jet deposits a large fraction of its total energy in the stellar envelope of the GRB progenitor. Analogous behaviour may have been seen in the “ultra-long” duration GRB 130925A [115], although in this case the event occurred in a low redshift (z = 0.35) galaxy, and apparently in a metal-rich environment [116].

Spectroscopy of such afterglows with 30 m-class telescopes (or JWST, if still operational) and ATHENA (probing photo-ionized gas closer to the burst) may reveal ultra-low metallicity if their line of sight intersects pristine gas in their host galaxy, and these signatures would represent a direct piece of evidence for the association between a GRB and a Pop-III star progenitor. Similarly, gas cloud pockets enriched by Pop-III star explosions and illuminated by Pop-II GRBs might provide us with another way to explore the metal abundance patterns characterizing the ISM enriched only by the very first generation of massive stars [117, 118]. Such abundance patterns may in fact reveal if the first heavy elements were produced by typical core-collapse explosions or by other mechanisms such as pair-instability supernovae, hence constraining the Initial Mass Function up to the earliest cosmic epochs [119, 120].

It has been suggested that 21 cm absorption features (the 21 cm forest) in the spectra of high-z radio sources could be used to gather information on the neutral hydrogen content in the IGM along the line of sight, similarly to what is presently done with the Ly α forest at lower redshift. Although this is an intriguing possibility, the lack of high-z sufficiently bright radio-loud sources makes it a purely theoretical one. If a GRB with a Pop-III star progenitor were detected, this could be an excellent candidate for follow-up searches of absorption features in its afterglow spectrum with the present or next generation of radio instruments such as SKA [121]. While the example in Fig. 10, refers to a source located at redshift 10, similar arguments apply to GRBs from Pop-III stars throughout the epoch of reionization.

Fig. 10
figure 10

Top: spectrum of a GRB with a Pop-III star progenitor located at zs = 10 with a flux density of 30 mJy. The lines refer to the intrinsic spectrum of the source (red dotted), the simulated spectrum for 21 cm absorption (blue dashed), and the absorption spectrum as it would be seen with a bandwidth of 10 kHz and an integration time of 1000 h (black solid). The left and right panels refer to observations with the LOFAR and SKA-1 telescopes. Bottom: corresponding S/N ratio. See [121] for further details

Full size image

Predicting how many of such events THESEUS will identify throughout the duration of its science operations is particularly challenging since the probe of such first light sources still remains a fully uncharted territory. Albeit rare and quickly transitioning to a more regular Pop-II regime [58] [122], Pop-III star formation at high redshift may represent a significant proportion of the massive stars that end their lives as long GRBs. Relying on Swift calibrations and numerical simulations of primordial structure formation, the fraction of Pop-III GRBs (born in metal-poor sites with metallicity Z < 10− 4Z⊙) may increase from low to high redshift up to 10% (40%) at z > 6 (z > 10) [54], in line with the trend in Fig. 6. Considering that THESEUS will observe around 40-80 GRB candidates at z > 6, there is an estimated probability of detecting at least 4-8 GRBs from Pop-III origin at z > 6. Such detections have the potential to be revolutionary, being the first primordial explosions ever observed.

In the era of 30 m class optical/nIR ground-based telescopes, THESEUS could be the only experiment enabling the discovery of these first very massive stellar explosions, hence providing exquisite targets for further follow-up with larger telescopes. Pop-III stars likely played a major role in the growth of the very first bound structures at early cosmic times, through chemical feedback and metal enrichment of the primordial IGM, and they may have also contributed a head start to the cosmic re-ionization process. To date, no direct evidence of the connection between Pop-III stars and GRBs has been observationally established. The identification by THESEUS of even a single GRB with metal abundance unveiling a Pop-III star progenitor or a Pop-III star enriched medium would put fundamental constraints on the unknown properties of the first stars and represent a major breakthrough in our understanding of first-light sources.

9 Probing the expansion history of the Universe and dark energy with GRBs

While the standardised candle correlations proposed for GRBs [123,124,125] exhibit larger scatter than those for SNIa, and are subject to potentially complex selection effects (e.g. [126]), GRBs can be observed to much higher redshift, and thus, in principle, provide much greater lever arm for testing cosmological world models.

By using the spectral peak energy – radiated energy (or luminosity) correlation Ep,i − Eiso [123] (or Liso; [124]), it has been argued that GRBs offer a promising tool to probe the expansion rate history of the Universe beyond the current limit of \(z\sim 2\) (Type-Ia SNe and Baryonic Acoustic Oscillations from quasar absorbers). With the present data set of GRBs, cosmological parameters consistent with the concordance cosmology can already be derived [127, 128]. The order of magnitude improvement provided by THESEUS on the size of the sample of GRBs above z ≈ 3, with more uniform selection (on-board measured redshifts in many cases and broad-band spectral parameters), will allow us to further refine the reliability of these methods and, possibly, to characterize the equation of state of Dark Energy. For example, the THESEUS sample may help assess whether the dark energy equation of state evolves with redshift, such as in extended dynamical dark energy models, like scalar field quintessence or alternative theories of gravity [129, 130]. Combined with the most recent inferrences from the Cosmic Microwave Background, this will offer the unique opportunity to constrain the geometry, and therefore the mass-energy content of the Universe back to \(z\sim 5\) or further, extending beyond the investigations of EUCLID and of next-generation large-scale structure surveys to the entire cosmic history.

Furthermore, at higher redshift, effects from alternative scenarios for dark matter nature could be detectable. They should impact the first collapsing phases of cosmic gas and cascade over the entire star formation history. In particular, warm dark matter (WDM) is an interesting possibility advocated to solve small-scale problems of the standard model based on cold dark matter (CDM). It relies on the assumption that dark matter is made by warm particles having a non-negligible (although small) thermal velocity. IGM based calibrations suggest a mass limit for WDM particles corresponding to about 3 keV which would produce a sharp decrease of power at relatively small scales. With this constraint, there is little room to disentangle between WDM and CDM in terms of low-redshift star formation, but it is possible to test WDM-induced delays in the gas molecular cooling and the onset of star formation at early times. In this respect, the high-redshift regime probed by THESEUS (z > 6) becomes an important new window that needs to be explored. During primordial epochs, WDM models predict a neat cut-off in the SFR and GRB history, departing from CDM expectations [131]. By tracing early GRBs with THESEUS, it will be possible to probe the nature of dark matter (WDM vs. CDM), independently from other techniques, and the lack of observed objects during the epoch of reionization will help pose serious constraints on WDM. This will be very interesting, as WDM effects on structure growth are thought to be more dramatic than the ones derived from dark energy [132] or from non-linear high-order corrections to cosmological perturbation theory [133].

Long GRBs detected by THESEUS in such early epochs could be useful also to constrain local non-Gaussianities in the matter distribution. Indeed, positive (negative) deviations from the Gaussian shape of the primordial density field are expected to enhance (inhibit) structure formation in early epochs. Assuming that GRBs are fair tracers of cosmic star formation, this means that positive non-Gaussianities might increase the expected GRB rate at high redshift [134]. Any excess of long GRB detections at z > 10 by THESEUS would impact the high-z GRB rate and help constrain possible non-Gaussian scenarios. Even one long GRB with a rate of at least 10− 6yr− 1sr− 1 in such early epochs would be sufficient to put additional constraints on the matter distribution. GRB detections in line with the cosmic star formation rate density, instead, would support the commonly adopted Gaussian assumption. In any case, once observed by THESEUS, long GRBs exploding during the epoch of reionization would provide new insights about the properties of dark matter.

10 Conclusions

Thanks to its unprecedented sensitivity and follow-up performance compared to past and current GRB-dedicated experiments, THESEUS will identify more than 40 events at z > 6, hence providing a much larger sample of high redshift GRBs than achieved so far. This will open a unique window on the early Universe, supplementing the picture of the first billion years of cosmic evolution to be painted in the coming decade by JWST, the ELTs and SKA. Spectroscopic follow-up of their afterglows will address a number of key issues such as the history of chemical enrichment in the ISM, the timeline of reionization and the escape fraction of ionizing radiation in early star-forming galaxies. The occurrence rate of high-z GRBs and the properties of their hosts will also further refine constraints on the cosmic history of star formation and the luminosity function at the very faint end, which may contribute solving the still open debate on the overall contribution of massive stars to reionization. THESEUS will finally open new routes for chasing explosion signatures from the still-elusive population III stars, hence offering exciting opportunities to probe primordial galaxies.