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

HAL
open science

# Modelling of long gamma-ray burst host galaxies at cosmic noon from damped Lyman-α absorption statistics

Jens-Kristian Krogager, A de Cia, K Heintz, J Fynbo, L Christensen, G

Björnsson, P Jakobsson, S Jeffreson, C Ledoux, P Møller, et al.

► To cite this version:

Jens-Kristian Krogager, A de Cia, K Heintz, J Fynbo, L Christensen, et al. Modelling of long gamma-ray burst host galaxies at cosmic noon from damped Lyman-α absorption statistics. Monthly Notices of the Royal Astronomical Society, 2024, 535 (1), pp.561-573. 10.1093/mnras/stae2348. hal-04802210

HAL Id: hal-04802210

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

Submitted on 30 Nov 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling of long gamma-ray burst host galaxies at cosmic noon from damped Lyman-α absorption statistics

J.-K. Krogager<sup>1,2*</sup> A. De Cia,<sup>1,3</sup> K. E. Heintz<sup>1,4</sup> J. P. U. Fynbo,<sup>4</sup> L. B. Christensen<sup>1,4</sup> G. Björnsson,<sup>5</sup> P. Jakobsson,<sup>5</sup> S. Jeffreson,<sup>6</sup> C. Ledoux,<sup>7</sup> P. Møller<sup>3,4</sup> P. Noterdaeme,<sup>8</sup> J. Palmerio,<sup>8</sup> S. D. Vergani<sup>8,9,10</sup> and D. Watson<sup>4</sup>

Accepted 2024 October 12. Received 2024 September 11; in original form 2023 November 10

## ABSTRACT

We study the properties of long gamma-ray burst (GRB) host galaxies using a statistical modelling framework derived to model damped Lyman-α absorbers (DLAs) in quasar spectra at high redshift. The distribution of $N_{H1}$ for GRB-DLAs is ~10 times higher than what is found for quasar-DLAs at similar impact parameters. We interpret this as a temporal selection effect due to the short-lived GRB progenitor probing its host at the onset of a starburst where the interstellar medium may exhibit multiple overdense regions. Owing to the larger $N_{H1}$, the dust extinction is larger with 29 per cent of GRB-DLAs exhibiting $A(V) > 1$ mag in agreement with the fraction of ‘dark bursts’. Despite the differences in $N_{H1}$ distributions, we find that high-redshift $2 < z < 3$ quasar- and GRB-DLAs trace the luminosity function of star-forming host galaxies in the same way. We propose that their differences may arise from the fact that the galaxies are sampled at different times in their star formation histories, and that the absorption sightlines probe the galaxy haloes differently. Quasar-DLAs sample the full H1 cross-section, whereas GRB-DLAs sample only regions hosting cold neutral medium. Previous studies have found that GRBs avoid high-metallicity galaxies (~0.5 $Z_⊙$). Since at these redshifts galaxies on average have lower metallicities, our sample is only weakly sensitive to such a threshold. Lastly, we find that the modest detection rate of cold gas ($H_2$ or C1) in GRB spectra can be explained mainly by a low volume filling factor of cold gas clouds and to a lesser degree by destruction from the GRB explosion itself.

**Key words:** galaxies: high-redshift – galaxies: ISM – galaxies: statistics – gamma-ray bursts.

## 1 INTRODUCTION

Studying galaxy evolution via flux-limited samples tends to bias our view towards the brightest galaxies at any cosmic epoch. One way of overcoming this selection effect is to identify galaxies based on methods that probe below conventional luminosity limits. For example, it is possible to study much fainter galaxies if they are magnified by a gravitational lens in the foreground (e.g. Pettini et al. 2000; Hainline et al. 2009; Christensen et al. 2012; Stark et al. 2013; Atek et al. 2018; Bouwens et al. 2022).

Alternatively, we can identify galaxies owing to their neutral gas absorption towards bright background sources (damped Lyman-α absorbers, DLAs, with $\log(N_{H1}/\text{cm}^{-2}) > 20.3$; Wolfe et al. 1986)

or by identifying luminous transients (such as gamma-ray bursts) that hint at the existence of these faint galaxies even if their emission is not detected (Tyson & Scalo 1988; Natarajan et al. 1997; Wijers et al. 1998; Hogg & Fruchter 1999; Vreeswijk et al. 2004; Jakobsson et al. 2005; Prochaska et al. 2007).

Gamma-ray bursts are powerful cosmic beacons that have been observed out to $z \sim 8$ (Salvaterra et al. 2009; Tanvir et al. 2009). A few host galaxies at $z \sim 6$ have been detected in emission (McGuire et al. 2016), and their afterglows can be bright enough to allow in-depth studies of the chemical abundances in their host galaxy at these high redshifts (Hartoog et al. 2015; Saccardi et al. 2023).

Long-duration gamma-ray bursts with a duration of prompt emission longer than 2 s (hereafter referred to simply as GRBs; Kouveliotou et al. 1993) are directly associated with the deaths of massive stars (Galama et al. 1998; Wijers et al. 1998; Bloom et al. 2002; Hjorth et al. 2003; Stanek et al. 2003) that trace regions

<footnote n="1"><sup>1</sup>Centre de Recherche Astrophysique de Lyon, Université de Lyon 1, ENS-Lyon, CNRS, UMR5574, 9 Av Charles André, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy, University of Geneva, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="3"><sup>3</sup>European Southern Observatory, Karl-Schwarzschildstrasse 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="4"><sup>4</sup>Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Jagtvej 128, DK-2200 Copenhagen N, Denmark</footnote>

<footnote n="5"><sup>5</sup>Centre for Astrophysics and Cosmology, Science Institute, University of Iceland, Dunhagi 5, 107 Reykjavík, Iceland</footnote>

<footnote n="6"><sup>6</sup>Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="7"><sup>7</sup>European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago, Chile</footnote>

<footnote n="8"><sup>8</sup>Institut d'Astrophysique de Paris, Université Paris 6-CNRS, UMR7095, 98bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="9"><sup>9</sup>GEPI, Observatoire de Paris, PSL Research University, CNRS, Univ. Paris Diderot, Sorbonne Paris Cité, Place Jules Janssen, F-92195 Meudon, France</footnote>

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

<footnote>*E-mail: jens-kristian.krogager@univ-lyon1.fr</footnote>
of active star formation in galaxies on time-scales of ~10 Myr. However, recent observations also indicate that some long-duration GRBs are associated with kilonovae which would indicate a binary progenitor in a merger scenario (Rastinejad et al. 2022). Other recent observations also seem to hint at a more complex population of progenitors in nearby galaxies (de Ugarte Postigo et al. 2024; Thöne et al. 2024). While some models of binary neutron star mergers predict a non-negligible fraction of mergers to occur on short time-scales (< 30 Myr; e.g. Beniamini & Piran 2024), the more typical lifetimes are expected to exceed 100 Myr. Thus, at the high redshifts ($z > 2$) studied in this work, the merger progenitors would contribute less to the GRB population as they require much longer time-scales.

Since the brightness of the GRB afterglow is not related to the luminosity of its host galaxy, we can identify and study fainter, more metal-poor galaxies that would otherwise go unseen in standard, flux-limited surveys (e.g. Fynbo et al. 2001, 2008). By obtaining spectra of the bright optical afterglow before it fades, we can furthermore study the interstellar medium (ISM) of the GRB host galaxy along the line of sight in absorption. In particular, we can characterize the column density of neutral hydrogen, $N_{\text{H}}$, and the abundance of metals and molecules in great detail (Jakobsson et al. 2006; Prochaska et al. 2007; Fynbo et al. 2009; De Cia et al. 2012; Cucchiara et al. 2015; Bolmer et al. 2019; Heintz et al. 2019b).

Despite their efficient use as cosmological probes of galaxies, GRBs may not trace the full population of star-forming galaxies in an unbiased manner; for example, an upper threshold in metallicity (or lower limit in specific star formation rate, SFR/$M_*$) has been invoked to explain the distribution of GRB host-galaxy properties (e.g. Christensen, Hjorth & Gorosabel 2004; Vergani et al. 2015; Japelj et al. 2016; Perley et al. 2016; Palmerio et al. 2019; Björnsson 2019; Metha & Trenti 2020). Understanding such selection effects is vital for the interpretation of how GRBs trace galaxies as well as how they trace the cosmic star formation rate (SFR) density at the highest redshifts (Kistler et al. 2009).

We can use simplified statistical models to investigate the link between the underlying galaxy population and the absorption properties that we observe. Such modelling has successfully been applied to high-redshift DLAs identified in random quasar sightlines (hereafter quasar-DLAs; Fynbo, Møller & Warren 1999; Fynbo, Burud & Møller 2000; Krogager et al. 2017, 2020; Rhodin et al. 2018; Krogager & Noterdaeme 2020) and to high-redshift GRB-DLAs (Fynbo et al. 2008).

If the probability of a GRB occurring in any given galaxy, and thus for that galaxy to be identified as a GRB host, is proportional to the SFR, as posited by Porciani & Madau (2001), then the ensemble of GRB host galaxies should be sampled uniformly from the overall population of star-forming galaxies weighted by their SFR (Fynbo et al. 2008).

Quasar-DLAs, on the other hand, sample the star-forming galaxy population weighted by their cross-section of DLA gas, $\sigma_{\text{DLA}}$, where $N_{\text{H}} > 2 \times 10^{20} \text{ cm}^{-2}$ (Wolfe et al. 1986). We assume that this projected DLA cross-section is proportional to luminosity (Krogager et al. 2020). At these high redshifts ($z \sim 2$), the passive galaxy population contributes no more than 20 per cent of the total stellar mass (e.g. Santini et al. 2022), and since the high-column density $N_{\text{H}}$ cross-section is expected to be suppressed for passive galaxies, we therefore neglect the contribution of passive galaxies in this model. Assuming that SFR scales directly with the ultraviolet (UV) luminosity as well (Kennicutt 1998, for unobscured star formation; see discussion in Section 4.2), GRB-DLAs and quasar-DLAs both sample the galaxy population weighted by luminosity. We would

therefore naively expect the host galaxies of high-redshift GRB-DLAs and quasar-DLAs to be drawn from the same underlying population of star-forming galaxies, i.e., they sample the luminosity function in the same way.

The aim of this paper is to test this simple scenario by comparing the observed properties of GRB-DLAs to the model by Krogager et al. (2020) and Krogager & Noterdaeme (2020). Using the two-phase model of the neutral medium implemented by Krogager & Noterdaeme (2020), we assume that GRBs only arise in parts of galaxies where the dense cold neutral medium (CNM) can be maintained in pressure balance, since the cold dense gas is needed for the formation of stars. Given the way quasar sightlines randomly probe the galactic environment, quasar-DLAs on the other hand probe predominantly the more extended warm neutral medium (WNM). We then compare the model predictions to the observed distributions of $N_{\text{H}}$, metallicity, $A(V)$, and impact parameters of GRB-DLAs. In this work, we only consider high-redshift bursts ($z \gtrsim 2$) for which $N_{\text{H}}$ can be measured directly from ground-based spectroscopy via the Ly$\alpha$ transition.

The paper is organized as follows: we first describe the compilation of data from the literature in Section 2. In Section 3, we present an overview of the statistical model used in our work and how this is applied to model GRB host galaxies before presenting our results of the model comparison in Section 3.3. In Section 4, we offer a discussion of how GRBs trace star-forming galaxies in light of our modelling results. Lastly, in Section 5, we provide a short summary of our findings.

Throughout this paper, we use the following cosmological parameters: a flat Lambda-cold dark matter cosmology with $H_0 = 68 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_M = 0.31$ (Planck Collaboration XIII 2016).

## 2 SAMPLE SELECTION

Since the model to which we compare our data has been tuned to absorption properties at redshift $z = 2.5$, we only consider bursts in the redshift range $2 \lesssim z \lesssim 3.5$.

### 2.1 Metallicity and dust extinction

We collect a sample of absorption metallicities and dust extinction measurements from GRB afterglow spectroscopy with intermediate (or higher) spectral resolution. We include measurements from VLT<fnref n="1" />/X-shooter, VLT/UVES, Keck/ESI, and Keck/HIRES in order to have robust metallicity measurements.

In order to constrain the gas-phase metallicity, the spectroscopic follow-up of the afterglow has to be of sufficient quality to constrain the column densities of the narrow metal lines. This therefore introduces a bias against dust-obscured bursts (see also Perley et al. 2016). This is directly evident in the distribution of $A(V)$ in our sample compared to the more complete study by Covino et al. (2013) who find values of $A(V)$ up to ~3 mag. This incompleteness is taken into account in our modelling, see Section 3.2.

The compilation of metallicities and dust extinction for the total of 24 GRBs is given in Table 1. This sample is referred to as the 'GRB metal sample' in the remainder of this paper.

<footnote n="1">¹Very Large Telescope of the European Southern Observatory</footnote>
<table><caption>**Table 1.** Compilation of GRB metal sample from afterglow absorption spectroscopy of bursts in the redshift range $2 < z < 3.5$.</caption><thead><tr><th>GRB</th><th>z<sub>abs</sub></th><th>[X/H]</th><th>X</th><th>A(V)</th><th>Ref.</th></tr></thead><tbody><tr><td>000926</td><td>2.0380</td><td>-0.11</td><td>Zn</td><td>0.38</td><td>(1, 9)</td></tr><tr><td>030226</td><td>1.9870</td><td>-1.05</td><td>Fe</td><td>0.00</td><td>(2)</td></tr><tr><td>050820A</td><td>2.6150</td><td>-0.39</td><td>Zn</td><td>0.27</td><td>(4, 2, 6, 9)</td></tr><tr><td>050922C</td><td>2.1990</td><td>-2.09</td><td>Si</td><td>0.00</td><td>(4, 3)</td></tr><tr><td>071031</td><td>2.6920</td><td>-1.76</td><td>Zn</td><td>0.00</td><td>(4)</td></tr><tr><td>080413A</td><td>2.4330</td><td>-1.63</td><td>Zn</td><td>0.00</td><td>(4)</td></tr><tr><td>081008</td><td>1.9685</td><td>-0.52</td><td>Zn</td><td>0.08</td><td>(5, 6)</td></tr><tr><td>090809A</td><td>2.7373</td><td>-0.86</td><td>Zn</td><td>0.11</td><td>(7)</td></tr><tr><td>090926A</td><td>2.1069</td><td>-1.97</td><td>Zn</td><td>0.03</td><td>(7)</td></tr><tr><td>111107A</td><td>2.8930</td><td>-0.74</td><td>Si</td><td>0.15</td><td>(7)</td></tr><tr><td>120327A</td><td>2.8143</td><td>-1.49</td><td>Zn</td><td>0.05</td><td>(7)</td></tr><tr><td>120716A</td><td>2.4874</td><td>-0.71</td><td>Zn</td><td>0.30</td><td>(7)</td></tr><tr><td>120815A</td><td>2.3582</td><td>-1.45</td><td>Zn</td><td>0.19</td><td>(7)</td></tr><tr><td>121024A</td><td>2.3005</td><td>-0.76</td><td>Zn</td><td>0.26</td><td>(7)</td></tr><tr><td>130408A</td><td>3.7579</td><td>-1.48</td><td>Zn</td><td>0.12</td><td>(7)</td></tr><tr><td>141028A</td><td>2.3333</td><td>-1.64</td><td>Si</td><td>0.13</td><td>(7)</td></tr><tr><td>141109A</td><td>2.9940</td><td>-1.63</td><td>Zn</td><td>0.16</td><td>(7)</td></tr><tr><td>150403A</td><td>2.0571</td><td>-1.04</td><td>Zn</td><td>0.12</td><td>(7)</td></tr><tr><td>151021A</td><td>2.3297</td><td>-0.98</td><td>Zn</td><td>0.20</td><td>(7)</td></tr><tr><td>160203A</td><td>3.5187</td><td>-1.31</td><td>S</td><td>0.10</td><td>(7)</td></tr><tr><td>161023A</td><td>2.7100</td><td>-1.23</td><td>S</td><td>0.09</td><td>(7)</td></tr><tr><td>170202A</td><td>3.6456</td><td>-1.28</td><td>S</td><td>0.08</td><td>(7)</td></tr><tr><td>181020A</td><td>2.9379</td><td>-1.50</td><td>Zn</td><td>0.27</td><td>(8)</td></tr><tr><td>190114A</td><td>3.3764</td><td>-1.16</td><td>Zn</td><td>0.36</td><td>(8)</td></tr></tbody></table>

Note. References: (1) Savaglio, Fall & Fiore (2003); (2) Prochaska et al. (2007); (3) Piranomonte et al. (2008); (4) Ledoux et al. (2009); (5) D'Elia et al. (2011); (6) Wiseman et al. (2017); (7) Bolmer et al. (2019); (8) Heintz et al. (2019b); and (9) Zafar & Møller (2019).

## 2.2 Neutral hydrogen column density

In order to improve the statistics on $N_{H1}$ measurements and to investigate possible biases related to the metallicity determinations, we collect measurements of $N_{H1}$ from the comprehensive work by Tanvir et al. (2019). Since we only consider GRB-DLAs in this work, we restrict the sample by Tanvir et al. (2019) to those absorbers with $log(N_{H1}/cm^{-2}) > 20.3$. When comparing the Tanvir et al. sample to our GRB metal sample there appears to be a slight lack of high-$N_{H1}$ absorbers in the metal sample, which could indicate a dust-obscuration bias as mentioned above. However, the difference is not statistically significant, as evidenced by the p-value of the two-sample Kolmogorov–Smirnov (KS) test $p = 0.72$.

## 2.3 Impact parameters

The sample of impact parameters, i.e., the projected distance from the GRB explosion site to the luminosity-weighted centre of the galaxy, used in this paper has been taken from the work by Lyman et al. (2017). These authors associate host galaxies to GRBs using near-infrared imaging data from the Hubble Space Telescope (i.e. rest-frame optical). We consider bursts down to a slightly lower redshift limit ($1.5 < z < 3.5$) than for the metallicities, since the sample is otherwise too small for a meaningful comparison. In total, there are nine bursts that meet our criteria. All impact parameters have been corrected for differences in assumed cosmology.

Lyman et al. (2017) report non-detections of impact parameters for 6 out of 15 bursts that meet our redshift criteria. Hence, $40 \pm 13$ per cent of GRB host galaxies at $z > 1.5$ are not detected down to a luminosity limit of $M \gtrsim -18$ (rest-frame V or R band depending

on the redshift; see their fig. 5 and section 4.2). This may introduce a bias in the distribution of impact parameters that will be discussed in more detail in Section 3.3.

We also compare to the work by Blanchard et al. (2016, see also Bloom et al. 2002; Fruchter et al. 2006) who study the associated host galaxies in a mix of near-infrared and optical photometric bands (rest-frame UV to optical). We again only consider bursts in the redshift range $1.5 < z < 3.5$, which gives a total sample of 32 bursts. Out of these, six bursts have no associated host. The fraction of non-detected bursts ($19 \pm 6$ per cent) is thus lower than what Lyman et al. (2017) infer but consistent within the rather large uncertainties due to the low number of bursts.

In some cases, the host galaxies associated by Blanchard et al. (2016) seem to be biased towards brighter and larger galaxies in the field. One example of this is seen for the burst 080319C at $z = 1.95$ which is studied by both teams. Blanchard et al. (2016) associate a bright spiral galaxy at a projected distance of 7 kpc as the host, whereas Lyman et al. (2017) assign a much fainter source as the host at a projected distance of 0.3 kpc. In this work, we compare our model to the measurements of both studies in order not to give preference to any particular sample.

## 3 STATISTICAL MODELLING

The model used in this work is based on the modelling approach by Fynbo et al. (2008), which has been developed further to include a statistical prescription for $N_{H1}$ (Krogager et al. 2020) as well as a two-phase description of the neutral gas as either a warm or cold phase (Krogager & Noterdaeme 2020). We here offer a short summary of the model framework. For details regarding the model implementation, see Krogager et al. (2020) and Krogager & Noterdaeme (2020).

The model has been designed to study the properties of quasar-DLAs in a simplified manner. These absorbers are selected based on their projected H1 cross-section with the probability density $P_{QSO-DLA}(L) \propto \sigma_{DLA}(L) \phi(L)$, where $\sigma_{DLA}(L)$ is the effective cross-section of DLAs and $\phi(L)$ is the UV luminosity function.

Krogager et al. (2020) use a simple circular projected cross-section, such that $\sigma_{DLA}(L) = \pi R_{DLA}(L)^2$. The typical radial extent of $\sigma_{DLA}$ is further assumed to scale with luminosity motivated by the tight mass–size relation for H1 (Stevens et al. 2019). We use a Holmberg relation of the form: $R_{DLA}(L) = R_{DLA}^*(L/L^*)^t$ where $R_{DLA}^*$ is the radial extent of the DLA cross-section for an $L^*$ galaxy. Krogager et al. (2020) obtain a value of $t = 0.5$, i.e., $\sigma_{DLA}(L) \propto L$. The absolute scaling of $R_{DLA}(L)$, which in turn determines the typical impact parameters of quasar DLAs, is obtained by normalizing the cross-section to the DLA incidence rate, $dn_{DLA}/dz$ (Zafar et al. 2013). Galaxies are then sampled from the luminosity function weighted by $\sigma_{DLA}$, and an impact parameter, $b$, for each randomly sampled galaxy is drawn with a probability $P(b) \propto b$ for $b \le R_{DLA}(L)$.

Each galaxy in the model population is assigned a global metallicity, $Z_0$, following a metallicity–luminosity relation: $log(Z_0) = log Z_0^* + \beta \times M_{UV}$ with $\beta = 0.2$. A radial metallicity gradient is assumed to obtain the absorption metallicity, $Z_{abs}$, at the radial position given by $b$. This gradient is taken to be luminosity dependent with a variable power-law index: $\gamma = \gamma^* L^{-0.5}$ following Boissier & Prantzos (2001) with $\gamma^* = -0.019 \text{ dex kpc}^{-1}$ (Krogager et al. 2020). However, at the small physical scales probed by the GRBs in this model, the metallicity gradient is negligible.

Together with the absorption metallicity, the model further assigns a total hydrogen column density, which is subsequently split into
<figure><img src="image_3.png" /></figure>

**Figure 1.** Average column density of total H, H₂, and HⅠ as a function of galactic radius for an L* galaxy in our model. The radius of R<sub>DLA</sub> (red dotted vertical line) scales with host galaxy luminosity as shown in Fig. 2 and denotes the outer extent of the randomly drawn impact parameters in our model. Within R<sub>CNM</sub> (blue dashed vertical line) defined by P > P<sub>min</sub>(Z) (Wolfire et al. 2003) the dense CNM clouds are stable in pressure equilibrium with the WNM giving rise to H₂ and CⅠ absorption. The R<sub>CNM</sub> from our model coincides with the transition from the HⅠ-dominated phase to the H₂-dominated phase. GRBs are expected to arise within R<sub>CNM</sub> (light blue shaded region) as star formation only occurs in these cold dense clouds. The central column density of total H is assumed to be constant for all galaxies motivated by local observations (Bigiel & Blitz 2012).

HⅠ and H₂ column densities (following Blitz & Rosolowsky 2006; Bigiel & Blitz 2012, see Fig. 1). N<sub>H1</sub> is then used to calculate the rest-frame optical dust extinction along the line of sight, A(V), assuming a constant dust-to-metal ratio (Krogager et al. 2019; Zafar & Møller 2019): log(A(V)) = log(Z<sub>abs</sub>) + log(N<sub>H1</sub> / cm<sup>-2</sup>) + κ, with κ = -21.4.

Krogager & Noterdaeme (2020) include a pressure-based two-phase description of the neutral gas following the canonical model by Field, Goldsmith & Habing (1969) describing the WNM and CNM in pressure equilibrium. The two phases are modelled by a metallicity-dependent pressure threshold, P<sub>min</sub>(Z), above which the neutral gas can exist in a stable cold phase (Wolfire et al. 1995, 2003). A radial pressure gradient (see Elmegreen & Parravano 1994) as a function of luminosity is included in the model, giving rise to a characteristic radius, R<sub>CNM</sub>, within which the pressure is large enough to sustain a stable CNM, i.e., P(r) > P<sub>min</sub>. The radial extent of these two phases, R<sub>DLA</sub>(L) and R<sub>CNM</sub>(L), are shown in Fig. 2 as a function of host galaxy luminosity. Since stars only form out of the dense and cold gas phase, we expect star formation to occur predominantly on spatial scales within R<sub>CNM</sub>. We therefore restrict our GRB model to only consider sightlines that probe this CNM region (Fig. 1). Quasar-DLAs, on the other hand, probe the full extent of neutral gas within R<sub>DLA</sub> and less frequently within R<sub>CNM</sub> due to the smaller projected cross-section. The relative cross-sections of DLA and CNM gas defined this way reproduces well the observed scarcity of H<sub>2</sub> detections in DLAs (e.g. Ledoux, Petitjean & Srianand 2003; Jorgenson, Wolfe & Prochaska 2010; Balashev & Noterdaeme 2018; Krogager & Noterdaeme 2020).

We highlight that our statistical model only considers the on-sky projected gas column density distribution. The column density and metallicity of a given absorption sightline is drawn from a one-

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

**Figure 2.** Radial extent of the projected DLA and CNM cross-section as a function of host galaxy luminosity in our model.

dimensional radial distribution function, such as shown in Fig. 1, based on the randomly assigned impact parameter.

## 3.1 The baseline GRB model

As a starting point for our model comparison, we make the simplifying assumption that GRBs directly trace SFR (e.g. Porciani & Madau 2001; Robertson & Ellis 2012), which is assumed to scale with the UV luminosity (Kennicutt 1998); Hence, P<sub>GRB-DLA</sub>(L) ∝ SFR φ(L) ∝ Lφ(L) ∝ P<sub>QSO-DLA</sub>(L). In other words, GRB-DLAs and quasar-DLAs in our model trace the luminosity function of high-redshift star-forming galaxies in a similar way. In our baseline model, their differences therefore arise from the fact that quasar-DLA sightlines are selected randomly within the full extent of neutral gas, i.e., within R<sub>DLA</sub>(L), whereas GRB-DLAs only trace the regions of recent star formation, i.e., within R<sub>CNM</sub>. This is similar to the cartoon illustration by Prochaska et al. (2007, see their fig. 1).

Since the GRB explosion arises within the host galaxy ISM, the absorption line of sight does not probe the full volume of the neutral gas in the same way as quasar sightlines do. We therefore calculate a statistical geometric correction to the GRB sightlines to take into account this diminished column density in front of a GRB as compared to a quasar that hits within R<sub>CNM</sub>. For uniformly distributed burst locations, the average value of this geometric correction is exactly 0.5, i.e., GRB-DLA sightlines should have half the total hydrogen column density compared to quasar-DLAs.

A full three-dimensional modelling of the complex density field is beyond the scope of this simplified model. We therefore scale down the column densities of GRB sightlines by a factor of q drawn randomly within the interval [0.1; 0.9] with <q> = 0.5. We still only consider sightlines in the model with log(N<sub>H1</sub> / cm<sup>-2</sup>) > 20.3 after the geometric correction, given our sample selection criterion. The statistical model explained above with the additional geometric correction will be referred to as our baseline GRB host galaxy model – or ‘baseline GRB model’.

## 3.2 GRB model corrections

A direct comparison of our GRB data to the baseline GRB model described above yields acceptable agreement for the metallicity and impact parameters, when assuming that GRBs only probe the
CNM regions and including the geometric correction. However, the distributions of $N_{H1}$ and $A(V)$ are both significantly underestimated for the same model. For reference, these baseline results are shown as blue lines in Figs 3 and 4. As the baseline model already takes into account the smaller average impact parameters of GRB-DLAs, the observed excess of $N_{H1}$ is therefore not caused by higher densities at smaller galactic radii (see also Section 4.4). Instead, we must introduce one extra degree of freedom in our model to reproduce all observations simultaneously. This additional parameter represents an overall scaling of the average $N_{H1}$ along GRB absorption sightlines. The interpretation of this additional parameter will be discussed in Section 4.1.

The density scaling of H I is included via the parameter $\Delta \log(N_{H1})$ which is the average overdensity of $N_{H1}$ within the CNM region of the host galaxy, i.e., a constant offset to $\log(N_{H1})$ for GRB sightlines. We note that the overdensity of H I might depend on other gas properties such as the metallicity. However, lacking any physically motivated scaling between the density and metallicity, we have ignored such higher order complexities in our analysis. However, we already note that a constant offset to $N_{H1}$ rather than $\log(N_{H1})$ would not provide a solution to the problem, as it would introduce a sharp edge to the distribution of $\log(N_{H1})$ around this value.

With an additional contribution to $N_{H1}$, the amount of dust extinction is also expected to increase for a fixed dust-to-gas ratio. We therefore similarly increase $A(V)$ for GRB sightlines by scaling the initial $A(V)$ to the new value of $N_{H1}$ for any given GRB sightline. Lastly, we change the implementation of the dust bias, which was modelled following Krogager et al. (2019) for quasar-DLAs. Since GRB detection is not susceptible to the complex colour and magnitude selection criteria that quasar-DLAs are, we simply include a binary selection criterion; GRBs with $A(V) > A(V)_{lim}$ are too obscured to yield metallicity measurements in absorption and will therefore be excluded in our model distributions (see Fynbo et al. 2009). The value of $A(V)_{lim}$ is freely variable in order to best reproduce the observed distribution of $A(V)$ for GRBs as mentioned in Section 2.1.

The two parameters, $A(V)_{lim}$ and $\Delta \log(N_{H1})$, are fitted to the observed distributions for the GRB metal sample. The best-fitting parameters are $\Delta \log(N_{H1}) = 1.06 \pm 0.09$ (i.e. the average H I density is 9–14 times higher) and $A(V)_{lim} = 0.44 \pm 0.05$ mag, see Fig. 5. This best-fitting model adaptation of the baseline model is referred to as the ‘GRB model’ in what follows.

To test the proposed metallicity thresholds for GRB progenitors, we have included calculations for two such thresholds: $Z < 0.7 Z_{\odot}$ (Palmerio et al. 2019) and $Z < 0.35 Z_{\odot}$ (Metha & Trenti 2020).

## 3.3 Model results

Fig. 3 shows the observed distributions of metallicity, $N_{H1}$ and $A(V)$ for the GRB metal sample in black. The additional GRB sample by Tanvir et al. (2019) used for comparison is shown in grey. The model predictions for ~15 000 CNM sightlines drawn randomly from the parent galaxy population are shown as the thick, red lines and the p-values based on a KS test comparing the model to the observed samples is indicated in the legend of each panel. For comparison, we also show the baseline model distribution for GRB-DLAs (without the additional density parameter) as thin, blue lines. The intrinsic model distribution without taking an observational dust bias into account is shown by the red, dashed lines. In the top panel of Fig. 3, we also show the distributions assuming the two proposed metallicity thresholds of 0.7 (Palmerio et al. 2019) and 0.35 (Metha & Trenti 2020) Solar as the red, dotted and purple, dash-dotted lines,

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

**Figure 3.** Cumulative distribution functions (top to bottom) of absorption metallicity from afterglow spectroscopy, H I column density from afterglow spectroscopy, and rest-frame optical dust extinction, $A(V)$. The thick red line shows the best model allowing for an increase in average density of the CNM and including GRB-specific dust bias. The red, dashed line shows the same model but without the dust bias. The thin blue line in all panels indicates the baseline model before the density and dust-bias tuning. The KS test p-values are calculated for the GRB model (red solid and dashed lines) compared to the metal sample. In the middle panel, we also give the KS p-value for the Tanvir et al. sample for comparison though the fit was performed using the GRB metal sample. In each panel, we further show a GRB model with an imposed $0.35 Z_{\odot}$ metallicity threshold (purple dash-dotted line). The higher metallicity threshold of $0.7 Z_{\odot}$ is not included in the middle and bottom panels, as it is indistinguishable from the GRB model.
<figure><img src="image_6.png" /></figure>

**Figure 4.** Cumulative distribution functions of observed impact parameters of GRBs from Lyman et al. (2017) in black and from Blanchard, Berger & Fong (2016) as the grey, dashed line. The thick red lines show the GRB model distribution, and the thin, blue lines show the baseline model results. Solid lines refer to the magnitude-limited sample (see the text) and the dotted lines indicate the distributions when considering all host galaxies. The purple dashed-dotted line includes a metallicity threshold.

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

**Figure 5.** Constraints on $\Delta \log(N_{H1})$ and $A(V)_{lim}$ parameters. $\Delta \log(N_{H1})$ parametrizes the overdensity of H I along GRB sightlines compared to quasar sightlines through the same regions of galaxies. $A(V)_{lim}$ quantifies the effects of a dust obscuration bias since we only consider GRB afterglows with absorption-derived metallicities.

respectively. Note that the red, dotted lines are not included in the other panels as they are indistinguishable from the thick red line.

We observe a good agreement between the absorption properties from the GRB metal sample when compared to the full GRB model for which we include one additional free parameter (the additional H I density). Moreover, the higher density of neutral gas needed to reproduce the observed $N_{H1}$ distribution leads to a high fraction of highly obscured bursts which would not enter the metal sample due to the need for good spectroscopic follow-up. This is illustrated in Fig. 6 where we compare the sample of $A(V)$ measurements by Covino et al. (2013) to the intrinsic distribution from the GRB model. The sample by Covino et al. is not restricted by the need to obtain detailed absorption measurements, and is thus more representative and less biased than the GRB metal sample. This is reflected in

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

**Figure 6.** Comparison of unbiased $A(V)$ for high-redshift GRBs by Covino et al. (2013). The grey shaded region indicates the range allowed by the upper limits given by these authors. The GRB model is shown as the dashed red line (same as bottom panel of Fig. 3) without the $A(V)_{lim}$ criterion in order to reproduce an intrinsic sample of GRB-DLAs without the additional requirement that the afterglow be bright enough for spectroscopic analysis. For comparison, we also show the GRB metal sample as the thin solid line.

the much larger $A(V)$ measurements obtained in their sample. We find good agreement between our model without the $A(V)$ limit and the unbiased sample by Covino et al. (2013) at moderate levels of extinction (for $A(V) \lesssim 1$ mag). Our model predicts a fraction of 29 per cent of so-called dark bursts when adopting the definition by these authors, i.e., $A(V) > 1$ mag.

In Fig. 4, we show the observed distribution of impact parameters from both samples described in Section 2 compared to the baseline model and the GRB model. In order to carry out a fair comparison, we restrict our models to similar host galaxy luminosities as the observational magnitude limit.<fnref n="2" /> This luminosity comparison of the host galaxy is not straightforward due to different wavelength ranges probed. The sample by Lyman et al. (2017) is restricted to the rest-frame optical (HST/F160W, ~4400 Å at $z = 2.5$) and the sample by Blanchard et al. (2016) is a mixture of rest-frame optical and UV, whereas our model only produces rest-frame luminosities at 1700 Å. To simplify the calculation, we assume a flat spectral shape (in $F_{\nu}$) between the UV and the optical. Yet, even a small amount of dust optical extinction can severely suppress the UV flux due to the steep extinction curves. We therefore need to correct the UV luminosities for dust-obscuration. Lacking a global dust attenuation model, we use the model $A(V)$ along the line of sight as a proxy of the overall attenuation of the host galaxy, though there is significant scatter between these two quantities (e.g. Krühler et al. 2011). The rest-frame $A(V)$ is converted to extinction at 1700 Å by multiplying by a factor of 4 assuming an average SMC extinction curve (Gordon et al. 2003). We then apply a luminosity cut to these dust-corrected model luminosities using a limit of -17.5 mag as inferred roughly from Lyman et al. (2017, their fig. 5). Furthermore, we include a cut on $A(V) < 3$ mag motivated by the maximum value observed by Covino et al. (2013) to take into account the bursts that have no optical/near-infrared afterglow identified and therefore cannot have an accurate impact parameter measurement.

<footnote n="2">²It is here assumed that all GRBs would intrinsically be bright enough to be detected if there were no dust along the line of sight.</footnote>
We show this luminosity restricted model in Fig. 4 as the solid, red line. The same luminosity limit of the host galaxies is applied to the baseline model (shown in blue). For comparison, we also show the distribution of impact parameters if no luminosity limit is applied (red, dashed line) and if a metallicity threshold is assumed (purple, dash-dotted line). Our model provides acceptable yet marginal agreement with both observed samples.

Lastly, we calculate the fraction of non-detected host galaxies given the adopted luminosity limit and redshift range (2 < z < 3.5). From our model, we find that 42 per cent of GRB host galaxies fall below the adopted luminosity limit for detection of the host galaxies. There is, however, a significant uncertainty ascribed to this estimate given the many assumptions mentioned above. We can gauge this uncertainty by varying the luminosity limit. A change of ±1 mag corresponds to a range of 35–49 per cent of undetected hosts.

# 4 DISCUSSION

The fair agreement between our baseline model and the observed impact parameters and metallicities of 2 < z < 3.5 GRB-DLAs suggests that GRBs sample the full luminosity function of star-forming galaxies weighted by their SFR. This agreement would indicate that GRB-DLAs and quasar-DLAs sample the luminosity function of star-forming galaxies in the same way. This is discussed further in Section 4.4. None the less, there are significant differences in the absorption properties between quasar- and GRB-DLAs. Most notably the distributions of $N_{H1}$ and $A(V)$, see middle and lower panels of Fig. 3. The higher average $N_{H1}$ is not straightforward to explain by our baseline model, and requires the inclusion of one additional free parameter, namely the overdensity of neutral gas, either from gas around the GRB progenitor or in the GRB host ISM in general. We will explore this scenario in more detail in Section 4.1. However, we emphasize that the increased gas density would not directly affect the distribution of impact parameters unless overdense star-forming regions only occur at specific galactic radii which we deem rather unlikely. Similarly for the metallicity, the increased density would not affect the metallicity of the gas. The added $N_{H1}$, however, affects the distributions indirectly through the higher fraction of GRBs that are missed due to optical extinction by dust. Our conclusion that GRB hosts should represent the underlying population of star-forming galaxies at z > 2 therefore still holds regardless of the additional density parameter introduced in this work. None the less, not all GRB hosts are detectable at current luminosity limits, see the discussion in Section 4.2.

While our model is consistent with the data without including a metallicity threshold, we cannot firmly exclude the possibility of a metallicity threshold mainly due to the low number statistics but also due to the redshift range in our study. In particular, the higher metallicity threshold of $Z < 0.7Z_⊙$ by Palmerio et al. (2019) leaves no significant impact on our model due to the low average metallicity of galaxies at these high redshifts. The lower and more restrictive threshold put forward by Metha & Trenti (2020) is not ruled out either but is disfavoured by the detections of higher metallicity GRBs as well as the few large impact parameters (see Fig. 4). We caution that the metallicities presented in our analysis only probe average line-of-sight measurements and are therefore not sensitive to small-scale variations that Metha & Trenti (2020) discuss, or like those observed in the Milky Way (De Cia et al. 2021). None the less, as the metallicity of galaxies globally increases with time (see compilation by De Cia et al. 2018), such metallicity thresholds will become more apparent towards lower redshifts than what we study in this work ($z < 2$; see also Perley et al. 2013, 2016; Palmerio et al. 2019), for

which we do not have absorption data in this work due to the need for space-based spectra to constrain $N_{H1}$ and thus the metallicity.

## 4.1 Overdensity of neutral hydrogen in GRBs

Our modelling indicates that the average CNM density is 9–14 times higher for GRB sightlines compared to quasar-DLAs with CNM tracers (C1 or $H_2$; see Section 3.2), which in our model arise at similar impact parameters. GRB sightlines (and thus GRB-DLAs) are mainly associated with the death of massive stars. We therefore propose two plausible explanations for the observed overdensity: (1) the excess $H1$ is cased by material locally around the massive progenitor, or (2) the host galaxy is overall exhibiting an increase in ISM density, which could be expected during a burst of star formation. In the following, we will explore both of these scenarios in more detail.

### 4.1.1 Are GRBs tracing specific overdense regions of the ISM?

Massive stars are only formed in very dense gas clouds. GRBs could therefore be biased tracers of massive, overdense regions of the ISM. Such overdensities around GRBs are also inferred from X-ray absorption (e.g. Watson et al. 2007; Schady et al. 2011), in which case the overdensity is mainly ascribed to an excess of ionized material along the line of sight. However, the X-ray absorption column density is found to correlate with the neutral gas absorption (Watson et al. 2013) which would further support the idea that GRBs could probe overdense regions of the ISM.

The GRB explosion itself may also affect the neutral gas around the progenitor star. However, photoionization modelling indicates that ionization from the burst itself has limited effects on the $H1$ column density beyond a few pc from the burst<fnref n="3" /> (Ledoux et al. 2009). We therefore do not consider this a significant effect to the full line-of-sight column density given the inferred distances between burst and bulk absorption of the order ~100 pc but may be up to 1.7 kpc (Vreeswijk et al. 2007, 2013; D'Elia et al. 2011; Hartoog et al. 2013). The effect of photoionization of the molecular gas phase will be discussed in Section 4.3.

The most massive stars, thought to give rise to GRBs, will explode very early in the star formation process. If the progenitor is more massive than ~30 $M_⊙$, the lifetime will be ~5–10 Myr (Woosley, Heger & Weaver 2002). Simulations show that the dense molecular clouds rarely survive more than a few Myr due to intense radiative feedback from the massive stars (Kimm et al. 2022). However, in a few cases, depending on geometry and density, the clouds may survive up to 5 Myr. Recent work looking at high-resolution simulations of a dwarf galaxy shows that such massive birth clouds may re-accrete after the initial dispersion by radiative feedback, extending the lifetimes of these clouds up to the order of 10 Myr (Jeffreson, Semenov & Krumholz 2024).

We estimate the time-dependent median gas volume density and $H_2$ gas volume density around massive, GRB-progenitor stars formed within dense molecular clouds, using the output from a high-resolution, chemodynamical simulation of an entire dwarf spiral galaxy (see Jeffreson et al. 2024). The mass resolution of the simulation is 859 $M_⊙$, and so all stellar populations are synthesized

<footnote n="3">³The effects may be significant out to ~100 pc if the initial column density is lower but this would argue against the higher column densities observed on average. Moreover, the models presented by Ledoux et al. (2009) do not consider the shielding from dust which would further reduce the ionizing effects of the burst.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Median total gas volume density (blue) and molecular hydrogen volume density (orange, at half the total density) as a function of time after stellar birth *t*, for each star particle of mass $\ge 30$ Myr (H/He-burning lifetime $\le 6$ Myr, see the text) in a simulation of a dwarf spiral galaxy.</figcaption></figure>

for star clusters (stellar particles) of median mass 859 $M_{\odot}$, drawn from a Chabrier (2003) initial mass function.

Because the runtime memory required to store the stellar population for each star particle in the high-resolution simulation is prohibitively large, we do not have access to the masses of individual stars that were formed during the simulation by Jeffreson et al. (2024). Instead, a supernova progenitor mass of $> 30 M_{\odot}$ is assumed corresponding to a lifetime of $< 6$ Myr (e.g. Chabrier 2003). We therefore select all stellar particles in the simulation that produce supernovae within 6 Myr of their formation, and count the number of supernovae $N_{SN}$ that occur in each particle before this time. We then compute the total gas density $\rho$ and molecular fraction $x_{H_2}$ within a sphere of radius 18 pc (similar to a giant molecular cloud) around each selected stellar particle, during each Myr before its death. In the simulation, there are $5 \times 10^4$ such massive progenitor stars over a period of 300 Myr. The $N_{SN}$-weighted median values of $\rho$ and $x_{H_2}$ across all selected star particles over this 6 Myr interval, along with their interquartile ranges, are shown in Fig. 7.

If we assume a standard average gas density of $\sim 1–10 \text{ cm}^{-3}$ (Ferrière 2001), our modelling suggests a gas density around GRBs of the order $10–150 \text{ cm}^{-3}$. This would imply an average lifetime of the GRB progenitor of $< 3$ Myr based on the results shown in Fig. 7. Such a short lifetime is still consistent with expected lifetimes of massive stars (Woosley et al. 2002); However, we find it unlikely that the re-accretion of only the progenitor birth cloud is enough to explain the overdensity that we infer for GRB sightlines compared to quasar-DLAs with CNM tracers. This is further supported by the rather large inferred distances between explosion site and the bulk of the gas giving rise to fine-structure excitation of metal lines (Vreeswijk et al. 2007, 2013; D'Elia et al. 2011). Hence, the bulk of the absorption most likely does not arise from the immediate surroundings of the GRB progenitor.

### 4.1.2 Are GRBs tracing the onset of a burst in star formation?

Another interpretation of the overdensity of H I would be that GRBs probe galaxies at specific times in their star formation histories. That is, if GRBs mainly occur at the onset of a starburst, or shortly thereafter, the host galaxy may contain more such overdense star-forming clouds on average. This is consistent with observations

of local starbursts whose overall gas surface densities are 1–2 orders of magnitude higher than regular galaxies (e.g. Kennicutt 1998; Bigiel et al. 2008). Such an overdensity would shift the radial column density profile used in our model (Fig. 1) towards higher log($N_{H1}$) for GRBs. Our modelling is consistent with the observations without any further changes to the slope of the radial column density distribution. The slight excess of GRB-DLAs with log($N_{H1}$/cm$^{-2}$) $\sim 21$ compared to our model (see Fig. 3) could be explained by a slightly flatter radial dependence. However, we caution that a full analysis of the radial distribution of log($N_{H1}$) would require more sophisticated models. In future work, we plan to investigate the radial log($N_{H1}$) profile in more detail assisted by numerical simulations.

In this temporal interpretation, the GRB sightlines are thus more likely to intercept high-density neutral gas compared to other times where such dense gas has been dispersed by radiative feedback from the embedded, young stars. A similar conclusion is reached by Hatsukade et al. (2020) who study the molecular gas mass of GRB host galaxies traced by CO emission. These authors find that the GRB host galaxies have higher molecular gas mass fractions than regular star-forming galaxies; However, this excess of molecular gas disappears when comparing to galaxies that have higher specific SFRs similar to GRB hosts.

The temporal interpretation of the overdensity is also consistent with the enhanced specific SFRs often observed in GRB host galaxies (Christensen et al. 2004; Perley et al. 2013; Björnsson 2019). If GRBs predominantly trace a specific time shortly after the onset of a starburst (see Hatsukade et al. 2020), the instantaneous SFR observed would be elevated compared to a random galaxy selected at a more representative phase in its star formation history such as quasar-DLAs with CNM tracers. A logarithmic offset in column density would then naturally follow from the power-law relation between SFR and gas density (Schmidt 1959; Kennicutt 1998). We therefore favour this interpretation of GRBs as tracers of a specific moment in the host galaxy’s star formation history.

## 4.2 The fraction of dark bursts

The fraction of so-called dark bursts is still a subject of debate, and several definitions exist based primarily on X-ray classification (e.g. Jakobsson et al. 2004; van der Horst et al. 2009) or optical dust obscuration (e.g. Krühler et al. 2011; Melandri et al. 2012). The fraction of dark bursts is estimated to be around 10–20 per cent (Jakobsson et al. 2004; van der Horst et al. 2009; Perley et al. 2013) but may be up to 30–40 per cent (Greiner et al. 2011; Melandri et al. 2012) for an adopted dark burst definition of $A(V) > 1$ mag. Using this definition, our GRB model predicts a dark burst fraction of 29 per cent, which is broadly consistent with the GRB observations. Yet, our model predicts a larger fraction of highly obscured bursts ($A(V) > 3$ mag) compared to the sample by Covino et al. (2013). We attribute this larger fraction of dark bursts in our model to the fact that some GRBs are expected to be so obscured that no optical or near-infrared afterglow is identified (Fynbo et al. 2009). These highly obscured bursts would therefore not be present in the sample by Covino et al. (2013) as these authors only include GRBs with a spectroscopic redshift determination. However, these highly obscured GRBs with no detected afterglow may also be caused by other factors, such as being at very high redshift or intrinsically faint bursts. It is therefore not straightforward to quantify the fraction of GRBs with no optical/near-infrared afterglow due to dust obscuration alone.
A related issue to the one of dust-obscured dark GRBs is the number of host galaxies that are not detected (see Hjorth et al. 2012, for a complete sample of GRBs). These undetected host galaxies may bias the associations of impact parameters (Blanchard et al. 2016; Lyman et al. 2017). In our model, we have tried to account for this undetected population of galaxies given the estimated limiting luminosity. We find that our estimated fraction of undetected host galaxies of 35–49 per cent agrees well with the 40 per cent of non-detected hosts by Lyman et al. (2017). However, the observations are based on inhomogeneous imaging data and varying detection limits in the two samples that we compare to. Hjorth et al. (2012) study a more complete sample and find a smaller fraction of about 20 per cent. A fair comparison to our GRB model is further complicated due to the differences in the modelled rest-frame UV luminosities as compared to the observed rest-frame optical range. Significant variations are observed between the line-of-sight dust extinction and the global attenuation of the host galaxy (Krühler et al. 2011; Friis et al. 2015; Heintz et al. 2017; Chrimes et al. 2019; Schroeder et al. 2022). These variations are most likely due to an inhomogeneous dust distribution and possibly very localized dust near the burst (Greiner et al. 2011; Krühler et al. 2011). However, on average there seems to be a tendency for more dusty host galaxies among the most dust-obscured GRB afterglows (Perley et al. 2013; Corre et al. 2018; Schroeder et al. 2022). These complications in dust corrections may explain the slight tension we observe in the predicted impact parameter distributions. From an observational point of view, it is also not straightforward to associate a host galaxy to a GRB as evidenced by the significant differences between the two samples of Blanchard et al. (2016) and Lyman et al. (2017).

## 4.3 CNM tracers in GRB sightlines

Regardless of whether GRBs arise in specific overdense pockets in the CNM or at specific times of starburst activity, one would naively expect to observe tracers of the CNM such as C I or H$_2$ in the majority of GRB sightlines. However, only around 30–40 per cent of GRBs have CNM tracers detected by either H$_2$ absorption (Bolmer et al. 2019) or C I absorption (Heintz et al. 2019a). Indeed, our statistical model predicts a detection rate of nearly 100 per cent based on the detection limit of $log(N_{H_2}/\text{cm}^{-2}) \gtrsim 17$ given by Bolmer et al. (2019).

One possible explanation for the modest detection rate of CNM tracers in GRB afterglow spectra is the efficient photodissociation of H$_2$ and ionization of C I by the GRB explosion itself (Draine & Hao 2002). The ionizing flux produced by the GRB explosion may destroy H$_2$ out to ~100 pc depending on the exact column density (Ledoux et al. 2009). Whalen et al. (2008) on the other hand find that the GRB explosion itself is not able to fully destroy all H$_2$ near the GRB. They argue that H$_2$ must already be suppressed before the GRB goes off. The simulations by Jeffreson et al. (2024) show that H$_2$ can indeed be significantly suppressed by radiative feedback around massive stars, thereby lowering the column density of H$_2$. If the column density is already suppressed before the burst, the GRB explosion will more effectively destroy the molecular gas phase (Ledoux et al. 2009). If we include this destruction of H$_2$ within 100 pc in our model, we find an H$_2$ detection fraction of 90 per cent. This is still inconsistent with the observations. The distance out to which H$_2$ is destroyed would need to be greater than 500 pc in order to match the moderate detection rate reported by Bolmer et al. (2019) and Heintz et al. (2019a).

Another explanation is related to the assumed azimuthally symmetric geometry. While this works well for the extended H$_1$ gas, the CNM bearing gas might realistically be confined to a more

flattened, disc-like structure. GRB sightlines arising from these CNM regions are then more likely to be perpendicular to the disc rather than piercing through it. This would effectively lower the volume of the CNM-bearing medium thereby lowering the probability of intercepting a CNM cloud. However, this would go against the excess of $N_{H_1}$ if the column density is not dominated by the environment close to the GRB, as discussed in Section 4.1. In this flattened CNM geometry, the probability of crossing other high-density star-forming regions would be lower, leading to the same $N_{H_1}$ distribution as for CNM-bearing DLAs, which we do not observe. Hence, unless the local environment of the GRB contributes significantly to the observed excess of $N_{H_1}$, or the H I distribution is different all together, we find this purely geometrical interpretation hard to reconcile with the data. Since our model is not based on a full three-dimensional geometry, a full test of the geometry is beyond the scope of this work. In future work, we will address such higher order geometrical considerations assisted by state-of-the-art numerical simulations.

Lastly, our model assumes a uniform covering fraction of the CNM gas for random background sources piercing the full volume of CNM-bearing gas. Since the volume-filling factor of CNM clouds may be quite low (e.g. Krogager et al. 2018), we do not expect sightlines arising from within the CNM-bearing volume to hit a cold gas cloud in all cases. Instead, the probability of piercing another cold gas cloud is proportional to the fraction of the volume through which the sightline travels (for uniformly distributed clouds). The GRB explosion may therefore locally destroy H$_2$ (or C I) but the absorption sightline still propagates through the rest of the host galaxy ISM. If the GRBs are randomly distributed within the volume where cold clouds are found, the average detection fraction of CNM tracers in our model would instead be ~50 per cent.

The projected covering fraction of CNM gas in quasar absorption systems may in fact be closer to ~80 per cent (Wiklind & Combes 1995; Wiklind, Combes & Kanekar 2018; Krogager et al. 2018; Boissé & Bergeron 2019). An overall lower projected covering fraction would further decrease the chance of a random GRB sightline piercing a cold cloud, bringing the expected H$_2$ detection rate into agreement with the observations (30–40 per cent, Bolmer et al. 2019; Heintz et al. 2019a). We therefore conclude that the most likely explanation for the low detection rate of H$_2$ and C I is due to a combination of local suppression of these species due to the prompt emission of the GRB and a low volume-filling factor of cold gas clouds in the host galaxy ISM.

## 4.4 How do quasar-DLAs and GRB-DLAs trace star-forming galaxies?

There is a fundamental relationship between the velocity width of quasar-DLA gas and the metallicity of the same gas (Ledoux et al. 2006); This relationship has a well-defined and clear evolution with cosmic time which has been traced back to $z = 5.1$ (Møller et al. 2013). If GRBs, at all redshifts, randomly select galaxies from the same population of galaxies as do quasar-DLAs, then the resulting sample should follow the same relation with the same evolution. Arabsalmani et al. (2015) have tested and confirmed that the two samples for GRB- and quasar-DLAs are indeed consistent with being drawn from the same relation (see the lower panel of their fig. 2), providing strong support for the validity of our combined model assumption in this work. Arabsalmani et al. (2015) also find that while following the same relation, GRB-DLAs preferentially populate the high-metallicity end of the relation suggesting that GRB selection is weighted either towards higher metallicity galaxies, towards higher metallicity sightlines in the same galaxies as quasar-
DLAs, or both. Christensen et al. (2014) find that DLA hosts have negative metallicity gradients from their centre and out, meaning that a GRB sightline (located close to the centre) in a DLA galaxy would indeed be expected to show higher metallicity than a random quasar-DLA sightline.

In agreement with the works mentioned above, our modelling supports the hypothesis that DLAs observed in both quasar and GRB sightlines trace the same underlying population of galaxies. On the one hand, quasar-DLAs trace galaxies weighted by their cross-section of neutral gas, which we assume is directly proportional to luminosity (Krogager et al. 2020). On the other hand, GRB-DLAs trace galaxies weighted by the SFR which is also assumed to scale with luminosity (Kennicutt 1998). This is slightly at odds with the previous work (Fynbo et al. 1999; Fynbo et al. 2008) who find that the quasar-DLA cross-section scales with a power of 0.8 instead of a power of 1 that we adopt in this work. This change in the scaling of DLA cross-section as a function of luminosity is related to the updated metallicity measurements used by Krogager et al. (2020) as well as the direct modelling of a dust-obscuration bias which alters the metallicity distribution. To illustrate this, we show a comparison of the metallicity distributions for quasar-DLAs and GRB-DLAs in Fig. 8. The quasar-DLAs extend to lower metallicities given the metallicity gradient and their larger impact parameters on average (see Fig. 9). We find that the median of the GRB-DLA metallicities is 0.2 dex higher than for the quasar-DLAs after applying the correction for the redshift evolution in the mass–metallicity relation (Møller et al. 2013), in agreement with the offset reported by Arabsalmani et al. (2015).

In Fig. 8, we also compare the $N_{H1}$ distributions of quasar-DLAs and GRB-DLAs. We interpret the significant offset as a result of GRBs tracing their host galaxies at the onset of a star formation event leading to an overall increase in gas surface density in the host galaxy. Quasar-DLAs, by selection, trace fully random sightlines through the DLA cross-section of the galaxies at no specific location nor time. We note that the $N_{H1}$ distribution of quasar CNM-DLAs compiled by Krogager & Noterdaeme (2020) shows a marginal excess of high $N_{H1}$ systems. We caution, however, that the sample is not homogeneously selected and is therefore not fully representative. Indeed, the sample includes a number of targets that have been selected based on having $log(N_{H1} / cm^{-2}) \gtrsim 21.5$ (Noterdaeme et al. 2014; Ranjan et al. 2020). Moreover, given the low number of CNM-DLAs, a KS test reveals that this excess is not significant given the p-value of 0.22.

Lastly, we compare the impact parameter distributions between quasar-DLAs and GRB-DLAs in Fig. 9. Our model has been restricted to only consider galaxies that would be bright enough for a detection in ground-based observational campaigns (see Section 3.3 for GRB-DLAs). For quasar-DLAs, we estimate the SFR from the model UV luminosity and consider only hosts with SFR > 0.2 M$_{\odot}$ yr$^{-1}$, see also Krogager et al. (2020). It is clearly seen that quasar-DLA sightlines on average probe the outer regions of their host galaxies. The impact parameters of quasar-DLAs are one order of magnitude larger than those of GRB-DLAs. We further observe that the impact parameters of quasar-DLAs are underpredicted by ~20 per cent in our model. This slight offset, however, is due to the complications of detecting the galaxies at small projected separations from the bright background quasar (less than ~ 5 kpc). Though we caution again that the impact parameter sample may further be biased by inhomogeneous observations often targeting metal-rich absorbers (Fynbo et al. 2010; Krogager et al. 2017) which in our model would also be larger and are thus more likely to have larger impact parameters.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Cumulative distributions of metallicity (top) and $N_{H1}$ (bottom) for the GRB metal sample shown by the solid black line and for quasar DLAs (Q16, Quiret et al. 2016) shown by the dashed grey line. We show the model distributions for GRBs in the thick solid, red line, and for quasar DLAs in the dashed, blue line. For comparison, we show the model distributions of quasar CNM-DLAs (i.e. with C1 or H2 absorption) as the thin, light blue line and the observed distribution of $N_{H1}$ for quasar CNM-DLAs (KN20, Krogager & Noterdaeme 2020) as the thick, light blue stepwise line. We do not show the metallicity distribution for the KN20 sample as this sample is biased towards high metallicity by being mostly selected based on C1 which favours high metallicity (Ledoux et al. 2015).</figcaption></figure>

## 5 SUMMARY

In this work, we study how GRB-DLAs probe their host galaxies using a statistical modelling approach. We compare observed distributions of $N_{H1}$, metallicity, $A(V)$, and impact parameters to the statistical model by Krogager & Noterdaeme (2020). We take into account a geometric correction since GRBs arise from within the volume of gas instead of piercing the full volume as a quasar sightline would. We find that this ‘baseline model’ provides an acceptable fit to the metallicity and impact parameter distributions of GRB-DLAs. These are the more fundamental quantities related to the host-galaxy luminosity (or mass). Yet, the ‘baseline model’ does not reproduce the observed $N_{H1}$ nor $A(V)$ distributions. We therefore include one additional free parameter to model the overdensity of neutral gas observed in GRB-DLAs (Heintz et al. 2019a). Lastly, we change the implemented dust bias prescription, which was tailored to model quasar-DLAs.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Cumulative distributions of impact parameters for the GRB sample by Lyman et al. (2017) shown by the solid black histogram and for quasar-DLAs (MC20, Møller & Christensen 2020) shown by the dashed grey histogram. We show the model distributions for GRBs in the solid red line (restricted in luminosity and non-obscured bursts, see Section 3.3). For quasar-DLAs, we show the model as the dashed blue line (restricted in SFR based on UV luminosity to mimic observational limitations at these high redshifts, see also Krogager et al. (2020).</figcaption></figure>

We fit the dust-bias parameter and the H I overdensity parameter simultaneously to a sample of GRBs detected in the redshift range $2 < z < 3.5$ (Section 2). From this fit, we find that bursts with $A(V) > 0.44 \pm 0.05$ mag are too heavily obscured to allow optical spectroscopy needed to determine the metallicity.<fnref n="4" /> The overdensity of neutral gas along the GRB sightline is found to be $\Delta \log(N_{H1}) = 1.06 \pm 0.09$ (a factor of 9–14), see also Fig. 8.

The overdensity of H I is interpreted as a temporal selection effect, namely that GRBs, due to their massive stellar progenitors, select galaxies in an early stage of a star formation event or a starburst. We speculate that galaxies in such a stage of a starburst may host many such massive and dense regions which would increase the column of gas along the line of sight (Hatsukade et al. 2020). Furthermore, local starburst galaxies exhibit an excess surface density of neutral gas compared to more regular galaxies (e.g. Kennicutt 1998; Bigiel et al. 2008). A similar effect could thus plausibly be expected for high-redshift galaxies as well. This interpretation is also qualitatively in agreement with the claim that GRB host galaxies have higher specific SFRs on average (Christensen et al. 2004; Perley et al. 2013; Björnsson 2019). We find it less convincing that the overdensity of neutral gas is caused by gas near the progenitor itself, unless the lifetime of the progenitor is < 3 Myr and the majority of the progenitors are able to re-accrete significant amounts of their birth clouds as suggested by Jeffreson et al. (2024).

Based on the acceptable agreement between metallicity and impact parameters from the model, both before and after fitting the $N_{H1}$ and $A(V)$ distributions, we conclude that GRB-DLAs and quasar-DLAs with CNM tracers sample the luminosity function of star-forming galaxies in the same way. We further test the metallicity thresholds put forward in the literature and find that at these high redshifts ($z > 2$), such metallicity thresholds do not significantly

affect the GRB population due to the overall lower metallicities at early times in the Universe. The host galaxies of GRB-DLAs are therefore representative of the underlying population in terms of metallicity and luminosity (or stellar mass), yet still susceptible to biases due to optical dust obscuration.

We quantify the effect of dust obscuration and compare the fraction of so-called dark bursts (with $A(V) > 1$ mag) to our model predictions. Observations infer a rather uncertain range of dark bursts ranging from 10 to 40 per cent (Jakobsson et al. 2004; Melandri et al. 2012). Based on our model, we expect 29 per cent of bursts to have $A(V) > 1$ mag.

Under the assumption that GRBs arise from the inner regions of galaxies where cold neutral gas should be present, we would naively expect that all GRB sightlines should show tracers of the CNM (C I or H$_2$). We discuss the detection rate of such CNM tracers and find that the observed fraction of 30–40 per cent (Bolmer et al. 2019; Heintz et al. 2019a) is consistent with our model, see Section 4.3.

Lastly, we compare the observed metallicity distribution of quasar- and GRB-DLAs to the expected distributions from our model framework. Quasar-DLAs have slightly lower metallicities on average due to the metallicity gradient and much larger impact parameters (Fig. 9). The GRB-DLAs have metallicities that are 0.2 dex higher on average, consistent with Arabsalmani et al. (2015), when taking into account the redshift evolution in the mass–metallicity relation (Møller et al. 2013). In contrast, the quasar-DLAs have significantly lower $N_{H1}$ on average, and even when comparing to the subset of quasar-DLAs with CNM gas that probe similar small impact parameters as GRB-DLAs, the median $N_{H1}$ differs by almost an order of magnitude. This bolsters our conclusion that the observed excess of $N_{H1}$ in GRB-DLAs is not simply due to the higher gas density at smaller galactic radii (see Figs 8 and 9).

## ACKNOWLEDGEMENTS

The authors would like to thank the referee for the very constructive and thorough review of this work. JKK would like to thank the ESO Visitor Programme for the financial support to visit ESO Garching in April 2023 which allowed us to focus on this project and make significant progress on the article. JKK and ADC acknowledge support from the Swiss National Science Foundation under grant 185692. KEH acknowledges support from the Carlsberg Foundation Reintegration Fellowship Grant CF21-0103. LBC was supported by the Independent Research Fund Denmark (DFF 2032-00071). SJ was supported by Harvard University via the ITC Fellowship. JPUF was supported by the Independent Research Fund Denmark (DFF-4090-00079) and thanks the Carlsberg Foundation for support. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant no. 140.

## DATA AVAILABILITY

All data analysed in this article are collected from literature and are referenced throughout the text. The main data underlying this analysis are available in Table 1. Additional data are available in the refereed works cited in this article.

## REFERENCES

* Arabsalmani M., Møller P., Fynbo J. P. U., Christensen L., Freudling W., Savaglio S., Zafar T., 2015, MNRAS, 446, 990

* Atek H., Richard J., Kneib J.-P., Schaerer D., 2018, MNRAS, 479, 5184

* Balashev S. A., Noterdaeme P., 2018, MNRAS, 478, L7

<footnote n="4">⁴The parametrization of the dust obscuration is of course overly simplified, as metallicity measurements could still be possible at higher levels of obscuration depending on the brightness of the GRB. None the less, the parameter is useful for the statistical comparison.</footnote>
Beniamini P., Piran T., 2024, ApJ, 966, 17

Bigiel F., Blitz L., 2012, ApJ, 756, 183

Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2002, AJ, 136, 2846

Björnsson G., 2019, ApJ, 887, 219

Blanchard P. K., Berger E., Fong W.-f., 2016, ApJ, 817, 144

Blitz L., Rosolowsky E., 2006, ApJ, 650, 933

Bloom J. S. et al., 2002, ApJ, 572, L45

Boissé P., Bergeron J., 2019, A&A, 622, A140

Boissier S., Prantzos N., 2001, MNRAS, 325, 321

Bolmer J. et al., 2019, A&A, 623, A43

Bouwens R. J., Illingworth G., Ellis R. S., Oesch P., Stefanon M., 2022, ApJ, 940, 55

Chabrier G., 2003, PASP, 115, 763

Chrimes A. A. et al., 2019, MNRAS, 486, 3105

Christensen L., Hjorth J., Gorosabel J., 2004, A&A, 425, 913

Christensen L. et al., 2012, MNRAS, 427, 1973

Christensen L., Møller P., Fynbo J. P. U., Zafar T., 2014, MNRAS, 445, 225

Corre D. et al., 2018, A&A, 617, A141

Covino S. et al., 2013, MNRAS, 432, 1231

Cucchiara A., Fumagalli M., Rafelski M., Kocevski D., Prochaska J. X., Cooke R. J., Becker G. D., 2015, ApJ, 804, 51

D'Elia V., Campana S., Covino S., D'Avanzo P., Piranomonte S., Tagliaferri G., 2011, MNRAS, 418, 680

De Cia A. et al., 2012, A&A, 545, A64

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

De Cia A., Jenkins E. B., Fox A. J., Ledoux C., Ramburuth-Hurt T., Konstantopoulou C., Petitjean P., Krogager J.-K., 2021, Nature, 597, 206

de Ugarte Postigo A. et al., 2024, preprint (arXiv:2406.16726)

Draine B. T., Hao L., 2002, ApJ, 569, 780

Elmegreen B. G., Parravano A., 1994, ApJ, 435, L121

Ferrière K. M., 2001, Rev. Mod. Phys., 73, 1031

Field G. B., Goldsmith D. W., Habing H. J., 1969, ApJ, 155, L149

Friis M. et al., 2015, MNRAS, 451, 167

Fruchter A. S. et al., 2006, Nature, 441, 463

Fynbo J. U., Møller P., Warren S. J., 1999, MNRAS, 305, 849

Fynbo J. U., Burud I., Møller P., 2000, A&A, 358, 88

Fynbo J. U. et al., 2001, A&A, 373, 796

Fynbo J. P. U., Prochaska J. X., Sommer-Larsen J., Dessauges-Zavadsky M., Møller P., 2008, ApJ, 683, 321

Fynbo J. P. U. et al., 2009, ApJS, 185, 526

Fynbo J. P. U. et al., 2010, MNRAS, 408, 2128

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

Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279

Greiner J. et al., 2011, A&A, 526, A30

Hainline K. N., Shapley A. E., Kornei K. A., Pettini M., Buckley-Geer E., Allam S. S., Tucker D. L., 2009, ApJ, 701, 52

Hartoog O. E. et al., 2013, MNRAS, 430, 2739

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

Hatsukade B., Ohta K., Hashimoto T., Kohno K., Nakanishi K., Niino Y., Tamura Y., 2020, ApJ, 892, 42

Heintz K. E. et al., 2017, A&A, 601, A83

Heintz K. E. et al., 2019a, A&A, 621, A20

Heintz K. E. et al., 2019b, A&A, 629, A131

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

Hjorth J. et al., 2012, ApJ, 756, 187

Hogg D. W., Fruchter A. S., 1999, ApJ, 520, 54

Jakobsson P., Hjorth J., Fynbo J. P. U., Watson D., Pedersen K., Björnsson G., Gorosabel J., 2004, ApJ, 617, L21

Jakobsson P. et al., 2005, MNRAS, 362, 245

Jakobsson P. et al., 2006, A&A, 460, L13

Japelj J. et al., 2016, A&A, 590, A129

Jeffreson S. M. R., Semenov V. A., Krumholz M. R., 2024, MNRAS, 527, 7093

Jorgenson R. A., Wolfe A. M., Prochaska J. X., 2010, ApJ, 722, 460

Kennicutt R. C., Jr, 1998, ARA&A, 36, 189

Kimm T., Bieri R., Geen S., Rosdahl J., Blaizot J., Michel-Dansac L., Garel T., 2022, ApJS, 259, 21

Kistler M. D., Yüksel H., Beacom J. F., Hopkins A. M., Wyithe J. S. B., 2009, ApJ, 705, L104

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

Krogager J. K. et al., 2018, A&A, 619, A142

Krogager J.-K., Noterdaeme P., 2020, A&A, 644, L6

Krogager J.-K., Møller P., Fynbo J. P. U., Noterdaeme P., 2017, MNRAS, 469, 2959

Krogager J.-K., Fynbo J. P. U., Møller P., Noterdaeme P., Heintz K. E., Pettini M., 2019, MNRAS, 486, 4377

Krogager J.-K., Møller P., Christensen L. B., Noterdaeme P., Fynbo J. P. U., Freuding W., 2020, MNRAS, 495, 3014

Krühler T. et al., 2011, A&A, 534, A108

Ledoux C., Petitjean P., Srianand R., 2003, MNRAS, 346, 209

Ledoux C., Petitjean P., Fynbo J. P. U., Møller P., Srianand R., 2006, A&A, 457, 71

Ledoux C., Vreeswijk P. M., Smette A., Fox A. J., Petitjean P., Ellison S. L., Fynbo J. P. U., Savaglio S., 2009, A&A, 506, 661

Ledoux C., Noterdaeme P., Petitjean P., Srianand R., 2015, A&A, 580, A8

Lyman J. D. et al., 2017, MNRAS, 467, 1795

McGuire J. T. W. et al., 2016, ApJ, 825, 135

Melandri A. et al., 2012, MNRAS, 421, 1265

Metha B., Trenti M., 2020, MNRAS, 495, 266

Møller P., Christensen L., 2020, MNRAS, 492, 4805

Møller P., Fynbo J. P. U., Ledoux C., Nilsson K. K., 2013, MNRAS, 430, 2680

Natarajan P. et al., 1997, New Astron., 2, 471

Noterdaeme P., Petitjean P., Pâris I., Cai Z., Finley H., Ge J., Pieri M. M., York D. G., 2014, A&A, 566, A24

Palmerio J. T. et al., 2019, A&A, 623, A26

Perley D. A. et al., 2013, ApJ, 778, 128

Perley D. A. et al., 2016, ApJ, 817, 8

Pettini M., Steidel C. C., Adelberger K. L., Dickinson M., Giavalisco M., 2000, ApJ, 528, 96

Piranomonte S. et al., 2008, A&A, 492, 775

Planck Collaboration XIII, 2016, A&A, 594, A13

Porciani C., Madau P., 2001, ApJ, 548, 522

Prochaska J. X., Chen H.-W., Dessauges-Zavadsky M., Bloom J. S., 2007, ApJ, 666, 267

Quiret S. et al., 2016, MNRAS, 458, 4074

Ranjan A., Noterdaeme P., Krogager J. K., Petitjean P., Srianand R., Balashev S. A., Gupta N., Ledoux C., 2020, A&A, 633, A125

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

Rhodin N. H. P., Christensen L., Møller P., Zafar T., Fynbo J. P. U., 2018, A&A, 618, A129

Robertson B. E., Ellis R. S., 2012, ApJ, 744, 95

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

Salvaterra R. et al., 2009, Nature, 461, 1258

Santini P. et al., 2022, ApJ, 940, 135

Savaglio S., Fall S. M., Fiore F., 2003, ApJ, 585, 638

Schady P., Savaglio S., Krühler T., Greiner J., Rau A., 2011, A&A, 525, A113

Schmidt M., 1959, ApJ, 129, 243

Schroeder G. et al., 2022, ApJ, 940, 53

Stanek K. Z. et al., 2003, ApJ, 591, L17

Stark D. P. et al., 2013, MNRAS, 436, 1040

Stevens A. R. H., Diemer B., Lagos C. d. P., Nelson D., Obreschkow D., Wang J., Marinacci F., 2019, MNRAS, 490, 96

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

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

Thöne C. C. et al., 2024, A&A, 690, A66

Tyson N. D., Scalo J. M., 1988, ApJ, 329, 618

van der Horst A. J., Kouveliotou C., Gehrels N., Rol E., Wijers R. A. M. J., Cannizzo J. K., Racusin J., Burrows D. N., 2009, ApJ, 699, 1087

Vergani S. D. et al., 2015, A&A, 581, A102

Vreeswijk P. M. et al., 2004, A&A, 419, 927

Vreeswijk P. M. et al., 2007, A&A, 468, 83
Vreeswijk P. M. et al., 2013, A&A, 549, A22

Watson D., Hjorth J., Fynbo J. P. U., Jakobsson P., Foley S., Sollerman J., Wijers R. A. M. J., 2007, ApJ, 660, L101

Watson D. et al., 2013, ApJ, 768, 23

Whalen D., Prochaska J. X., Heger A., Tumlinson J., 2008, ApJ, 682, 1114

Wijers R. A. M. J., Bloom J. S., Bagla J. S., Natarajan P., 1998, MNRAS, 294, L13

Wiklind T., Combes F., 1995, A&A, 299, 382

Wiklind T., Combes F., Kanekar N., 2018, ApJ, 864, 73

Wiseman P., Schady P., Bolmer J., Krühler T., Yates R. M., Greiner J., Fynbo J. P. U., 2017, A&A, 599, A24

Wolfe A. M., Turnshek D. A., Smith H. E., Cohen R. D., 1986, ApJS, 61, 249

Wolfire M. G., Hollenbach D., McKee C. F., Tielens A. G. G. M., Bakes E. L. O., 1995, ApJ, 443, 152

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278

Woosley S. E., Heger A., Weaver T. A., 2002, Rev. Mod. Phys., 74, 1015

Zafar T., Møller P., 2019, MNRAS, 482, 2731

Zafar T., Péroux C., Popping A., Milliard B., Deharveng J. M., Frank S., 2013, A&A, 556, A141

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