# MUSEQuBES: Connecting H i Absorption with Ly α Emitters at $z \approx 3.3$

Eshita Banerjee, Sowgat Muzahid, Joop Schaye, Jérémy Blaizot, Nicolas

Bouché, Sebastiano Cantalupo, Sean Johnson, Jorryt Matthee, Anne

Verhamme

► To cite this version:

Eshita Banerjee, Sowgat Muzahid, Joop Schaye, Jérémy Blaizot, Nicolas Bouché, et al. MUSEQuBES: Connecting H i Absorption with Ly α Emitters at $z \approx 3.3$. The Astrophysical Journal, 2025, 980 (2), pp.171. 10.3847/1538-4357/ada7e9. [hal-05371869](https://hal.science/hal-05371869v1)

HAL Id: hal-05371869

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

Submitted on 19 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

MUSEQuBES: Connecting H I Absorption with Ly$\alpha$ Emitters at $z \approx 3.3$

Eshita Banerjee<fnref n="1" />, Sowgat Muzahid<fnref n="1" />, Joop Schaye<fnref n="2" />, Jérémy Blaizot<fnref n="3" />, Nicolas Bouché<fnref n="3" />, Sebastiano Cantalupo<fnref n="4" />, Sean D. Johnson<fnref n="5" />, Jorryt Matthee<fnref n="6" />, and Anne Verhamme<fnref n="7" />

Received 2024 July 23; revised 2025 January 6; accepted 2025 January 7; published 2025 February 12

Abstract

We present a comprehensive analysis of H I absorption around 96 Ly$\alpha$ emitters (LAEs) at $z \approx 3.3$ (median Ly$\alpha$ luminosity $\approx 10^{42}$ erg s$^{-1}$). These LAEs were identified within eight MUSE fields, each $1' \times 1'$ on the sky and centered on a bright background quasar, as part of the MUSEQuBES survey. Using Voigt profile fitting for all H I absorbers detected within $\pm 500$ km s$^{-1}$ of these LAEs, we compiled a catalog of 800 H I absorption components. Our analysis shows that H I absorption is enhanced near the LAEs compared to the intergalactic medium. However, no trend is found between the column densities of H I absorbers and their impact parameters from the LAEs (spanning $\approx 54-260$ pkpc). Additionally, all galaxies associated with Lyman-limit systems have impact parameters $> 50$ pkpc from the quasar sightlines, suggesting that true absorber hosts may be too faint to detect. The LAEs show an overall H I covering fraction ($f_c(H I)$) of $\approx 88\%$ for a threshold $\log_{10}(N/\text{cm}^{-2})(\text{H I}) = 15$. Notably, at the same threshold, the LAEs in pairs/groups exhibit a $100\%$ H I covering fraction out to $\approx 250$ pkpc. In contrast, isolated LAEs consistently show a lower $f_c(H I)$ of $\approx 80\%$. This environmental influence on $f_c(H I)$ is also evident up to $\approx 300$ km s$^{-1}$ in differential bins of line-of-sight velocity. We find an anticorrelation between $f_c(H I)$ and the equivalent width of rest-frame Ly$\alpha$ emission (EW$_0$). Based on the Ly$\alpha$ shell model, this could imply that gas-rich galaxies tend to reside in gas-rich environments or that the LAEs with higher EW$_0$ are more efficient at ionizing their surrounding medium.

*Unified Astronomy Thesaurus concepts: Circumgalactic medium (1879); Galaxy evolution (594); Ly$\alpha$ galaxies (978); Quasar absorption line spectroscopy (1317)*

*Materials only available in the online version of record: machine-readable table*

1. Introduction

The circumgalactic medium (CGM) is a dynamic, complex, multiphase gaseous region, extending out to a few hundred kiloparsecs and enveloping galaxies (see the review by J. Tumlinson et al. 2017). It contains the imprints of crucial processes such as galactic-scale winds and the accretion of intergalactic baryons via the “hot mode”/“cold mode” along cosmic filaments (D. Kereš et al. 2005; F. van de Voort et al. 2011b). This cyclic interplay of baryons in the CGM, often referred to as the “baryon cycle,” is a key factor governing galaxy evolution (e.g., C. Péroux & J. C. Howk 2020). It is widely believed that these processes played a crucial role in determining the evolution of the cosmic star formation rate density (SFRD), which shows a peak at a redshift of $z \approx 2-3$ and declines by almost an order of magnitude at both higher and lower redshifts (e.g., F. van de Voort et al. 2011a; P. Madau & M. Dickinson 2014). Given that hydrogen accounts for the majority of the baryonic matter in the Universe, it has become a vital tool for probing the cool gas that acts as the fuel for star formation.

Probing the cool gas in the CGM in emission is challenging owing to its low density ($n_H \sim 10^{-3}$ cm$^{-3}$ at $z \approx 3$). State-of-the-art integral field spectrographs (IFSs) such as VLT/Multi-Unit Spectroscopic Explorer (MUSE; R. Bacon et al. 2010) and Keck Cosmic Web Imager (KCWI; P. Morrissey et al. 2018) have enabled the exploration of the diffuse medium surrounding galaxies in Ly$\alpha$ emission within a few tens of kiloparsecs from normal galaxies (e.g., D. K. Erb et al. 2018; L. Wisotzki et al. 2018; Y. Guo et al. 2024a, 2024b) and a few hundred physical kiloparsecs around active galaxies (e.g., S. Cantalupo et al. 2014; E. Borisova et al. 2016; F. Arrigoni Battaia et al. 2019). Extended Ly$\alpha$ and metal line emission is also detected in several individual cases and in stacks (e.g., C. C. Steidel et al. 2011; J. Zabl et al. 2021; S. D. Johnson et al. 2022, 2024; F. Leclercq et al. 2022; R. Dutta et al. 2023; Y. Guo et al. 2023; Z. Liu et al. 2024). However, the absorption spectra from bright background sources, such as quasars, remain the most sensitive tool to probe this elusive medium.

Historically, the relation between galaxies and the neutral gas surrounding them was studied by identifying host galaxies of strong H I absorbers, such as Lyman-limit systems (LLSs; $N_{H I} > 10^{17.2}$ cm$^{-2}$; e.g., E. K. Lofthouse et al. 2023), and damped Ly$\alpha$ absorbers (DLAs; $N_{H I} \gtrsim 10^{20.3}$ cm$^{-2}$; e.g., C. Péroux et al. 2012; J. K. Krogager et al. 2017; R. Mackenzie et al. 2019), seen in quasar spectra. Several investigations have primarily focused on high-metallicity DLAs, driven by the

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

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

<footnote n="1">¹IUCAA, Post Bag 04, Ganeshkhind, Pune, 411007, India; eshitaban18@iucaa.in, sowgat@iucaa.in</footnote>

<footnote n="2">²Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 AA Leiden, The Netherlands</footnote>

<footnote n="3">³CNRS, Centre de Recherche Astrophysique de Lyon (CRAL), Université Lyon 1, Saint-Genis-Laval, France</footnote>

<footnote n="4">⁴Department of Physics, University of Milan Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="5">⁵Department of Astronomy, University of Michigan, 1085 S. University, Ann Arbor, MI 48109, USA</footnote>

<footnote n="6">⁶Institute of Science and Technology Austria (ISTA), Am Campus 1, Klosterneuburg, Austria</footnote>

<footnote n="7">⁷Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, CH-1290 Versoix, Switzerland</footnote>
presumption that these DLAs possess more luminous galaxy counterparts, owing to the underlying mass–metallicity relation (e.g., M. Neeleman et al. 2013). These studies utilized facilities such as VLT/X-Shooter (e.g., J. K. Krogager et al. 2017), VLT/SINFONI (e.g., C. Péroux et al. 2012), and the Atacama Large Millimeter/submillimeter Array (e.g., M. Neeleman et al. 2017, 2019), revealing close associations (often within a few tens of kiloparsecs) between DLAs and their host galaxies. However, recent studies with state-of-the-art IFSs such as the MUSE with a field of view (FOV) of $1' \times 1'$ revealed that galaxies associated with the DLAs/LLSs are often located at larger impact parameters (>50 pkpc; see R. Mackenzie et al. 2019; E. K. Lofthouse et al. 2023). This phenomenon could potentially be attributed to the fact that the true host of the absorber remains undetected owing to its relative faintness (e.g., A. Rahmati & J. Schaye 2014).

Recently, using MUSE observations, E. K. Lofthouse et al. (2023) reported a detection rate of ≈80% for Lyα-emitting galaxies (LAEs) around absorbers with $N_{H\,I} > 10^{16.5} \text{ cm}^{-2}$, with their impact parameter ranging from 18 to 275 pkpc. Based on the higher number density of these LAEs compared to field galaxies, they concluded that LAEs are significantly clustered around absorbers with high H I column density. Besides, they found that LAEs in groups exhibit a three times higher covering factor of optically thick gas than isolated systems. The absence of a substantial correlation between emission and absorption characteristics led them to conclude that the CGM is patchy and inhomogeneous. They argued that very large samples of quasar–galaxy pairs are required to reveal the signatures of the baryon cycle. While H I absorbers with high column density serve as reliable indicators of galaxy overdensities, it is crucial to investigate the connection(s) between galaxies and gas in an “absorption-blind” manner.

G. C. Rudie et al. (2012), as a part of the Keck Baryonic Structure Survey (KBSS; see also K. L. Adelberger et al. 2005; C. C. Steidel et al. 2010; O. Rakic et al. 2012; G. C. Rudie et al. 2013, 2019; M. L. Turner et al. 2014; Y. Chen et al. 2020), performed a complete Voigt profile decomposition of all the H I absorbers down to $N_{H\,I} \sim 10^{12} \text{ cm}^{-2}$ around a sample of Lyman-break galaxies (LBGs) and found a strong enhancement of the H I absorption relative to randomly located regions. This enhancement extends out to ≈300 pkpc in the transverse direction and ±700 km s⁻¹ along the line of sight (LOS). The halo masses of the LBGs in their sample are ~$10^{12} M_{\odot}$. Clearly, these are massive structures and represent only the “tip of the iceberg” of the high-z galaxy population. It is important to investigate the connections between “typical” high-z galaxies and neutral gas around them. The Lyα-emitting galaxies with halo masses of ~$10^{11} M_{\odot}$ are a good candidate for low-mass, high-z galaxies (Y. Herrero Alonso et al. 2021). Such low-mass galaxies are more susceptible to galactic-scale outflows, owing to their shallow gravitational potentials. Mapping the distribution of gas and metals around LAEs is therefore crucial to exploring the galaxy–CGM connections in the unexamined low-mass regime.

The MUSE Quasar-field Blind Emitters Survey (MUSE-QuBES; S. Muzahid et al. 2020, 2021; E. Banerjee et al. 2023) is dedicated to investigating the correlation between gas properties and low-mass galaxies at redshifts $z \approx 3–4$. The MUSEQuBES galaxies, identified as LAEs, have a median stellar mass $M_* \approx 10^{8.9} M_{\odot}$. The MUSEQuBES galaxy sample is obtained without any prior knowledge of the presence or

absence of gas (it is “absorption-blind”). S. Muzahid et al. (2021) carried out a spectral stacking analysis of the CGM absorption around $96 z \approx 3.3$ LAEs and found that the median rest-frame equivalent width (W) of H I absorption correlates with the star formation rate (SFR) and environment of the LAEs. For the same LAE sample, E. Banerjee et al. (2023) studied the connection between galaxies and C IV absorbers at $z \approx 3.3$. In continuation of the MUSEQuBES survey, here we present a thorough Voigt profile decomposition of H I absorbers detected within ±500 km s⁻¹ LOS velocity windows around each of the 96 MUSEQuBES LAEs.

There is no consensus on the extent of the CGM from its host galaxy. Some low-z surveys suggest that it extends out to the virial radius $R_{vir}$ of galaxies (e.g., J. Tumlinson et al. 2017; S. Dutta et al. 2024). In our current work, we are probing a large distance of up to ≈320 pkpc from these galaxies, much larger than their typical $R_{vir}$ (≈42 pkpc). Therefore, it is more reasonable to say that we are probing the extended gaseous medium associated with these LAEs. However, for simplicity, here we will call it the CGM.

This paper is organized as follows: in Section 2, we provide a concise overview of the quasar and galaxy data employed in this paper, with particular emphasis on the preparation of the H I absorption line catalog. In Section 3, we present our analysis and results. In Section 4, we discuss and interpret our results. Finally, a brief summary is presented in Section 5.

Throughout this study, we adopt a flat $\Lambda$CDM cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_{\Lambda} = 0.7$, and all distances given are in physical kiloparsecs (hereafter, pkpc) unless specified otherwise.

## 2. Data and Measurements

The quasar and galaxy data presented in this study are drawn from the MUSEQuBES survey. Our survey contains excellent-quality optical spectra of eight UV-bright background quasars, obtained using the VLT/UVES and Keck/HIRES spectrographs with a spectral resolution of $R \approx 45,000$. All spectra exhibit a signal-to-noise ratio (SNR) $\gtrsim 30$ per pixel within the Lyα forest, except for quasar BRI1108-07, which has an SNR of ≈15. Further details on the data and quasar continuum fitting can be found in section 2 of S. Muzahid et al. (2021). The subsequent subsections contain brief discussions on the galaxy data and details of absorption line measurements.

### 2.1. LAE Sample

A total of 50 hr guaranteed-time MUSE observations (PI: J. Schaye) of eight quasar fields allowed us to identify 96 LAEs at $2.9 < z < 3.8$, spanning an impact parameter ($\rho$) range of 16 pkpc $< \rho < 320$ pkpc. These LAEs were detected solely based on their Lyα emission lines, which typically exhibit offsets of hundreds of km s⁻¹ from the galaxies’ systemic redshifts (see C. C. Steidel et al. 2010; O. Rakic et al. 2011; T. Shibuya et al. 2014; A. Verhamme et al. 2018). An empirical relation between this offset and the FWHM of the Lyα emission peak was provided for our sample by S. Muzahid et al. (2020), i.e., $V_{offset} = 0.89 \times \text{FWHM} - 58 \text{ km s}^{-1}$. The systemic redshifts of the LAEs in this study are derived using this relation.

The median stellar mass of our sample is $M_* = 10^{8.9} M_{\odot}$ and median SFR is ≈1.2 $M_{\odot}$ yr⁻¹. Note that these SFR values are not corrected for dust and are derived from the measured values
<table><caption>**Table 1**
Best-fitting Parameters Obtained from Voigt Profile Fitting of the H I Absorbers.</caption><thead><tr><th>QSO</th><th>z</th><th>Error in z</th><th>b<br>(km s<sup>-1</sup>)</th><th>Error in b<br>(km s<sup>-1</sup>)</th><th>log<sub>10</sub><br>(N/cm<sup>-2</sup>)</th><th>Error in log<sub>10</sub><br>(N/cm<sup>-2</sup>)</th><th>Detection ID</th></tr></thead><tbody><tr><td>Q1422+23</td><td>2.958128</td><td>2.8e-06</td><td>23.00</td><td>0.35</td><td>13.52</td><td>0.01</td><td>0</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>BRI1108-07</td><td>3.584494</td><td>8.8e-06</td><td>28.95</td><td>3.46</td><td>17.25</td><td>0.80</td><td>1</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Note. Detection ID = 1: lower limit on $N_{H\alpha}$; ID = 0: well-constrained $N_{H\alpha}$.

(This table is available in its entirety in machine-readable form in the online article.)

of UV luminosity density using the local calibration relation of R. C. J. Kennicutt (1998) corrected to the G. Chabrier (2003) initial mass function. The stellar masses of these LAEs are estimated from the measured SFR values, assuming the star-forming main-sequence relation of P. Behroozi et al. (2019). Although there are several assumptions involved in estimating these values, they are very close to the stellar mass and SFR of typical LAE populations at these redshifts, i.e., $M_* = 10^8-10^{8.5} M_{\odot}$ and an SFR of $\approx 2.3 M_{\odot} \text{yr}^{-1}$, as estimated based on broader spectral energy distribution fitting and rest-optical emission line studies (see Y. Ono et al. 2010; R. F. Trainor et al. 2015). Further, using the halo abundance matching relation from B. P. Moster et al. (2013), we estimated the halo mass ($M_{\text{vir}}$) and corresponding virial radius ($R_{\text{vir}}$) for a subset of 39 out of 96 LAEs, where the SFR could be measured from their UV continuum. For the remaining LAEs, the UV continuum was either undetected with $>5\sigma$ confidence or blended with foreground sources (see S. Muzahid et al. 2020). The halo masses (virial radii) range from $10^{11.0}$ to $10^{11.7} M_{\odot}$ (33–59 pkpc) with a median value of $10^{11.3} M_{\odot}$ (42 pkpc). The 5th and 95th percentiles of the $Ly\alpha$ luminosity for the MUSEQuBES LAEs are $\log_{10}(L_{Ly\alpha}/\text{erg s}^{-1}) = 41.5$ and 42.7; the corresponding $Ly\alpha$ line fluxes are $3 \times 10^{-18}$ and $4 \times 10^{-17} \text{ erg cm}^{-2} \text{ s}^{-1}$. Here we focus on the H I absorbers located within $\pm 500 \text{ km s}^{-1}$ of the corrected LAE redshifts. Further details regarding the MUSE data analysis and galaxy properties can be found in S. Muzahid et al. (2021).

## 2.2. H I Absorption Line Measurements

The accurate measurement of H I column densities in the spectra of high-z quasars is challenging due to (a) the blending of H I absorption with other H I or metal lines present in the $Ly\alpha$ forest and (b) the saturation of strong H I absorption components. To overcome these challenges, it is imperative to simultaneously fit all the higher-order Lyman-series lines. This approach not only facilitates the accurate measurement of column density but also aids in effectively identifying contamination, thereby simplifying its removal.

In our study, we conducted a detailed Voigt profile decomposition of H I absorption detected within $\pm 500 \text{ km s}^{-1}$ of the 96 MUSEQuBES LAEs. Note that a given absorber may be associated with multiple LAEs. For Voigt profile decomposition, we used the Vpfit code (R. F. Carswell & J. K. Webb 2014), which works based on $\chi^2$-minimization. Vpfit allowed us to simultaneously fit the available Lyman-series lines as well as the contaminating metal lines. We obtained the best-fitting values for the redshift, Doppler parameter ($b$), and column density ($N$) of the 800 associated individual H I components in the redshift range $2.914 \le z \le 3.823$ (median $z = 3.364$). The best-fitting

parameters are summarized in Table 1. To ensure the reliability of our measurements, we manually examined the fitted profiles and repeated the fitting process until a satisfactory fit was achieved, in addition to values of $\chi^2$ and the Akaike information criterion, corrected for small sample sizes (AICc), of the model.

Due to blending, the presence of higher-order Lyman-series lines does not always guarantee a well-resolved structure for heavily saturated absorption profiles. In such instances, we had to rely mostly on the saturated $Ly\alpha$ line. Previous studies (e.g., J. M. Shull et al. 2000; A. Songaila 2001; C. W. Danforth et al. 2010) have indicated that, when using a single-line measurement of H I, the fitted profile tends to be $\approx 1.5$ times broader than what is estimated through curve-of-growth measurements derived from higher-order Lyman-series lines. An overestimate of the Doppler width leads to an underestimate of the column density. As a result, the estimate of column density for these components represents a lower limit. Within a velocity range of $\pm 500 \text{ km s}^{-1}$ from the LAEs, we have identified 66 components (comprising 8% of the total number of components) that are treated as lower limits in our analysis. These components either suffer from saturation without the availability of other Lyman-series lines or, if available, the true velocity component structure cannot be resolved due to blending. They are classified as detection ID = 1 absorbers in Table 1.

Our absorber catalog contains 48 weak, unsaturated $Ly\alpha$ absorption components, comprising 6% of the total, for which no higher-order lines are available. Since the lines are deemed to be unsaturated, obtaining accurate column densities is not problematic for these components (see also T.-S. Kim et al. 2021). Hence, these absorbers, along with the rest of the absorbers for which we can constrain $N_{H\alpha}$ using the unsaturated higher-order Lyman-series lines, are classified as ID = 0 components.

Figure 1 shows an example of our Voigt profile decomposition. Each panel displays the spectra (in black) in LOS velocity space along with the best-fitting profiles (in red) for an LAE toward the Q1422+23 sightline. A velocity of $0 \text{ km s}^{-1}$ corresponds to the systemic redshift of the LAE. The centroid of each H I component within $\pm 500 \text{ km s}^{-1}$ of the LAE (region marked by orange dashed vertical lines) is indicated by blue ticks. All components, except number 13, are classified as ID = 0. Owing to the saturation of all the identified higher-order lines, only a lower limit on $N_{H\alpha}$ can be derived for the component labeled 13, which is classified as an ID = 1 absorber.

Figure 2 shows a scatter plot of $b(H\alpha)$ versus $N(H\alpha)$ of all 800 H I components in our catalog with arrows indicating lower limits on $N(H\alpha)$. The distributions of these two parameters are shown in the side panels. The absorber sample spans the range
<figure><img src="image_4.png" /></figure>

**Figure 1.** H I absorption profiles associated with one of the LAEs toward Q1422+23. The normalized quasar spectra are represented in black, with the best-fitting Voigt profiles shown in red. The LOS velocity of 0 km s⁻¹ denotes the systemic redshift of the corresponding LAE. The ±500 km s⁻¹ velocity range is depicted by the orange dashed vertical lines. The centroids of individual Voigt profile components within this velocity window are indicated by blue ticks. The Lyman-series lines are plotted from top to bottom in descending order of wavelength, as mentioned in each panel. The column density is well constrained for the unsaturated Lyα components (i.e., 0–3, 5–9, and 14–17) and for the components with unsaturated high-order lines (i.e., 10–12). All of these components are classified as ID = 0 absorbers in Table 1. On the other hand, all the higher-order lines corresponding to component 13 are saturated. We therefore treat the column density returned by Vpfit for this component as a lower limit (i.e., ID = 1). The total N(H I) within ±500 km s⁻¹ of this LAE is log₁₀(N/cm⁻²) = 17.1, which is also treated as a lower limit. The unfitted absorptions are contaminating lines originating from other redshifts.

log₁₀(N/cm⁻²) = 11.7–20.3 with a mean (median) of 13.9 (13.7). For components with log₁₀(N/cm⁻²) < 16, the median error in column density is 0.05 dex, while for components with log₁₀(N/cm⁻²) > 16, it is 0.04 dex. The mean (median) value of the Doppler parameter is 34 km s⁻¹ (27 km s⁻¹). Only a small fraction of the components show very narrow (b ≈ 6.6 km s⁻¹; ≈0.6%) or very broad (b > 100 km s⁻¹; ≈3%) profiles. We verified that, except for a few cases, the components with large b-value arise in blends and are susceptible to continuum fitting uncertainties. We find that some of the narrow components can be attributed to unidentified low-z metal lines. Hence, they were excluded from our catalog.

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

**Figure 2.** Scatter plot of b vs. N for all the 800 H I absorption components detected within ±500 km s⁻¹ of the LAE systemic redshifts. The points marked with a rightward arrow indicate the components for which N(H I) cannot be measured due to blending and/or saturation. The side panels present the distributions of the N and b values, with the black dashed lines indicating their respective median values.

### 3. Analysis and Results

We will now proceed with the analysis of the overall properties of the H I absorbers, which are located within ±500 km s⁻¹ of the MUSEQuBES LAEs. This analysis will involve comparing the characteristics of these absorbers to the typical absorption features observed in the intergalactic medium (IGM) at a similar redshift.

#### 3.1. H I Column Density versus Impact Parameter

First, we focus on the distribution of H I-absorbing gas surrounding the LAEs. Figure 3 (left panel) shows the distribution of the total column density (ΣN_HI) calculated within a velocity range of ±500 km s⁻¹ from the LAEs as a function of the impact parameter (ρ) of the LAEs. The empty circles represent lower limits on N_HI for the LAEs that possess at least one ID = 1 component within ±500 km s⁻¹. The plot reveals a lack of any noticeable correlation between the total H I column density and ρ. This observation is further confirmed by a Spearman correlation test with rcorr = -0.01 and a p-value of 0.95. When we exclude the lower limits, the results of the Spearman correlation test yield rcorr = -0.07 and p-value = 0.54, indicating no significant correlation. In right panel, we show the same plot, but here the impact parameter is normalized by R_vir, considering only those LAEs with a continuum detected at >5σ confidence. Here, also, we do not see any significant correlation.

The mean value of total N_HI within ±500 km s⁻¹ of the LAEs is ≈10¹⁶.⁰ cm⁻² (blue dashed line in the left panel of Figure 3), which may be an underestimate due to the presence of several lower limits on N_HI. This is quite enhanced compared to the IGM, where the typical overdensity is ≈1. In the same plot, the green (red) dashed line represents N_HI = 10¹³.⁹ (10¹⁵.⁰) cm⁻², corresponding to overdensities of ≈1 (10) at z ≈ 3 (see J. Schaye 2001). Our observed values are generally higher, indicating that H I absorption is elevated near the LAEs compared to random regions in the IGM at similar redshifts. We obtained consistent results when we used a
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left: the total $N_{H1}$ ($\Sigma N_{H1}$) within $\pm 500$ km s$^{-1}$ of the LAEs against the impact parameter ($\rho$). The data points represented by empty circles with upward arrows denote lower limits on $N_{H1}$, indicating that these data points have some contribution from unresolved/blended H1 component(s). The average value of the total $N_{H1}$ within $\pm 500$ km s$^{-1}$ of the LAEs is $\log_{10}(N/cm^{-2}) > 16.0$ (indicated by the blue dashed line), where the inequality is due to the presence of unresolved lower limits in our data. The Spearman correlation results, mentioned in the plot, reveal the absence of any significant trend between $\Sigma N_{H1}$ and $\rho$. The green and red dashed lines in the plot represent $N_{H1}$ corresponding to overdensities of $\delta \approx 1$ and $\approx 10$ respectively, suggesting that there is an enhancement of neutral hydrogen in the vicinity of the LAEs. Right: similar to the left panel, but here the impact parameter is normalized by $R_{200}$ or $R_{vir}$. Only LAEs with a continuum detected at $>5\sigma$ are shown.</figcaption></figure>

smaller velocity window of $\pm 250$ km s$^{-1}$ to associate LAEs and H1 absorption.

Finally, another intriguing aspect of Figure 3 is that the majority of the absorbers with $\log_{10}(N/cm^{-2}) > 17.2$ are associated with LAEs with large impact parameters, in contrast to conventional wisdom. Considering that excessive noise from quasars can result in the nondetection of LAEs at smaller impact parameters, note that a total of 96 LAEs are detected within the 320 pkpc $\times$ 320 pkpc FOV. Assuming a random distribution of galaxies, we expect approximately two LAEs within 42 pkpc $\times$ 42 pkpc, where 42 pkpc is the typical virial radius of the MUSEQuBES LAEs. In fact, we detect three LAEs within 42 pkpc of the quasar sightline. This suggests that the nondetection of LAEs at smaller impact parameters is not because of some issue with subtraction of the point-spread function of the quasar. We discuss these systems further in Section 4.

### 3.2. H1 Covering Fraction around LAEs

The CGM is known to be an inhomogeneous and patchy medium, exhibiting variations in its physical and chemical properties across different regions. The covering fraction of atomic and ionic species within the CGM serves as a quantitative measure of this patchiness, providing valuable insight into the underlying processes that govern galaxy evolution. In this section, we focus on the H1 covering fraction, $f_c(H1)$, near the LAEs. For a given threshold column density, $N_{Th}$, it is defined as$^8$

$$f_c(\text{H1}) = \frac{n_{\text{Hit}}(N \ge N_{\text{Th}})}{n_{\text{Total}}}, \quad (1)$$

where $n_{Hit}$ is the number of LAEs hosting a total H1 column density $N_{H1} \ge N_{Th}$ in their CGM, measured within a specific

LOS velocity window (here, $\pm 500$ km s$^{-1}$) around the LAE redshifts, and $n_{Total}$ is the total number of LAEs for which the quasar spectra have sufficient SNR to detect absorption down to the threshold column density $N_{Th}$.

To determine the $3\sigma$ limiting column densities of our quasar spectra, we used the relation provided by U. Hellsten et al. (1998), following a similar approach to that described in Section 3.1.1 of E. Banerjee et al. (2023). These limits are calculated using the median SNR in the $\text{Ly}\alpha$ forest region for each spectrum, as reported in Table 2 of S. Muzahid et al. (2021). The undetected line is assumed to have a $b$ parameter of 35 km s$^{-1}$, which is the median value for the detected components. It is worth noting that the limiting column density is $\approx 10^{12.7}$ cm$^{-2}$ for all eight sightlines, which is lower than the threshold value we have used to calculate $f_c(H1)$ in this paper. Additionally, we excluded the LAEs for which the lower limit on the detected H1 is lower than $N_{Th}$, since it is not guaranteed that the total $N_{H1}$ for such cases is higher than the threshold value.

Figure 4 presents the H1 covering fraction for different threshold column densities. The filled data points represent the average $f_c(H1)$ measured within the impact parameter range $16 < \rho/\text{pkpc} < 320$ and an LOS velocity window of $\pm 500$ km s$^{-1}$. The error bars indicate the 68% Wilson-score confidence intervals. $f_c(H1)$ for low-column-density H1 absorbers (approximately $\log_{10}(N/cm^{-2}) \approx 14.5$) is close to unity. However, as the threshold column density increases, $f_c(H1)$ declines.

To compare these results with random regions, we calculated the H1 covering fractions for the IGM at similar redshifts using the following equation from M. C. Wilde et al. (2021):

$$f_{\text{random}} = 1 - \exp\left(-\left\langle \frac{dn}{dz} \right\rangle \Delta z\right) \quad (2)$$

where $\langle dn/dz \rangle$ is the average number of H1 absorbers with $N_{H1} > N_{Th}$ per unit redshift range, calculated from the best-fitting

<footnote>$^8$ In this study, we computed $f_c$ within a defined transverse distance, i.e., $f_c(\langle \rho \rangle)$, unless we are doing differential binning of impact parameter.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The average H I covering fraction for $\rho < 320$ pkpc for different threshold $N(HI)$ values. The red squares represent $f_c$ estimated using H I absorbers within $\pm 500$ km s$^{-1}$ of the LAEs, while the empty squares depict $f_c$ for the IGM (see text). The error bars on the red squares indicate the 68% Wilson-score confidence interval. The numbers adjacent to the squares correspond to the respective values of $n_{\text{Hit}}$ and $n_{\text{Total}}$ as in Equation (1). For all threshold column densities, $f_c(\text{HI})$ is consistently significantly higher near LAEs than in random locations. The line and the right y-axis represent the fractional difference in $f_c(\text{HI})$ with respect to random regions, which exhibits a gradual increase with the threshold column density.</figcaption></figure>

model to the column density distribution function (CDDF) from T.-S. Kim et al. (2021). Here, $\Delta z$ is the redshift range corresponding to a velocity window of $\pm 500$ km s$^{-1}$ at $z \approx 3.3$. The open squares in Figure 4 represent the covering fractions for the IGM. The H I covering fraction near LAEs is consistently higher than in random IGM locations. We obtained consistent results for the other CDDF models mentioned in the literature (e.g., P. Noterdaeme et al. 2009; G. Altay et al. 2011; G. C. Rudie et al. 2013).

The light blue line in the plot illustrates the fractional difference in $f_c(\text{HI})$, i.e., $(f_c(\text{HI}) - f_{\text{random}})/f_{\text{random}}$, with its values depicted on the opposite y-axis. This ratio gradually increases for higher threshold column densities, reflecting the greater overabundance of higher-column-density absorbers near the LAEs than in random regions.

We observe no notable changes in the H I covering fraction around LAEs across different redshifts for any threshold $N_{\text{HI}}$. This is expected given the small redshift range ($\Delta z/(1 + z_{\text{median}}) \approx 0.2$) for our sample. Consequently, we do not account for any redshift evolution of $f_c(\text{HI})$ near LAEs in our subsequent analysis.

### 3.3. Influence of Galaxy Environments on $f_c(HI)$

In previous sections, we explored the distribution of H I absorbers and noted their significant incidence around LAEs compared to the IGM. In this section, we will investigate how the environment of the LAEs affects the distribution of the neutral gas around them.

To examine this, we first identify LAEs that have at least one other companion LAE within $\pm 500$ km s$^{-1}$ in LOS velocity space within the MUSE field of view. We refer to them as galaxy “pairs/groups.” The remaining LAEs are considered “isolated.” Based on this simple classification, we find that 33 out of the total of 96 LAEs are part of a pair/group.

The left panel of Figure 5 shows the differential H I covering fraction against the impact parameter for a threshold

$\log_{10}(N/cm^{-2}) = 15.0$ and a velocity window of $\pm 500$ km s$^{-1}$. The dashed-dotted line represents $f_c(\text{HI}) (\approx 0.49)$ at the same threshold column density as obtained for the IGM, similar to the empty squares in Figure 4. The first two bins cover the ranges 0–100 pkpc and 100–200 pkpc, respectively, and the third bin covers 200–320 pkpc. Although there is no obvious trend between $f_c$ and $\rho$, it is evident that $f_c(\text{HI})$ is consistently higher ($\approx 1$) for the pair/group LAEs than for the isolated ones ($\approx 0.8$). We obtain similar trends if we use a smaller velocity window of $\pm 250$ km s$^{-1}$ for associating LAEs and H I absorption.

The right panel of Figure 5 shows the differential $f_c(\text{HI})$ in bins of LOS velocity for $\rho \le 320$ pkpc (median $\rho \approx 160$ pkpc). As in the left panel, a threshold $\log_{10}(N/cm^{-2}) = 15.0$ is used here. The dashed-dotted line represents $f_c(\text{HI})$ at the same threshold column density in the IGM, but for a velocity window of 100 km s$^{-1}$, consistent with the dv-bin size used in the plot. The pairs/groups are again seen to have higher $f_c(\text{HI})$ than the isolated LAEs, although at large LOS velocities they become consistent within the $1\sigma$ allowed uncertainties.

Finally, we confirm that using a smaller velocity window of $\pm 250$ km s$^{-1}$ for defining the galaxy groups yielded a consistent trend as our fiducial value of 500 km s$^{-1}$. See the tables in Tables 2, 3, and 4 in the Appendix for the $f_c(\text{HI})$ values calculated in different impact parameter bins when varying the threshold H I column densities and with different LOS velocities for associating H I absorbers with the LAEs.

### 3.4. Influence of EW$_0$ on $f_c(HI)$

The rest-frame equivalent width of Ly$\alpha$ emission ($EW_0$) is defined as

$$ EW_0 = \frac{F_{Ly\alpha}^{\text{line}}}{f_{Ly\alpha}^{\text{cont}} \times (1 + z)} \quad (3) $$

Here, $F_{Ly\alpha}^{\text{line}}$ is the Ly$\alpha$ emission line flux, and $f_{Ly\alpha}^{\text{cont}}$ is the UV continuum flux density at the observed Ly$\alpha$ wavelength. The continuum flux density is estimated from the extrapolation of measured continuum at rest-frame 1500 Å assuming a UV continuum slope ($\beta_{UV}$) of -2.0 (R. J. Bouwens et al. 2014). In our sample, for some of the LAEs, the UV continuum is detected with $>5\sigma$ significance. In the case of a nondetection, we place a lower limit on $EW_0$ (see S. Muzahid et al. 2020 for details).

Figure 6 shows $f_c(\text{HI})$ within $\pm 500$ km s$^{-1}$ as a function of $EW_0$ for two different threshold $N_{\text{HI}}$ values, as indicated in the legend. To divide the LAEs into two $EW_0$ bins, we first exclude 15 LAEs for which the continuum is blended with foreground sources. Next, we separate the remaining 81 nonblended LAEs based on their median $EW_0$ value (51.8 Å). Subsequently, we exclude the lower limits on $EW_0$ from the lower $EW_0$ bin, as their values may not fall within that specific bin. However, this is not the case for LAEs with lower limits on $EW_0$ falling in the upper $EW_0$ bin, thus they are considered. Although the $f_c(\text{HI})$ values for the two $EW_0$ bins are consistent within the error bars for the lower threshold $\log_{10}(N/cm^{-2}) = 15.0$, a significant decrease in the covering factor with $EW_0$ is evident for the threshold $\log_{10}(N/cm^{-2}) = 16.0$. We find a similar trend using a smaller velocity window of $\pm 250$ km s$^{-1}$ for associating absorbers with the LAEs. See Table 5 in the Appendix for the tabulated values of $f_c(\text{HI})$. Additionally, we noted that the
<figure><img src="image_8.png" /></figure>

**Figure 5.** Left: the differential H I covering fraction for a threshold $log_{10}(N/cm^{-2}) = 15.0$ and velocity window $\pm 500$ km s$^{-1}$ for the isolated (red squares) and pairs/group (blue circles) LAEs against impact parameter. The error bars along the y-axis indicate the 68% Wilson-score confidence intervals. Along the x-axis, the data points and error bars, respectively, represent the median values and the 16th–84th percentile range of the $\rho$ distribution within each bin. The black dashed-dotted line shows $f_c(H I)$ for the IGM at the same threshold column density. Although no visible trend is observed with respect to $\rho$, the pairs/group LAEs consistently exhibit a higher $f_c(H I)$ across all three bins. Right: $f_c(H I)$ with a differential LOS velocity window used to associate the H I absorbers with the LAEs. The x-axis denotes the midpoint of each velocity bin, with error bars representing the span of that bin. As in the left panel, the threshold column density used here is $log_{10}(N/cm^{-2}) = 15.0$. The dashed-dotted line represents $f_c(H I)$ at that threshold in the IGM for a velocity window of 100 km s$^{-1}$. Other details are as in the left panel.

median value of the total $N_{H I}$ within $\pm 500$ km s$^{-1}$ of the LAEs contributing to the lower $EW_0$ bin is $log_{10}(N/cm^{-2}) = 16$, whereas for those contributing to the higher $EW_0$ bin it is $log_{10}(N/cm^{-2}) = 15.5$.

We verified that almost 50% of the LAEs in each $EW_0$ bin reside in a pair/group environment. Therefore, it is unlikely that the observed trend is driven by environmental differences.

# 4. Discussion

## 4.1. The Host Galaxies of Strong H I Absorbers

The strong H I absorbers are categorized as LLSs, sub-damped Ly$\alpha$ absorbers (sub-DLAs), and DLAs, depending on their column densities ($N_{H I}$) being $>10^{17.2}$, $>10^{19.0}$, and $>10^{20.3}$ cm$^{-2}$, respectively. Among the 11 optically thick (i.e., $N_{H I} > 10^{17.2}$ cm$^{-2}$) absorbers in our absorber catalog, only one is classified as a DLA and four are sub-DLAs. In Figure 7, we present a comparison of our estimated $N_{H I}$ values with those from the literature for the same systems, when available (M. Neeleman et al. 2013; J. X. Prochaska et al. 2015; E. K. Lofthouse et al. 2023). The plot indicates a good agreement between our measurements and the reported values in the majority of cases. Although one absorber deviates significantly from the 1:1 line, it is actually consistent due to it being a lower limit.

High-column-density H I absorbers, particularly DLAs, are considered effective probes for studying gas within a few tens of kiloparsecs of galaxies (e.g., R. Cooke et al. 2010; A. Rahmati & J. Schaye 2014; K. H. R. Rubin et al. 2015). Pre-MUSE surveys utilizing instruments with a small FOV, such as VLT-SINFONI (FOV of $10'' \times 10''$; C. Péroux et al. 2012) or long-slit observations using X-Shooter at multiple position angles (J. P. U. Fynbo et al. 2010; H. Rahmani et al. 2016; J. K. Krogager et al. 2017), have successfully detected DLA host galaxies, often within a few tens of kiloparsecs. Several of these earlier studies

have targeted high-metallicity DLAs, as their host galaxies are likely to adhere to a mass-metallicity relation, followed by the observed metallicity-velocity relation for DLAs (C. Ledoux et al. 2006; M. Neeleman et al. 2013). Assuming that luminosity scales with stellar mass, it is expected that high-metallicity DLAs would have brighter galaxy counterparts, making them easier to detect.

Of the 96 LAEs in our “absorption-blind” sample, only one is associated with a DLA, with an impact parameter of 230 pkpc. The LAEs associated with the four sub-DLAs have impact parameters ranging from 160 to 240 pkpc. These findings are in line with the MUSE-based, absorption-centric study by R. Mackenzie et al. (2019). They identified 14 LAEs associated with six metal-poor ($log_{10}(Z/Z_\odot) \lesssim -1$) DLAs. All but two LAEs were found to be located at $\rho > 50$ pkpc. Explaining such high column densities of neutral gas at these large distances remains challenging, particularly in light of studies estimating the radius of DLAs to be $\lesssim 10$ kpc (e.g., R. Cooke et al. 2010; J. Cooke & J. M. O’Meara 2015; K. H. R. Rubin et al. 2015). The large DLA-LAE separation suggests that they may not be directly linked to the absorbing gas.

All the LAEs associated with LLSs in MUSEQuBES have impact parameters of $>50$ kpc with no significant trend between $N_{H I}$ and $\rho$ (see Figure 3). Recently, as part of the MUSE Analysis of Gas around Galaxies (MAGG) survey, E. K. Lofthouse et al. (2023) also reported no observed trend between $N_{H I}$ of their LLS sample and the impact parameter of LAEs detected within $\pm 1000$ km s$^{-1}$ out to $\approx 300$ pkpc. They presented a detailed discussion on the connection between optically thick gas and LAEs in their sample (see their Section 7.2). They suggested that LLSs originate primarily either from the outer CGM (defined as $3-4 R_{vir}$) or in the IGM but in proximity to galaxies. When multiple LAEs are detected near an LLS, they observe a preferential alignment of the LAEs, which is suggestive of the fact that LAEs and optically thick
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Covering fraction as a function of the rest-frame equivalent width (EW₀) of the Lyα emission. Different colors represent different threshold column densities, as indicated by the legend. The error bars along the y-axis indicate the 68% Wilson-score confidence interval. The numbers adjacent to the data points correspond to the respective values of *n*<sub>Hit</sub> and *n*<sub>Total</sub> (see Equation (1)). Along the x-axis, the points represent the median values and the error bars are the 16th–84th percentile ranges of EW₀ in each bin. While *f*<sub>c(H I)</sub> for the two EW₀ bins are consistent within the error bars for the lower threshold log₁₀(N/cm⁻²) = 15.0, a significant anticorrelation between *f*<sub>c(H I)</sub> and EW₀ is present for the threshold log₁₀(N/cm⁻²) = 16.0.</figcaption></figure>

absorbers may lie within filamentary structures. Such a
filamentary structure is also evident in the MUSEQuBES
survey for the system with the highest LAE overdensity
detected toward Q1317-0507 (E. Banerjee et al. 2025, in
preparation).

Exploiting ≈140 hr of MUSE observations, R. Bacon et al. (2021) found that extended diffuse Lyα emission around LAEs at z ≈ 3.1–4.5 originates from gas beyond their CGM, likely driven by a large population of ultra-low-luminosity LAEs with L<sub>Lyα</sub> < 10<sup>40</sup> erg s<sup>-1</sup>. In comparison, the median L<sub>Lyα</sub> for our sample, i.e., ≈10<sup>42</sup> erg s<sup>-1</sup>, is two orders of magnitude higher. Those faint LAEs will remain undetected in MUSEQuBES with on-source exposure times of ≲10 hr per field, but could potentially give rise to strong H I absorption in background quasar spectra. This association of strong H I absorbers with ultra-low-luminosity galaxies agrees with predictions from cosmological simulations (A. Rahmati & J. Schaye 2014) and provides a plausible explanation for the presence of high-N<sub>H I</sub> absorbers at larger impact parameters in MUSE surveys.

Finally, H. Kusakabe et al. (2020) showed that at z ≈ 3–6, only up to 30% of star-forming galaxies are LAEs. This indicates that surveys focusing solely on LAEs may not capture the full galaxy population. Consequently, the nondetection of the true hosts of the absorbers can result in a scatter in the distribution of N<sub>H I</sub> with impact parameters, as shown in Figure 3, leading to the observed lack of trend.

4.2. H I Covering Fractions

The characteristics of H I absorption are primarily influenced
by three main factors: (a) the Hubble expansion; (b) changes in
the ionizing UV background radiation field resulting from star-
forming galaxies and active galactic nuclei; and (c) the growth
of large-scale structure. R. Davé et al. (1999) demonstrated
that, as a consequence of these factors, the overdensity

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison of the H I column densities measured in this work with those reported in the literature for the same systems. The dashed line indicates the 1:1 relationship. For the majority of the cases, they are in good agreement with each other.</figcaption></figure>

associated with a fixed N<sub>H I</sub> decreases with increasing redshift
(see also T. Theuns et al. 1998; J. Schaye 2001). Consequently,
absorbers with the same N<sub>H I</sub> at different redshifts correspond to
different types of cosmic structures. Using Equation (10) from
J. Schaye (2001), it can be inferred that log<sub>10</sub>(N/cm<sup>-2</sup>) = 14.0
corresponds to an overdensity of δ ≈ 15 at z ≈ 0, whereas at
z ≈ 3 it corresponds to an overdensity of δ ≈ 1.4
(I + δ = ρ/ρ̄, where ρ̄ is the cosmic mean matter density).
These findings result in a variation in dn/dz of H I absorbers at
different cosmic epochs (e.g., T.-S. Kim et al. 2021). This
evolution of the H I number density is expected to result in a
higher H I covering fraction at higher redshift for a fixed
comoving impact parameter and an even stronger effect for a
fixed proper impact parameter provided collisional ionization is
unimportant.

In Figure 8 we compare $f_c(\text{H I})$ at threshold $\log_{10}(N/\text{cm}^{-2}) = 14.0$ that we measure for MUSEQuBES with a few surveys from the literature: the KBSS-LBGs at $z \approx 2.3$ from G. C. Rudie et al. (2012, see their Table 6), the $L_B \gtrsim 0.25 L_B^*$ galaxies from H.-W. Chen et al. (2001b, $f_c(\text{H I}) = 34/47$; $0.1 < z < 0.9$), and $L > 0.1 L^*$ galaxies from J. X. Prochaska et al. (2011, $f_c(\text{H I}) = 18/25$; $0.005 < z < 0.4$) and B. P. Wakker & B. D. Savage (2009, $f_c(\text{H I}) = 6/14$; local Universe). Along the x-axis, we indicate the LOS velocity windows and impact parameter ranges used for the respective $f_c(\text{H I})$ calculations. The MUSEQuBES LAEs have the highest $f_c(\text{H I})$ in all cases. Notably, it becomes comparable to that of the $z \approx 2.3$ LBGs for $\rho < 200$ pkpc and $\Delta v = 400$ km s$^{-1}$. If instead of the proper impact parameter we calculate $f_c(\text{H I})$ of the MUSEQuBES sample for the same comoving impact parameter as these low-$z$ surveys, it reaches 100% for all three different cuts explored here.

Note that the overdensity (≈15) traced by an absorber of N(H I) ≈ 10<sup>14.0</sup> cm<sup>-2</sup> at z ≈ 0 would be traced by an absorber of N(H I) ≈ 10<sup>15.6</sup> cm<sup>-2</sup> at z ≈ 3.3. *f*<sub>c(H I)</sub> we obtained for the threshold N(H I) of 10<sup>15.6</sup> cm<sup>-2</sup> is 0.38 within 200 pkpc and 0.62 within 320 kpc from our survey. It is intriguing to note that these values are consistent with low-z measurements of B. P. Wakker & B. D. Savage (2009) for <200 pkpc and H.-W. Chen et al. (2001a) and J. X. Prochaska et al. (2011) for
<figure><img src="image_11.png" /><figcaption>**Figure 8.** A comparison of H I covering fraction for the MUSEQuBES sample (represented by red circles) and several samples from the literature (indicated by the legends) at different redshifts, calculated at the same threshold $log_{10}(N/cm^{-2}) = 14.0$. The literature samples consist of the high-z KBSS-LBGs from G. C. Rudie et al. (2012), low-z $L_B \gtrsim 0.25 L_B^*$ galaxies from H.-W. Chen et al. (2001b), and low-z $L > 0.1 L_B^*$ galaxies from J. X. Prochaska et al. (2011) and B. P. Wakker & B. D. Savage (2009). Along the x-axis, we indicate the respective LOS velocity windows and physical impact parameter ranges used for the $f_c(HI)$ measurements. The error bars shown for the KBSS-LBG data are directly taken from G. C. Rudie et al. (2012). For the remaining cases, the error bars represent the 68% Wilson-score confidence intervals. For the MUSEQuBES sample, the error bars are smaller than the circle size. MUSEQuBES consistently exhibits the highest $f_c(HI)$.</figcaption></figure>

≤330 pkpc for threshold $log_{10}(N/cm^{-2}) = 14.0$, as can be seen from Figure 8.

S. D. Johnson et al. (2015) found $f_c(HI) \approx 90\%$ at threshold $log_{10}(N/cm^{-2}) = 14$ within the virial radii of $L > 0.1 L_B^*$ galaxies at $z < 0.4$, which reduces to only ≈16% at larger distances. The COS-Halos survey (J. Tumlinson et al. 2013) reported $f_c(HI) \approx 91\%$ for $log_{10}(N/cm^{-2}) > 14$ within 150 pkpc and $|\Delta v| \approx 200 \text{ km s}^{-1}$ for ≈L galaxies at $z \approx 0.2$. In our sample, for the same cuts on LOS velocity and impact parameter, $f_c(HI) \approx 92\%$ (38/41). Thus, it appears that the outskirts (≈3$R_{vir}$) of the LAEs in our high-z sample are as H I-rich as the inner CGM ($< 0.5 R_{vir}$) of ≈L galaxies at low z.

The MAGG survey measured $f_c(HI) \approx 0.2$ within ±500 km s⁻¹ and 300 pkpc from the $z \approx 3$ LAEs in their sample for a threshold column density $log_{10}(N/cm^{-2}) = 17.2$ (E. K. Lofthouse et al. 2023). We obtain $f_c(HI) = 0.16_{-0.04}^{+0.04}$ for MUSEQuBES for the same threshold. Additionally, E. K. Lofthouse et al. (2023) reported $f_c(HI) = 0.29 \pm 0.03$ for group galaxies and $f_c(HI) = 0.08 \pm 0.02$ for isolated LAEs for the same threshold $N(HI)$. In contrast, MUSEQuBES yields $f_c(HI) = 0.10_{-0.06}^{+0.04}$ (3/29) for group galaxies and $f_c(HI) = 0.19_{-0.05}^{+0.04}$ (9/46) for isolated LAEs. The measurements of $f_c(HI)$ in both surveys use similar ranges of $\Delta v$ (±500 km s⁻¹) and $\rho$ (≈300 pkpc). Recall that only 11 absorbers associated with the LAEs in our sample have $N(HI) > 10^{17.2} cm^{-2}$. We could not confirm the enhancement of $f_c(HI)$ for group galaxies as reported by E. K. Lofthouse et al. (2023) owing to the scarcity of LLSs in our study (their survey had more than 50 absorbers with $log_{10}(N/cm^{-2}) \ge 17.2$). Nevertheless, we do see a strong environmental effect on $f_c(HI)$ at a smaller threshold $log_{10}(N/cm^{-2}) = 15.0$ out to ≈250 pkpc (see Figure 5).

<figure><img src="image_12.png" /><figcaption>**Figure 9.** H I covering fraction (within 320 pkpc and ±500 km s⁻¹ LOS velocity of the LAEs) calculated at two different threshold column densities as indicated by the legend for two bins of SFR. The error bars along the y-axis indicate the 68% Wilson-score confidence intervals. The numbers adjacent to the blocks correspond to the respective values of $n_{\text{Hit}}$ and $n_{\text{Total}}$. Along the x-axis, the points represent the median values and the error bars are the 16th–84th percentile ranges of SFR in each bin. No significant trend is found between $f_c$ and SFR of the LAEs.</figcaption></figure>

S. Muzahid et al. (2021) reported a correlation between the rest-frame equivalent width of stacked circumgalactic Lyα absorption and the SFR (determined via the UV continuum emission) of the MUSEQuBES LAEs. They split the LAEs into three subsamples: two subsamples with high and low SFRs, and an additional subsample encompassing LAEs for which the UV continuum was not detected with >5σ significance, i.e., the nondetections. Their study revealed that the higher-SFR subsample has a mean equivalent width of H I absorption ≈1.6 times larger than that of the lower-SFR subsample. However, from our individual Voigt profile analysis, we did not find any significant correlation between the total $N_{\text{H I}}$ and SFR for these LAEs. The median total $N_{\text{H I}}$ for the higher-SFR subsample ($10^{15.7} cm^{-2}$) is similar to that of the lower-SFR subsample ($10^{15.5} cm^{-2}$).

This lack of a trend is also evident in Figure 9, where we plot $f_c(HI)$ against $log_{10}(SFR)$. The data are binned based on the median value of the SFR of the LAEs. We excluded the LAEs with SFR nondetections (i.e., upper limits) from the higher SFR bin but considered the LAEs when their upper limits fall in the lower SFR bin. $f_c(HI)$ for LAEs in both the SFR bins, calculated for two different threshold column densities, are consistent within the error bars, indicating a lack of any significant difference in the H I distribution in the CGM for the subsamples with higher and lower SFR. Note that the median SFRs of these two bins only differ by a factor of ~3 from each other.

Additionally, we calculated the column density-weighted LOS velocity of H I absorbers from the respective host LAEs ($\Delta v_{\text{wgt}}$), obtaining a median value of 4.8 km s⁻¹ and $\sigma = 175.4$ km s⁻¹. However, we do not find any correlation between the impact parameter and $\Delta v_{\text{wgt}}$ as indicated by the Spearman correlation test results: $r_{\text{corr}} = 0.01$ and $p$-value = 0.91. We performed this correlation test by creating smaller subsamples of galaxies based on their environments
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The mean excess Lyα transmission (⟨δ_T⟩) as a function of 3D distance. The ⟨δ_T⟩ values from this work (in orange) are compared with those at z ≈ 4.3 from J. Matthee et al. (2024; in blue). The error bars for our sample represent the 1σ bootstrap error, calculated by resampling the LAE sample 1000 times. The black dotted-dashed line indicates ⟨δ_T⟩ = 0.</figcaption></figure>

(i.e., pairs/group or isolated LAEs) and did not find any significant correlations in any of the cases.

### 4.3. The Relation between Neutral Gas in the ISM and CGM

In Figure 6, we observed an anticorrelation between $f_c(H1)$ and EW$_0$. This trend is particularly evident for the threshold $log_{10}(N/cm^{-2}) = 16.0$, where the difference in $f_c(H1)$ between the two EW$_0$ bins is significant at the ≈2σ level. Since ≈50% of LAEs in both the EW$_0$ bins are associated with the pair/group subsample, it is unlikely that environmental effects are driving the trend.

Simple shell models of Lyα radiative transfer predict a strong anticorrelation between $N(H1)$ and EW$_0$ (A. Verhamme et al. 2015, 2017). This is because as $N(H1)$ increases, Lyα photons experience more scattering, unlike continuum photons. Hence, Lyα photons cover a longer distance in order to shift their frequencies off-resonance to escape the interstellar medium (ISM) of the galaxies, which increases their probability of being destroyed by dust grains. Thus, the observed anticorrelation between $f_c(H1)$ and EW$_0$ could imply a positive trend between the covering fraction of neutral gas in the CGM and $N(H1)$ in the ISM of the LAEs. In other words, galaxies in whose environments cool gas is abundant are also inherently gas-rich themselves.

The shell model (A. Verhamme et al. 2006) also predicts that galaxies with lower $N_{H1}$ exhibit narrower line profiles, and vice versa. However, we did not observe a significant difference in the distribution of FWHM of Lyα emission lines between the high- and low-EW$_0$ samples ($p_{KS} = 0.2$). Both samples have a median FWHM of ≈250 km s⁻¹.

An alternative explanation for the observed trend could be enhanced photoionization from the local radiation fields of LAEs with relatively higher EW$_0$, which eventually leads to a lower H1 covering fraction (e.g., J. Matthee et al. 2022).

### 4.4. Correlation between Lyα Transmission and LAEs

In this section, we discuss the influence of LAEs on the Lyα transmission as measured in the quasar spectra. We define the excess Lyα transmission ($\delta_T$) with respect to the average

transmission at a given redshift (i.e., ⟨$T_z$⟩) as

$$ \delta_T = \frac{T}{\langle T_z \rangle} - 1. $$

Here, T represents the normalized spectral flux.

In Figure 10, we present the mean excess Lyα transmission (⟨$\delta_T$⟩: shown in orange) as a function of the 3D distance, $\Delta r_{3D}$, where $r_{3D} = \sqrt{n^2 + r_{\perp}^2}$. Here, $r_{\parallel}$ and $r_{\perp}$ represent the LOS and transverse distances with respect to the LAEs, respectively. To calculate this, we summed the total flux from the normalized quasar spectra within the LOS velocity window $\Delta v$ around the systemic redshift of each LAE and averaged it over the total number of pixels contributing to the respective $\Delta v$ bin considering all LAEs. We converted $\Delta v$ to $r_{\parallel}$ (in comoving Mpc or cMpc) using ($\Delta v/H(z)) \times (1+z$), assuming a redshift of $z=3.3$, the median for these LAEs. Additionally, we assumed $r_{\perp} \approx 0.7$ cMpc, corresponding to the median impact parameter of 160 pkpc for our sample. To model the redshift evolution of ⟨$T_z$⟩ in our data, we followed a method similar to J. Matthee et al. (2024, see also G. D. Becker et al. 2015). We used the following polynomial function that best fits our data:

$$ \langle T_z \rangle = -1.30 + 1.34z - 0.23z^2. $$

We have also shown ⟨$\delta_T$⟩ near the LAEs at a slightly higher redshift of ≈4.3, as found by J. Matthee et al. (2024; in blue). The data points are mostly within each other's error bars.

## 5. Summary and Conclusions

This work is a part of the MUSEQuBES survey at high z (S. Muzahid et al. 2020, 2021), which focuses on the CGM of high-z LAEs. The MUSEQuBES sample consists of 96 LAEs (median Lyα luminosity ≈10⁴² erg s⁻¹ and a median SFR of ≈1.3 M$_\odot$ yr⁻¹, without accounting for dust extinction) within the redshift range 2.9–3.8 that were detected in eight MUSE fields centered on eight bright background quasars. The high-resolution optical spectra of these eight quasars obtained with VLT/UVES and/or Keck/HIRES enabled us to explore the neutral hydrogen (H1) absorption arising near these 96 galaxies.

We estimated the systemic redshifts of the LAEs using the empirical relation given by S. Muzahid et al. (2020). Subsequently, we conducted a Voigt profile decomposition of the Lyman-series lines detected within a velocity range of ±500 km s⁻¹ around the LAEs. The column densities of the 800 H1 absorption components observed in these spectra span a range from $log_{10}(N/cm^{-2}) = 11.7$ to 20.3, with a median value of $log_{10}(N/cm^{-2}) = 13.7$. In cases where we were unable to precisely constrain the column density due to blending or the absence of higher-order unsaturated Lyman-series lines, we reported the lower limit on the column density. Our main findings are:

1. We do not find any significant correlation between the total $N_{H1}$ and impact parameter within 16–315 pkpc (median 165 pkpc; Figure 3). However, the total $N_{H1}$ within ±500 km s⁻¹ of the LAEs is generally higher than the typical $N_{H1}$ expected from the IGM at similar redshifts. This trend also holds for absorbers within ±250 km s⁻¹ of the LAE redshifts, suggesting that H1 absorption is elevated near the LAEs compared to the IGM.
2. All the LAEs associated with HI absorbers of $log_{10}(N/cm^{-2}) > 17.2$ are detected at impact parameters >50 pkpc (Figure 3). This is likely due to the nondetection of the actual absorber hosts owing to their faintness, as discussed in Section 4.1.

3. The covering fraction, $f_c(HI)$, is consistently higher near the LAEs than in random regions for thresholds $log_{10}(N/cm^{-2}) = 14-16$, with $\rho \le 320$ pkpc. The relative enhancement in $f_c(HI)$ with respect to random regions increases with the threshold column density (Figure 4).

4. The pairs/group galaxies exhibit a higher $f_c(HI)$ than isolated LAEs for a threshold $log_{10}(N/cm^{-2}) = 15$ along transverse directions (Figure 5, left). The pairs/group LAEs display a 100% HI covering fraction, extending up to ≈250 pkpc.

5. $f_c(HI)$ decreases with LOS velocity between the LAEs and absorbers (Figure 5, right). The trend of pairs/groups having higher $f_c(HI)$ than isolated LAEs is also visible along the LOS, although only up to ≈±300 km s⁻¹, beyond which $f_c(HI)$ values are comparable within $1\sigma$.

6. For a threshold $log_{10}(N/cm^{-2}) = 16$, $f_c(HI)$ anticorrelates with the rest-frame equivalent width (EW₀) of Lyα emission (Figure 6). According to the Lyα shell model, a lower EW₀ results from a large reservoir of neutral gas present within the galaxies. This could imply that galaxies residing in cool-gas-rich environments are also inherently gas-rich. The other reason might be that high-EW₀ LAEs are efficient in ionizing their surrounding medium, which decreases $f_c(HI)$ (see Section 4.3).

7. Finally, we compared the excess Lyα transmission, $\langle \delta_T \rangle$, from our observation and that of the LAEs at a slightly

higher redshift of 4.3 and found they are mostly consistent with each other.

This study presents the connection between $z \approx 3$ LAEs and cool neutral gas in their surroundings, for the first time, in an “absorption-blind” manner using detailed Voigt profile decomposition of HI absorbers. In the future, we will present ionization models to derive the density and metallicity of these absorbers and investigate their relation with the LAEs.

## Acknowledgments

We would like to thank the anonymous referee for useful comments. We thank Marijke Segers, Lorrie Straka, and Monica Turner for their early contributions to the MUSE-QuBES project. We thank Raghunathan Srianand for useful suggestions. E.B. thanks Labanya Kumar Guha and Yucheng Guo for helpful discussions. S.C. gratefully acknowledges the fund support from the European Research Council (ERC).

*Software:* NumPy (C. R. Harris et al. 2020), SciPy (P. Virtanen et al. 2020), Matplotlib (J. D. Hunter 2007), and AstroPy (Astropy Collaboration et al. 2013, 2018).

## Appendix

In this section we have provided Tables 2, 3, 4, and 5, containing the covering fractions of HI in different column density thresholds, impact parameters and LOS velocities as mentioned in the titles of the tables.

<table><caption>**Table 2**
$f_c(HI)$ at Threshold $log(N(HI)/cm^{-2}) = 15$ for Different Impact Parameters ($\rho$)</caption><thead><tr><th rowspan="2">Impact Parameter (pkpc)</th><th colspan="2">Δv = 500 km s<sup>-1</sup></th><th colspan="2">Δv = 250 km s<sup>-1</sup></th></tr><tr><th>Pairs/Group</th><th>Isolated</th><th>Pairs/Group</th><th>Isolated</th></tr></thead><tbody><tr><td>0–100</td><td>1.0<sup>+0.18</sup> (4/4)</td><td>0.81<sup>+0.08</sup><sub>-0.13</sub> (9/11)</td><td>1.0<sup>+0.18</sup> (4/4)</td><td>0.67<sup>+0.11</sup><sub>-0.13</sub> (8/12)</td></tr><tr><td>100–200</td><td>1.0<sup>+0.05</sup> (17/17)</td><td>0.83<sup>+0.06</sup><sub>-0.08</sub> (24/29)</td><td>0.82<sup>+0.07</sup><sub>-0.1</sub> (14/17)</td><td>0.73<sup>+0.07</sup><sub>-0.08</sub> (22/30)</td></tr><tr><td>200–320</td><td>1.0<sup>+0.07</sup> (12/12)</td><td>0.79<sup>+0.07</sup><sub>-0.10</sub> (15/19)</td><td>0.83<sup>+0.08</sup><sub>-0.12</sub> (10/12)</td><td>0.44<sup>+0.11</sup><sub>-0.10</sub> (8/18)</td></tr><tr><td>≤320</td><td>1.0<sup>+0.02</sup> (33/33)</td><td>-0.81<sup>+0.04</sup><sub>-0.05</sub> (48/59)</td><td>0.84<sup>+0.05</sup><sub>-0.07</sub> (28/33)</td><td>0.63<sup>+0.06</sup><sub>-0.06</sub> (38/60)</td></tr></tbody></table>

Note. We have included $n_{\text{Hit}}(N \ge N_{\text{Th}})/n_{\text{Total}}$ values in parentheses beside each entry.

<table><caption>**Table 3**
$f_c(HI)$ at Threshold $log(N(HI)/cm^{-2}) = 16$ for Different Impact Parameters ($\rho$)</caption><thead><tr><th rowspan="2">Impact Parameter (pkpc)</th><th colspan="2">Δv = 500 km s<sup>-1</sup></th><th colspan="2">Δv = 250 km s<sup>-1</sup></th></tr><tr><th>Pairs/Group</th><th>Isolated</th><th>Pairs/Group</th><th>Isolated</th></tr></thead><tbody><tr><td>0–100</td><td>0.5<sup>+0.21</sup><sub>-0.21</sub> (2/4)</td><td>0.33<sup>+0.16</sup><sub>-0.12</sub> (3/9)</td><td>0.5<sup>+0.21</sup><sub>-0.21</sub> (2/4)</td><td>0.3<sup>+0.15</sup><sub>-0.11</sub> (3/10)</td></tr><tr><td>100–200</td><td>0.36<sup>+0.13</sup><sub>-0.10</sub> (5/14)</td><td>0.39<sup>+0.09</sup><sub>-0.08</sub> (11/28)</td><td>0.28<sup>+0.12</sup><sub>-0.1</sub> (4/14)</td><td>0.37<sup>+0.10</sup><sub>-0.08</sub> (11/29)</td></tr><tr><td>200–320</td><td>0.45<sup>+0.14</sup><sub>-0.13</sub> (5/11)</td><td>0.35<sup>+0.11</sup><sub>-0.10</sub> (6/17)</td><td>0.36<sup>+0.14</sup><sub>-0.12</sub> (4/11)</td><td>0.22<sup>+0.10</sup><sub>-0.08</sub> (4/18)</td></tr><tr><td>≤320</td><td>0.41<sup>+0.09</sup><sub>-0.08</sub> (12/29)</td><td>0.37<sup>+0.06</sup><sub>-0.06</sub> (20/54)</td><td>0.34<sup>+0.09</sup><sub>-0.07</sub> (10/29)</td><td>0.31<sup>+0.06</sup><sub>-0.05</sub> (18/57)</td></tr></tbody></table>

Note. We have included $n_{\text{Hit}}(N \ge N_{\text{Th}})/n_{\text{Total}}$ values in parentheses beside each entry.
<table><caption>**Table 4**
$f_c(H I)$ for Different LOS Velocity Windows ($\rho \le 320$ pkpc)</caption><thead><tr><th rowspan="2">LOS Velocity<br>(km s<sup>-1</sup>)</th><th colspan="2">Threshold log(N(H I)/cm<sup>-2</sup>) = 15</th><th colspan="2">Threshold log(N(H I)/cm<sup>-2</sup>) = 16</th></tr><tr><th>Pairs/Group</th><th>Isolated</th><th>Pairs/Group</th><th>Isolated</th></tr></thead><tbody><tr><td>0–100</td><td>0.61<sup>+0.08</sup><sub>-0.08</sub> (20/33)</td><td>0.39<sup>+0.06</sup><sub>-0.06</sub> (23/59)</td><td>0.26<sup>+0.08</sup><sub>-0.07</sub> (8/30)</td><td>0.12<sup>+0.05</sup><sub>-0.03</sub> (7/57)</td></tr><tr><td>100–200</td><td>0.45<sup>+0.08</sup><sub>-0.08</sub> (15/33)</td><td>0.28<sup>+0.06</sup><sub>-0.05</sub> (17/60)</td><td>0.06<sup>+0.05</sup><sub>-0.03</sub> (2/32)</td><td>0.17<sup>+0.05</sup><sub>-0.04</sub> (10/59)</td></tr><tr><td>200–300</td><td>0.24<sup>+0.08</sup><sub>-0.06</sub> (8/33)</td><td>0.12<sup>+0.05</sup><sub>-0.03</sub> (7/57)</td><td>0.0<sup>+0.03</sup><sub>-0.00</sub> (0/31)</td><td>0.07<sup>+0.04</sup><sub>-0.02</sub> (4/57)</td></tr><tr><td>300–400</td><td>0.24<sup>+0.08</sup><sub>-0.06</sub> (8/33)</td><td>0.16<sup>+0.05</sup><sub>-0.04</sub> (9/56)</td><td>0.15<sup>+0.07</sup><sub>-0.05</sub> (5/33)</td><td>0.02<sup>+0.03</sup><sub>-0.01</sub> (1/55)</td></tr><tr><td>400–500</td><td>0.18<sup>+0.07</sup><sub>-0.06</sub> (6/32)</td><td>0.14<sup>+0.05</sup><sub>-0.04</sub> (9/61)</td><td>0.06<sup>+0.05</sup><sub>-0.03</sub> (2/31)</td><td>0.03<sup>+0.03</sup><sub>-0.02</sub> (2/60)</td></tr></tbody></table>

Note. We have included $n_{\text{Hit}}(N \ge N_{\text{Th}})/n_{\text{Total}}$ values in parentheses beside each entry.

<table><caption>**Table 5**
$f_c(H I)$ for Different EW$_0$</caption><thead><tr><th rowspan="2">EW<sub>0</sub><br>(Å)</th><th colspan="2">Threshold log N(H I)/cm<sup>-2</sup> = 15</th><th colspan="2">Threshold log N(H I)/cm<sup>-2</sup> = 16</th></tr><tr><th>Δν = 500 km s<sup>-1</sup></th><th>Δν = 250 km s<sup>-1</sup></th><th>Δν = 500 km s<sup>-1</sup></th><th>Δν = 250 km s<sup>-1</sup></th></tr></thead><tbody><tr><td>&lt;52</td><td>1.0<sup>+0.04</sup><sub>-0.04</sub> (21/21)</td><td>0.86<sup>+0.05</sup><sub>-0.08</sub> (19/22)</td><td>0.52<sup>+0.10</sup><sub>-0.10</sub> (10/19)</td><td>0.4<sup>+0.10</sup><sub>-0.10</sub> (8/20)</td></tr><tr><td>&gt;52</td><td>0.87<sup>+0.04</sup><sub>-0.06</sub> (36/41)</td><td>0.70<sup>+0.06</sup><sub>-0.07</sub> (29/41)</td><td>0.21<sup>+0.07</sup><sub>-0.05</sub> (8/38)</td><td>0.21<sup>+0.07</sup><sub>-0.05</sub> (8/38)</td></tr></tbody></table>

Note. 52 Å is the median EW$_0$ for our sample. Also, we have included $n_{\text{Hit}}(N \ge N_{\text{Th}})/n_{\text{Total}}$ values in parentheses beside each entry.

## ORCID iDs

Eshita Banerjee https://orcid.org/0009-0002-7382-3078

Sowgat Muzahid https://orcid.org/0000-0003-3938-8762

Joop Schaye https://orcid.org/0000-0002-0668-5560

Nicolas Bouché https://orcid.org/0000-0003-0068-9920

Sebastiano Cantalupo https://orcid.org/0000-0001-5804-1428

Sean D. Johnson https://orcid.org/0000-0001-9487-8583

Jorryt Mathee https://orcid.org/0000-0003-2871-127X

Anne Verhamme https://orcid.org/0000-0002-2201-1865

## References

Adelberger, K. L., Shapley, A. E., Steidel, C. C., et al. 2005, ApJ, 629, 636

Altay, G., Theuns, T., Schaye, J., Crighton, N. H. M., & Dalla Vecchia, C. 2011, ApJL, 737, L37

Arrigoni Battaia, F., Hennawi, J. F., Prochaska, J. X., et al. 2019, MNRAS, 482, 3162

Astropy Collaboration, et al. 2013, A&A, 558, A33

Astropy Collaboration, et al. 2018, AJ, 156, 123

Bacon, R., Accardo, M., Adjali, L., et al. 2010, Proc. SPIE, 647, 773508

Bacon, R., Mary, D., Garel, T., et al. 2021, A&A, 647, A107

Banerjee, E., Muzahid, S., Schaye, J., Johnson, S. D., & Cantalupo, S. 2023, MNRAS, 524, 5148

Becker, G. D., Bolton, J. S., & Lidz, A. 2015, PASA, 32, e045

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Borisova, E., Cantalupo, S., Lilly, S., et al. 2016, ApJ, 831, 39

Bouwens, R. J., Lllingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115

Cantalupo, S., Arrigoni-Battaia, F., Prochaska, J. X., Hennawi, J. F., & Madau, P. 2014, Natur, 506, 63

Carswell, R. F., & Webb, J. K., 2014 VPFIT: Voigt Profile Fitting Program, Astrophysics Source Code Library, ascl: 1408.015

Chabrier, G. 2003, PASP, 115, 763

Chen, H.-W., Lanzetta, K. M., & Webb, J. K. 2001a, ApJ, 556, 158

Chen, H.-W., Lanzetta, K. M., Webb, J. K., & Barcons, X. 2001b, ApJ, 559, 654

Chen, Y., Steidel, C., Hummels, C., et al. 2020, MNRAS, 499, 1721

Cooke, J., & O'Meara, J. M. 2015, ApJL, 812, L27

Cooke, R., Pettini, M., Steidel, C. C., et al. 2010, MNRAS, 409, 679

Danforth, C. W., Stocke, J. T., & Shull, J. M. 2010, ApJ, 710, 613

Davé, R., Hernquist, L., Katz, N., & Weinberg, D. H. 1999, ApJ, 511, 521

Dutta, R., Fossati, M., Fumagalli, M., et al. 2023, MNRAS, 522, 535

Dutta, S., Muzahid, S., Schaye, J., et al. 2024, MNRAS, 528, 3745

Erb, D. K., Steidel, C. C., & Chen, Y. 2018, ApJL, 862, L10

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

Guo, Y., Bacon, R., Bouche, N., et al. 2023, Natur, 624, 53

Guo, Y., Bacon, R., Wisotzki, L., et al. 2024a, A&A, 688, A37

Guo, Y., Bacon, R., Wisotzki, L., et al. 2024b, A&A, 691, A66

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hellsten, U., Hernquist, L., Katz, N., & Weinberg, D. H. 1998, ApJ, 499, 172

Herrero Alonso, Y., Krumpe, M., Wisotzki, L., et al. 2021, A&A, 653, A136

Hunter, J. D. 2007, CSE, 9, 90

Johnson, S. D., Chen, H.-W., & Mulchaey, J. S. 2015, MNRAS, 449, 3263

Johnson, S. D., Schaye, J., Walth, G. L., et al. 2022, ApJL, 940, L40

Johnson, S. D., Liu, Z., Li, J., et al. 2024, ApJ, 966, 218

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

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Kim, T.-S., Wakker, B P, Nasir, F, et al. 2021, MNRAS, 501, 5811

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

Kusakabe, H., Blaizot, J, Thubault, G, et al. 2020, A&A, 638, A12

Leclercq, F., Verhamme, A, Epinat, B., et al. 2022, A&A, 663, A11

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

Liu, Z., Johnson, S. D., Li, J., et al. 2024, MNRAS, 527, 5429

Lofhouse, E. K., Fumagalli, M., Fossati, M., et al. 2023, MNRAS, 518, 305

Mackenzie, R., Fumagalli, M., Theuns, T., et al. 2019, MNRAS, 487, 5070

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Matthee, J., Golling, C., Mackenzie, R., et al. 2024, MNRAS, 529, 2794

Matthee, J., Naidu, R. P., Pezzulli, G., et al. 2022, MNRAS, 512, 5960

Morrissey, P., Matuszewski, M., & Martin, D. C. 2018, ApJ, 864, 93

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Muzahid, S., Schaye, J., Cantalupo, S., et al. 2021, MNRAS, 508, 5612

Muzahid, S., Schaye, J., Marino, R. A., et al. 2020, MNRAS, 496, 1013

Neeleman, M., Kanekar, N., Prochaska, J. X., et al. 2017, Sci, 355, 1285

Neeleman, M., Kanekar, N., Prochaska, J. X., Rafelski, M. A., & Carilli, C. L. 2019, ApJL, 870, L19

Neeleman, M., Wolfe, A. M., Prochaska, J. X., & Rafelski, M. 2013, ApJ, 769, 54

Noterdaeme, P., Petitjean, P., Ledoux, C., & Srianand, R. 2009, A&A, 505, 1087

Ono, Y., Ouchi, M, Shimasaku, K, et al. 2010, MNRAS, 402, 1580

Péroux, C., Bouché, N., Kulkarni, V. P., York, D. G., & Vladilo, G. 2012, MNRAS, 419, 3060

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363
Prochaska, J. X., O'Meara, J. M., Fumagalli, M., Bernstein, R. A., & Burles, S. M. 2015, ApJS, 221, 2

Prochaska, J. X., Weiner, B., Chen, H.-W., Mulchaey, J., & Cooksey, K. 2011, ApJ, 740, 91

Rahmani, H., Peroux, C., Turnshek, D. A., et al. 2016, MNRAS, 463, 980

Rahmati, A., & Schaye, J. 2014, MNRAS, 438, 529

Rakic, O., Schaye, J., Steidel, C. C., & Rudie, G. C. 2011, MNRAS, 414, 3265

Rakic, O., Schaye, J., Steidel, C. C., & Rudie, G. C. 2012, ApJ, 751, 94

Rubin, K. H. R., Hennawi, J. F., Prochaska, J. X., et al. 2015, ApJ, 808, 38

Rudie, G. C., Steidel, C. C., Pettini, M., et al. 2019, ApJ, 885, 61

Rudie, G. C., Steidel, C. C., Shapley, A. E., & Pettini, M. 2013, ApJ, 769, 146

Rudie, G. C., Steidel, C C, Trainor, R F, et al. 2012, ApJ, 750, 67

Schaye, J. 2001, ApJ, 559, 507

Shibuya, T., Ouchi, M, Nakajima, K., et al. 2014, ApJ, 788, 74

Shull, J. M., Giroux, M. L., Penton, S V, et al. 2000, ApJL, 538, L13

Songaila, A. 2001, ApJL, 561, L153

Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2011, ApJ, 736, 160

Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289

Theuns, T., Leonard, A., & Efstathiou, G. 1998, MNRAS, 297, L49

Trainor, R. F., Steidel, C. C., Strom, A. L., & Rudie, G. C. 2015, ApJ, 809, 89

Tumlinson, J., Peeples, M., & Werk, J. 2017, ARA&A, 55, 389

Tumlinson, J., Thom, C., Werk, J., et al. 2013, ApJ, 777, 59

Turner, M. L., Schaye, J., Steidel, C. C., Rudie, G. C., & Strom, A. L. 2014, MNRAS, 445, 794

van de Voort, F., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011a, MNRAS, 415, 2782

van de Voort, F., Schaye, J., Booth, C. M., Haas, M. R., & Dalla Vecchia, C. 2011b, MNRAS, 414, 2458

Verhamme, A., Garel, T., Ventou, E., et al. 2018, MNRAS, 478, L60

Verhamme, A., Orlitová, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7

Verhamme, A., Orlitová, I., Schaerer, D., et al. 2017, A&A, 597, A13

Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397

Virtanen, P., Gommers, R., Oliphant, T., et al. 2020, NatMe, 17, 261

Wakker, B. P., & Savage, B. D. 2009, ApJS, 182, 378

Wilde, M. C., Werk, J. K., Burchett, J. N., et al. 2021, ApJ, 912, 9

Wisotzki, L., Bacon, R., Brinchman, J., et al. 2018, Natur, 563, E31

Zabl, J., Bouche, N. F., Wisotzki, L., et al. 2021, MNRAS, 507, 4294