The H i covering fraction of Lyman Limit Systems in FIRE haloes

Abstract

Atomic hydrogen (H i) serves a crucial role in connecting galactic-scale properties such as star formation with the large-scale structure of the Universe. While recent numerical simulations have successfully matched the observed covering fraction of H i near Lyman Break Galaxies (LBGs) and in the foreground of luminous quasars at redshifts $z\lesssim 3$, the low-mass end remains as-of-yet unexplored in observational and computational surveys. We employ a cosmological, hydrodynamical simulation (FIREbox) supplemented with zoom-in simulations (MassiveFIRE) from the Feedback In Realistic Environments (FIRE) project to investigate the H i covering fraction of Lyman Limit Systems ($N_{\text{H\,{i}}}\gtrsim 10^{17.2}$ cm^-2) across a wide range of redshifts ($z=0-6$) and halo masses ($10^{8}-10^{13}\,M_{\sun}$ at $z=0$, $10^{8}-10^{11}\,M_{\sun}$ at $z=6$) in the absence of feedback from active galactic nuclei. We find that the covering fraction inside haloes exhibits a strong increase with redshift, with only a weak dependence on halo mass for higher-mass haloes. For massive haloes ($M_{\mathrm{vir}}\sim 10^{11}-10^{12}M_{\sun}$), the radial profiles showcase scale-invariance and remain independent of mass. The radial dependence is well-captured by a fitting function. The covering fractions in our simulations are in good agreement with measurements of the covering fraction in LBGs. Our comprehensive analysis unveils a complex dependence with redshift and halo mass for haloes with $M_{\mathrm{vir}}\lesssim 10^{10}M_{\sun}$ that future observations aim to constrain, providing key insights into the physics of structure formation and gas assembly.

1 Introduction

In the concordance $\Lambda$CDM cosmology, the large-scale structure of the Universe results from the collapse of a collisionless, cold dark matter and subsequent clustering of baryons. Cosmological simulations show that initial dark matter overdensities hierarchically assemble into virialized structures called dark matter haloes which are interconnected into a filamentary structure known as the cosmic web (see Peebles, 1980; Efstathiou & Silk, 1983). Baryons later fall into the potential wells mapped out by these haloes and become gravitationally bound to them. The dark matter haloes hence form the building blocks for large-scale structure formation in which galaxies and clusters of galaxies are born (e.g., White & Rees, 1978; White & Frenk, 1991; Guo et al., 2010; Wechsler & Tinker, 2018)

Galaxies continue interacting with gas from the intergalactic medium (IGM) during their lifetime, within an intricate combination of accretion and feedback processes. To sustain their growth, they need to obtain fresh gas from the IGM (e.g. Kereš et al., 2005; Bauermeister et al., 2010). The nature of this supply is strongly dependent on redshift and halo mass. For instance, the gas is shock-heated in massive haloes and takes a long time before settling in the galactic disk (e.g., Binney, 1977; Rees & Ostriker, 1977; Birnboim & Dekel, 2003). In less massive haloes, much of the accretion occurs instead via the cold mode (e.g., Dekel et al., 2009; Faucher-Giguère et al., 2011; Ho et al., 2019; Stern et al., 2020). This gas can collapse into molecular clouds which then become stellar nurseries (e.g., Hayashi & Nakano, 1965; Bromm et al., 2002; Dessauges-Zavadsky et al., 2019), or fall into the center of galaxies and ignite active galactic nuclei (AGN; e.g. Antonucci, 1993; Urry & Padovani, 1995; Padovani et al., 2017). Feedback processes regulate star formation and gas accretion by launching powerful outflows into the surrounding regions of these galaxies, affecting their dynamics and morphology (see Tumlinson et al., 2017; Hopkins et al., 2018; Biernacki & Teyssier, 2018; Valentini et al., 2019; Faucher-Giguère & Oh, 2023). The complex distribution and physics of gas around galaxies thus contains the fingerprint of galactic properties. Studying the absorption features of elements such as hydrogen and heavier metals in the spectra of bright background sources allows us to understand the principal physical processes governing galactic properties.

Advances in both simulations and observational campaigns have led to significant advances in comprehending stellar properties and molecular gas within galaxies (e.g. Guglielmo et al., 2015; Aravena et al., 2016; Le Fèvre et al., 2020; Tacconi et al., 2020; Feldmann, 2020). Much is still unknown about atomic hydrogen H i, specifically at higher redshifts, due to its lack of allowed transitions at cold temperatures resulting in a challenging detection process. As the most abundant element in the Universe, mapping its distribution through cosmic time promises to offer crucial constraints on galactic evolution and cosmology (see e.g., Padmanabhan, 2017; Dutta, 2019). In the low redshift Universe ($z\lesssim 1$), the (highly forbidden) 21cm line can be used to directly map the H i distribution (e.g., Kirby et al., 2012; Reeves et al., 2015). Current and future observing campaigns with improved sensitivity, such as MeerKAT (Booth et al., 2009; Jonas & MeerKAT Team, 2016) and SKA (Weltman et al., 2020), will use the 21cm emission to also map the distribution of hydrogen at higher redshifts. At these redshifts, absorption lines in the spectra of bright background sources are generally used to investigate the presence of atomic hydrogen along specific sight-lines (see e.g., Altay et al., 2011; Glowacki et al., 2019). These studies showed that H i can be found across a large range in column densities, including Lyman Limit Systems (LLSs; with column density $N_{\text{H\,{i}}}>10^{17.2}$ cm^-2), see e.g. Noterdaeme et al. (2012); Crighton et al. (2015); Padmanabhan (2017).

Current simulations predict that a substantial fraction of accreted material that enters haloes is relatively cold and therefore could contain considerable amounts of neutral gas (Fumagalli et al., 2011; Faucher-Giguère & Kereš, 2011; van de Voort et al., 2012; Nelson et al., 2013; Fumagalli et al., 2014). Furthermore, the number and column density of H i absorbers is predicted to increase closer to galaxy centers, suggesting that absorbers with high H i column density are better probes of gas in the proximity of galaxies (e.g., Rahmati et al., 2015; Diemer et al., 2019; Stern et al., 2021). However, strong H i absorbers such as LLSs are predicted to be often close to galaxies that may be too faint to be detected in actual surveys (Rahmati & Schaye, 2014). The study of strong H i absorbers around bright galaxies residing in massive haloes ($\geq 10^{12}M_{\sun}$) can help to overcome this problem. In fact, many modern observations and simulations make use of this galaxy-centered approach to measure covering fractions of neutral hydrogen clouds (e.g. Chen & Mulchaey, 2009; Rakic et al., 2012; Tumlinson et al., 2013; Prochaska et al., 2013b; Turner et al., 2014; Rubin et al., 2015; Prochaska et al., 2017).

The observational constraints have motivated several groups to investigate the distribution of H i surrounding galaxies via the use of simulations (for instance, Fumagalli et al., 2011; Fumagalli et al., 2014; Shen et al., 2013; Faucher-Giguère et al., 2015; Meiksin et al., 2015, 2017; Gutcke et al., 2017; Suresh et al., 2019; Nelson et al., 2020; Garratt-Smithson et al., 2021; Stern et al., 2021; Weng et al., 2023). Historically, the high observed covering fractions reported by Rudie et al. (2012) around Lyman Break Galaxies (LBGs) and by Prochaska et al. (2013a); Prochaska et al. (2013b) in the vicinity of quasars (QSOs) have been a challenge to reproduce in cosmological simulations (see e.g. Faucher-Giguère & Kereš, 2011; Fumagalli et al., 2014; Faucher-Giguère et al., 2015). These simulations are based on zoom-ins that focus on one galaxy (Faucher-Giguère & Kereš, 2011; Shen et al., 2013) or include only a few galaxies with a limited range of masses and redshifts (Fumagalli et al., 2014; Faucher-Giguère et al., 2015). There are several caveats to these approaches. For instance one could expect, given the diversity of the observed objects, that a large sample of simulated galaxies is required to accurately compare with observations of the H i distribution (Rahmati et al., 2015). Furthermore, other constraints such as the cosmic distribution of H i should be satisfied. Finally, constraints on current instrumentation limit observations to massive haloes of over $10^{12}M_{\sun}$, hence requiring a large volume or a high number of zoom-ins to be able to obtain a meaningful statistical distribution of these systems in simulations (Altay et al., 2011; Barnes et al., 2020).

Past works have demonstrated that accurately replicating the properties of the circumgalactic medium (CGM) and the resulting gas covering fractions of galaxies is complicated by the effects of resolution and the details and implementation of both star formation and feedback mechanisms (e.g. Faucher-Giguère et al., 2015; Suresh et al., 2015; Rahmati et al., 2015; van de Voort et al., 2019; Sorini et al., 2020). Resolution is particularly critical, as the quantity of cold gas in the CGM is not converged (see e.g. Faucher-Giguère & Oh, 2023; Ramesh & Nelson, 2023). Although the precise details vary, simulations with a generally stronger stellar feedback implementation have produced consistently higher values of covering fraction (Faucher-Giguère et al., 2015, 2016; Rahmati et al., 2015), and hence reconciled the results found by Rudie et al. (2012) around LBGs, as opposed to works with weaker feedback (Faucher-Giguère & Kereš, 2011; Fumagalli et al., 2014). Observational data from the COS-Halos survey (Tumlinson et al., 2013; Prochaska et al., 2017) and simulations (Faucher-Giguère et al., 2015; Rahmati et al., 2015) show little correlation between the instantaneous star-forming activity or specific star-formation rate (sSFR) and column density of H i. On the other hand, the precise role of AGN feedback in shaping the CGM cool gas distribution remains open to discussion, with the existing literature citing either insignificant or important effects when including them (Rahmati et al., 2015; Weng et al., 2023; Khrykin et al., 2024). While AGNs are ubiquitous in reality and their various feedback mechanisms are necessary to construct a full picture of galaxy formation, how to best model them in cosmological simulations remains uncertain, and the current state-of-the-art models necessarily introduce modelling degeneracies which significantly hinder predictive power.

The most recent studies on the covering fraction of H i have produced a robust agreement with results from both the LBGs and QSOs observations, mainly by remedying the hindering factors mentioned above and studying a great number of simulated objects in large cosmological volumes and concluding on the average distribution of H i. Rahmati et al. (2015) surveyed the covering fraction of atomic hydrogen over a large range of masses ($M_{200}\in[10^{11}-10^{13.7}]M_{\odot}$) in the full-scale cosmological simulations EAGLE (Crain et al., 2015; Schaye et al., 2015), and found an agreement with the radial distribution of H i around QSOs reported by (Prochaska et al., 2013b). Additionally, Faucher-Giguère et al. (2016) complemented their previous works with multiple higher resolution zoom-ins and found their results for the most massive haloes in their simulations ($M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\sun}$) to also be consistent with the covering fractions observed around quasars.

In this work, we offer a description of covering fractions of Lyman Limit Systems using the high-resolution cosmological FIREbox simulation (Feldmann et al., 2023) and a series of zoom-ins from the MassiveFIRE suite (Feldmann et al., 2016, 2017) rerun with FIRE-2 physics (Anglés-Alcázar et al., 2017). This allows us to resolve haloes with very-low mass (from $M_{\mathrm{vir}}=10^{7.75}M_{\sun}$) and investigate the profiles of covering fraction from a much lower mass range than that of other simulations. The analysis is done for redshifts $0\leq z\leq 6$, offering a more consistent investigation of the redshift dependence of covering fractions. The properties of these simulations mean we have the high resolution necessary to resolve the small-scale gas structure around haloes, while also ensuring that we have a large enough sample of haloes to obtain systematic information about them.

The paper is structured as follows. In section  2, we introduce the cosmological simulations and the methodology for our analysis. In section  3 we present the relevant results and discuss our findings. We discuss our results in the context of other works and observations in section  4. We finally conclude in section  5.

2 METHODOLOGY

2.1 Simulations

FIREbox is a high-resolution hydrodynamic cosmological volume simulation with a box size of 22.1 cMpc (Feldmann et al., 2023). FIREbox is part of the Feedback In Realistic Environments (FIRE) project (Hopkins et al., 2014; Hopkins et al., 2018, 2023). In the next paragraphs, we briefly discuss the details of the simulation.

The simulation volume contains 1024³ dark matter particles and 1024³ gas particles at the initial redshift ($z=120$). Dark matter and baryon masses are $m_{\mathrm{DM}}=3.35\times 10^{5}M_{\sun}$ and $m_{\mathrm{b}}=6.26\times 10^{4}M_{\sun}$ respectively. Dark matter (star) particles have a fixed softening length of 80 pc (12 pc). Gas softening is adaptive with a minimum softening length of 1.5 pc. Initial conditions were generated with the MUSIC (MUlti-Scale Initial Conditions) code (Hahn & Abel, 2011) and with 2015 Planck cosmological parameters (Ade et al., 2016): $H_{0}=67.74$ km/s/Mpc (or $h=0.6774$), $\Omega_{\mathrm{m}}=0.3089$, $\Omega_{\Lambda}=0.6911$ , $\Omega_{\mathrm{b}}=0.0486$ , $\sigma_{8}=0.8159$ and $n_{s}=0.9667$.

FIREbox is run with the gravity-hydrodynamics solver GIZMO (Hopkins, 2015) using the FIRE-2 physics model (Hopkins et al., 2018). The simulation incorporates multiple gas-cooling processes (such as: free-free, photo-ionization, recombination, Compton, photoelectric, metal-line, molecular, fine-structure, dust collisional, and cosmic ray physics) following an implicit algorithm described in Hopkins et al. (2018). FIREbox includes radiative feedback in the form of photo-ionization and photoelectric heating. Radiative transfer effects are accounted for in our simulations via the Locally Extincted Background Radiation in Optically thin Networks approximation (LEBRON; Hopkins et al., 2012; Hopkins et al., 2014; Hopkins et al., 2018; Hopkins & Grudić, 2019). To compute the H i fraction in each resolution element, the model assumes equilibrium between recombination, collisional ionization and photoionization by local stellar sources. The process of self-shielding is taken into account using a local Sobolev/Jeans-length approximation which is calibrated from radiative transfer experiments (Faucher-Giguère et al., 2010; Rahmati et al., 2013). Approximating the attenuation of incident flux in this manner was shown to reproduce neutral fractions calculated with full radiative transfer codes in post-processing (see Faucher-Giguère et al., 2010, 2015; Rahmati et al., 2013; Stern et al., 2021, for further details). Relevant metal ionization states are tabulated from the CLOUDY simulations (Ferland et al., 1998). Feedback from AGNs is not included. We refer the interested reader to Hopkins et al. (2018) for the full description of the FIRE-2 simulations.

To offer a more complete comparison of our results with observing campaigns, we supplement our analysis of the higher mass end of haloes with four zoom-in simulations (A1, A2, A4, A8) from the MassiveFIRE (Feldmann et al., 2016, 2017) suite, re-simulated with FIRE-2 physics (Anglés-Alcázar et al., 2017). We briefly outline the selection criteria of the zoom-ins. Isolated haloes with $M_{\mathrm{vir}}\,(z=2)\approx 10^{12.5}M_{\sun}$ were selected from low-resolution, dark matter-only runs. The environmental density of each halo was then measured by computing the total mass in a sphere of radius 1.8 Mpc. The final haloes of the zoom-in suite were selected from the 5th, 25th, 50th, 75th and 95th percentile of the distribution of environmental density. This choice results in the haloes having varied accretion histories (see Feldmann et al., 2017, Section 2, for details). The masses of particles are $m_{\mathrm{DM}}=1.76\times 10^{5}M_{\sun}$ and $m_{\mathrm{b}}=3.29\times 10^{4}M_{\sun}$, and dark matter (gas) particles have softening lengths of 143 pc (9 pc).

2.2 Covering fractions

Our goal is to compute the H i covering fraction of the simulated haloes. We define this quantity and outline the procedure to obtain column density maps from simulation snapshots in the next paragraphs.

The covering fraction is used to quantify the distribution of atomic hydrogen around haloes. In our study, we mainly consider the covering fraction of strong H i absorbers, specifically LLSs (with $N_{\text{H\,{i}}}>10^{17.2}$cm^-2) without further separating them from Damped Lyman-$\alpha$ systems (with $N_{\text{H\,{i}}}>10^{20.3}$cm^-2). We use two separate measures to quantify the H i distribution around haloes: the cumulative and differential covering fractions, respectively denoted by $f_{\mathrm{cov}}(<R)$ and $f_{\mathrm{cov}}(r)$. They are related according to the following formula:

$$f_{\mathrm{cov}}(<R)=\frac{\int_{0}^{R}f_{\mathrm{cov}}(r)\cdot 2\pi r\cdot dr}{\pi\cdot R^{2}}\equiv\frac{A_{\mathrm{abs}}\,(<R)}{\pi\cdot R^{2}}$$ (1)

where $A_{\mathrm{abs}}\,(<R_{\mathrm{vir}})$ is the area covered by LLSs within the field of view defined by the virial radius of a given halo. The cumulative covering fraction of LLSs essentially measures the probability of finding such systems in line-of-sights within some radius $R$.
The differential covering fraction is used to quantify the spatial distribution of LLSs around haloes and complements the previous measure. Effectively, it is computed according to:

$$f_{\mathrm{cov}}(r)=\frac{A_{\mathrm{abs}}\,(r_{i}<r<r_{i+1})}{\pi\cdot(r_{i+1}^{2}-r_{i}^{2})}$$ (2)

with $A_{\mathrm{abs}}\,(r_{i}<R<r_{i+1})$ being the area covered by LLSs within an annulus of inner and outer impact parameters $r_{i}$ and $r_{i+1}$ respectively. The $r_{i}$’s are consistently chosen such that $r_{0}=0$, and we decided to use more annuli closer to the center of the halo. This choice is motivated by previous works and ensures we can capture details of the shape of the differential profile.

2.3 Data generation

We use the snapshots corresponding to $z=0,1,2,2.5,3,4,5,6$ from the FIREbox simulation for our study. In order to obtain covering fractions, we compute column density maps by depositing particle densities from the simulated cube onto a uniform grid using the smooth and tipgrid algorithms. For each particle, smooth calculates a smoothing length specified as half the distance to its $n^{\mathrm{th}}$ neighbour particle. We set $n=10$ for this work. We note that lower values of $n$ translate into higher particle noise but allow us to better resolve small structures. Following this, we subdivide our simulated box into 20 equally spaced slabs (of thickness $\sim$1.1 cMpc), along each of the 3 spatial directions. The tipgrid algorithm then interpolates particles within the same slab onto a two-dimensional grid, by distributing the mass using a spherically-symmetric kernel according to the aforementioned smoothing lengths. We chose a grid resolution of 32768² pixels, such that individual pixels resolve roughly 0.68 ckpc. This choice ensures we can observe the finer details in the structures formed in the simulation.

We identify simulated haloes and recover their main properties using the AMIGA Halo Finder (AHF; Gill et al., 2004; Knollmann & Knebe, 2009). We include all haloes with at least 168 dark matter particles in our analysis. This means that all AHF haloes with roughly $M_{\mathrm{vir}}\gtrsim 10^{7.75}M_{\sun}$ are included, which enables us to study the behaviour for much smaller structures and masses than done in previous works. The virial mass ($M_{\mathrm{vir}}$) and virial radius ($R_{\mathrm{vir}}$) of dark matter haloes are computed following the virial overdensity definition outlined in Bryan & Norman (1998). Using the AHF particle files, we identify haloes by their center and virial radius and distinguish between ‘main haloes’ and ‘sub-haloes’. In this work, the term ‘sub-haloes’ refers to dark matter haloes that are nested within another dark matter halo. All other haloes identified by AHF are ‘main haloes’. For sub-haloes, the virial radius $R_{\mathrm{vir}}$ reported by AHF and used throughout this work refers to the smaller of the virial and tidal radius (see, e.g., van den Bosch et al., 2018). In particular, the tidal radius is significantly smaller than the virial radius for sub-haloes close to the center of the host main halo.

There are roughly 96000 haloes at redshift $z=6$, growing to a peak of 160000 haloes at redshift $z=2$ and finally 130000 at $z=0$. The main-to-sub proportion evolves from 95%-5% at $z=6$ to 80%-20% at $z=0$. This selection offers a statistically sound sample of covering fractions for all halo masses below some redshift-dependent higher-mass of $M_{\mathrm{vir}}\sim 10^{11}M_{\sun}$ at redshift $z=6$, up to $M_{\mathrm{vir}}\sim 10^{13}M_{\sun}$ at redshift $z=0$.

For our study, we always evaluate the covering fraction in Eq. (1) for $R=R_{\mathrm{vir}}$. It is computed for each halo according to the following straightforward procedure. We place a circular mask around haloes so that only pixels within the virial radius are considered. We then place another circular mask so that only those pixels with column density above the threshold are counted. The ratio of the two gives the covering fraction. The differential covering fraction is obtained in the same way, with annuli being used instead of disks. We repeat this procedure for all redshifts and for all 3 orientations of the simulated cube, and compute statistics using each orientation as an independent data point for our final results.

3 RESULTS

3.1 Cumulative covering fraction of LLSs

Figure 1: H i column density distribution in surroundings of randomly selected main haloes, at redshift $z=3$. Top, middle and bottom rows show haloes with $M_{\mathrm{vir}}\approx 10^{12},10^{11.5}$ and $10^{11}\,M_{\sun}$ respectively. The columns in each row show the same halo from three different orthogonal projections. The panels show a region of size 1 $\times$ 1 cMpc², with the same depth of projection. The full and dotted circles in white are centred on the haloes and display their virial radius $R_{\mathrm{vir}}$ and 10% of their virial radius, respectively. Black contours highlight LLS sightlines (with $N_{\text{H\,{i}}}>10^{17.2}\mathrm{cm}^{-2}$), and the covering fraction of LLSs within $R_{\mathrm{vir}}$ is indicated in the top-right of each panel. The covering fraction does not evolve strongly with mass for $M_{\mathrm{vir}}\geq 10^{11}M_{\sun}$ haloes, but can significantly differ from one orthogonal projection to another of the same halo.

Figure 2: H i column density distribution in surroundings of randomly selected main haloes. The rows show haloes with $M_{\mathrm{vir}}\approx 10^{11.5}M_{\sun}$ at redshifts $z=6,4,2$ and $0$ (respectively, from top to bottom). The columns in each row show the same halo from three different orthogonal projections. The panels show a region of size 1 $\times$ 1 cMpc², with the same depth of projection. The full and dotted circles in white are centred on the haloes and display their virial radius $R_{\mathrm{vir}}$ and 10% of their virial radius, respectively. Black contours highlight LLS sightlines (with $N_{\text{H\,{i}}}>10^{17.2}\mathrm{cm}^{-2}$), and the covering fraction of LLSs within $R_{\mathrm{vir}}$ is indicated in the top-right of each panel. The filamentary structure of H i in and around FIREbox haloes evolves very strongly with redshift, from almost complete coverage of the virial cross-section at redshift $z=6$ to being highly centrally contracted within haloes at redshift $z=0$. The covering fraction thus increases significantly with increasing redshift, and we find that this is the predominant parameter affecting values of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$.

In the following paragraphs, we begin investigating the covering fraction of Lyman Limit Systems in FIREbox haloes by directly taking a look at some snapshots from the simulations and performing a visual inspection of the pictures. Figures  1 and  2 show examples of the distribution of H i around randomly-chosen main haloes, with either changing mass or redshift.

Fig. 1 shows haloes of different mass from different orthogonal projections at redshift $z=3$. We see that the gas distribution is highly inhomogeneous, with a filamentary structure that can extend beyond the virial radius of the haloes. The shape and extent of this structure also change significantly when viewing the same halo from different angles.

We note that a significant fraction of the projected area within one virial radius of the center of the haloes is covered by LLSs. The covering fraction of the haloes does not seem to be strongly correlated with virial mass: the gas around $M_{\mathrm{vir}}=10^{12}M_{\sun}$ haloes covers a much greater area than that of $M_{\mathrm{vir}}=10^{11}M_{\sun}$ haloes and the virial radius of $10^{12}M_{\sun}$ haloes is also much bigger than that of $10^{11}M_{\sun}$ haloes; the two effects are comparable in magnitude so that $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ is about the same for these haloes. The inhomogeneity of the gas distribution noted above results in different values of covering fraction for different projections of the same halo. Specifically, the individual value of H i covering fraction of a halo can vary up to a factor of $\sim 2$ from one angle of viewing to another.

The pictures can also be used to visually highlight any redshift dependence of the distribution of neutral hydrogen around haloes. Fig.  2 shows haloes with mass $M_{\mathrm{vir}}\approx 10^{11.5}M_{\sun}$ at 4 different redshifts used in the study for three different orthogonal projections each. We find that the filamentary structure described in the previous figure evolves very strongly with redshift. At $z=6$ we observe large filaments and clumps of atomic gas that extend far beyond the virial radius of the studied haloes, interlinking huge regions of space. Such structure diminishes in extent and compactness with redshift, until they have almost disappeared by $z=0$. This is found to be consistent with the reports that at $z\lesssim 2$, essentially all the H i is found inside dark matter haloes (e.g., Villaescusa-Navarro et al., 2018; Feldmann et al., 2023). Consequentially, the covering fraction increases very strongly with increasing redshift: only around 3% of a halo’s virial radius is covered in LLSs at $z=0$, whereas almost 90% is covered at $z=6$. We find this redshift relation to be the predominant parameter in determining the covering fraction of H i gas for a halo.

Figure 3: Cumulative covering fraction of LLSs within impact parameter $R_{\mathrm{vir}}$ as a function of halo mass $M_{\mathrm{vir}}$ (left column) and specific star formation rate (sSFR) $\dot{M}_{\star}/M_{\star}$ (right column). We distinguish between main- (top row) and sub- (bottom row) haloes. Different colours are used to distinguish between redshift from $z=0$ to $z=6$ as per the legend. Solid lines indicate the median covering fraction and shaded areas show the 5th-95th percentile error on the median obtained from bootstrapping. All haloes with $M_{\mathrm{vir}}\geq 10^{7.75}M_{\sun}$ identified in the simulation were used to create these plots. The behaviour of main- and sub-haloes is very different and details are discussed in the main text. The covering fraction of LLSs increases strongly with redshift and (for low-mass haloes) with halo mass. The covering fraction is (weakly) anticorrelated with sSFR in main haloes.

We proceed to a systematic statistical study of covering fractions over a wider range of masses and redshifts. Figure  3 shows the measured $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ for all sampled haloes at $z=0-6$. The covering fraction is studied as a function of halo virial mass $M_{\mathrm{vir}}$ and specific star-formation rate $\dot{M}_{\star}\,/\,M_{\star}$, for all resolved haloes with $M_{\mathrm{vir}}>10^{7.75}M_{\sun}$. We distinguish between main- and sub-haloes. The haloes are grouped in equidistant mass bins, and for each one we show the median covering fraction and highlight the 5th-95th percentile error on the median found via bootstrapping. We find that the results significantly change between main- and sub-haloes.

For the main haloes (top row of Fig.  3), we conclude that the median covering fraction of LLSs strongly increases with increasing redshift, at all halo masses. We note that there is a particularly noticeable rise from $z=2$ to $z=3$ and all redshifts thereafter, as compared to the evolution between redshifts $z=0,1$ and 2. Our results are consistent with the idea that haloes contain higher gas fractions at high redshift, due to increased accretion of gas and a higher mean density of the Universe (see e.g. Rahmati et al., 2015). This finding is of particular importance for comparisons with observations. Indeed, since observed samples contain galaxies with a wide range of redshifts, the observed probability of finding an LLS sightline within a given impact parameter is not directly equal to the covering fraction of LLSs at the mean redshift of the sample, because higher-redshift galaxies contribute more to the covering fraction than lower-redshift ones (Rahmati et al., 2015). The predicted scatter of covering fractions at all redshifts further accentuates this issue and magnifies errors, underscoring that even minor inaccuracies in redshift estimation can lead to significant misjudgements of covering fractions from observations. Finally, we should mention that the virial radius of haloes cannot be directly observed and predominantly serves as a parameter in theoretical and numerical studies, meaning that $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ as shown in Fig.  3 is inherently challenging to compare with observations.

The complex halo-mass dependence of the covering fraction of main-haloes can also be studied with Fig.  3. At lower masses, the covering fraction is close to zero until some threshold mass between $10^{8}M_{\sun}$ and $10^{10}M_{\sun}$ depending on redshift. The covering fraction then rapidly increases with halo mass, at all redshifts. At higher masses, the covering fraction does not evolve strongly with halo mass and eventually plateaus at some value, which is higher for higher redshifts. We conclude that the covering fraction of massive ($M_{\mathrm{vir}}\approx 10^{11}-10^{13}M_{\sun}$, depending on redshift) haloes is nearly independent of halo mass, at any given redshift. We conduct a deeper analysis of this mass-dependence when studying the H i differential covering fraction of haloes in the following subsection.

As the top-right panel of Fig.  3 highlights, the covering fraction is weakly anticorrelated with sSFR in main haloes. Together with results from Faucher-Giguère et al. (2015), this suggests that the covering fraction is not affected by the details of star formation, namely the instantaneous star formation rate, in haloes. On longer timescales, as highlighted previously, stellar feedback is very important for shaping the CGM and its properties such as the covering fraction of atomic gas (see e.g. Faucher-Giguère et al., 2016). This argument more easily enables comparisons of predicted covering fractions because they do not need to be compared with galaxies in the exact same stages of star formation. This immediate result also strengthens our conclusion that covering fraction strongly increases with redshift and that this is the main parameter with which it varies.

For the sub-haloes (bottom row of Fig.  3), it is also found that $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ increases with redshift. We note that the covering fraction is generally higher for sub-haloes than main haloes, albeit with significant scatter, at the same virial mass. This is likely due to a geometric effect, whereby the sub-haloes can be entirely covered by gas found within the virial radius of more massive main haloes. As such, some sub-haloes can be absorbed ‘into’ or pushed ‘out’ of the field of view of the main halo when viewed from different orientations. This usually results in a notable increase of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$, but also produces significant scatter as is seen by the large shaded areas around the median in Fig.  3. The scatter can be explained by a combination of both the randomness of projection effects discussed above and the smaller sample size ($\sim 4-20$ times fewer sub- than main-haloes). The bottom-right panel of Fig.  3 suggests that the H i covering fraction of sub-haloes is dependent on their specific star formation rate. We can explain this apparent dependence by noting that for sub-haloes $R_{\mathrm{vir}}$ refers to the tidal radius, which depends on the distance $d_{\mathrm{host}}$ to the center of the main halo. Specifically, for sub-haloes located far from the host, the tidal radius approaches the definition of the virial radius of main haloes, whereas for those close to the host, the tidal radius is significantly smaller. Consequently, sub-haloes with $d_{\mathrm{host}}>0.7\,R_{\mathrm{vir,host}}$ exhibit trends of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ as a function of sSFR that are very similar to those of main haloes. While for sub-haloes closer to the host, the much smaller tidal radius results in median covering fractions that approach unity. The characteristic shape seen in the bottom-right panel of Fig.  3 emerges from the combination of these two trends for all sub-haloes.

3.2 Differential covering fraction of LLSs

Figure 4: Differential covering fraction profiles of LLSs as a function of normalised impact parameter for four different mass bins. The different lines indicate the median $f_{\mathrm{cov}}(r)$ at each impact parameter for different redshifts, and shaded areas show the 5th-95th percentile obtained from bootstrapping. The central mass of the bins is shown in the top right of each panel. All haloes with $M_{\mathrm{vir}}\pm 0.5\,(0.3)$ dex are included in the $M_{\mathrm{vir}}=10^{11};10^{12}\,(10^{10})\,M_{\odot}$ bins, while only haloes with $M_{\mathrm{vir}}>10^{12.5}M_{\sun}$ are chosen for the bottom-right panel. The curves systematically shift to the left for decreasing redshift at all mass bins, indicating that the covering fraction of LLSs drops at all impact parameters for decreasing $z$.

Figure 5: Differential covering fraction profiles of LLSs as a function of normalised impact parameters for redshifts $z\geq 2$. The solid lines indicate the median $f_{\mathrm{cov}}(r)$ at each impact parameter for the different mass bins and the shaded areas show the 5th-95th percentile error on the median obtained from bootstrapping. The redshift is indicated at the bottom left of each panel. The reference profile of $M_{\mathrm{vir}}=10^{11}M_{\sun}$ haloes at redshift $z=3$ is shown as a black-dashed line in every panel to aid in analysing the general trends. We see a strong increase in covering fraction at all impact parameters with increasing redshift. At these redshifts, we note that the differential covering fraction does not depend strongly on mass for haloes with $M_{\mathrm{vir}}=10^{11}-10^{12}M_{\sun}$, which hints at scale-invariance of the profiles. This does not hold for lower-mass ($M_{\mathrm{vir}}\leq 10^{10}M_{\sun}$) haloes, for which we note that the differential covering fraction profiles are generally less extended and evolve with halo mass.

In the previous section, we investigated the cumulative covering fractions of haloes. We now look at how the spatial distribution of LLSs around halo centers influences the covering fraction, by studying the differential covering fraction. We show the profiles of $f_{\mathrm{cov}}(r)$, disentangling the effects of mass and redshift, in Figures  4 &  5 respectively. We choose to depict the differential covering fraction as a function of normalised impact parameter $r\ /\ R_{\mathrm{vir}}$ guided by the previous study in Rahmati et al. (2015).

Fig.  4 shows the predicted median differential covering fraction for four mass bins, distinguishing between different redshifts. Each bin is centered at the indicated value and includes all haloes within $M_{\mathrm{vir}}\,\pm\,0.5$ dex for $M_{\mathrm{vir}}=10^{11};10^{12}M_{\odot}$ haloes and within $M_{\mathrm{vir}}\,\pm\,0.3$ dex for $M_{\mathrm{vir}}=10^{10}M_{\odot}$ haloes. For the last bin ($M_{\mathrm{vir}}>10^{12.5}M_{\sun}$), only haloes with mass above the threshold are chosen. We note that the profiles show higher values of differential covering fraction for all impact parameters with increasing redshift, for all mass bins. This result is consistent with the previous section and Feldmann et al. (2023), namely that covering fractions of atomic hydrogen increase with increasing redshift. For instance, the pictures of the haloes show that H i clouds extend far outside the virial radius of haloes for higher redshifts (see Fig.  2, wherein we observe the filamentary structure around haloes to be more extended at higher redshifts). The differential covering fraction is thus expected to increase with redshift. It is noted that the radial distribution of H i around haloes has a sigmoidal shape, with some asymmetries that are more pronounced at higher redshift. As evident in the top-right panel of Fig.  4, at redshift $z=0$, the median $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ rapidly goes from 1 to 0 as the normalised impact parameter increases, with steep turning points. On the other hand, we see that for redshift $z=6$ the median $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ steeply decreases from 1, but goes down more slowly towards 0 for large impact parameters away from the center.

In Fig.  5 we show the predicted differential covering fraction for individual redshifts $z=2,3,4$ and 6, and classify the haloes in each panel by mass bins as in Fig.  4. We include two lower-mass bins for this analysis, namely all haloes within $M_{\mathrm{vir}}\pm\,0.1$ dex for $M_{\mathrm{vir}}=10^{9}\ M_{\sun}$ haloes and within $M_{\mathrm{vir}}\,\pm\,0.2$ dex for $M_{\mathrm{vir}}=10^{9.5}\ M_{\sun}$ haloes. At these redshifts, we find that the profiles of haloes with $M_{\mathrm{vir}}\simeq 10^{11}-10^{12}M_{\sun}$ take roughly the same values at all impact parameters. The curves are almost superimposed, hinting at the existence of some characteristic length scale similar to the virial radius for this mass of haloes (Rahmati & Schaye, 2014; Rahmati et al., 2015). This explains the weaker dependence (i.e. flattening) of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ noted in haloes with $M_{\mathrm{vir}}\gtrsim 10^{11}M_{\sun}$ (Fig.  3). Although the specific geometry and physical scale of individual haloes are different, we find that the gas is distributed to a similar relative extent away from their center. This hints at scale invariance, which is studied in the following subsection. We note that a similar trend is found for redshifts $z=0$ and 1 (not shown in Fig.  5), although with lesser agreement.

The scale-invariance observed in more massive haloes does not seem to hold for $M_{\mathrm{vir}}\leq 10^{10}M_{\sun}$ haloes. This explains the mass-dependence of the cumulative covering fraction of $M_{\mathrm{vir}}\lesssim 10^{11}M_{\sun}$ haloes (see Fig.  3). Indeed, the LLS differential covering fraction profiles of $10^{9}-10^{10}M_{\sun}$ haloes are shifted to smaller radii, indicating that neutral hydrogen is less extended in lower mass haloes and hence the cumulative covering fraction is also smaller. As redshift increases, we find that the radial profiles of $10^{9}-10^{10}M_{\sun}$ haloes evolve from being zero at all radii to gradually shifting to higher radii, until they nearly converge with those of massive haloes at $z=6$. The evolution of the cumulative covering fraction of LLSs with halo mass is thus explained by the evolving radial concentration of H i with halo mass. A further investigation reveals that the H i column density in lower mass haloes is distributed similarly to that in more massive haloes in the outskirts ($\gtrsim R_{\mathrm{vir}}$) but is significantly lower in the central regions, see Fig.  11 in the appendix. Consequently, the density threshold of Lyman Limit Systems is reached at smaller radii in lower mass haloes and the resulting differential covering fraction profiles are thus less extended. This suggests that gas which is located in the center of more massive haloes is expulsed from the halo and redistributed via heating by the UV background, stellar feedback and tidal interactions in lower-mass haloes (Jaiswal & Omar, 2020; Ayromlou et al., 2023; Zheng et al., 2024).

There is tentative evidence that the scale-invariance is also not exhibited in the most massive haloes ($M_{\mathrm{vir}}\geq 10^{12.5}M_{\sun}$), particularly at redshifts $z=0,1$ (not shown here). It was found in other studies that, at this mass, the cooling time of gas in the inner parts of massive haloes is long, such that there is a low fraction of neutral gas there, which could explain the lower differential covering fraction of H i (Stern et al., 2021). However, there are only a few objects in our simulation that reach these masses at redshifts $z=0,1$ and 2 (respectively: 12, 5 and 1), and a more robust sample is needed to draw conclusions.

3.3 Fitting function for differential covering fraction

Halo mass [$M_{\sun}$]	Parameter	$a$	$b$	$c$	$d$
$M_{\mathrm{vir}}=10^{11}$	$A$	$-0.0065$	$0.092$	$-0.153$	$0.012$
	$B$	$0.0017$	$-0.016$	$0.026$	$0.998$
	$L_{z}$	$0.0008$	$0.005$	$0.084$	$0.13$
	$\alpha$	$-0.06$	$0.57$	$-1.58$	$5.86$
$M_{\mathrm{vir}}=10^{12}$	$A$	$-0.0013$	$0.054$	$-0.061$	$0.034$
	$B$	$0.0003$	$-0.0070$	$0.0086$	$0.989$
	$L_{z}$	$-0.007$	$0.045$	$0.041$	$0.085$
	$\alpha$	$0.099$	$-0.39$	$-0.017$	$4.33$

Table 1: Best-fit values for the free parameters of the differential covering fraction fitting function for LLSs of haloes with $M_{\mathrm{vir}}=10^{11}M_{\sun}$ or $10^{12}M_{\sun}$.

Figure 6: Comparison of differential covering fraction profiles obtained from the simulations (dashed lines) and the fitting function of equation (3) (solid lines). The different colours distinguish between different redshifts, and the mass bin of the haloes considered is shown at the top right of each panel. We note that the fit performs very well in recovering the features of the radial profiles and that we can attribute some physical length scale to certain parameters (see the main text for details).

Our study of the differential covering fraction of LLSs showed that haloes of a certain mass range share very similar, possibly scale-invariant, profiles. This can be further investigated by way of a generalised fitting function for the radial profiles of atomic hydrogen around haloes. Let us denote the normalised impact parameter as $x=r\ /\ R_{\mathrm{vir}}$. The differential covering fraction of LLSs around haloes with $M_{\mathrm{vir}}\simeq 10^{11}-10^{12}M_{\sun}$ at a given redshift $z$ can be fitted via:

$$f_{\mathrm{cov}}(x,z)=1+\left(\frac{A(z)}{1+x}-1\right)\cdot\frac{1}{B(z)+\left(L_{z}(z)/x\right)^{\alpha(z)}}$$ (3)

where the four free parameters $A,B,L_{z}$ and $\alpha$ are fitted by way of a 3-rd degree polynomial in redshift $z$. This fitting function is a revision of a similar characterization of differential covering fraction profiles of LLSs proposed in Rahmati et al. (2015) (see Eq. (5) therein). It satisfies important properties that a physical distribution should respect in the appropriate limits. In particular, it approaches unity in the two limits $x\to 0$ and $z\to\infty$, and approaches some asymptotic value $1-1/B$ at large impact parameters, which depends on redshift.

The parameters determining our empirical fit are physically meaningful. For instance, $L_{z}$ can be interpreted as some typical projected distance between galaxies and H i absorbers (similarly to Rahmati et al., 2015). $x\sim L_{z}$ corresponds roughly to where $f_{\mathrm{cov}}(r)\sim 0.5$, and hence the parameter $L_{z}$ can be used to estimate the distance between haloes-absorbers. $A$ dictates the first turning point where $f_{\mathrm{cov}}(r)$ decreases from 1; $B$ determines the asymptotic value of the differential covering fraction of haloes at large impact parameters; $\alpha$ roughly describes the slope of $f_{\mathrm{cov}}(r)$ (i.e., how quickly it goes from 1 to 0 as a function of $x$). All these parameters $P$ are described via $P(z)=a\cdot z^{3}+b\cdot z^{2}+b\cdot z+d$. The values for the best-fit coefficients for each parameter used in equation (3) are summarised in Table  1 and shown in Figure  6.

At all redshifts, the fitting function reproduces accurately the characteristic behaviour we described in section  3.2. We remark that it is not as precise at redshifts $z=0$ and 1, but this is expected as it corresponds to the redshifts for which the actual profiles are the least superimposed. The fitting function highlights the asymmetries of the differential covering fraction profiles observed in the previous section. Using the values listed in Table  1, we find that at $z\sim 3$ the expected projected distance between LLSs and their host haloes should be around $L_{z}\approx 0.42-0.45\ R_{\mathrm{vir}}$ for $M_{\mathrm{vir}}=10^{12}-10^{11}M_{\odot}$ haloes respectively. These values are smaller than findings from the OWLS simulations (Rahmati & Schaye, 2014) and from EAGLE (Rahmati et al., 2015), wherein it was found that such distance $L_{z}\approx R_{\mathrm{vir}}$. This suggests that H i gas in FIREbox is concentrated closer to the center of haloes than in the EAGLE simulations.

4 DISCUSSION

4.1 Cumulative covering fraction

The most recent and statistically significant constraints from observations of the H i covering fraction of Lyman Limit Systems come from the Keck Baryonic Structure Survey (Rudie et al., 2012; Steidel et al., 2014; Strom et al., 2017) and the Quasars Probing Quasars project (QPQ; see Findlay et al., 2018, and references therein). On the one hand, Rudie et al. (2012) report a H i covering fraction $f_{\mathrm{cov}}(<R_{\mathrm{vir}})=0.30\pm 0.14$ around Lyman Break Galaxies (LBGs) at $z\sim 2-2.5$, residing in haloes with $M_{\mathrm{vir}}\approx 10^{12}M_{\sun}$. On the other hand, Prochaska et al. (2013a) predict a covering fraction $f_{\mathrm{cov}}(<R_{\mathrm{vir}})=0.64^{+0.06}_{-0.07}$ around $z\sim 2-2.5$ quasars (QSOs) residing in haloes with characteristic halo mass of $M_{\mathrm{vir}}\approx 10^{12.5}M_{\sun}$ (White et al., 2012). These results have historically been a challenge to reproduce in simulations, and different suites and physical implementations lead to varying predictions.

The most recent numerical works carried out on this topic (Fumagalli et al., 2014; Rahmati et al., 2015; Faucher-Giguère et al., 2015, 2016; Meiksin et al., 2015, 2017; Gutcke et al., 2017; Suresh et al., 2019) have all been able to broadly reproduce covering fractions found around LBGs by Rudie et al. (2012), while using a vast range of numerical solvers and sub-grid physics. The high covering fractions observed around QSOs by Prochaska et al. (2013a) have posed a greater challenge to reproduce (see e.g., Fumagalli et al., 2014; Faucher-Giguère et al., 2015). Nonetheless, further work conducted by Faucher-Giguère et al. (2016) was able to replicate values for the QSOs using higher resolution zoom-ins and the same (stellar feedback driven) physics, arguing that the high resolution enabling more finely-resolved stellar feedback from satellites is a key ingredient in matching the observations. Additionally, Rahmati et al. (2015) have succeeded in matching both observations using the EAGLE suite of cosmological simulations (Crain et al., 2015; Schaye et al., 2015), via implementation of both stellar and AGN feedback at a lower numerical resolution.

We now discuss the meaning of our results and compare them with those of the previous works described above. We found that the covering fraction of LLSs in FIREbox haloes increases with increasing redshift, and is roughly independent of mass in massive ($M_{\mathrm{vir}}\gtrsim 10^{11}-10^{12}M_{\sun}$) haloes. These conclusions are largely in agreement with other simulations (see Fumagalli et al., 2014; Faucher-Giguère et al., 2015; Rahmati et al., 2015; Gutcke et al., 2017; Stern et al., 2021). Our values of H i covering fraction for the very-massive ($M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\sun}$) haloes are lower than found in other cosmological-size suites as in Rahmati et al. (2015). They report that the covering fraction in $M_{\mathrm{vir}}\approx 10^{12.5}M_{\sun}$ haloes is roughly $f_{\mathrm{cov}}(<R_{\mathrm{vir}})\approx 0.27$ at $z=2$, whereas FIREbox results at the same redshift are about $f_{\mathrm{cov}}(<R_{\mathrm{vir}})\sim 0.15-0.2$. This overall trend is the same at redshifts $z=1,3$ and 4. We note, however, that there are significantly fewer of these very-massive haloes in our simulation than in Rahmati et al. (2015). Specifically, there is only one very-massive halo in FIREbox at $z=2$, whereas EAGLE has 39 very-massive haloes at redshift $z=3$ and 116 at redshift $z=2$.

We proceed to compare our results with observations. Given the above discussions, we complement our FIREbox results for the very-massive haloes by adding 4 ‘zoom-in’ simulations of $M_{\mathrm{vir}}\,(z=2)\sim 10^{12.5}M_{\sun}$ haloes run with the same FIRE-2 physics (see section  2 for details). The results are shown in Figure  7. We highlight the LLS covering fractions predicted by FIREbox for redshifts $z=2,2.5,3,4$, the results for our zoom-ins, and the data obtained by Rudie et al. (2012); Prochaska et al. (2013a). Our results for redshifts $z=2,2.5$ are well within the confidence interval of the LBGs observations, and we predict that the sample used to obtain these covering fractions is matched similarly well by $10^{12}M_{\sun}$ haloes at redshifts $z=2,2.5$. This conclusion is consistent with other numerical works of similar scope (see e.g. Faucher-Giguère et al., 2015, 2016; Rahmati et al., 2015).

Our FIRE-2 simulations do not reproduce the high covering fraction observed in QSOs by Prochaska et al. (2013a). The zoom-in haloes presented here tend to have a larger range of covering fractions and span varied accretion histories. One halo, in particular, shows a higher average covering fraction of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})\approx 0.35$ at $z=2$, but the zoom-ins do not constitute any major improvement against the QSOs observations. This is to be contrasted with results from Faucher-Giguère et al. (2016), wherein the large H i covering fraction in very-massive haloes was attributed to the enhanced resolution of low-mass satellite galaxies and their associated winds interacting with filaments of cosmic origin in the zoom-ins.

Our investigation reveals that, on average, the covering fraction of our haloes is lower than the mean covering fraction in the sample introduced in Faucher-Giguère et al. (2016). Several factors contribute to this disparity. On the one hand, they analyzed a more extensive ensemble of 15 haloes to our limited sample of 4, granting them more robustness in analyzing mean values of covering fraction. Specifically, at redshift $z=2$ and for $M_{\mathrm{vir}}\approx 10^{12.5}M_{\odot}$ haloes, their analysis revealed a broad halo-to-halo scatter of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})\in[0.3,0.6]$ while the average exceeded $f_{\mathrm{cov}}(<R_{\mathrm{vir}})\gtrsim 0.5$, bringing them considerably closer to the observed value for QSOs. It is hence plausible that the values we obtain from our FIRE-2 simulations are on the low end of a distribution, and that selecting more haloes to simulate in zoom-ins might raise the average covering fraction at the high-mass end. Furthermore, their analysis includes very-massive haloes at higher redshift, with $M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\odot}$ at $z=2.5$. This is significant because of the notable surge in the typical covering fraction from $z=2$ to $z=2.5$, naturally producing a higher mean covering fraction when including such haloes, which we have not done in this work.

Our analysis of the cumulative covering fraction in $M_{\mathrm{vir}}<10^{12.5}M_{\odot}$ haloes demonstrates good agreement with observational data for $M_{\mathrm{vir}}=10^{12}M_{\odot}$ LBGs. This work robustly extends the statistics down to very-low masses of $M_{\mathrm{vir}}\sim 10^{8}M_{\odot}$, unveiling a complex halo-mass dependence of the H i covering fraction. Given the very low number of $M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\odot}$ haloes in our sample from FIREbox supplemented with 4 MassiveFIRE (FIRE-2) haloes, we cannot robustly compare our covering fractions with the QSOs sample and leave this for future work.

Figure 7: Cumulative covering fraction of LLSs in FIRE simulations and in observations at redshifts $z=2-4$. Solid lines: Median cumulative covering fraction from FIREbox for different mass bins. Coloured circles: Average (over three orthogonal projections) of the cumulative covering fraction of the 5 most massive haloes at each redshift in FIREbox. Coloured stars: Average (over three orthogonal projections) of the cumulative covering fraction of the 4 most massive haloes at each redshift in the MassiveFIRE (FIRE-2) runs (see main text for details). The lines and points are colour-coded by redshift as per the rest of the paper and shown in the top-left of the figure. We compare the simulations with observations, indicated by a square (Rudie et al., 2012) and a triangle (Prochaska et al., 2013a). We note that the covering fractions from the LBGs sample match well with the expectations for $M_{\mathrm{vir}}\approx 10^{12}M_{\odot}$ FIREbox haloes at redshift $z=2-2.5$. However, the observed covering fractions from QSOs are neither reproduced by FIREbox nor MassiveFIRE (FIRE-2) simulations.

4.2 Differential covering fraction

Figure 8: Simulated and observed differential H i covering fraction for redshift $z\sim 2$. Solid lines: Median differential covering fraction as a function of physical projected impact parameter away from the center of $M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\sun}$ haloes in FIREbox (1 halo at $z=2$). Dashed line: Median differential covering fraction as a function of physical projected impact parameter away from the center of the four most massive haloes in the MassiveFIRE (FIRE-2) zoom-ins (4 haloes at $z=2$). In both cases, the shaded areas display the 5th-95th percentile error on the median obtained from bootstrapping. Note that there are no haloes of this mass at redshift $z=2.5$ in either the FIREbox or MassiveFIRE (FIRE-2) simulations. Data points: the diamond data points show the observations of Prochaska et al. (2013b) for a sample of quasars at $z\sim 2-2.5$. The square data points show the simulations of Fumagalli et al. (2014) for five $M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\sun}$ galaxies at $z=2$. In both, the scatter in the abscissa indicates the lowest and highest impact parameter of each bin for which the value of $f_{\mathrm{cov}}(r)$ is calculated. We find good agreement with observations close to the halo centers ($\leq 50$ kpc) but systematically underestimate the covering fraction at large impact parameters ($\geq 100$ kpc).

We continue our comparisons with previous works, turning now to the radial profiles of covering fraction of LLSs. Rahmati et al. (2015) were able to reconcile observations around quasars by comparing results in physical units rather than fractions of the virial radius. We follow this convention in this section and discuss its implications in the text.

In Figure  8, we present our results for the differential covering fraction of very-massive ($M_{\mathrm{vir}}\gtrsim 10^{12.45}M_{\odot}$) haloes in both FIREbox and MassiveFIRE (FIRE-2), with data points for previous simulations and observations. The simulations by Fumagalli et al. (2014) reported very low covering fractions of LLSs which did not match the observations at any impact parameter. Rahmati et al. (2015) (not shown in Fig.  8) report excellent agreement between observed and simulated covering fractions of LLSs with the EAGLE simulations, at all impact parameters. Both FIREbox and MassiveFIRE (FIRE-2) results agree with observations of radial profiles of H i in the vicinity of very-massive haloes, but systematically underestimate the covering fraction further away from their center. In particular, we predict that the covering fraction of LLSs tends to 0 as $r$ increases, effectively becoming null for gas extending beyond $\sim 800$ kpc from the center, whereas the quasar sample from Prochaska et al. (2013a) suggests that it stagnates at a non-zero value of $f_{\mathrm{cov}}(r)\approx 0.2$.

FIREbox underestimates radial profiles of the observed covering fraction of LLSs, particularly in the outer regions of haloes. This particular outcome is not an exception: it was noted in most numerical works which could not reproduce the QSOs values (e.g., Meiksin et al., 2015, 2017; Gutcke et al., 2017; Suresh et al., 2019), and in related observing campaigns using the QPQ data (Rubin et al., 2015). This result can be expected for FIREbox, which slightly underestimates the overall covering fraction of LLSs around very-massive haloes, seemingly due to a smaller amount of H i found at large radii away from the center of haloes.

Figure 9: H i cumulative covering fraction for FIREbox main haloes within impact parameter $R_{\mathrm{vir}}$ as a function of halo mass $M_{\mathrm{vir}}$ for diverse density cuts. From top-left to bottom-right, we show $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ for H i density cuts of $N_{\text{H\,{i}}}>15.5$ cm^-2 (sLLSs), $N_{\text{H\,{i}}}>17.2$ cm^-2 (LLSs; shown here for reference), $N_{\text{H\,{i}}}>19.0$ cm^-2 (sDLAs) and $N_{\text{H\,{i}}}>20.3$ cm^-2 (DLAs). Different colours are used to distinguish between redshift from $z=0$ to $z=6$ as per the legend. Solid lines and shaded areas respectively show the median cumulative covering fraction and the 5th-95th percentile error on the median obtained from bootstrapping. All haloes with $M_{\mathrm{vir}}\geq 10^{7.75}M_{\sun}$ identified in the simulation were used to create these plots. Square data points correspond to measurements around $z\sim 2,2.5$ LBGs at the corresponding density cuts from Rudie et al. (2012) and hexagonal data points indicate the median and 16th-84th percentile of the covering fraction for a sample of DLAs from the EAGLE simulations reported in Garratt-Smithson et al. (2021). We find good agreement with the covering fractions from the LBGs sample for LLSs, sDLAs and DLAs. For the sLLS threshold ($N_{\text{H\,{i}}}>15.5$ cm^-2), we find that FIREbox underestimates the cumulative covering fraction by $\sim$30%. FIREbox shows good agreement with the EAGLE simulations at low redshifts and extends the relation down to low masses and to higher redshifts.

In the EAGLE simulations, Rahmati et al. (2015) find agreement with observations of differential covering fraction by Prochaska et al. (2013b) for $M_{\mathrm{vir}}\gtrsim 10^{12.5}M_{\odot}$ haloes. We stress however that Rahmati et al. (2015) cannot reproduce the observed cumulative covering fraction $f_{\mathrm{cov}}(<R_{\mathrm{vir}})=0.64^{+0.06}_{-0.07}$, underestimating it by a factor $>2$. The agreement here comes from considering the observational biases present in the QPQ sample (for details, see discussions in Rahmati et al. (2015) and Faucher-Giguère et al. (2016)). Essentially, they argue that Prochaska et al. (2013b) overestimate the covering fraction because of using a fixed virial radius typical of $M_{\mathrm{vir}}\approx 10^{12.5}M_{\sun}$ haloes, whereas they probe quasars of higher mass than predicted, closer to $M_{\mathrm{vir}}\gtrsim 10^{13}M_{\sun}$. They further argue that most of the sight lines at high impact parameters actually come from objects with redshift $z>2$. Both effects essentially lead to an overestimation of the cumulative covering fraction in the QPQ sample, particularly at high impact parameters away from the center, and Rahmati et al. (2015) find agreement with observations by correcting for these effects. Given the absence of haloes exceeding $M_{\mathrm{vir}}\gtrsim 10^{13}M_{\sun}$ at redshifts $z=2-2.5$ in our FIREbox and MassiveFIRE (FIRE-2) simulations, we were unable to employ this method to our results, and hence cannot conclude on its effectiveness in recovering the observed radial profiles of covering fraction.

Our discussion underscores the lack of consensus and challenges in comparing observations and simulations of H i covering fraction, as well as the large range of predictions produced by different models.

4.3 Cumulative covering fraction for additional $N_{\text{H\,{i}}}$ cuts (sLLS, LLS, sDLA, DLA)

To offer a more complete overview of strong H i absorbers in FIRE, we extend our analysis of the covering fraction of Lyman Limit Systems by providing the H i covering fraction of haloes in FIREbox at additional density cuts. The methods for computing these covering fractions are identical to that introduced in section  2. Figure  9 shows the covering fraction of atomic hydrogen $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ of main haloes in FIREbox as a function of halo virial mass, for density cuts corresponding to the class of absorbers known as sub-Lyman Limit Systems (sLLSs; $N_{\text{H\,{i}}}>15.5$ cm^-2), sub-Damped Lyman-$\alpha$ Absorbers (sDLAs; $N_{\text{H\,{i}}}>19.0$ cm^-2) and Damped Lyman-$\alpha$ Absorbers (DLAs; $N_{\text{H\,{i}}}>20.3$ cm^-2). We also compare with measurements by Rudie et al. (2012) of the H i covering fraction at these density cuts in each panel and show results for the cumulative covering fraction of DLAs in the EAGLE simulations from Garratt-Smithson et al. (2021).

The covering fractions measured by Rudie et al. (2012) for sDLAs and DLAs are in good agreement with results from FIREbox, as was the case for LLSs (Fig.  7). For the lower column density threshold of $N_{\text{H\,{i}}}>15.5$ cm^-2, the simulations slightly underestimate the cumulative covering fractions of $M_{\mathrm{vir}}=10^{12}M_{\sun}$ LBGs at redshifts $z\sim 2-2.5$. In the zoomed view of the bottom-left panel of Fig.  9, we show that the covering fractions of DLAs in FIREbox are similar to those of EAGLE massive haloes reported by Garratt-Smithson et al. (2021). Our analysis robustly extends the redshift and mass dependence for these systems, offering greater possibilities to compare results from future simulations and observations.

5 SUMMARY AND CONCLUSIONS

In this work, we have used simulations run with the FIRE-2 physics prescription (Hopkins et al., 2018) at high numerical resolution ($m_{\mathrm{b}}=3.6-6.24\times 10^{4}M_{\sun}$) to comprehensively investigate the atomic hydrogen covering fraction of Lyman Limit Systems in haloes spanning the mass range from $M_{\mathrm{vir}}\sim 10^{12.5}M_{\sun}$ down to very-low mass haloes with $M_{\mathrm{vir}}\sim 10^{8}M_{\sun}$, across redshifts $z=0-6$. Our analysis includes haloes from the (22.1 cMpc)³ cosmological volume simulation FIREbox (Feldmann et al., 2023), which currently constitutes the highest-resolution cosmological volume simulation with the largest dynamical range of its kind, supplemented with zoom-ins of four massive haloes from the MassiveFIRE (FIRE-2) sample (Feldmann et al., 2016, 2017; Anglés-Alcázar et al., 2017). This comprehensive approach allows for a more rigorous comparison with both observations and prior numerical investigations on the subject.

Our main results can be summarized as follows:

• The H i cumulative covering fraction of LLSs in FIREbox exhibits a pronounced dependence on redshift, showing significant increase at all halo masses from redshift $z=0$ to $z=6$ (Fig.  3). The complex halo mass dependence of $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ can be divided into two regimes. Notably, the high resolution of our study enables an investigation of the dependence for low-mass haloes, revealing that the cumulative covering fraction steeply increases from zero to some maximal value, which increases with increasing redshift, from $M_{\mathrm{vir}}=10^{7.75}M_{\odot}$ to a threshold mass which decreases with increasing redshift. For instance, at $z=0$, this threshold resides at $\ M_{\mathrm{vir}}\sim 10^{10.5}M_{\odot}$ while at $z=6$, it is approximately $\ M_{\mathrm{vir}}\sim 10^{9.2}M_{\odot}$. Beyond this threshold, the covering fraction plateaus and remains nearly independent of mass for higher-mass haloes at all redshifts.

• The H i differential covering fraction of LLSs in FIREbox is also highly dependent on redshift, showing a similar increase at all halo masses from redshift $z=0$ to $z=6$, see Fig.  4. The radial profiles $f_{\mathrm{cov}}(r)$ as a function of projected impact parameter $r$ from the center are found to resemble an inverse-sigmoid. It takes the value of $f_{\mathrm{cov}}(r)=1$ for inner impact parameters close to the center, before steeply decreasing toward 0 further away from the center.

• The turning points in the radial profiles happen at lower impact parameters with decreasing redshift, indicating that a greater fraction of H i is found closer to the center of haloes with decreasing redshift (Fig.  5). In particular, we note that almost all of the H i is found within the virial radius of haloes at lower redshifts ($z\lesssim 2$), in agreement with Villaescusa-Navarro et al. (2018); Feldmann et al. (2023).

• The radial profiles of strong H i absorbers in massive ($M_{\mathrm{vir}}\simeq 10^{11}-10^{12}M_{\odot}$) haloes in FIREbox are very similar to each other and show scale-invariance, see Fig.  5.

• Lower mass ($M_{\mathrm{vir}}\leq 10^{10}M_{\odot}$) haloes have differential covering fraction profiles with a similar shape but are significantly less extended than more massive haloes. The profiles get more extended with increasing mass and redshift until they overlap at all halo masses. This explains the evolution of the cumulative covering fraction with mass and redshift (Fig.  3 and Fig.  5).

• We presented a fitting function which accurately captures the radial profiles of our simulations (see Eq. (3); Fig.  6). The free parameter $L_{z}$ in our fitting function can be thought of as the typical projected radial extent of the H i halo with respect to the center of haloes. We found that the H i radial profiles of LLSs in FIREbox are much less extended than those studied in the OWLS and EAGLE simulations (Rahmati & Schaye, 2014; Rahmati et al., 2015).

• Our FIRE-2 simulations predict cumulative covering fractions for $M_{\mathrm{vir}}=10^{12}M_{\odot}$ haloes that are in agreement with those observed in $M_{\mathrm{vir}}\approx 10^{12}M_{\odot}$ LBGs at redshifts $z=2-2.5$ by Rudie et al. (2012), see Fig.  7.

• However, our simulations do not reproduce the observations of Prochaska et al. (2013a) for QSOs hosted in $M_{\mathrm{vir}}\approx 10^{12.5}M_{\odot}$ haloes at redshifts $z=2-2.5$. In particular, the average covering fraction of the observed sample is almost twice as large as the typical value found in our simulations (Fig.  7).

• Comparing the radial plots of H i covering fraction from FIREbox with the observations indicates that the simulations agree with the distribution inside haloes, but underestimate the extent of H i in the outer regions (Fig.  8). This seems to point at missing H i at large radii in the simulations. In particular, Prochaska et al. (2013b) find that the H i radial covering fraction stagnates at around 20% at all radii, whereas FIRE-2 simulations predict a decay to zero outside the virial radius of $M_{\mathrm{vir}}\sim 10^{12.5}M_{\sun}$ haloes.

• We compute the H i covering fraction of strong absorbers at additional density cuts of $N_{\text{H\,{i}}}>15.5$ cm^-2, $N_{\text{H\,{i}}}>19.0$ cm^-2 and $N_{\text{H\,{i}}}>20.3$ cm^-2, and show that they are in good agreement with observational measurements (Rudie et al., 2012) and the EAGLE simulations (Garratt-Smithson et al., 2021) (see Fig.  9). FIREbox’s high dynamical range allows us to extend the study of such systems to much lower halo masses and in a broader range of redshifts, providing avenues for future comparison with new simulations and observational campaigns.

Further work investigating the covering fraction in haloes hosting massive quasars needs to be conducted. This can be realized by examining much larger datasets of very-massive haloes in simulations and exploring the role of different feedback origins. The expanded samples of very-massive FIRE haloes, with and without AGN feedback, introduced in Wellons et al. (2023); Byrne et al. (2023), present a substantial opportunity to both extend the study of the H i covering fraction in FIRE to more massive haloes and to assess the effective contribution of super-massive black holes in shaping the distribution of cool gas in the CGM. Likewise, including feedback from AGNs in future iterations of FIREbox and MassiveFIRE simulations will be an essential complement to the findings presented in this work, albeit at the cost of increased uncertainty in the form of modelling degeneracies.

Measurements of the covering fraction of H i in lower mass haloes are also needed to more comprehensively compare with the dependence for $M_{\mathrm{vir}}\lesssim 10^{11}M_{\sun}$ haloes presented in our study. Although this remains a challenge, future campaigns promise advances in this regard thanks to considerable strides in instrumental capabilities and analytical methodologies over recent years. For instance, instruments such as MUSE (Multi Unit Spectroscopic Explorer; Bacon et al., 2010) and KCWI (Keck Cosmic Web Imager; Morrissey et al., 2018) are poised to explore absorption lines of gas in the circumgalactic medium and extend the analysis to lower galaxy masses than previously achievable (Dutta et al., 2020). Additionally, the planned instruments MOSAIC (Multi-Object Spectrograph; Evans et al., 2015) and ANDES (ArmazoNes high Dispersion Echelle Spectrograph; Marconi et al., 2022) on the Extremely Large Telescope (Gilmozzi & Spyromilio, 2007; Neichel et al., 2018) are anticipated to offer the highest precision attainable across a great range of objects and redshifts in the coming decades. These advancements open exciting new avenues for comparisons with our research and other studies in the field.

Data Availability

The data underlying this article are available on reasonable request to the corresponding author. A public version of the GIZMO code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.

Appendix A Resolution effects on H i covering fraction of Lyman Limit Systems

In this section, we discuss the numerical convergence of the H i covering fraction of LLSs FIREbox. It is well-known today that cold gas in the CGM is not a converged quantity (Faucher-Giguère & Oh, 2023; Ramesh & Nelson, 2023). Recent work shows that increased resolution in the CGM significantly boosts the simulated neutral hydrogen column density and resulting covering fractions (e.g., Faucher-Giguère et al., 2016; van de Voort et al., 2019).

We use the lower resolution runs FB512 and FB256 from the FIREbox suite of simulations to investigate resolution effects in our study. The dark matter and baryon masses are 8 and 64 times lower, respectively, and the softening lengths are adjusted accordingly (for details, see Feldmann et al., 2023, section 2). We show the cumulative covering fraction within the virial radius $f_{\mathrm{cov}}(<R_{\mathrm{vir}})$ as a function of virial mass $M_{\mathrm{vir}}$ for the different runs in Fig. 10. We find that the increased resolution of the fiducial N1024 run produces slightly higher H i covering fractions, at all redshifts and for all masses. This increase is not significant, consistent with findings from Gensior et al. (2023), wherein it was found that other H i properties in FIREbox are converged (see Fig. A1 in their online supplementary material). The neutral hydrogen maps used throughout this work (N1024) thus appear effectively converged with numerical resolution.

Figure 10: Cumulative covering fraction of LLSs within $R_{\mathrm{vir}}$ as a function of halo mass $M_{\mathrm{vir}}$ for different numerical resolutions. Only main haloes were used for this analysis. Different colours indicate different redshifts from $z=0$ to $z=6$ as per the legend. Solid, dashed and dotted lines show results for the N1024 run (fiducial FIREbox run used in this work), the N512 run (8 times lower resolution than N1024) and the N256 run (64 times lower resolution than N1024) respectively. The dashed (dotted) black vertical line indicates the N1024-equivalent lowest halo mass resolved for the N512 (N256) run, corresponding to haloes with at least 168 dark matter particles. We find that the covering fraction slightly increases with resolution, at all halo masses and redshifts. However, this effect is small which suggests that the covering fractions studied throughout this work are close to being converged with resolution.

Appendix B H i column density distribution

We show the H i column density distribution for FIREbox haloes of different masses and for redshifts $z=2,3,4$ and 6 in Fig.  11. We indicate the neutral hydrogen column density threshold of Lyman Limit Systems ($N_{\mathrm{HI}}=10^{17.2}$ cm^-2) by a black dotted line.

We find that the profiles of neutral hydrogen column density are very similar in the outer edges ($\gtrsim R_{\mathrm{vir}}$) of FIREbox haloes, for all halo masses and at all redshifts. In the inner regions of haloes, however, we see that the H i column density strongly increases with increasing halo mass until it reaches a maximum threshold. For massive haloes with $M_{\mathrm{vir}}\geq 10^{11}M_{\sun}$, the density threshold of Lyman Limit Systems is reached at the same normalised radius $r/\,R_{\mathrm{vir}}$, further hinting at the scale-invariance discussed in sections  3.2 and  3.3. While in lower mass haloes ($M_{\mathrm{vir}}\lesssim 10^{10}M_{\sun}$), the density threshold is reached at lower radii or not reached at all. The LLS differential covering fraction profiles of lower mass haloes are therefore less extended than for more massive haloes.

Figure 11: H i column density profiles of FIREbox haloes as a function of normalised impact parameter for redshifts $z\geq 2$. Solid lines indicate the median $N_{\text{H\,{i}}}(r/\,R_{\mathrm{vir}})$ at each impact parameter and the shaded areas show the 5th-95th percentile error on the median obtained from bootstrapping. The mass bins are described in section  3.2. The redshift depicted in each panel is shown at the bottom left. We indicate the LLS column density cut ($N_{\mathrm{HI}}=10^{17.2}$ cm^-2) with a black dashed line. The column density profiles of H i are similar in the outskirts ($\gtrsim R_{\mathrm{vir}}$) for all halo masses. In the inner regions, $N_{\text{H\,{i}}}$ increases with increasing halo mass. LLSs are thus less extended in lower mass haloes because the density threshold is reached at lower radii than for more massive haloes.