KODIAQ-Z: Metals and Baryons in the Cool Intergalactic and Circumgalactic Gas at 2.2

We present the KODIAQ-Z survey aimed to characterize the cool, photoionized gas at 2.2<z<3.6 in 202 HI-selected absorbers with 14.6<log N(HI)<20, i.e., the gaseous interface between galaxies and the intergalactic medium (IGM). We find that the 14.6<log N(HI)<20 gas at 2.2<z<3.6 can be metal-rich gas (-1.6<[X/H]<-0.2) as seen in damped Ly-alpha absorbers (DLAs); it can also be very metal-poor ([X/H]<-2.4) or even pristine gas ([X/H]<-3.8) not observed in DLAs, but commonly observed in the IGM. For 16<log N(HI)<20 absorbers, the frequency of pristine absorbers is about 1%-10%, while for 14.6<log N(HI)<16 absorbers it is 10%-20%, similar to the diffuse IGM. Supersolar gas is extremely rare (<1%) in this gas. The factor of several thousand spread from the lowest to highest metallicities and large metallicity variations (a factor of a few to>100) between absorbers separated by less than 500 km/s imply that the metals are poorly mixed in 14.6<log N(HI)<20 gas. We show that these photoionized absorbers contribute to about 10% of the cosmic baryons and 30% of the cosmic metals at 2.2<z<3.6. We find the mean metallicity increases with N(HI), consistent with what is found in z<1 gas. The metallicity of gas in this column density regime has increased by a factor ~8 from 2.2<z<3.6 to z<1, but the contribution of the 14.6<log N(HI)<19 absorbers to the total metal budget of the universe at z<1 is half that at 2.2<z<3.6, indicating a substantial shift in the reservoirs of metals between these two epochs. We compare the KODIAQ-Z results to FOGGIE cosmological zoom simulations. The simulations show an evolution of [X/H] with N(HI) similar to our observational results. Very metal-poor absorbers with [X/H]<-2.4 at z~2-3 in these simulations are excellent tracers of inflows, while higher metallicity absorbers are a mixture of inflows and outflows.


Abstract
We present the KODIAQ-Z survey aimed to characterize the cool, photoionized gas at 2.2  z  3.6 in 202 H Iselected absorbers with 14.6 N log H I < 20 that probe the interface between galaxies and the intergalactic medium (IGM). We find that gas with N 14.6 log 20 H I <  at 2.2  z  3.6 can be metal-rich (−1.6  [X/H]  − 0.2) as seen in damped Lyα absorbers (DLAs); it can also be very metal-poor ([X/H] < − 2.4) or even pristine ([X/H] < − 3.8), which is not observed in DLAs but is common in the IGM. For N 16 log 20 H I < < absorbers, the frequency of pristine absorbers is about 1%-10%, while for N 14.6 log 16 absorbers it is 10%-20%, similar to the diffuse IGM. Supersolar gas is extremely rare (<1%) at these redshifts. The factor of several thousand spread from the lowest to highest metallicities and large metallicity variations (a factor of a few to >100) between absorbers separated by less than Δv < 500 km s −1 imply that the metals are poorly mixed in N 14.6 log 20 H I <  gas. We show that these photoionized absorbers contribute to about 14% of the cosmic baryons and 45% of the cosmic metals at 2.2  z  3.6. We find that the mean metallicity increases with N H I , consistent with what is found in z < 1 gas. The metallicity of gas in this column density regime has increased by a factor ∼8 from 2.2  z  3.6 to z < 1, but the contribution of the N 14.6 log 19 H I <  absorbers to the total metal budget of the universe at z < 1 is a quarter of that at 2.2  z  3.6. We show that FOGGIE cosmological zoom-in simulations have a similar evolution of [X/H] with N H I , which is not observed in lower-resolution simulations. In these simulations, very metal-poor absorbers with [X/H] < − 2.4 at z ∼ 2-3 are tracers of inflows, while highermetallicity absorbers are a mixture of inflows and outflows.

Introduction
The intergalactic medium (IGM) and the circumgalactic medium (CGM) are massive reservoirs of baryons (e.g., Fukugita et al. 1998;McQuinn 2016;Macquart et al. 2020) and major fuel sources for star formation in galaxies (e.g., Tumlinson et al. 2017). They play a vital role in the formation and evolution of galactic and large-scale structures in the universe. The empirical and theoretical characterizations of the CGM and IGM before, during, and after the peak of cosmic star formation (z ∼ 2) are therefore of prime importance to understanding their role in the evolution of galaxies and large-scale structures. Gathering data that can probe both the low-and high-redshift universe is critical for robustly constraining the evolution of these cosmic structures. This is also the best route to fully test always-improving cosmological simulations. Modern numerical simulations are now reaching new milestones in fidelity. Notably, a new generation of simulations boasts unprecedentedly high spatial resolution, even in the more diffuse regions of the simulation boxes (e.g., Hummels et al. 2019;Mandelker et al. 2019;Peeples et al. 2019;Suresh et al. 2019;). This allows for far more realistic modeling of the CGM and IGM and their interrelationships with galaxies.
Over the past few years, our group has engaged in several surveys at both low and high redshifts to target absorbers in the H I column density range N 15 log 19 H I   . This column density interval spans a range of physical environments from the diffuse IGM (Lyα forest, hereafter LYAF) to the denser portions of the CGM and the edges of galaxy disks. Following Lehner et al. (2018) ). With overdensities between the diffuse IGM (probed by the LYAF) and galaxies (probed by DLAs), the SLFSs, pLLSs, and LLSs spanning N 15 log 19 H I   should be at the heart of the exchange of matter between galaxies, their CGM, and the diffuse IGM.
One of the main goals of our H I absorption surveys is to provide a census of the chemical enrichment of the absorbers with N 14.6 log 20 H I   over cosmic time (e.g., the COS CGM Compendium (CCC) at low redshift, Lehner et al. 2018Lehner et al. , 2019Wotta et al. 2019; see also Lehner et al. 2013;Wotta et al. 2016; the HD-LLS survey, Prochaska et al. 2015;Fumagalli et al. 2016b; the KODIAQ and KODIAQ-Z surveys, Lehner et al. 2014Lehner et al. , 2016; this paper). The metallicity of the absorbers is a direct measure of their metal enrichment and a key diagnostic of their origins. The metallicity provides direct information on how efficient (or not) galaxies are at enriching their immediate surroundings and beyond, as well as on the level of metal mixing in the diffuse regions of the universe. The selection of our surveys provides a relatively unbiased way to obtain a census of the metallicity, since selecting on hydrogen is sensitive to both very low and very high metallicities.
While pLLSs and LLSs have been found often in the CGM of galaxies at both low and high redshift (e.g., Lehner et al. 2009Lehner et al. , 2013Ribaudo et al. 2011;Kacprzak et al. 2012;Prochaska et al. 2017b;Cooper et al. 2021;Berg et al. 2022), there is growing evidence that at least some of the pLLSs and LLSs with very low metallicities may probe the denser IGM (very low metallicity being [X/H]  −1.4 at z  1, Berg et al. 2022; and [X/H]  −3 at 2  z  3.5, Fumagalli et al. 2011aFumagalli et al. , 2016aCrighton et al. 2016;Lofthouse et al. 2020). 13 Recent high-resolution zoom-in cosmological simulations also show that metal-poor ([X/H] < − 3) LLS-like absorbers can be found in the IGM far from any galaxies (Mandelker et al. , 2021. Determining the metallicity distributions of the SLFSs, pLLSs, and LLSs at low and high z directly informs us on the metallicity enrichment (or lack thereof) relative to galaxies and their inner regions probed by DLAs and SLLSs. This technique can also assess the amount of pristine gas in these overdense regions of the universe and how it compares to the more diffuse IGM probed by LYAF absorbers (Fumagalli et al. 2011a;Mandelker et al. 2019).
In previous papers, our group presented a sample of 261 z < 1 absorbers with N 15 log 19 H I   (Lehner et al. 2018Wotta et al. 2019, hereafter CCCI, CCCII, CCCIII, respectively). This sample includes an unexpectedly large fraction of metal-poor absorbers, many of which have metallicities [X/H] < − 2, implying little chemical enrichment since z ∼ 2-3. The finding that metal-poor absorbers with [X/H] − 1.4 represents about half of the population of z  1 absorbers with N 15 log 19 H I   was unanticipated. This is because neither the LYAF (due, in retrospect, to a lack of sensitivity) nor the SLLSs/DLAs hinted at the presence and importance of such low-metallicity gas at low redshift. Although not fully appreciated at the time, there was, however, already some evidence of very low metallicity gas at low redshift based on individual studies of LLSs (Cooksey et al. 2008;Ribaudo et al. 2011; see also Tripp et al. 2005 for a rare example of a primitive SLLS at z ∼ 0.06). CCC also reveals that, at z < 1, there is no evidence of pristine gas, implying that >99% of the gas with N 15 log 19 H I   has been polluted by z  1, even if it is only at a very low level. Low-redshift photoionized gas with N 15 log 19 has a wide range of metallicities observed (−3  [X/H]  + 0.4) and can show large metallicity variations (factor up to 25) over small redshift/velocity separation (Δv < 300 km s −1 ), implying that these absorbers sample a highly inhomogeneous medium probed by these absorbers. The CCC results raise the following questions: Are these properties unique to z  1? How do they evolve from earlier times, across the peak of cosmic star formation at z ∼ 2.5?
To address these questions, we present here the Keck Database of Ionized Absorbers toward Quasars (KODIAQ-Z) survey, a survey designed to identify H I absorbers with N log 14.5 H I  at z > 2 in Keck HIRES spectra available in the KODIAQ database (O'Meara et al. , 2017 and determine their metallicities. Due to the higher sensitivity of the Keck HIRES data, compared to the Hubble Space Telescope (HST) Cosmic Origins Spectrograph (COS) data, we can probe low-N H I absorbers at high redshift and still be sensitive to lowmetallicity gas. In fact, due to the new sensitivity gained with 8-10 m ground-based telescopes, early studies on the metal enrichment at high redshift largely focused mostly on IGM absorbers with N log 14.5 H I  (e.g., Songaila 1998;Ellison et al. 2000;Schaye et al. 2003;Aguirre et al. 2004;Simcoe et al. 2004). Our sample of H I-selected absorbers consists of 155 SLFSs, 24 pLLSs, 16 LLSs, and 7 SLLSs at 2.2  z  3.6, totaling 202 absorbers, 195 of them with N 14.5 log 19 H I <  . On account of the nature of the Keck HIRES data, the sample of pLLSs and LLSs has only marginally increased by a factor 1.3 compared to our earlier pilot survey . However, in this paper we undertake a systematic search and characterization of H I-selected SLFSs. 14 In our analysis, we include the H I-selected absorbers from the HD-LLS survey (Prochaska et al. 2015;Fumagalli et al. 2016b), which adds another 46 LLSs, increasing the total sample to 241 absorbers with N 14.5 log 19 H I   . The HD-LLS survey is directly complementary to KODIAQ-Z because most of the HD-LLS absorbers have N log 17.5 H I  (including 73 SLLSs), whereas KODIAQ-Z largely samples gas with lower N H I .
The main goals of the present paper are to determine the metallicity distributions of these absorbers, assess how their metallicity changes with N H I and z, and provide a robust limit on the amount of pristine gas at high redshift. We also present and discuss the physical properties of these absorbers and use these results to estimate the metal and baryon budgets of the cool photoionized CGM at z ; 2.2-3.6. Our paper is organized as follows. In Section 2, we describe the survey design of 13 For the metallicity, we use conventional square bracket notation , where X is an α-element, unless otherwise stated.
KODIAQ-Z, the search for the H I-selected absorbers, and the sample of SLFSs, pLLSs, and LLSs. In Section 3, we present the methods to derive the column densities of the metals and H I for each absorber, while in Section 4 we present basic statistical properties of our sample (z, N H I , and the Doppler parameter b). In Section 5, we present ionization modeling of the KODIAQ-Z absorbers and assumptions for estimating their metallicity. Our main results on the metallicities of the H Iselected absorbers are presented in Section 6. In Section 7, we derive the physical properties of the absorbers, including their density (n H ), total column density (N H ), temperature, and linear scale (l ≡ N H /n H ), and we model the neutral fraction as a function of N H I in the range N 14.6 log 20 H I   . In Section 8, we estimate the cosmic metal and baryon budgets of these absorbers. In Section 9, we compare the KODIAQ-Z results to the properties of the gas in cosmological zoom-in simulations from the Figuring Out Gas & Galaxies in Enzo (FOGGIE) project (Peeples et al. 2019;Corlies et al. 2020) and use these simulations to gain insight into the connection of the N 14.5 log 20 H I   gas depending on its metallicity with galaxies. In Section 10, we discuss the results, and in Section 11, we summarize our main conclusions.

Sample Selection Criteria
Our survey is based on a search for H I absorption in the KODIAQ DR2 database (O'Meara et al. 2017) of 300 QSO spectra observed with HIRES. The QSOs useful for our survey lie at 2  z em  4.5, limiting the sample size to 235 QSOs. We also searched in the spectra of higher-redshift (z em > 4.5) QSOs, but the LYAF is too dense to allow us to reliably determine the H I and metal-line column densities at these redshifts. The lower redshift threshold z  2 is required to have wavelength coverage of at least Lyα and Lyβ. In each of the QSO spectra, we search for absorbers with H I column density N log 14.5 H I  . This contrasts with our previous KODIAQ-Z survey , where our search for H I absorbers was limited to absorbers with N log 16.2 H I  (to match our initial lower z survey, Lehner et al. 2013) and was by no means a systematic and complete examination of the KODIAQ DR1 database .
A major difference (besides the redshift probed) with the low-redshift survey CCC (CCCI) is that all the HIRES spectra are flux normalized for the co-addition of the individual exposures and echelle orders (O'Meara et al. , 2017. This makes determination of large-scale flux decrements (e.g., at the Lyman limit break or damping wings of DLAs) unreliable. Therefore, to estimate the H I column density, we can rely only on the absorption observed in the Lyman series transitions. This limits especially the sample size of absorbers with N 17.2 log 18 H I   because it requires having enough weak Lyman transitions that are not too contaminated to robustly derive N H I (see Section 3.2). For absorbers with N log 18 H I  , damping wings start to be observable in the Lyα absorption, and therefore this transition can be used to derive N H I depending on the level of contamination of the absorption from the LYAF. On the other hand, while CCC is limited to absorbers with N log 15.3 H I  , due to the need to be sensitive to metallicities [X/H]  − 1, the typically higher signal-tonoise ratios (S/Ns) and higher spectral resolution of the Keck HIRES spectra allow us to probe lower N H I , down to the LYAF threshold of about N log 14.5 H I  , and still be sensitive to metallicity with [X/H]  − 2. Although with very high S/N Keck HIRES spectra a strong limit ([X/H]  − 2) on the metallicities can be placed down to N log 13.6 H I  (e.g., Simcoe et al. 2004), our lower limit on N H I was selected to correspond to the division between the diffuse IGM and denser IGM/diffuse CGM. Following Schaye (2001b), a 14.5 dex limit on the H I column density corresponds to an overdensity of δ  5 at z ∼ 2.8 (about the average redshift of the absorbers in our sample; see Section 4.1).

Search for H I Absorbers in the KODIAQ Database
To identify the strong H I absorption, we developed an automated search tool to help identify the Lyman series transitions from 1215 Å (Lyα) down to 917 Å. In Figure 1, the colored profiles show single-component Voigt models for absorbers with a Doppler parameter b = 25 km s −1 and N log 15, 16, 17, 18 H I = .
That b-value is selected based on previous profile fitting (PF) of strong H I absorbers at similar redshifts (e.g., Lehner et al. 2016); as we will see below, that choice is close a posteriori to the mean b-value in our sample. We use these H I models to set the following criteria for our automatic search: (1) Lyα must have flux F λ 0 over 0.4 Å in contiguous pixels; (2) at the average velocity/redshift (v c ) of the identified Lyα, the flux in the Lyβ absorption must also reach F λ 0; and (3) at v c , we also require that there is substantial absorption in at least Lyγ (and Lyδ if there is wavelength coverage of that transition). Additional conditions were added to ensure that when a DLA or an SLLS is present the search program would not continue to look for additional absorbers near the redshift of the DLA or SLLS. For this search, we therefore require for each QSO wavelength coverage of at least Lyα, Lyβ, and Lyγ. However, since the damping wings of Lyα can be used to determine N H I if N log 18.5 H I  , we also undertook a separate search that relies solely on Lyα where the flux f λ ; 0 over a wavelength width of  2 Å. A posteriori, only a few of these very strong H I absorbers ended up in our sample, due to the fact that the KODIAQ spectra are normalized prior to co-addition, and the normalization was often not reliable enough to model the damping wings for these absorbers, especially the strong SLLSs and DLAs (DLAs not being an issue, since they are not part of our survey).
For each candidate absorber, our search tool created a stack of H I transitions as shown on Figure 1. Over 3000 candidates were created, many of these being false positives owing to the large contamination from the LYAF and other metal lines. That number is still small enough that a quick look at the stack of H I profiles could easily weed out most of the false positives, using the Voigt models as a guide. This vetting led the sample to be reduced by more than a third, to 949 candidate absorbers, the redshift distributions of which we show in Figure 2. In this figure, we also show the redshift distribution of the QSOs. About 83% of the candidate absorbers are at 2.2 z abs 3.6; candidate absorbers at z abs 2 are detected solely based on Lyα.
To assemble our final sample for metallicity determinations, we considered both the ability to estimate N H I and the coverage of metal ions (in order to be able to estimate the metallicity). For each candidate absorber, we therefore created another stack of velocity profiles with Lyα, Lyβ, Lyδ, and some key metal transitions used in the estimation of the metallicity (including Si II, Si III, Si IV, C II, C III, C IV, Al II, Al III, Fe II, O I, and O VI) to assess whether there is enough metal-ion coverage to determine the metallicity. After all the measurements are done (detailed in Section 3), the final sample size consists of 202 absorbers, i.e., ≈7% of the candidates from the initial automatic search were absorbers where we could place a robust estimate on N H I and metal-line column densities (including strong upper limits). In Table 1, we summarize the sample of QSOs with information on their redshifts, the wavelength coverage, S/N, and spectral resolution of the Keck HIRES observations. The final KODIAQ-Z sample includes 155 SLFSs, 24 pLLSs, 16 LLSs, and 7 SLLSs. We note that the redshift of each absorber was determined from the peak optical depth of the strongest H I component in our automatic search process. In most cases (89%), the velocities in the rest frame of the absorber redshift are |v| < 20 km s −1 , but in some cases they can be somewhat larger.

Comparison Samples
One of the goals of our survey is to determine how the metallicity changes with N H I at similar redshifts and between low and high z. For the low redshift, we use results from our CCC survey, including the compilation of H I absorbers at N log 19 H I  (see CCCI; CCCII; CCCIII and references therein). At z  2, we use two surveys: (1) the HD-LLS survey (Prochaska et al. 2015;Fumagalli et al. 2016b (which is directly complementary to KODIAQ-Z; HD-LLS also follows the same H I selection); and (2) the DLA survey by Rafelski et al. (2012), hereafter the These are used to help assess the level of contamination in each H I transition. Ultimately, for this candidate absorber, two absorbers were identified at z = 2.55125 and 2.55216 (corresponding to v ; −50 and +17 km s −1 in this figure). R12-DLA survey. Our use of these surveys was mostly dictated by the ease of access to the data while working in the early stage of KODIAQ-Z, including the metal and H I column densities for the HD-LLS. We note that the Magellan DLA survey by Jorgenson et al. (2013) has essentially the same metallicity distribution as R12-DLA in our redshift of interest 2.2  z  3.6 (see below). The recent XQ-100 survey of SLLSs by Berg et al. (2021) also shows the same metallicity distributions as in HD-LLS.
For the R12-DLA, we adopt the metallicities from the αelements (or, if not available, Zn). For the HD-LLS, we reestimate all the metallicities using the same extreme ultraviolet background (EUVB) adopted in CCC, in order to avoid mismatches in the comparisons of the metallicities-although, as we will show in Section 5.3, the systematic error arising from the EUVB in the photoionization modeling of the absorbers is much smaller (in fact, it is negligible) at z > 2 than at lower redshift, especially for the strong LLSs and SLLSs that are part of the HD-LLS survey.

Estimation of the Column Densities
Our analysis follows closely that undertaken in CCC to estimate the column densities of H I and metal ions. Throughout we assume the atomic parameters (central wavelengths and f-values) listed in Morton (2003). For each absorber, we first attempt to estimate the H I column density since that step is necessary to estimate the metallicity of the absorbers. For H I, we use a combination of different methods-apparent optical depth (AOD) method (Savage & Sembach 1991), curve-of-growth (COG) method, and Voigt PF method-while for the metal ions, we use solely the AOD method. We employ several methods to estimate the H I column densities, on account of the contamination, blending, and saturation effects being more severe than for the metal lines (see Figure 1).

Metal Column Densities
Although we only estimate the column densities of the metal species once the H I transitions have been modeled (and hence the candidate absorber has become a selected absorber), we first describe the metal lines, as we use only the AOD method to extract the properties of the absorption. In this method, the absorption profiles are converted into apparent optical depth per where F c (v) and F obs (v) are the modeled continuum and observed fluxes as a function of velocity, respectively. The AOD, τ a (v), is related to the apparent column density per unit velocity, N a (v), through the relation where f is the oscillator strength of the transition and λ is the wavelength in Å. The total column density is obtained by integrating the profile over the predefined velocity interval, ]are the boundaries of the absorption. We computed the average line centroids through the first moment of the AOD v a = ∫vτ a (v)dv/∫τ a (v)dv km s −1 . The velocity range over which the line was integrated was determined from weak H I transitions and metal line.
We emphasize that, as much as possible, the integration range v v , min max [ ] corresponds to a single absorbing H I complex, as shown, e.g., in Figure 3, but this does not necessarily mean that the metal lines consist solely of a single component (often they are, but not always). For the absorbers considered in this survey, the main H I components can typically be separated when the central velocities between two H I components are at least separated by |Δv|  40 km s −1 and N log 18 H I  . As we will see below, the mean Doppler parameter of the H I components is around 27 km s −1 , corresponding to a temperature of T < 4 × 10 4 K, typical for photoionized gas. This means that the thermal broadening is important and hence the H I lines are broader than the metal lines, i.e., there are H I transitions that can appear as a single component (at the observed spectral resolution) despite the metal ions showing multiple components over a velocity width of ∼50 km s −1 . In that case, we integrate the metal ions over a similar velocity width to the H I when that information is available from weak H I transitions; if the latter is not available (e.g., SLLSs where weak H I transitions are not covered, which is the case for 5% of the sample), we guide our analysis using the low ions or O I (if it is detected), which are better tracers of denser regions probed by the strongest H I absorbers. In the latter case, we integrate over the entire profiles of the low ions that may have a single or several components; this velocity integration range dictates the integration range of the higher ions, even though they may extend well beyond the absorption seen in the low ions (see, e.g., Lehner et al. 2014). As illustrated in Figures 1 and 3, even deep in the LYAF, the original continuum placement (see O'Meara et al. 2015O'Meara et al. , 2017 is satisfactory and can be used in many cases to estimate the physical parameters of the absorption lines accurately. We therefore took full advantage of the Keck HIRES spectra already being normalized to use the AOD method in a fully automated way. For each identified absorber, we typically made two iterations to estimate the physical parameters of the lines (v a , N a , and equivalent width W λ ). In the first iteration, the velocity range is estimated based on the stack of velocity profiles used to preview the data. Our program created figures that show for each transition the normalized and apparent column density profiles (such as Figure 3 in CCCI), allowing us to refine the integration velocity range and determine whether the continuum needed any local refinement. In the second iteration, the refined velocity range is applied to the original or renormalized spectra. If we determine that the continuum model needs to be reestimated, we either use the automated Legendre polynomial fitting described in CCCI or manually select the continuum region to be fitted by low-order Legendre polynomials.
We integrate the equivalent widths using the same v v , ] integration range adopted for the N a (v) profiles. The main use of the equivalent width for the metal lines in our survey is to determine whether the absorption is detected at the 2σ level. If it is not, we quote a 2σ upper limit on the column density, which is simply defined as twice the 1σ error derived for the column density assuming that the absorption line lies on the linear part of the curve of growth. The 1σ error is determined by integrating the spectrum over a velocity interval similar to that of a detected ion or very weak H I transition, or otherwise ±20 km s −1 , based on the typical smallest velocity intervals in other absorbers with weak detection of metals.
When an absorption feature is detected at 2σ significance, the next step is to assess whether there are any contamination or saturation issues. The main ions detected in this work are C II, C III, C IV, Si II, Si III, Si IV, and O VI (see Figure 3 and the Appendix). We also consider other ions, such as Fe II, Fe III, Al II, and Al III, or atoms such as O I, but in most cases these Figure 3. Example of normalized spectra against the rest-frame velocity of the absorber at z = 2.55125 toward J010311+131617 that show detection of H I, C III, C IV, and O VI (the integration ranges of the profiles are shown in red). Another absorber separated by 65 km s −1 is found at z = 2.55216; in this case, the column densities can be estimated reliably in each component/absorber (see also Figure 4).
are not detected ( Figure 3 and the Appendix). For C IV λλ1548, 1550, Si IV λλ1393, 1402, O VI λλ1031, 1037, and S II λλ1193, 1260, 1304, and 1526, these are either doublets or ions with multiple transitions, and both contamination and saturation can be readily checked by comparing their N a (v) profiles. Ions with λ  1215 Å are more often likely to be contaminated than ions with transitions at longer wavelengths, due to the dense LYAF at these redshifts. For example, for the absorber shown in Figure 3, the C IV transitions are not contaminated, but the weak transition of O VI is. In contrast, the weakest H I absorber at about −55 km s −1 seen in Figure 3 (which is not analyzed as part of this work because its H I column density is N log 14.5 H I < ) shows absorption in both transitions of each doublet of C IV and O VI. The C II ion has two transitions at 1334 and 1036 Å, and we always check the latter transition if λ1334 is detected; more often than not, it is, however, contaminated by the LYAF.
On the other hand, ions such as C III λ977 and Si III λ1206 are in the LYAF and have only a single transition. In this case, our strategy for assessing whether a given transition is contaminated by another unrelated absorber is based on a combination of visual inspection and comparison of the relevant quantities (e.g., mean velocity and apparent optical depth) with other transitions. Strong contamination can often be diagnosed visually, by comparing both the central velocities of the peak optical depth and the shape of the velocity structure of the absorption profiles across transitions from the same or similar ions. 15 The mean velocities of the absorption profiles estimated from the integration of the absorption profiles can often point to strong contamination if the mean velocity from one transition is significantly different from others that are aligned. In the case of obvious contamination, the absorption feature is discarded. The only exceptions are as follows. If the contamination is very mild, only occurring in the wing of a relatively strong absorption, we slightly change the velocity integration range to avoid the contaminated portion for that absorption feature. In cases where we are not sure that the absorption arises solely from the understudied metal ion (e.g., detection of C III but no coverage or detection of C IV or other ions), but the feature could very likely be produced by that ion based on the same velocity as that of H I and narrow absorption (as expected from a metal, low ion), we treat the column density from that ion as an upper limit in our ionization modeling. 16 For most of the metal ions in our survey, saturation is not a major issue, due to the lower column H I densities than those of DLAs and even SLLSs (there are only seven SLLSs in our sample and no DLA), the lower metallicities on average than at z < 1, and the high spectral resolution of the Keck data. To determine whether there is saturation, we proceed similarly to how we did for the contamination assessment. Where there is no obvious sign of contamination and the absorption from the metal ion reaches zero flux, the absorption is automatically marked as a lower limit. The comparison of the AOD profiles for doublet ions or ions with multiple transitions can readily help determine whether there is saturation if the profiles do not match at the peak optical depth and if the apparent column density of the stronger transitions is systematically smaller than that of the weaker transitions.
For atomic or ionic transitions that have f log 0.3  ( ) l D and the difference in apparent column densities of the weak and strong transitions is Wotta et al. 2016), we are able to correct for the mild saturation using the method described in Savage & Sembach (1991) and Wotta et al. (2016). In this case, we report the apparent column densities for each transition and then the adopted column density, corrected for saturation. In cases where we were not able to correct for saturation, we report the column density as a lower limit determined by the AOD method. 17 For multiple transitions where there are no saturation or contamination issues, we give velocities and column densities as the weighted average. All the column densities, velocity integration ranges, average velocities, and redshifts for the metal ions and H I (see Section 3.2) are summarized in Table 2. We also provide all the results in a machine-readable format (see the Appendix).

Determination of the H I Properties
To derive the metallicity, we can rely in all the cases on having column densities for several metal ions (detections and/ or nondetections). However, to have information on the amount of hydrogen atoms, we rely solely on H I. The advantage with H I is that we can access, in most cases, the entire Lyman series from 1215 Å all the way down to 913 Å (depending on wavelength coverage, redshift of the absorber, presence of an SLLS/DLA along the sight line, and strength of the absorption), spanning a difference in fλ of over 2700. However, the entire Lyman series is in the LYAF (and Lyβ forest), and therefore many transitions can be and often are heavily contaminated, necessitating a close look at all the H I absorption features to assess the level of contamination. Therefore, our first task is to systematically look simultaneously at several H I profiles (from 1215 to 915 Å) to assess whether there is enough information (i.e., portion of the profiles of H I transitions that are uncontaminated) to derive N H I . If it is judged that the contamination can be dealt with, we proceed by estimating the apparent column densities from the AOD profiles (and equivalent widths) and creating a mask of the H I profiles that are contaminated (i.e., flagging the pixels that are clearly contaminated), a mask that is then used in the PF (see below).
Due to the LYAF contamination, each transition is handled manually, and we often also reassess the continuum that can be more uncertain than in the redder part by modeling it with loworder Legendre polynomials. As mentioned above, we treated as much as possible an absorber as a single H I component. For the example shown in Figure 3, there are two components separated by 65 km s −1 , which correspond to two absorbers at z abs = 2.55125 and 2.55216. We therefore integrate the AOD profiles over each component (for the absorber at z abs = 2.55125, this corresponds to the velocity interval of about [−30, 40] km s −1 ). For absorbers with N log 18 H I  , this is not always possible, especially if there is no coverage of the weaker Lyman series transitions (for these absorbers, the damping wings of Lyα and Lyβ can be used to determine N H I ; 15 Of course, if the investigated ion is not in the same gas phase as the ion that is uncontaminated (e.g., O VI in some cases), these criteria cannot be used. 16 This upper limit can be differentiated from a nondetection upper limit, as the error bars of the former are not the standard nondetection upper limit errors +0.18, −0.30 dex. 17 Contamination and saturation can be a priori both present at the same time and deceitful. However, this is unlikely to occur frequently over the same pixels, and therefore each can be disentangled reliably in most cases. see below). However, these absorbers are rare in our sample, 5% (see below). For each transition, the AOD results are checked for contamination and/or saturation. For the transitions where there is no sign of contamination and/or saturation, we use a weighted average of the integrated apparent column densities, which provides the N H I estimates from the AOD method.
The equivalent widths, estimated along with the apparent column densities and average velocities, are used for the COG analysis, which employs a χ 2 minimization and χ 2 error derivation approach outlined by Sembach & Savage (1992). A single-component COG is assumed. Our procedure solves for N log and b independently in estimating the errors. We start with all the H I transitions for which we have estimated the equivalent widths and that do not appear a priori to be contaminated based on the AOD analysis. We check the results from the COG model and then remove only the transitions clearly departing from the model. We then run our COG program again with the new set of transitions, to obtain the final estimate of N H I with the COG method.
Using the normalized and masked H I spectra, we also perform a Voigt profile fit of the H I transitions using an inhouse PF program described in Lehner et al. (2011Lehner et al. ( , 2014Lehner et al. ( , 2018, which was updated and refined from Fitzpatrick & Spitzer (1997). A major difference from the previous methods is that the model profiles are convolved with the HIRES instrumental line-spread function, assumed to be Gaussian. The three parameters for each component i (typically i = 1, but in some cases i = 2 or 3, or in rare cases i > 3)-column density (N i ), Doppler parameter (b i ), and central velocity (v i )-are input as initial guesses and then subsequently varied to minimize the χ 2 of the fit to all the fitted H I transitions. As for the COG, we start with all the H I transitions that do not appear a priori contaminated and iterated, removing any transitions that a posteriori appear contaminated (or sometimes revisiting the masked pixels). In Figure 4, we show an example of PF that consists of a two-component fit corresponding to the absorbers at z abs = 2.55125 and 2.55216. As can be seen in each panel of Figure 4, there is some evidence of contamination especially in the component at 65 km s −1 (z abs = 2.55216), but even in that component there is enough information to derive robustly N H I when using the different H I transitions simultaneously.
In   require a more hands-on approach owing to having mostly saturated absorption in all the H I transitions and only weak damping features in Lyα. These were all analyzed as part of Lehner et al. (2014) and Lehner et al. (2016), and we refer the reader to these papers for a full description, noting that very conservative errors were adopted with a minimum error on the N H I for any LLS of 0.15 N log s = using this methodology. However, the metal lines associated with these absorbers were all reanalyzed as part of this work. The adopted N H I results are summarized in Table 2, and we also report in Table 3 the PF results, i.e., the velocities, Doppler parameters, and column densities, as well as errors associated with these quantities. In that table, the listed redshift corresponds to the adopted redshift (see above), and the third column indicates the component number for each absorber.

Empirical Properties of the KODIAQ-Z Survey
Before describing the ionization modeling of the absorbers, it is helpful to summarize the basic properties of the KODIAQ-Z survey, in particular in terms of N H I and redshift distributions, as we will compare our sample of absorbers to other surveys (such as surveys of stronger H I absorbers at similar redshift or absorbers with similar N H I but at lower redshift). As we show below, the Doppler parameters also Note. A 999.0 value corresponds to a nondetection, and the upper limit on the column density is quoted at the 2σ level. The detection flag has the following definition: 0 = detection; − 1 = upper limit (nondetection); − 2 = detection but a lower limit due to saturation of the line. Depending on whether a single transition or several were detected, a reliability flag is assigned as follows: 1 = most reliable, i.e., results based on several transitions or a nondetection (a nondetection is always reliable because it is estimated in an uncontaminated region of the spectrum); 2 = reliable, i.e., results based on at least two transitions, but where only one transition is detected at the 2σ level and the upper limits agree with that detection; 3 = less reliable, i.e., results based only on a single transition or several saturated transitions.
(This table is available in its entirety in machine-readable form.) motivate the use of photoionization as the main ionization mechanism in the ionization modeling to determine the metallicity of the absorbers.

H I Column Density and Redshift Distributions
In Figure 5, we show the distributions of the absorber redshifts and H I column densities, as well as the scatter plot between these two quantities. ) and 〈z abs 〉 = 2.79 ± 0.31 (the median being z abs = 2.75), respectively. There is an evident peak in the redshift distribution at 2.4  z abs  2.6 where there are about 65 absorbers; otherwise, between z abs ; 2.6 and 3.35 the z-distribution is relatively flat, with about 30-35 absorbers in each 0.2 redshift bin. This peak in the redshift distribution corresponds to redshifts where the optimal conditions are assembled to derive N H I (access to the entire Lyman series transitions and smaller contamination from the LYAF): as z decreases below 2.4, the number of H I transitions diminishes rapidly owing to the available wavelength coverage; as z increases above 3.6, the level of contamination from the LYAF increases (and indeed at 3.5  z abs  4.4, only for five absorbers were we able to robustly derive N H I ). As alluded to before, we use two other main surveys at z > 2 for comparison with KODIAQ-Z: the HD-LLS survey (Prochaska et al. 2015;Fumagalli et al. 2016b) and the R12-DLA survey (Rafelski et al. 2012). Based on the redshift distribution of KODIAQ-Z, we limit these surveys to the redshift range 2.2 z  3.6 when comparing the absorbers as a function of N H I over the same redshift interval.
Due to the selection of the absorbers, the N H I distribution should follow to some extent the column density distribution of H I (Prochaska et al. 2014), which is why the N H I distribution is not flat. From approximately 10 15 to 10 16.5 cm −2 , we see in Figure 5 the expected drop in the number of absorbers as N H I increases. At N log 15 H I < , our sample is not complete, as there is a sharp drop in the number of absorbers. This is because as N H I decreases, we rely more and more on high-S/N spectra in order to be sensitive to low metallicities, significantly reducing the sample.

 
should be less frequent than the previous category, there are more of these strong LLSs/ SLLSs in our sample owing to the fact that we can use the damping of Lyα and Lyβ to derive N H I (see also Section 3.2).  Target

H I Column Density and Doppler Parameter Distributions
In Figure 6, we compare the Doppler parameter, b, and N H I derived from the PF in the individual components of H I. There is no strong trend between b and N H I , although broad components with b > 40 km s −1 are more frequent at N log 16 H I < (b = 40 km s −1 implies T < 10 5 K, which is the threshold b-value between broad and narrow H I absorbers as defined by Lehner et al. 2007). 18 We, however, note that the majority of these broad components are found in more complex velocity profiles (two or more components), lower-S/N spectra, and/or strong H I absorbers where the b-value is not well constrained. It is quite possible that in most of these cases narrow components are blended, mimicking a broad component.
Considering the entire sample, the mean and standard deviation of b are 〈 b〉 = 30 ± 11 km s −1 , and the median is 28 km s −1 . About 90% of the components have 13.3 b 40 km s −1 , where 〈b〉 = 27 ± 6 km s −1 and the median is 27 km s −1 . The latter value corresponds to a temperature of the gas of T < 4 × 10 4 K, consistent with the gas being primarily photoionized. At z < 1 for the CCC absorbers, the results are remarkably similar since 91% of the profile-fitted absorbers have also b < 40 km s −1 and a mean of 28 ± 8 km s −1 (CCCI and Figure 7 therein). This contrasts with the weaker LYAF absorbers with N log 14 H I  , where the median and mean b-values are systematically larger by 15%-30% at z  0.5 than at 1.5  z  3.6, implying that a larger fraction of the low-z IGM is hotter and/or more kinematically disturbed than the high-z IGM (Lehner et al. 2007).

Motivations for Photoionization Modeling
Several empirical studies at both low and high redshifts have shown that SLFSs, pLLSs, and LLSs have properties that are consistent with the gas being predominantly photoionized (Prochaska 1999;Lehner et al. 2013;2016;Prochaska et al. 2017b;CCCI;CCCII;Crighton et al. 2013;Fumagalli et al. 2016b;Cooper et al. 2021;Zahedy et al. 2021). In Section 4.2, we show that the KODIAQ-Z absorbers have temperatures on average T  4 × 10 4 K, consistent with the gas being photoionized (the same result was found for the low-redshift absorbers with similar N H I ; CCCI). We also observe that when detected metal ions, such as the suite of carbon ions (C II, C III, C IV) or silicon ions (Si II, Si III, Si IV), align well with H I (see velocity profiles provided in the Appendix), this is consistent with these metal ions and H I probing a single ionization gasphase (we, however, emphasize again that sometimes metal ions, due to their smaller b, can have more than one component within the breadth of the H I profiles even though the H I absorption shows only a single componentthat also implies that the thermal broadening is a large contributor to the broadening).
In contrast to the low redshift universe, for the ionizing radiation field and typical lower metallicities at z ∼ 2-3 (about 0.1% solar or [X/H] = −2; see below; see also Fumagalli et al. 2016b;Lehner et al. 2016), even strong transitions like C II λ1334 and Si II λ1260 are often not detected, so we have to rely more on intermediate (Si III, C III, S IV) and high (C IV) ions to determine the metallicity. Because the current survey probes lower-N H I absorbers than our previous surveys at similar redshifts (Lehner et al. 2014, we find that the O VI absorption can often be narrow and align with the absorption seen in H I (and other lower metal ions if present). This contrasts with the very broad and strong O VI absorption frequently found in strong LLSs, SLLSs, and DLAs, which cannot be explained by photoionization (Lehner et al. 2014); it is quite plausible that this highly, collisionally ionized gas has also a different metallicity than the cooler, photoionized gas as discussed in Lehner et al. (2014). When O VI is narrow, it is used to constrain the photoionization models, and we find a posteriori that O VI is consistent with being produced by photoionization in many of these cases. The higher intensity of the EUVB at high z than at low z explains that O VI absorption associated with H I absorbers with N log 17.2 H I  can be more often easily photoionized in low-density gas with densities that are reasonable (i.e., densities that do not imply extremely large path lengths for these absorbers, and hence a Hubble flow broadening that would be inconsistent with the observed broadening of the absorption lines).

Photoionization Modeling
To determine the metallicity for each absorber, we follow the same methodology and make similar assumptions to those described and discussed in CCCII and CCCIII. As laid out in the previous section, due to the higher-redshift absorbers and the change in the EUVB, we include high ions in the photoionization modeling. We model the photoionization using Cloudy (version C13.02; see Ferland et al. 2013), assuming a uniform slab geometry in thermal and ionization equilibrium. In all cases the slab is illuminated with a Haardt-Madau EUVB radiation field from quasars and galaxies (HM05 and HM12; see Haardt & Madau 1996. We adopt HM05 (Haardt & Madau 2001, as implemented in Cloudy as the fiducial radiation field for KODIAQ-Z to allow for a direct comparison with the CCC results at lower redshift. However, we also use the HM12 EUVB to explore systematics associated with uncertainties in the radiation field. In Figure 7, we show the HM05 and HM12 EUVB at the average redshift of our survey (z = 2.8) along with the ionizing energy ranges for key ions used in KODIAQ-Z to constrain the metallicity of the absorbers. We also show in this figure the recent Khaire & Srianand (2019) EUVB models for their three favored slopes of the EUVB. While there are some differences between these EUVBs, these are not as large as at low z (Gibson et al. 2022). And indeed, we show a posteriori that the systematic errors arising from using a different EUVB are not as large as observed at z  1. We note that for stronger absorbers Fumagalli et al. (2016b) show that local ionizing sources did not affect much the metallicity estimates, and therefore we do not expect that local sources of (galactic) ionizing radiation would drastically change the metallicities in KODIAQ-Z. We adopt the same grid of Cloudy models as summarized in Table 3 of CCCII.
The main variables for the photoionization models are the ionization parameter-U ≡ n γ /n H , the ionizing photon density/total hydrogen number density (neutral + ionized)-and the metallicity, [X/H]. We assume solar relative| heavyelement abundances from Asplund et al. (2009) but allow for possible variation between between C and α-elements (where α can be O or Si) since this ratio can be sensitive to nucleosynthesis effects (see, e.g., Cescutti et al. 2009;Mattsson 2010, for more details). [C/α] is therefore not necessarily solar in the gas probed by absorbers studied here (see Aguirre et al. 2004; see also Lehner et al. 2013; CCCII, at z < 1).
To derive the metallicities, we compare the column densities for each absorber with a grid of photoionization models in a Bayesian context. We make use of the publicly available codes described in Fumagalli et al. (2016b; see also Prochaska et al. 2017b) and CCCII. 19 We perform the Markov Chain Monte Carlo (MCMC) analysis on an absorber-by-absorber basis. For each absorber we apply Gaussian priors on N H I and redshift. 20 For all the absorbers, we first assume a flat prior on U. After examination of the model outputs (in particular, convergence of walkers, comparison of predicted and observed column densities, and "corner plots" showing the posterior probability density function (pdf's)), we can assess which absorbers can be reliably modeled with this flat prior on U and which may require a Gaussian prior on log U to help the models converge (see below). Similarly, during this first iteration, we also determine which absorbers can be modeled with a flat prior on [C/α] (i.e., for which we could estimate [C/α]), which absorbers require [C/α] = 0 (e.g., those with no detections for all the ions or no carbon ions), and which absorbers may require a Gaussian prior on [C/α] (absorbers with poor constraints on carbon ions and/or α-elements).
When the constraints from the metal ions are not optimal, we adopt a Gaussian prior on log U (CCCII). These concern only SLFSs and pLLSs; the 16 LLSs and 7 SLLSs have enough metal-ion constraints for the ionization models to converge. A total of 58% of the SLFSs and pLLSs do not require a Gaussian prior on log U. In Figure 8, we show the log U median values as a function of N H I for the absorbers that can be modeled assuming a flat prior on log U. As observed at low z (CCCIII; CCCII), a strong anticorrelation is observed between the ionization parameter and N H I , with a slope steeper by a factor 2.3 compared to the low-redshift absorbers. To characterize the distribution of log U appropriate for use as a Bayesian prior, we split the N H I range into two bins, N 14.4 log 16.2 H I <  and N 16.2 log 17.2 H I <  . In the former, 52% (75/143) can be modeled with a flat log U prior, while for the latter that number is 81% (30/37). In both cases, we find that each pdf is well-fit with a normal distribution with U log 1.43 0.59 á ñ = - and −2.06 ± 0.53. We adopt these results for our Gaussian prior on log U in the N 14.4 log 16.2 For the Gaussian prior on [C/α], we similarly use absorbers with relatively good constraints from carbon and α-element ions to determine a mean and standard deviation of We apply these for the Gaussian prior on the [C/α] ratio for absorbers with poorly constrained [C/α]. We note that for a few (13) absorbers, [C/α] is clustered around − 1  [C/α]  − 0.9. Not including those in the mean would of course increase the mean [C/α] value, but it would not change the posterior pdf of [C/α] for these absorbers. We discuss in Sections 5.3.4 and 6.5 in more detail the [C/α] ratio. We summarize the results of our photoionization modeling in Table 4 for all the absorbers in the KODIAQ-Z sample. For each absorber, we list the sight-line name, absorber redshift (z abs ), log N H I , metallicity, ionization parameter ( U log ), and [C/α] ratio (when estimated). Each quantity is reported with the 68% CI and median values, except in the cases where we derive an upper or lower limit, in which case we report the 80% CI and median. We emphasize that these values are attempts to summarize the posterior pdf's provided by the MCMC sampling of the CLOUDY models. While those pdf's can be well behaved and nearly Gaussian, some are not as well behaved (see corner plots in the Appendix and supplemental materials).
For the HD-LLS absorbers used to compare with KODIAQ-Z, we apply the same overall methodology and the same HM05 EUVB in the photoionization models in order to reduce any systematic errors between the two surveys. However, we note that the absorption in the absorber in the HD-LLS and KODIAQ-Z is not treated in the same way. In the HD-LLS survey, a velocity interval of ±500 km s −1 defines an absorber, which is driven by the methods to derive N H I that use either the break at the Lyman limit or the damped wings from Lyα. For the metal lines, they sum the absorption components identified within that interval even though it is rarely larger than 200 km s −1 . As detailed in Section 3, we define an absorber using the Lyman series, and as much as possible we consider the absorption component by component, so in a window of ±200-500 km s −1 KODIAQ-Z can have several absorbers. However, as N H I increases and overlaps with the HD-LLS survey, i.e., for N log 18 H I  , we lose our ability to derive the column densities in individual components, especially when we have information only from Lyα and Lyβ. In these cases, we rely on the low ions to define the integration range in the intermediate and low ions. In our approach, we make the physically motivated assumption that the optical depths of the low ions follow that of the stronger H I absorption while the higher ions follow better weaker H I absorbers that may not contribute much to the total H I column-but the extra absorption in the higher metal ions could dominate their column density (Section 3). However, a posteriori we will show that this different treatment does not lead to large different results.

Uncertainties on the Metallicities
All the errors for the output quantities (e.g., metallicity) reported in this work are statistical errors from the Bayesian MCMC ionization modeling. There are additional systematic and other uncertainties that we discuss here.

Ionizing Background
The first one is caused by the uncertainty in the ionizing EUVB as alluded to previously. At low redshift, Wotta et al. (2019) show that the metallicity is systematically shifted on average by about + 0.4 dex for the SLFSs/pLLSs and + 0.2 dex for the LLSs when the ionizing EUVB changes from HM05 to HM12 (see also Gibson et al. 2022 for a study with other EUVBs). We have done a similar experiment with the KODIAQ-Z and HD-LLS absorbers, showing that the impact of the EUVB is overall smaller at high z than observed at low z.
In Figure 9, we show the comparison between the metallicities using HM12 and HM05 for the KODIAQ-Z absorbers where we did not apply a Gaussian prior on U log (which essentially removes all the limits, although including those would not change the results). The immediate conclusion from this figure is that the systematic differences between the metallicities derived using HM05 and HM12 are quite small, and much smaller than observed at low z. On average, considering all the absorbers, we   As illustrated in Figure 10, where the differences of the median metallicities between HM12 and HM05 versus N H I are shown for both the KODIAQ-Z and HD-LLS absorbers, that difference decreases with increasing N H I as observed for the low-redshift absorbers. In fact, for N log , the HM12 and HM05 EUVBs essentially give the same metallicities and any small differences are essentially random. For the SLFSs and pLLSs, 〈[X/H] HM12 − [X/H] HM05 〉 ; + 0.12 ± 0.15. These results are not entirely surprising when considering the different EUVBs in Figure 7, with their difference being smaller than at low z over similar ionizing energies (Gibson et al. 2022). Based on these results and Figure 7, the metallicities would not change much if we had used the more recent EUVBs from Khaire & Srianand (2019). For comparison with CCC, we adopt hereafter the results from the HM05 EUVB.

Multiple Ionization Phases
Throughout we also assume a single ionization phase. In cases where a higher ion (O VI or sometimes C IV in our survey) is clearly underproduced by the photoionizing model, this given high ion is not included in the final model. In the latter case, this would be evidence for an absorber probing multiple gas phases. As shown by Lehner et al. (2014), this can especially happen in the strong LLSs, where O VI can be very strong and can rarely be modeled with a single gas-phase ionization model with the lower ions. In these cases, the velocity profiles of the O VI profiles are also very different from those of the lower ions (including C IV), strongly hinting at multiple gas phases.
As we explore lower-N H I absorbers in this present survey, we find in many cases that a single ionization gas phase can in fact match well the observed metal-ion column densities, including those of O VI. In these cases and in contrast to stronger H I absorbers, the absorption is also often quite simple (typically a single component), and the velocity profiles of the different metal ions align well with each other (including O VI when present) and with the H I velocity profiles. This, of course, does not necessarily mean that a more complex ionization structure might not be present, but a single gasphase model is often sufficient to reproduce the observed column densities, especially in the SLFS and pLLS regimes.
We adopt here therefore the simplest approach and assumption, a methodology also employed in CCC. Nevertheless, we note that, at least at low redshift, the effects on the metallicities between single versus multiple gas-phase modelings typically lead to small changes in metallicities of the order of < 0.1-0.2 dex (Howk et al. 2009;Haislmaier et al. 2021). Therefore, we do not expect that considering a more complex ionization structure for the absorbers in our sample would drastically change the metallicities.

Ions Included in the Models
Depending on wavelength coverage, contamination, and H I column density, the ions used to constrain the photoionization Note. The lower to upper bounds for each quantity represent the 68% CI for detections and the 80% CI for the upper or lower limits; the middle number is the median metallicity. When a colon is present after the median value of log U, this implies that a log U Gaussian prior was used to determine the properties of the absorber.
(This table is available in its entirety in machine-readable form.) Figure 9. Comparison of the median metallicities (with 68% CI) of KODIAQ-Z derived using the HM12 and HM05 EUVBs. The dashed line is the 1:1 relationship. The dotted line shows the mean of the differences between the metallicities derived using HM12 and HM05. models can be different. This may introduce some additional systematic uncertainties when comparing the metallicities across N H I and z. For example, the suite of carbon ions available in the redshift probed by our sample consists of C II, C III, and C IV. Sometimes all these ions are available, but it is not always the case. In particular, the inclusion of C IV can a priori change the density (ionization parameter) to lower (higher) values owing to the 48 eV ionization energy required to ionize C III to C IV. This change in log U could possibly affect also the metallicities. At z < 1, Lehner et al. (2019) investigated the effect on the metallicity and other quantities (e.g., log U) if a high ion (in that case O IV) is included or not in the photoionization models. They estimated that this effect would only change the metallicities on average by 0.06 dex, with the largest metallicity difference being around 0.25 dex. However, and importantly, for the physical conditions of the absorbers such as U, n H , and N H , that effect can modify the values for these quantities on average by 0.5 dex.
Using our sample, we specifically undertook a similar experiment but using C IV instead of O IV. To select this sample, we require that the absorbers have coverage of C IV and C II and/or C III (and often there is also coverage of other ions, such Si II, Si III, and Si IV). For 38 absorbers with N 14.8 log 20 H I   , we systematically removed C IV in the photoionization models and ran the same photoionization models. In Figure 11, we show the comparison of [X/H] and log U as a function of N H I for these 38 absorbers with and without C IV included in the photoionization models. There is no trend in the difference between these quantities and N H I . The errors are smaller when C IV is included in the models, which is expected. The overall impression from this experiment is that neither the metallicities nor the log U values (or quantities not shown here, such as n H , N H ) drastically change with the inclusion or not of C IV in the models. The mean difference in metallicities between absorbers with and without C IV is −0.06 ± 0.21 dex. Only 8% of the 38 absorbers have metallicity differences larger than 0.4 dex (0.4−0.7 dex). For log U, the mean difference is, however, much smaller than at low redshift with −0.04 ± 0.29 dex. In that case, 10% of the absorbers have log U differences larger than 0.4 dex (0.4-1.0 dex). Therefore, the change in the ionic constraints between the different absorbers at 2.2  z  3.6 is unlikely to produce large systematic uncertainties on the metallicities and other physical parameters such as the density or ionization parameter.

Nonsolar Relative Abundances
Other uncertainties may arise from dust or nonsolar relative abundances. Since carbon, silicon, and oxygen ions are typically used to constrain the ionization modeling at high redshift, we explicitly take account of possible nonsolar [C/α] by letting this ratio vary in the range − 1 [C/α] + 1. In Figure 12, we show [C/α] as a function of N H I (left panel) and the neutral fraction f H I ≡ N H I /(N H I + N H II ) (right panel) for the KODIAQ-Z and HD-LLS absorbers. We also show in this figure the results from CCC at z  1 (CCCIII). For [C/α]−N H I , a Spearman rank-order test shows no strong evidence for a correlation between [C/α] and N H I with a correlation coefficient r S = 0.14 and a p-value = 0.04. On the other hand, the same statistical test shows a positive monotonic correlation between [C/α] and f H I with a correlation coefficient r S = 0.44 and a p-value = 0.05%. However, removing the 13 KODIAQ-Z absorbers with low [C/α]  − 0.9 eliminates any evidence of correlation between [C/α] and f H I with a correlation coefficient r S = − 0.05 and a p-value = 0.88.  So far CCC, HD-LLS, and KODIAQ-Z have been treated as a single sample. Considering the samples individually would not change the conclusions for the lack of correlation between [C/α] and N H I . However, while neither CCC nor HD-LLS shows any correlation between [C/α] and f H I , for KODIAQ-Z the Spearman rank-order test shows a strong positive monotonic correlation between [C/α] and f H I with a correlation coefficient r S = 0.72 and a p-value = 0.05%. In that case, removing the 13 absorbers with low [C/α]  − 0.9 still yields r S = 0.55 and a p-value = 0.05%. Therefore, in the KODIAQ-Z sample, [C/α] is more affected by ionization correction than in the other samples, most likely because larger ionization corrections are applied (for CCC, C II and C III are used to constrain [C/α] instead of C IV, C III, and C II for KODIAQ-Z, while for HD-LLS the ionization corrections are typically smaller). In the CCC and HD-LLS samples, we always find [C/α] − 0.6. Hence, [C/α] values with [C/α] < − 0.6 are most likely predominantly caused by ionization rather than nucleosynthesis effects. Excluding the 13 absorbers with low [C/α]  − 0.9, we find that 〈[C/α]〉 = − 0.20 ± 0.31 in the KODIAQ-Z sample compared to 〈[C/α]〉 = − 0.05 ± 0.30 in the CCC survey and 〈[C/α]〉 = − 0.05 ± 0.24 in the HD-LLS survey. The similarity between CCC and HD-LLS is quite remarkable since different z, N H I , and [X/H] are probed. The 1σ dispersion is about the same in the three surveys, but on average the KODIAQ-Z sample has a systematic effect of −0.15 dex on [C/α], which is most likely caused by ionization correction. We can test this further using the sample of 38 absorbers in Section 5.3.3, where C IV was removed in the ionization models. From this sample, there are 26 absorbers for which there is a good constraint on [C/α]. We find that the mean of the [C/α] median values increases from −0.38 ± 0.35 dex when C IV is included in the ionization models to −0.20 ± 0.35 when C IV has been removed. There is therefore a systematic ionization effect of about 0.18 dex on [C/α]. In conclusion, the variable C/α is important in that we can marginalize over that as a "nuisance" variable in order to get a better sense of the uncertainties in the metallicity, but it may not provide deep insights on the nucleosynthesis effects that may affect C/α (see also Section 6.5).
Finally, we note that for stronger LLSs and SLLSs such as those in the HD-LLS survey (Fumagalli et al. 2016b), dust depletion is found to be a small effect for strong LLSs and even SLLSs at high z and hence not very likely to have any appreciable effect on the metallicity estimates, especially considering that the typical species-carbon, silicon, oxygen -used to constrain the ionization models are known not to be strongly depleted onto dust (e.g., Savage & Sembach 1996;Jenkins 2009;Jenkins & Wallerstein 2017).

Metallicity Distributions
In Figures 13-15, we show the posterior pdf's of the metallicities of the SLFSs, pLLSs, and LLSs for KODIAQ-Z. In Figure 15, we also overplot the metallicity pdf of the LLSs from the HD-LLS survey, showing similar pdf's between the two surveys, and hence implying that we can combine both samples to improve the statistics for the LLSs. These posterior pdf's are constructed by combining the normalized metallicity pdf's of all of the absorbers within a given H I column density regime. In Table 5, we summarize the medians, means, standard deviations, and interquartile ranges (IQRs) for the various absorbers, as well as a combination of these. We also include in this table the results for the HD-LLS, KODIAQ-Z +HD-LLS, and R12-DLA samples.   Considering first the SLFSs, the metallicity pdf is not unimodal. It is negatively skewed with a long and prominent tail extending well below [X/H] < − 3. Its main peak is around the median value of −2.4 dex. The metallicity pdf of the SLFSs dips around −1.1 dex, with a second weaker peak around −0.6 dex. The IQR is nearly 2 dex from −3.6 to −1.8 dex. As N H I increases, the overall median metallicity increases and the IQR decreases. For pLLSs+LLSs, the median metallicity is −2.2 and the IQR is only ∼1 dex from −2.8 to −1.7 dex. These features are distinctly illustrated in Figure 16, where we compare the pdf's of the SLFSs and pLLSs+LLSs: there are more frequently pLLSs+LLSs in the metallicity range − 3.2  [X/H]  − 1.2 than SLFSs, but SLFSs become more frequent at extremely low metallicity [X/H]  − 3.5. Compared to the SLFS and pLLS pdf's, the LLS pdf is not as negatively skewed, with about the same values for the mean and median.
For completeness, and if one ignores the H I column density dependence of the metallicity pdf's between the SLFSs, pLLSs, and LLSs, we show in Figure 17 the resulting metallicity pdf of the absorbers with N 14.5 log 19 H I < < at 2.2  z  3.6. Since the SLFSs are about 2/3 of the sample, the metallicity pdf of the N 14.5 log 19 H I < < absorbers is not too different from that of the SLFSs, with a similarly skewed distribution toward low metallicities but also a second weaker peak around −0.6 dex.

Metallicity versus N H I
As shown in the previous section, there is already evidence for an overall increase of the metallicities with N H I in the range     H I   . We now explore further how the metallicity varies with N H I considering the higher-N H I absorbers (the SLLSs and DLAs). From Table 5, the trend described above of increasing median/mean metallicity with increasing N H I continues to apply in the SLLS and DLA regimes. The IQR and standard deviation for the SLLSs are about similar to those of pLLSs+LLSs, but those for the DLAs are substantially smaller (the IQR is 0.65 dex smaller than that of the SLFSs). The negative skewness is also less pronounced for the SLLSs and absent for the DLAs. In fact, the DLA metallicity pdf, contrary to the pdf's of the other weaker absorbers, is fully consistent with a Gaussian distribution according to the Kolmogorov−Smirnov test with a p-value =0.1% (see also Rafelski et al. 2012).
These conclusions are further strengthened considering the cumulative distribution functions (CDFs) of the metallicity pdf's (see Wotta et al. 2019 for the description of the estimation of the CDFs including upper and lower limits). We show in Figure 18 the metallicity CDFs for the SLFSs, pLLSs, LLSs, SLLSs, and DLAs. With only 24 pLLSs, this regime remains the least sampled. Nevertheless, the CDF confirms the overall evolution of the metallicities as N H I increases: very metal-poor (VMP) absorbers with [X/H] − 2.4 are more likely to be found in the SLFS, pLLS, and LLS regimes than in the SLLS or DLA regime ([X/H] VMP = − 2.4 is the 2σ lower bound of the DLA metallicities). Notably the median metallicity of the SLFS is right at that threshold of [X/H] VMP = − 2.4, implying that 50% of these absorbers are VMP while only <5% of the DLAs may be VMP by definition.
Thus far, we have separately analyzed the absorbers in different N H I categories. However, we can also simply plot the metallicities of the absorbers as a function of N H I , which is shown in Figure 19. This has the obvious advantage that one does not have to make an a priori differentiation between the different absorbers in classes of H I column densities. For the SLFSs, pLLSs, LLSs, and SLLSs with well-constrained metallicities (i.e., not including the lower and upper limits), the central values represent the median of the posterior pdf's, and the error bars represent the 68% CI. For the upper and lower limits, the downward-and upward-pointing triangles give the 90th and 10th quantiles, while the vertical bar gives the 80% CI. For the DLAs, the best estimates with their 1σ error bars are shown. The  . While the VMP definition implies that this is VMP gas for DLAs, for low-N H I absorbers the VMP definition corresponds to the median value for SLFSs and is about 0.3 dex below the median of pLLSs and LLSs. As already pointed out above, the frequency of VMP absorbers also increases as N H I decreases, which is quantitatively summarized in Table 6. We therefore also define extremely metal-poor absorbers as absorbers with [X/H] < − 3, which are mostly observed at some level in    (Fumagalli et al. 2016b), but we reestimated the metallicities of these absorbers following the same EUVB and methodology as in KODIAQ-Z. The purple data (DLAs) are from Rafelski et al. (2012).  Frebel et al. 2007;Crighton et al. 2015).
We finally note that the metallicity dip around −1.1 dex observed especially in the SLFS metallicity pdf is also evident in the nonbinned distribution shown in Figure 19. In that figure, there are a lack of data points at − 1. . This is reminiscent of the metallicity dip at a similar value of [X/ H] ; − 1 observed in the pLLS+LLS pdf at 0.45  z  1 ; see also Section 6.3 for the cosmic evolution of the metallicity).

Cosmic Evolution of the Metallicity
In Lehner et al. (2016), for the first time, we studied the redshift evolution of the metallicity absorbers with N 16.2 log 19 H I <  , but with a much smaller sample at both low and high z. The samples of absorbers at all surveyed z have increased at least by a factor 8 and provide far more robust results. Furthermore, we also now sample at low and high z the N H I regimes below N log 16.2 H I  . In Figure 20, we show the metallicity pdf's at low (z < 1) and high (2.2 z 3.6) redshift for H I-selected absorbers with 15  N H I < 19. There is some overlap between the two distributions between − 3  [X/H]  − 1.4, but evidently there are more absorbers with [X/H]  − 2 at high z and many more with [X/H] > − 1.2 at z < 1. In fact, the latter value corresponds to the median metallicity at z < 1 for N 16.2 log 19 H I <  absorbers, and only 10% of the high-redshift absorbers have metallicities [X/H] − 1.2. The metallicity median and mean at z < 1 are about +1 dex higher than at 2.2 z 3.6. Both metallicity distributions are negatively skewed with similar IQRs ; 1.5 dex at low and high z.
In Figure 21, we compare the metallicities of the absorbers (following the same methodology used to produce Figure 19) at z < 1 and 2.2 z 3.6 over the H I column density range 14.5  N H I < 22. For the low-redshift absorbers, we shifted the metallicity by −0.9 dex, which is about the average difference between the mean/median metallicity of the absorbers at low and high z (−0.8 dex for DLAs, −1 dex for lower-N H I absorbers). The notable feature from this figure is that despite different metallicity behaviors at low and high H I column densities, the offset metallicities and the overall dispersions at low and high z show a striking overlap at any given N H I (where there is a good sampling in both data sets). A closer look reveals a small excess of higher-metallicity absorbers in the SLFS, SLLS, and DLA regimes at high z (i.e., a slightly larger dispersion of the metallicity distribution). However, overall the major change between 2.2 z 3.6 and z < 1 absorbers is an increase of their metallicity by a factor ∼8.
In Figure 22, we now show the metallicities as a function of the redshift for the absorbers included in CCC, HD-LLS, and KODIAQ-Z (for clarity we do not show the R12-DLA sample in this figure). For the HD-LLS sample, we differentiate the LLSs and SLLSs. This figure strengthens the previous conclusions, showing the overall increase in metallicities as z decreases. At z < 1, there is a substantial fraction of N log 19 H I < absorbers with − 0.15  [X/H]  + 0.15, while these are mostly absent at higher redshift. On the other hand, at z < 1 there is no evidence of population of absorbers with [X/ H]  − 3, while there is at 2.2 z 3.6. As already noted by Lehner et al. (2016) with a much smaller sample, while the VMP threshold increases by 1 dex over 2.2 z 3.6 to z < 1 (the 5% quantile of DLA metallicities from −2.4 to −1.4 dex), a similar fraction of about 50% of N log 18 H I  VMP absorbers at low and high z lie in the respective VMP regime.
In Figure 22, there is also some evidence that the metallicity distribution changes to a certain degree around z ∼ 2.85. Considering the KODIAQ-Z and HD-LLS sample with   Figure 19 and Figure 9 of CCCIII for the low-redshift sample). Median metallicities are plotted, but we treat upper and lower limits as in Figure 19, i.e., the lower and upper values represent the 10th and 90th percentiles, respectively.
N log 19 H I < absorbers at z > 2.85, the fraction of absorbers with [X/H] − 1 is very small (4/112 ; 3.6%), while at z 2.85 it is a factor ∼4 larger (17/130 ; 13.1%). To explore that potential of the change of the metallicity with z, in Figure 23 we show the comparison between the metallicity CDFs of the absorbers with N 14.5 log 19 H I <  from KODIAQ-Z and HD-LLS above and below several redshift thresholds, z th = 2.8, 2.9, and 3.0. The large number of absorbers with [X/H]  − 1 at z  2.9 compared to higher z is further demonstrated in this figure. However, this figure also shows that at z  2.9, and even more so at z 3.0, there is a much larger fraction of metal-poor absorbers with [X/H]  − 2.8 than at lower redshifts, demonstrating that within the redshift interval 2.2  z  3.6 there is a change in the metallicity pdf's with an overall increase of the metallicity at 2.2  z < 2.8 compared to 2.8  z  3.6.

Metallicity Variation over Small Velocity Scales
Due to the treatment of the absorbers where we separate individual H I components as much as possible (see Section 3), our KODIAQ-Z survey also provides information on the metallicity variation over small redshift/velocity separation between absorbers. Following Lehner et al. (2019), we define here paired absorbers as absorbers that are closely separated in redshift space so that v z z z c 1 abs 2 abs 1 abs 1 |( ) ( ) | D º -+ < 500 km s −1 . We have a sample of 36 such paired absorbers.
In Figure 24, we show the median or upper limit on the metallicities of the 36 paired absorbers as a function of their absolute velocity separation, where the two absorbers of each pair are connected by a dotted vertical line. In Figure 25, we show in the left panel the absolute metallicity difference (Δ[X/H]) between the higher and lower median metallicities of the paired absorbers against their absolute velocity separation and in the right panel the histogram distribution of Δ[X/H]. From these two figures, at any velocity separation, there is a large scatter from about 0.2 to >2 dex. At Δv  150 km s −1 , the paired absorbers are more clustered above [X/H] ; − 2.4 and their metallicity difference is smaller on average than for larger Δv. Using the survival method to account for the limits (Feigelson & Nelson 1985 Ignoring the velocity separation, we also observe that there are two main clusters separated at Δ[X/H] ; 1 with about 50% of the sample in each bin and with a mean 〈Δ[X/H]〉 = 0.39 ± 0.05 and 2.15 ± 0.13 below and above this limit, respectively. Therefore, for about half of the paired absorbers, there is evidence for metallicity variations over Δv  500 km s −1 of a factor 2-3 and for the other half of   a factor >140. For absorbers with Δv < 150 km s −1 , the metallicity variation is on average around a factor 5 but grows to ∼50 for absorbers with Δv > 150 km s −1 . Due to the large metallicity variation over Δv  500 km s −1 , a posteriori it is justified to treat these paired absorbers as individual absorbers. The other direct consequence is that, when possible, components in absorbers should be treated individually to derive the true metallicity of the gas rather than a column-density-weighted average metallicity between very different types of gaseous components. This was a consequence first inferred by Prochter et al. (2010) at similar redshifts for one LLS and then in the CCC survey at z < 1 for 30 paired absorbers Wotta et al. 2019; see also Kacprzak et al. 2019;Zahedy et al. 2021). However, we also find that for the smaller Δv 150 km s −1 separation the metallicity variation is substantially smaller than for the paired absorbers with larger Δv. We discuss the implications of these findings in Section 10.4 in more detail.

The [C/α] Ratio as a Function of the Metallicity
In the photoionization modeling used to determine the metallicities of the absorbers, the [C/α] ratio is allowed to vary in the range − 1  [C/α]  + 1 to take into account possible nucleosynthesis effects that likely arise from the time lag in the production between α-elements and carbon (e.g., Cescutti et al. 2009;Mattsson 2010) With KODIAQ-Z and our remodeling of the HD-LLS absorbers, we can now revisit these conclusions with a much larger sample at high redshift, especially in the low-metallicity regime. In Figure 26, we show [C/α] as a function of [X/H] for all the absorbers in KODIAQ-Z, HD-LLS, and CCC, where the relative abundance and metallicity are simultaneously constrained. In that figure, we use different symbols to differentiate LLSs from SLLSs in the HD-LLS survey (there is only one SLLS in KODIAQ-Z with information on [C/α]; see Figure 12; we did not highlight this in Figure 26). Excluding the 13 absorbers with [C/α]  − 0.9 (because those are more largely affected by ionization; see Section 5.3 and below), the addition of the KODIAQ-Z and HD-LLS data confirms the first two characteristics described above and summarized in Lehner et al. (2019). However, at [X/H]  − 2, there is no clear upturn in [C/α] for the absorbers with N log 19 H I < , although the scatter in the [C/α] distribution becomes smaller. On the other hand, considering the SLLSs from HD-LLS (diamond symbols in Figure 26), the upturn in [C/α] is apparent when the metallicity becomes smaller than at [X/H]  − 2.2.

. The general trend between [C/α] and
Since the derived ratio of [C/α] is controlled by both ionization corrections and nucleosynthesis effects, it is plausible that the ionization corrections wash away some of the trends that would appear if only nucleosynthesis effects were at play (such as for the DLAs). In Section 5.3.4, we show that [C/α] in the KODIAQ-Z sample is more affected by ionization than in the two other samples in view of the correlation between [C/α] and f H I . However, no such correlation is observed in CCC and HD-LLS, and the overall Figure 25. Left: scatter plot of the difference between the higher and lower median metallicities for the paired absorbers as a function of their absolute velocity difference. For the lower limits (blue triangles), we use the lower bound of the 80% CI to be the most conservative. Right: distribution of the metallicity differences between paired absorbers. 7. Physical Properties of the SLFSs, pLLSs, LLSs, and SLLSs at 2.2  z  3.6 The photoionization models discussed in Section 5 also predict the physical properties of the absorbers. In particular, these include the neutral fraction of the gas ( f HI ), the density of the gas (n H ), the total H column density (N H ≡ N H I + N H II ), the length scale (l ≡ n H /N H ), and the temperature of the gas (T). These can in turn be used to constrain the cosmological baryon and metal budgets for the absorbers (Fumagalli et al. 2016b;Lehner et al. 2019). An important caveat to keep in mind and discussed at length in Fumagalli et al. (2016b) is that there is a degeneracy between the density and intensity of the radiation field. This dependency produces much less robust inference on the density and sizes of absorbers than for the metallicity, which depends only on the shape of the radiation field. We also remind the reader that only photoionization by the EUVB has been used in this work, and it is quite plausible that more than one ionization process is at play, which would also affect more the physical quantities discussed in this section.

Densities, H Column Densities, Temperatures, and Length Scales
We summarize the physical quantities (n H , N H , l, and T) derived from our modeling for the KODIAQ-Z sample in Table 7. Values of n log H with a trailing colon represent absorbers for which we adopted a Gaussian prior on log U and therefore on n log H . In Table 8, we summarize the statistics of the physical properties derived for the entire sample (median, mean, and dispersion of N H I , n H , N H , l, and T) and a restricted sample where only a flat prior on U log was used. Comparing the results from Table 8 between the entire and restricted samples, it is apparent that the results are not statistically different. We emphasize that most of the pLLSs, LLSs, and SLLSs did not require a Gaussian prior on U log . However, for about 45% of the pLLSs we had to use a Gaussian prior on log U, and yet the results are essentially the same in both samples. Hereafter, we therefore consider the results from the entire sample. In Figure 27, we show the posterior pdf's of N H , n H , and l for the SLLSs and LLSs (using results from KODIAQ-Z and HD-LLS) in the top two rows and the pLLSs and SLFSs (from KODIAQ-Z) in the bottom two rows. These pdf's contrast remarkably from the metallicity pdf's with much narrower distributions.
Temperatures of ionized gas predicted by CLOUDY are based on the balance of atomic heating (photoionization) and cooling (recombination and collisional effects). As such, temperatures are dependent on the radiation field shape and intensity, density, and metallicity of the models. The distribution of predicted temperatures has overall a small dispersion, and there is a clear trend of decreasing temperatures of about 0.1-0.2 dex between each absorber category from the SLFSs, pLLSs, LLSs, to SLLSs (see Table 8), i.e., the temperature of the gas decreases as the H I column density increases. For the SLFSs, the average temperature is around 3 × 10 4 K, while it is a factor 2 less on average for the SLLSs. In Section 4.2, we show that for about 90% of the H I components 〈b〉 = 27 ± 6 km s −1 . This places an upper limit on the temperature from the observations of T < 4 × 10 4 K, consistent with the ionization modeling. It also implies that for the SLFSs, on average, the thermal broadening component is somewhat larger than the nonthermal component, while for the other absorbers nonthermal motions take over the thermal ones. Compared to Lehner et al. (2019), the predicted temperatures of the z < 1 SLFSs, pLLSs, and LLSs are about 1.4 times smaller than at 2.2  z  3.6.
Considering the SLFSs, pLLSs, and LLSs, the total H column densities decrease with decreasing N H I , but their distributions overlap within 1σ dispersion (which are large for Note. The lower and upper bounds for each quantity represent the 68% CI, while the middle value is the median value. When a colon is present after the median value of n log H , this implies that a log U Gaussian prior was used to determine the properties of the absorber (see Section 5.2).
(This table is available in its entirety in machine-readable form.) these quantities; see Table 8 and Figure 27). The differences of the mean log N H I for SLFSs compared with the pLLSs and of the pLLSs compared with the LLSs are a factor 4 and 28, respectively. In contrast, the differences in total H are significantly smaller with differences of a factor ∼4 and 6.5 between the SLFSs and pLLSs and between the pLLSs and LLSs, respectively. Thus, despite their lower N H I values, the contributions of pLLSs and SLFSs to the total mass and baryon budgets are higher compared with the LLSs than their H I columns would suggest, especially since they are more numerous (e.g., Prochaska et al. 2014; see also Section 8). The N H column densities for SLLSs are similar to those for LLSs; N H is about the same even though N H I is on average a factor 34 smaller for the LLSs than for SLLSs. Given the higher incidence of LLSs, they make a larger contribution to the cosmic baryon budget than the SLLSs. Compared to Lehner et al. (2019), N H at the z < 1 SLFSs, pLLSs, and LLSs are a factor ∼10-15 smaller than at 2.2  z  3.6.
The estimated hydrogen densities increase as N H I increases. The SLLSs have the largest dispersion on n H , which is caused by a long tail at low and high n H values (see Figure 27). While there is an increase in n H with increasing N H I , the bulk of the hydrogen density for the pLLSs, LLS, and LLSs is within about n log 3.4 H and −1.8, while for SLFSs it is between n log 4.2 H and −2.5. At z ∼ 2.8 (the mean redshift of KODIAQ-Z), the mean cosmic density is about n 10 H 5  cm −3 (Schaye 2001a). The pLLSs, LLSs, and SLLSs therefore probe overdensities of n n 1 40 H H d º -~-1600, while for the lower column density SLFSs δ ∼ 5-320. A large part of the observed pLLSs, LLSs, and SLLSs largely probe virialized (mini)halos (δ > 100; Schaye 2001a). For SLFSs, there is more a mixture of systems probing overdensities below and above 60, at the transition between the IGM and virialized halos (Schaye 2001a). Compared to the CCCIII, n H at the z < 1 SLFSs, pLLSs, LLSs, SLLSs are a factor ∼1.5-2.5 smaller than at 2.2  z  3.6. However, due to the redshift dependence of the mean cosmic density that is ∝ ((1 + z)/4) 3 (Schaye 2001a), the density contrast at z < 1 for the SLFSs and other stronger H I absorbers is always δ > 100, even for the lower densities probed by SLFSs (although we note that at z < 1 only N log 15 H I > SLFSs could be probed).
The length scales, which are directly related to n H and N H , have the largest dispersion among all these physical quantities. For the SLFSs and pLLSs, the length-scale IQR is about 4  l  150 kpc. For the LLSs, the IQR for l is 18  l  288 kpc, and for the SLLSs it is 8  l  80 kpc. In contrast, at z < 1, the length-scale IQR of these absorbers is 0.3  l  3.2 kpc ). At 2.2  z  3.6, the SLFSs, pLLSs, LLSs, and even the SLLSs may probe a variety of structures from gas embedded in virialized halos to largescale structures. Recent MUSE galaxy surveys show strong hints that at least LLSs may probe a wide variety of environments at z > 2, including gaseous filaments, CGM of galaxies, or intragroup gas (e.g., Fumagalli et al. 2016aFumagalli et al. , 2017Lofthouse et al. 2020).

Neutral Fractions
The behavior of the neutral fraction with N H I can be used to constrain the cosmic baryon and metal budgets (see Section 8). In Figure 28, we show the neutral fraction of the absorbers estimated from the CLOUDY models as a function of N H I from N log 14.6 H I  to 20 using KODIAQ-Z and HD-LLS (see Section 5.2; see also the recent work for the SLLSs from Berg et al. 2021). There is a clear correlation between f H I and N H I . This is confirmed by the Spearman rank-order test that shows a strong positive monotonic correlation between f H I and N H I with a correlation coefficient r S = 0.87 and a p-value = 0.05%. However, the slope of the correlation is shallower between the SLFS, pLLS, and LLS regimes compared to the LLS and SLLS regimes. The dashed red line in Figure 28 , which is consistent with the linear fit derived in the HD-LLS survey for the LLSs +SLLSs (Fumagalli et al. 2016b). The Spearman rank-order test also produces a similar conclusion: for the pLLSs+LLSs, the correlation coefficient is r S = 0.65 and the p-value is = 0.05%, while for the SLFSs, it drops to r S = 0.27 and a p-value = 0.011. There is a similar effect if we consider N H  instead of f H I . This trend can be understood as follows: as the gas becomes more optically thick, N H I increases more rapidly relative to the total H column density.
To model this overall trend between f H I and N H I , we use a Gaussian process (GP) model, providing a generic supervised learning method designed to solve a regression. We refer the reader to Appendix F in Lehner et al. (2020) for more details, but in short, the advantage of this method is that the prediction interpolates the observations in a nonparametric way and is probabilistic so that empirical confidence intervals can be computed. We use the Python GAUSSIAN PROCESS REGRES-SION package within SCIKIT-LEARN (Pedregosa et al. 2011;Buitinck et al. 2013) to model the data with a squaredexponential kernel (the use of a more complex kernel like a Matern kernel would yield similar results when using similar bounds). The overall scatter in the data is typically much larger than the errors on a single data point, and we use the 1σ dispersion derived for SLFSs, pLLSs, LLSs, and SLLSs to inform the GP model of the intrinsic scatter in the data (this can be understood as a prior factor to smooth out the scatter of the data). The black curve in Figure 28 is the GP model, with the dark area showing the standard deviation determined by the GP model. There is a large overlap with the two linear fits, but the GP fit has the advantage of providing a global model to the data without a priori selecting manually the inflection point in the curve. The ). Finally, we show in Figure 28 with the orange line a GP model if instead we had used the results obtained with the HM12 EUVB. This model is about within the standard deviation of the GP model using the HM05 data, demonstrating that a change of the EUVB from HM05 to HM12 would not have a large impact on the relationship between f H I and N H I (and therefore on the estimates of the cosmological baryon and metal budgets).

Cosmological Baryon Budget
Observations of the LYAF show that about 80%-90% of the baryons at 2 < z < 4 are found in the cool photoionized phase of IGM (e.g., Rauch et al. 1997;Weinberg et al. 1997;Kim et al. 2001). In Lehner et al. (2014), we show that the photoionized and collisionally ionized gas found in τ LL > 1 systems is very likely the second-largest contributor to the baryon budget at z ∼ 2-3. However, this result could be affected by small number statistics, as the Lehner et al. (2014) sample included only seven LLSs and eight SLLSs. Our earlier study also did not consider the baryon contributions from the SLFSs and pLLSs. Here we review the impact of the cool, photoionized gas in these absorbers on the cosmic baryon budget at 2.2  z  3.6.
The mass gas density relative to the critical density, Ω g , can be estimated via (e.g., Tytler 1987;Lehner et al. 2007;O'Meara et al. 2007;Prochaska et al. 2010;CCCIII): , we adopt the spline function model from Prochaska et al. (2014) at z ≈ 2.5 (see their Table 2 and Figure 7). For the neutral fraction, f H I , we adopt the GP model from Section 7.2. Integrating the above equation over the different N H I regimes of the SLFSs, pLLSs, LLSs, and SLLSs, we find that Ω g ; (3.3, 1.1, 1.1, 1.2) × 10 −3 , respectively. The contributions of Ω g to the total baryon density are Ω g /Ω b ; 6.8%, 2.3%, 2.3%, and 2.5% for the SLFSs, pLLSs, LLSs, and SLLSs, respectively (Ω b = 0.0486; Planck Collaboration et al. 2016). The values for the LLSs and SLLSs are consistent with the ranges in Lehner et al. (2014, their Table 6, column with Si IV). 23 Importantly, these estimates do not include a contribution from the highly ionized, hotter gas that is not photoionized. Lehner et al. (2014) found that the O VI collisionally ionized gas in LLSs and SLLSs could contribute 0.4%-10% to the baryonic budget. As noted in Section 5.1, for the pLLSs and SLFSs we often find that the O VI absorption is narrow and aligns with the absorption seen in H I (and other lower metal ions if they are detected), i.e., the properties of the O VI are quite different as N H I decreases to N log 17 H I < compared to stronger H I absorbers. The observed column densities for these narrow O VI components are also well reproduced by singlephase photoionization models such as those described in Section 5 where the source of the ionizing photons is EUVB. Providing a full description of the O VI associated with the  KODIAQ-Z absorbers is beyond the scope of this paper, but based on a visual inspection of the profiles and our ionization modeling, very broad, strong O VI absorbers associated with SLFSs and pLLSs are rare, whereas they are common in LLSs, SLLSs, and even DLAs. Therefore, the contribution of broad, strong O VI associated with SLFSs and pLLSs to the cosmic baryon budget is very likely small.
The photoionized gas associated with pLLSs, LLSs, and SLLSs combined therefore contributes Ω g /Ω b ; 7% and the SLFSs contribute another 7% to the total baryon budget. 24 In contrast, the neutral gas that is mostly found in DLAs yields only Ω g /Ω b ; 2% at z ∼ 2-5 (e.g., Prochaska et al. 2005;Wolfe et al. 2005;Lehner et al. 2014;Péroux & Howk 2020). At z < 1, Lehner et al. (2019) find fractional contributions to the baryonic budget for the SLFSs, pLLSs, and LLSs that are about 5 times smaller. At z < 1, the SLFSs, pLLSs, and LLSs are probing much larger overdensities, and, at these redshifts, only ∼5% of the baryons are thought to be in the CGM of galaxies (Shull et al. 2012).

Cosmological Metal Budget
The metal density of the absorbers can be estimated with the baryon density via Ω m = Ω g Z, where Z = 10 [X/H] Z e is the metallicity of the gas in mass units and Z e = 0.0142 is the bulk solar metallicity in mass units from Asplund et al. (2009). For the metallicity, we use the linear mean metallicities in an about 0.5 dex bin of N H I from N log 14.5 H I = to N log 20 H I = . We note that the linear mean metallicities for SLFSs, pLLSs, and LLSs are 0.034, 0.032, and 0.022 Z e , respectively. We find Ω m = (15.9, 5.0, 3.6) × 10 −7 for the SLFSs, pLLSs, and LLSs, respectively. The latter value is consistent with that derived by Fumagalli et al. (2016b), Ω m = 5.1 × 10 −7 . Since we do not have SLLSs with N log 20 H I > in our sample, we use the result from the HD-LLS survey where Fumagalli et al. (2016b) derived Ω m = 1.6 × 10 −6 for the SLLSs. For reasons (sample variance, using the arithmetic mean) discussed in detail in Fumagalli et al. (2016b), these values are uncertain by a factor ∼3 (see also discussion in Berg et al. 2021). We note that the comoving metal-mass density can be estimated via ρ m = Ω m ρ c (e.g., Péroux & Howk 2020). Therefore, for the combined photoionized SLFSs, pLLSs, and LLSs, we find ρ m = 3.1 × 10 5 M e cMpc −3 , while combined with the SLLSs, this increases to 5.1 × 10 5 M e cMpc −3 .
The total expected cosmic metal density from the output of stars at z ∼ 2.8 is about 9 10 m exp 6  W´or, expressed in terms of cosmic metal mass density, 10 m exp 6  r M e cMpc −3 (Péroux & Howk 2020). Péroux & Howk (2020) also recently estimated for the DLAs, i.e., for the neutral gas, that Ω m ; 3.5 × 10 −6 and for the dust Ω m ; 1 × 10 −6 at 2.2  z  3.6, i.e., about 40% of the metals are in the neutral gas and about 10% are in the form of dust at 2.2  z  3.6. The SLLSs account for 17.7%, while the SLFSs, pLLSs, and LLSs account for 17.7%, 5.6%, and 4.0%, respectively, and all combined they account for about 45% of the metals at 2.2  z  3.6. We note that the recent estimates from Péroux & Howk (2020) for the SLLSs and LLSs are about within 1σ errors from those derived here and Fumagalli et al. (2016b). As for the baryons, another contribution is from the hot O VI gas associated with LLSs, SLLSs, and DLAs that can contribute another 7% (Lehner et al. 2014).
While most of the baryons at 2.2  z  3.6 are in the LYAF, this is not the case for the metals since the DLAs, dust, SLLSs, LLSs, pLLSs, and SLFSs contribute to about  95% of the cosmic metal budget. Direct estimates of the metal budgets in the LYAF are uncertain owing to ionization correction and the fact that some of the metals are not solely associated with N log 14.5 H I  absorbers. Using the recent results on Ω C IV from D' Odorico et al. (2013Odorico et al. ( , 2010; but see also, e.g., Schaye et al. 2003;Simcoe et al. 2004;Songaila 2005) and correcting their estimate from C IV to carbon and carbon to the total metal content yields Ω m  3.4 × 10 −7 at z ∼ 2.8 for the "IGM" or  4% of the total cosmic metal budget (although we note that this may include contributions from higher H I column density absorbers than LYAF absorbers). Overall the photoionized universe probed by absorbers N log 20 H I < contains about 50% of the metals at 2.2  z  3.6. Combining these results with those from Péroux & Howk (2020), the majority of metals are in the neutral and ionized gas at z ∼ 2.8.
At z < 1, Ω m = (4.6, 6.5, 25.9) × 10 −7 for the SLFSs, pLLSs, and LLSs, respectively . 25 In terms of comoving mass density, the combined SLFSs, pLLSs, and LLSs at z < 1 have ρ m = 4.7 × 10 5 M e cMpc −3 . Therefore, the amount of metals in SLFSs, pLLSs, and LLSs has increased by a factor of 1.5 from z ; 2.2-3.6 to z  1. However, the total expected cosmic metal density at z ∼ 0.5 is about 6 10 m exp 5  W´-, corresponding to a cosmic metal mass density 8 10 m exp 6  r´M e cMpc −3 (Péroux & Howk 2020), a factor 7 larger than at z ∼ 2.8. This implies that the (photoionized) metals associated with the SLFSs, pLLSs, and LLSs account for about 6% of the cosmic metal budget at low redshift, a factor of 4.4 lower than the contribution of the same absorbers at 2.2  z  3.6. As shown by Péroux & Howk (2020), the majority of the metals at low redshift remain in galaxies (including stars), in their vicinity, and in the hot intracluster medium, demonstrating a major shift in the distribution of the metals from 2.3 to 9 Gyr after the big bang.

Motivation and Methodology
In this section, we compare the KODIAQ-Z metallicity measurements as a function of H I column density to the simulated CGM from the FOGGIE project (Peeples et al. 2019, and https://foggie.science/). This project is composed of cosmological zoom-in simulations of six Milky Way-like halos. While KODIAQ-Z selects for absorbers without knowing their galaxy associations, absorbers with the highest column densities probe overdensities that are consistent with being virialized (mini)halos of galaxies (see Sections 7.1 and 9.3). The FOGGIE simulations let us address whether or not the metallicity-H I column density distribution seen in KODIAQ-Z could be explained by only the CGM and nearby IGM of progenitors to Milky Way-mass galaxies.
Many groups, including FOGGIE, have recently shown that better spatial resolution in the simulated CGM leads to a natural development of small clouds Peeples et al. 2019;Rhodin et al. 2019;Suresh et al. 2019;, largely because once the cooling length of the gas is resolved, the gas can cool and condense into smaller structures (Hummels et al. 2019;Zheng et al. 2020). Importantly for comparisons to observables such as metallicity, improved spatial resolution in the low-density gas also reduces the artificial mixing of metals, thereby allowing gas to remain at "high" and "low" metallicities even in the CGM, that is, while the median metallicities do not change much, the scatter about that median increases with improved resolution (Corlies et al. 2020). In the FOGGIE cosmological zoom-in simulations analyzed here the CGM resolution at z = 2 is 91 proper pc in the cooler gas, with a coarsest resolution of 365 proper pc. (The resolutions are correspondingly a factor of 0.85 and 0.75 smaller at z = 2.5 and z = 3, due to the spatial resolution being at a fixed comoving size.) This spatial resolution profile both enables the development of small, relatively dense and cool cloudlets and allows for the hotter, more rarefied, generally more metal-rich gas to be well resolved. Our analysis uses the z = 2, 2.5, and 3 outputs for each of the six FOGGIE halos. As presented in Simons et al. (2020), these halos are selected to be "Milky Way-like" in that they are roughly ∼ 10 12 M e at z = 0 and have no mergers of mass ratios more than 1:10 after z = 2 (see also Zheng et al. 2020 andLochhaas et al. 2021 for the full histories and z = 0 properties of the Tempest, Squall, and Maelstrom halos). While the KODIAQ-Z absorbers are not selected to be associated with Milky Way-like halos, galaxy halos of this size account for most of the stellar mass in the universe and therefore most of the chemical enrichment.
We extracted H I absorbers from each sight line using the Synthetic Absorption Line Surveyor Application (SALSA; Boyd et al. 2020) and its Simple Procedure for Iterative Cloud Extraction (SPICE) method. Across all six FOGGIE galaxies, we found 8825 H I absorbers at z = 3, 8207 at z = 2.5, and 7694 at z = 2, for a total of 24,726 H I absorbers. SPICE identifies absorbers based on H I number density, without the need for synthetic spectra. First, a number density cutoff is determined such that 80% of the total column density lies above the cutoff. 26 An interval is defined encompassing the sight line's cells that lie above this cutoff. This interval is then masked, and a new number density cutoff is found. The intervals from these two generations are then combined if they overlap along the line of sight and their average line-of-sight velocities are within 10 km s −1 of each other. These steps are repeated until the total column density of the remaining, unmasked cells falls below N log 12 H I = . SALSA then calculates the total column density and column-densityweighted average metallicity of each identified absorber.
While KODIAQ-Z is an H I-selected survey, we also extracted C IV and Si IV absorbers from these same sight lines in order to assess how our sample of simulated absorbers would differ if we also required the presence of metal absorption. Carbon and silicon are not tracked by the FOGGIE simulations, so they are inferred from the metallicity field assuming solar abundances. Their ionization states are estimated using the Trident code (Hummels et al. 2017), based on a table of CLOUDY equilibrium models (Ferland et al. 2013) using the Haardt & Madau (2012) UV background with self-shielding (Emerick et al. 2019). We pair our H I and metal absorbers based on their position within the simulations and their velocities. To be paired, the two absorbers must be within the same sight line, within 100 pc of each other, and have lineof-sight velocities within 10 km s −1 of each other. Pairing absorbers in this way is robust to minor changes in distance and velocity. Of our 24,736 H I absorbers at redshifts 2, 2.5, and 3, 1189 have an associated C IV absorber and 358 are associated with an Si IV absorber.

Comparison of the Metallicities
The column density and average metallicity of our extracted FOGGIE H I absorbers are shown in Figure 29. The plot bounds and dotted and dashed lines are the same as in Figure 19. The left panel shows all H I absorbers in cyan, with the observational samples overplotted. Grayscale contours have been added to illustrate the structure of the scatter plot where the density of points is highest. For the rightmost panel, no metallicity limit or metal-line requirement has been imposed on the H I absorbers. The middle panel only shows H I absorbers that are close to a C IV or Si IV absorber in position and velocity space, which-due to column density limits on our selectionis analogous to requiring a C IV or Si IV detection with column densities > 10 12 cm −2 . This limits the metallicity of the H I absorbers to [X/H]  − 2.4, which is the limit defining the metal-poor regime for KODIAQ-Z. Imposing requirements on metal-line detection can therefore limit the observed metallicity range (reinforcing the idea that an H I selection is required to probe the entire metallicity distribution). It should be noted that the gradual increase in the absorber metallicity as N H I decreases seen in the middle panel is caused by the column density requirement of N log 12 X > placed on the simulated C IV and Si IV absorbers. The far right panel in this figure compares probability densities for the full, C IV-paired, and Si IV-paired H I absorber populations in FOGGIE.
There is a broad agreement between the FOGGIE H I absorbers and those seen by the KODIAQ-Z, HD-LLS, and R12-DLA surveys. Both the observed and simulated absorbers follow similar trends in metallicity with N H I . Table 9 gives the median, mean, standard deviation, and IQR range on the metallicity for different classes of H I absorber, which can be compared to the observed populations in Table 5 (although note that the observations probe a larger redshift range). Both the mean and median metallicities increase with N H I , with the spread tending to decrease as N H I increases. This is remarkable because neither EAGLE nor FIRE zoom-in simulations could reproduce the similar trend seen at low redshift Lehner et al. 2019). In Figure 30, we compare the metallicity CDFs for the various absorbers at 2  z  3 in the FOGGIE simulations and observations, further demonstrating the similar trends of increasing metallicity with increasing N H I . 27 However, this figure also shows that the metallicity of the FOGGIE absorbers is higher than that seen in the observations (even though their IQR or standard deviations are similar), especially for the LLSs and DLAs. 26 Adjusting this cutoff has a slight effect on the total number of absorbers but does not significantly impact the resulting distribution of properties. Using a 70% threshold results in a 6% drop in the number of H I absorbers, and a 90% threshold results in an 11% increase. 27 The FOGGIE H I absorbers can have metallicity as low as Z Z log 8  ( )= -(the FOGGIE metallicity floor). There are absorbers with Z Z log 5  ( )<classified as SLFSs and pLLSs. These absorbers are not excluded from the CDFs in Figure 30, so these CDFs do not start at zero.
The FOGGIE absorbers therefore mimic the general trend seen in the observations, where higher H I column density regimes are dominated by higher metallicities. In FOGGIE, however, the metallicity at which a particular H I regime is most prevalent is higher than seen in the observational data. FOGGIE's SLFSs have a higher metallicity than both the observed SLFSs and the observed LLSs. Moreover, the FOGGIE LLSs more closely match the CDF of the observed SLLSs. Similarly, the FOGGIE SLLS CDF is closer in metallicity to the observed DLAs. The FOGGIE DLAs have a substantially higher metallicity (the FOGGIE median metallicity is 0.5 dex higher, or a factor 3 higher, than in the R12-DLA sample).
The metallicities used for the FOGGIE absorbers are conceptually slightly different from the metallicities derived in KODIAQ-Z and other spectral surveys, in that the metallicities reported for FOGGIE are the column-density-weighted average of the gas cells identified as the absorber rather than the result of a fit to single-phase ionization models. However, Marra et al. (2021) show that these two approaches, when applied to the same physical system, do not lead to large discrepancies in the average metallicity of the gas. The precise metallicities of FOGGIE's absorbers likely depend strongly on its implementation of supernova feedback, and these would be expected to differ slightly (e.g., be generally lower or higher) for different implementations. Yet despite FOGGIE's overall higher metallicity, the general similarity between the observed and simulated absorbers in Figure 29 seems to be a good indication that the trends observed in the KODIAQ-Z, HD-LLS, and DLA surveys can in principle be explained by CGM and nearby IGM gas around typical galaxies at these redshifts.

Insights from FOGGIE Simulations
Given some similarities in the metallicity trends with N H I between our samples of observed and simulated absorbers, we can use the FOGGIE simulations to shed light on the origin(s) of these absorbers and their connection to galaxies. From the empirical results from KODIAQ-Z, the overdensities (at least Figure 29. Metallicities of the H I absorbers extracted from the FOGGIE cosmological zoom-in simulations (Peeples et al. 2019;Simons et al. 2020). On the left, all H I absorbers (cyan) are compared to the KODIAQ-Z, HD-LLS, and R12-DLA surveys (where we limit these surveys to 2 z 3 to match the FOGGIE redshift range). Contours have been added where the density of cyan scatter plot points is highest. In the middle panel, we only show H I absorbers that have been paired with a C IV (blue filled circles) or Si IV (orange open circles) absorber, analogous to requiring detection of C IV or Si IV with N log 12 X > . Probability densities of both the full FOGGIE H I absorber population (cyan) and the subset with paired metal absorbers (blue and orange) are shown on the right. The dashed and dotted lines are the solar and VMP metallicity levels, respectively.  for the pLLSs, LLSs, and SLLSs) imply that these absorbers should largely probe gas within virialized (mini)halos (see Section 7.1). The Keck Baryonic Structure Survey (KBSS) also shows that at z ∼ 2-3 there is a strong incidence of absorbers with N log 14.5 H I > with galaxies at transverse physical distance  300 kpc and velocity separation between the absorber and galaxy redshifts  300 km s −1 (Rudie et al. 2012). This connection is not observed for lower-N H I absorbers, also implying some connection between N log 14.5 H I > absorbers and galaxies. On the other hand, recent Very Large Telescope (VLT)/MUSE IFU observations of high-z LLSs also show that at least some of these absorbers may probe IGM filaments (Fumagalli et al. 2016a;Lofthouse et al. 2020). At z < 1, similar findings are found using Keck/KCWI and VLT/MUSE observations, where VMP pLLSs are found either in the CGM of galaxies or well outside (likely in the IGM), while metalenriched are always found in the CGM of galaxies (Berg et al. 2022). With the FOGGIE simulations, we can explore the origins of these absorbers without observational constraints due to, e.g., the brightness limit of the galaxies. It is possible, however, that biases due to the selection criteria of the FOGGIE galaxies as Milky Way-like analogs at z = 0 could bias the results in some way. Importantly, we can also explore whether and how the origins may depend on the metal enrichment of the absorbers.
Each of the central FOGGIE galaxies is surrounded by a number of satellites at z = 2-3. These satellite galaxies likely contribute to the enrichment seen in our absorbers, and it is worth determining whether a given absorber is more likely to be associated with such a satellite than with the central galaxy. Including only those satellites with at least 1% of the central galaxy's stellar mass, we identify which galaxy each absorber is physically closest to. In Figure 31, we compare the absorber metallicities to their distance from the central galaxy and the distance to their closest galaxies, be that a satellite or the central galaxy. Generally, distance to an absorber's closest galaxy increases with distance from the central galaxy, which is the behavior we expect for absorbers whose closest galaxy is the central galaxy. Metallicity also tends to decrease with increasing distance from the central galaxy. Yet there are important exceptions to these two trends: some absorbers are far from the central galaxy yet close to a satellite, and these absorbers span the range of extracted metallicities. This includes some absorbers distant from the central galaxy that have metallicity close to and above solar. Therefore, not all high-metallicity absorbers should be assumed to be associated with the central galaxy.
In Figures 32 and 33, we select FOGGIE absorbers that are inflowing toward or outflowing from the central galaxy. This separation is made based on the angle between an absorber's 3D velocity vector and its position vector, both of which are calculated in the rest frame of the relevant central FOGGIE galaxy. Outflowing gas must have an angle θ > 60°, while inflowing absorbers must have θ < 120°. This is equivalent to v r cos 0.5 ·ˆq = = . This cut prevents gas that is neither strongly inflowing nor outflowing from being selected. Figure 32 shows the metallicity and H I column densities of our inflowing (blue) and outflowing (red) FOGGIE absorbers. Figure 31. Metallicities of the FOGGIE H I absorbers vs. their distance from the central galaxy, colored by the distance to their closest galaxy. Possible "closest galaxies" include the central satellite and all satellites with more than 1% of the central's stellar mass. High-metallicity absorbers far from the central galaxy tend to be close to a satellite galaxy. Dashed and dotted lines are the same as in Figure 19. The pdf's are also included for these two quantities, with both metallicity pdf's repeated in the top and bottom panels for comparison. The inflowing absorbers dominate by number, but both categories have similar H I column distributions, with absorbers becoming more numerous at lower H I column densities. The outflowing absorbers skew toward the high-metallicity end of the distribution. This is not surprising, as in FOGGIE winds are launched by the same supernova feedback that also enriches the surrounding gas. There is a plateau in the outflowing metallicity distribution around Z Z 2 log 1  --  , which may be tracing older outflows that have slowed and mixed with the ambient CGM. These absorbers number significantly less than the inflowing absorbers at the same metallicity, however, so it would be difficult for observations to disentangle these less enriched outflows. The inflowing absorbers follow a more symmetric metallicity distribution, with a peak about 0.5 dex above the metalpoor limit for KODIAQ-Z; however, as stated above, the FOGGIE absorbers generally tend to have higher metallicities than seen in the random sight lines observed by KODIAQ-Z. The inflowing absorbers also extend to the lowest metallicities. Even though Figure 31 reveals that some absorbers may be closer to a satellite than to the central galaxy, the trends in Figure 32 are unchanged when the absorbers are limited to those whose closest galaxy is the central. The biggest difference is that there are no outflowing absorbers with Z Z log 2.4  < -. Although the outflowing absorbers dominate the high end of the metallicity pdf, they are far outnumbered by the inflowing absorbers even at the highest metallicities. This is largely due to the selection on H I and the underlying physics included in the FOGGIE simulations. There is very little cool gas that would form H I within the outflows, as these are driven only by thermal pressure. Simulations with additional nonthermal pressure sources such as magnetic fields and cosmic rays (Butsky & Quinn 2018;Ji et al. 2020) demonstrate an increase in H I at large radii within the CGM. It is therefore quite likely that the number of H I absorbers classified as either inflowing or outflowing following our definition would change if these physical processes were included in the FOGGIE simulations, though it is likely that the numbers would still favor inflowing gas.
With Figure 33, we look at the metallicity and velocity of the FOGGIE absorbers as a function of their distance from the central simulated galaxy. We also compare the absorbers to the underlying gas distributions. This distribution compiles the gas in all 18 FOGGIE outputs we sampled. Probability densities are shown for the metallicity and radial distributions, and the mean metallicity of the gas is shown with a black line. Colors indicate the mean radial velocity at that metallicity and radius; we note that, unlike with the absorbers, the gas radial velocity has no selection based on the angle toward or away from the galaxy. First, we consider outflowing material: outflowing absorbers with the highest metallicities tend to be at or within the virial radius of these central halos, which at z = 2-3 is about 30-50 kpc, and have the highest velocities. These velocities are generally higher than what is seen in the gas at similar radii and metallicities, and they appear consistent with fast-moving enriched material launched by stellar feedback. On the other hand, high-metallicity gas can be found even beyond the virial radius, implying that galactic outflows become more diffuse as they travel farther away from the galaxy. The outflowing absorbers with the lowest metallicity are found well outside the virial radius, and their velocities are comparable to the gas in this regime. This supports the interpretation that these absorbers come from previously ejected gas.
Inflowing absorbers can be found at all radii, with their metallicity generally increasing the closer they are to their central galaxy, similar to the outflowing gas. The velocity of these inflowing absorbers is also higher near the galaxy. The FOGGIE zoom-in simulations contain several satellite and premerger galaxies that produce their own enriched outflows, which explains why high-metallicity absorbers can be found at large distances from the central galaxy. These absorbers do not appear to reflect the fastest inflowing gas, which has very high metallicity and can again be found at all radii.
Together, Figures 32 and 33 imply that, given a set of H Iselected absorbers that are associated with a galaxy system, most of these absorbers are likely to be inflowing toward a central galaxy. The precise number of inflowing and outflowing absorbers predicted by simulation will depend on the included physics, especially nonthermal pressure sources, though inflowing absorbers will likely still dominate. This dominance is true across column densities and metallicities, as satellites and nearby companions enrich gas flowing toward the central galaxy. Whether inflowing or outflowing, more enriched absorbers are more likely to be located physically close to the central galaxy.

Bridging the LYAF and DLAs
The history of metallicity in the universe provides an important constraint on models of galaxy formation and evolution. The metals in the universe represent a fossil record of star formation. Characterizing how the metallicities change with H I column densities (and hence overdensities or densities) yields information on the transport of metals and efficiency of that transport from the densest to the most diffuse regions of the universe. At high redshift, much of the effort has focused on the metal enrichment of the most diffuse regions probed by the LYAF (e.g., Cowie et al. 1995;Songaila 1998;Ellison et al. 2000;Schaye et al. 2000Schaye et al. , 2003Aguirre et al. 2004;Simcoe et al. 2004;Simcoe 2011) and the densest regions of the universe probed by DLAS (e.g., Pettini et al. 1997Pettini et al. , 1999Prochaska 1999;Prochaska et al. 2003;Rafelski et al. 2012;Jorgenson et al. 2013). The HD-LLS survey (Prochaska et al. 2015;Fumagalli et al. 2016b; see also Berg et al. 2021 for the XQ-100 survey of SLLSs) has surveyed H I column density absorbers mostly in the range N 17.8 log 20.3 H I <  . With KODIAQ-Z (this work, and earlier KODIAQ surveys; Lehner et al. 2014Lehner et al. , 2016, we bridge the gap between the DLA/SLLS regimes and the LYAF. 28 In Section 6.2, we show that the scatter of absorber metallicities at [X/H] − 2.4 is quite similar for the DLAs and the N 14.6 log 20 H I   regime, i.e., metal-enriched gas not only is observed in the densest regions (in or near galaxies) but also has spread to the more diffuse regions down to overdensities δ  10 (see Figure 19). However, as N H I decreases below N log 14.4 H I  (diffuse IGM), most of the absorbers have metallicities [X/H]  − 1.6 (see Figure 5 in Simcoe et al. 2004), while gas with − 1.6  [X/H]  − 0.2 is commonly observed over the entire range N 14.6 log 22 H I   . Therefore, the diffuse IGM probed by LYAF absorbers has not been enriched at the same level as stronger H I absorbers.
In contrast, for metallicities [X/H]  − 2.4, blind surveys of DLAs do not reveal a significant population of VMP DLAs (e.g., Rafelski et al. 2012;Jorgenson et al. 2013). Yet VMP absorbers are regularly observed at N log 20 H I  , and their frequency increases with decreasing N H I (see Table 6 and Figure 19). These VMP absorbers are, of course, also observed in the LYAF. The median metallicity of the LYAF at z ∼ 2.5 is around [X/H] = − 2.8 (Simcoe et al. 2004;Simcoe 2011; see also, e.g., Ellison et al. 2000;Schaye et al. 2003;Aguirre et al. 2004), which is a factor 2.5 lower than the median metallicity of the SLFSs at 〈z〉 ; 2.8 (Table 4). However, due to sensitivity issues, the lowest metallicities found for LYAF absorbers are around [X/H] ∼ − 3.5, while our survey reveals a population of absorbers (most with N log 18 H I < ) with [X/H] < − 3.5 that is nearly completely metal-free (see Section 10.2; see also Fumagalli et al. 2011a;Lehner et al. 2016). Simcoe et al. (2004) estimated that about 70%-80% of the LYAF has been enriched to [X/H]  − 3, but 20%-30% might be chemically pristine gas at low densities, which is about the same as the fraction of SLFSs with [X/H] < − 3.
Therefore, gas in the intermediate overdensity regime between the DLAs and LYAF has metallicities that are found in both the most diffuse and densest regions of the universe.

Pristine Gas
From the KODIAQ-Z (2.2  z  3.6) survey and CCC at low redshift (z < 1), we find that there is an abundance of VMP absorbers in relatively high overdensity regions (except for the DLAs by definition). While the metallicities have increased by about a factor 8-10 from 2.2  z  3.6 to z < 1, the fractions of 28 We show in Section 6.3 that similar trends of the metallicities with N H I are observed at low and high redshifts in the range N 15 log 22 H I   . However, the current UV observations at low z do not have the sensitivity to probe much lower metallicities than [X/H] < − 1 in the LYAF at z < 1 (this would need to await a future ∼6 m UV space-based telescope), and we cannot yet assess the full range of metallicities for absorbers with N log 15 H I < at z < 1. Therefore, we keep our discussion focused on the high-redshift universe. VMP absorbers are strikingly similar at low and high redshifts in N log 19 H I < absorbers. On the other hand, owing to the overall increase of the metallicity from 2.2  z  3.6 to z < 1, it is not surprising that extremely metal-poor gas with [X/H] < − 3.5 or even < − 3.0 is scant at z < 1 . With KODIAQ-Z, we determine that the fractions of SLFSs and pLLSs with [X/H] < − 3.5 and < − 3.0 are about 3%-10% and 15%-25% (90% confidence intervals; see Table 6). Combining KODIAQ-Z and HD-LLS, similar numbers are derived for the LLSs. For the SLLSs, the fractions of absorbers with these metallicities are smaller by about a factor 2.
The pristine LLSs reported by Fumagalli et al. (2011a) at z ∼ 3.5 have [X/H] < − 3.8 and < − 4.2, but at the time of this discovery it was impossible to say how common this population was. With KODIAQ-Z, we can revisit this question at 2.2  z  3.6 (see also Lehner et al. 2016). The lowest metallicity where some metals are detected is around [X/H]  − 3.8 in KODIAQ-Z (see Figure 19). We use that value of [X/H] ; − 3.8 to separate pristine (no metal) from metal-enriched (even at very low levels) absorbers. Considering Figure 19, Simcoe et al. 2004). The fraction of pristine gas in overdensities δ  50-100 at 2.2  z  3.6 is therefore at the level of a few percent, while for smaller overdensities it is likely in the range 10%-20%.
The frequency of mock absorbers with [X/H] < − 3.8 in the FOGGIE simulations is very similar to the frequencies found in KODIAQ-Z over similar redshifts, despite that the FOGGIE absorbers have on average higher metallicities than observed. In contrast, cruder resolution simulations show a tendency to have overenriched pLLSs and LLSs (see discussion in Fumagalli et al. 2011a). In the EAGLE HiRes cosmological hydrodynamical simulations at z ∼ 2-4, Rahmati & Oppenheimer (2018) find a fraction of pristine gas with [X/H] < − 3.8 of about 10% in the pLLS and LLS regimes. However, contrary to FOGGIE and the observations reported here or at low redshift Lehner et al. 2019), the EAGLE simulations show little change of the metallicity with N H I (see also Lehner et al. 2019;Wotta et al. 2019); the FIRE simulations also show similar results at z < 1 (Hafen et al. 2017). Consistent with the lack of N H I dependence, these simulations find a significant population of [X/H] < − 3 DLAs that is not observed.
This highlights some issues in simulations that may be related to the implementation of feedback physics but is probably also related to insufficient resolution in the more diffuse gas. As discussed in Section 9.2, low-metallicity inflowing gas is able to naturally get much closer to the galaxy when the CGM structures are more resolved than in the standard-resolution simulations because the metals are not forced to overmix in the regions around galaxies (see also Hummels et al. 2019;Peeples et al. 2019;Suresh et al. 2019;. As discussed in Section 9.3 and shown in Figures 31-33, the pristine and the VMP absorbers predominantly probe inflowing cool gas in the context of these simulations (see also Hafen et al. 2017;Rahmati & Oppenheimer 2018;Suresh et al. 2019), consistent with the idea that these absorbers are ideal candidates for the cold-mode accretion (e.g., Kereš et al. 2005;Dekel et al. 2009;Fumagalli et al. 2011b).
We finally note that achieving high spatial resolution not only in the CGM but also in the IGM is critical to fully capture the various structures probed by strong H I absorbers surveyed in KODIAQ-Z or CCC at low redshift. In recent highresolution simulations of an intergalactic sheet between two massive halos at z ∼ 3-5, Mandelker et al. (2019Mandelker et al. ( , 2021 show that pristine LLSs with [X/H] < − 3 such as those discussed in this work can be found in the IGM well outside the virial radius of any galaxy. Increasing the resolution in the IGM results in more H I in smaller, dense clouds, which are essentially absent in lower resolution. In their highest-resolution simulations (which have not converged yet), Mandelker et al. (2021) find that the covering fraction of LLSs in the IGM with [X/ H] < − 3 is about 3%. Therefore, the IGM itself could also provide an alternative origin for the observed pristine strong H I absorbers. Observationally, galaxy counterpart searches for two pristine LLSs with [X/H] < − 3.3 at z > 2 show no associated galaxies for one and a group of galaxies for the other one (Fumagalli et al. 2016a). At z < 1, the recent BASIC galaxy survey of 19 pLLSs shows two populations for the VMP pLLSs, one associated with galaxy halos and another one well outside the virial radius of any galaxy, likely originating in the IGM (Berg et al. 2022). It is therefore likely based on these observations and high-resolution simulations that a combination of these dual origins (IGM and cold inflow) can most likely explain the pristine and VMP absorbers.

Extremely High Metallicities and Metal-rich Gas
The metal-enriched diffuse gas constitutes a record of the transport of metals from the densest to the most diffuse regions of the universe. In Section 10.1, we discuss that while the LYAF shows evidence for metal enrichment, it is not enriched to levels seen in higher H I column absorbers. Metallicities in the range − 1.6  [X/H]  − 0.2 are commonly observed over the entire range N 14.6 log 22 H I   , but not in the LYAF. Therefore, at 2.2  z  3.6, galaxies have not yet polluted the diffuse IGM with metals at the same level as denser regions of the universe.
In Lehner et al. (2016), we reported the discovery of one supersolar pLLS at z ; 2.5 with [X/H] ; + 0.2. 29 Besides its metallicity, this pLLS is unique on several other levels, including the detection of O I, its small physical size (0.35 pc), its relatively high density (n H ; 0.2 cm −3 ), and its multiphase nature, with C IV having a very different velocity profile compared to the low ions (see Figure 14 in Lehner et al. 2016). We argued in that paper that its high metallicity and multiphase nature strongly suggest that it directly probes an active outflow from a protogalaxy at z ; 2.5. At z ∼ 1.8, Prochaska et al. (2006) uncovered two supersolar SLLSs, suggesting that these metal-rich absorbers may represent a significant metal reservoir in the young universe.
However, the HD-LLS survey does not report any supersolar SLLSs (see Figure 19), and the XQ-100 survey only reports one supersolar SLLS at z ∼ 2.5 (see Figure 2 in Berg et al. 2021; we exclude poorly constrained upper limits on the metallicity that are above solar in that survey). Adding that supersolar pLLS to the KODIAQ-Z+HD-LLS survey would imply that only 1 out of 242 absorbers has a supersolar metallicity at 2.2  z  3.6, i.e., 0.41% 0.32% . Supersolar-metallicity absorbers are therefore very rare at 2.2  z  3.6 for any N H I absorber.
At low redshift, the situation is quite different, with at least 2%-6% of supersolar absorbers at z < 1 (see Figure 22 and Lehner et al. 2019). This number could increase by a factor ∼3 if a harder EUVB was used instead of the HM05 EUVB to model the absorbers (at low redshift the effect of the EUVB is more important; see Wotta et al. 2019 and Section 5). This increase of supersolar absorbers is consistent with the overall increase of the metallicities by a factor 8-10 from 2.2  z  3.6 to z < 1 (see Figure 21).
The larger fraction of metal-enriched and even supermetallicity absorbers at low redshift indicates that a much larger volume of the universe has been exposed to metal pollution than at 2.2  z  3.6. Ferrara et al. (2000) developed a model that predicts that at z < 1 essentially all absorbers should have associated metal absorption, and the spread in metallicity should be less than 1 dex. This is not the case. Even at z < 1, the enrichment of the most diffuse gas in the universe is still very inhomogeneous in view of the nearly 3 dex spread in metallicities at z < 1 for absorbers with N 15 log 19 H I   (see Figure 22 and next section).

Inhomogeneous Metal Mixing
There is ample evidence for inhomogeneous abundances in H I column density N 14.6 log 20 H I   gas where the metallicities range from pristine ([X/H] < − 4) to about [X/H] ; − 0.2 dex, a factor of 6000 spread from the lowest to highest metallicities. Similar results are found at z < 1, but a factor somewhat smaller with 1000 variation between the highest and lowest metallicities  and Figure 22). However, our results also show that there is inhomogeneous metal mixing in single halos as directly evidenced by the large variation in metallicity in closely redshift-separated absorbers (see Figure 25). Such large metallicity variation over small redshift separations was initially reported by Prochter et al. (2010) between an SLFS, a pLLS, and an LLS at z ∼ 3.5 separated in velocities along the same QSO from about 130 to 180 km s −1 , where the metallicity differences were a factor 3, 40, and 158. Now with a sample of 36 paired absorbers, we find similarly large metallicity variations between absorbers separated by less than Δv < 500 km s −1 along a given QSO sight line where about half of paired absorbers have metallicity variations of a factor 2-3 and for the other half of a factor >140. However, on average the absorbers separated by Δv < 150 km s −1 have a metallicity variation that is on average around a factor 5, substantially smaller than in paired absorbers with larger Δv, where on the average the difference in metallicities is about 10 times larger. This would be consistent with the smaller variations being related to the variation within a halo and the larger variations being related to gas in interhalos.
Thus, within the overdensities of a few to several hundreds, the gas probed by absorbers with N 14.6 log 20 (see Section 10.1). It is also much smaller for the LYAF (e.g., Simcoe et al. 2004), implying overall a more chemically homogeneous gas in the LYAF and DLA regimes. Therefore, in and quite near galaxies, in the gas traced by DLAs, the metal enrichment from supernovae is quite efficient, but beyond the immediate vicinity of the galaxies, the volume filling factor for metals must be low, which is consistent with models of metal ejection from supernovae that is not expected to be efficient and confined to small regions (e.g., Ferrara et al. 2000;Scannapieco 2005). at 2.2 z 3.6 comparable in size to our companion survey, CCC, at z < 1. The H I selection and the H I column density range both ensure that no bias is introduced in the metallicity distribution of these absorbers, i.e., we are sensitive to any absorbers with [X/H]  − 3.5. The H I selection also provides a clean separation of the absorbers based on H I column density, which is largely used to separate the LYAF, absorbers with N log 14.5 H I  (IGM), from the denser regions of the universe probed by stronger H I absorbers. By definition, the SLFSs all have overdensities δ > 3 in our sample (z = 3.6; the maximum redshift in our statistical sample corresponds to this δ value using the analytical expression in Schaye 2001a). In contrast, a metal selection using C IV or O VI can include a wide range of H I column density absorbers from the LYAF to the pLLS regime.

Summary and Concluding Remarks
In the KODIAQ-Z survey, we have 155 SLFSs, 24 pLLSs, 16 LLSs, and also 7 SLLSs, for a total of 202 absorbers. To increase the number of LLSs and SLLSs in our sample, we use the H I-selected absorbers from the HD-LLS survey that have an H I column density range  Fumagalli et al. 2016b). We also use the results from the survey from Rafelski et al. (2012) to study the metallicity changes with N H I in the DLA regime. For both of these surveys, we restrict them to absorbers at 2.2 z 3.6, unless otherwise stated.
For all the absorbers with N log 20 H I  in the KODIAQ-Z and HD-LLS, we derive the posterior pdf's of the metallicities and other physical quantities using a Bayesian formalism that employs an MCMC sampling of a grid of CLOUDY photoionization models (where we assume that photons from the HM05 EUVB provide the source of photoionization). This follows directly from the methodology used in HD-LLS (Fumagalli et al. 2016b) and CCC Lehner et al. 2019). For absorbers with less-than-ideal constraints (about half of the SLFS sample and a few pLLSs), we adopt the "low-resolution" method developed at low redshift by Wotta et al. (2016) and refined in Wotta et al. (2019) to estimate the metallicities. Using SLFSs and pLLSs with reliable constraints from the metal ions, we find that log U can be reasonably well modeled by a Gaussian for the SLFSs and pLLSs, which can then be used reliably as a prior in the Bayesian MCMC modeling. We explore the effects of changing the EUVB from HM05 and HM12 to estimate the metallicities (Figures 9, 10), and we find that the changes in the metallicities between these two EUVBs are negligible for the absorbers with N log 17.2 H I  . For the SLFSs and pLLSs, the effect is only an increase of +0.12 ± 0.15 dex. The effect from changing the EUVBs on the metallicities is therefore much smaller than at z < 1.
Our main findings on the properties of our statistical sample of H I-selected absorbers with N 14.6 log 20 H I   at 2.2 z 3.6 can be summarized as follows.
1. From the comparison of the absorption profiles and the width of the profiles, we conclude that singly, doubly, and triply ionized species (e.g., Si II, Si III, Si IV, C II, C IV) and H I often trace the same gas in absorbers with N log 19 H I < . For SLFSs and pLLSs, we also find that when O VI absorption is present, it has often a similar velocity structure to lower ions, i.e., the O VI absorption is commonly narrow in these absorbers. This contrasts from the O VI that is frequently strong and broad when detected in absorbers with N log 17.8 H I  . 2. From the PF of the H I transitions, we show that 90% of the components have 13.3 b 40 km s −1 with a mean 〈b 〉 = 27 ± 6 km s −1 . This implies a temperature of the gas of T < 4 × 10 4 K, consistent with the gas being primarily photoionized. 3. We find that a single-phase photoionization model is appropriate to match the column densities of the low ions to high ions (including O VI) for the majority of the SLFSs and pLLSs. For the LLSs and SLLSs, when O VI is detected, a single-phase photoionization model cannot commonly reproduce the observed O VI column density, implying that as N H I increases, the multiple gas-phase nature becomes more important. 4. In our ionization models, [C/α] is allowed to vary in the range [−1, + 1] to accommodate for nonsolar relative abundances between carbon and α-elements caused by nucleosynthesis effects. This approach allows us to have a better sense of the uncertainties in the metallicity caused by the variability of C/α. Overall, we find that [C/α] is in the range − 0.6  [C/α]  + 0.5 as observed in other environments (DLAs, stars, H II regions), but robustly studying the variation of [C/α] with [X/H] for these absorbers is hindered by ionization corrections that add enough noise in the [C/α] distribution to conceal any subtle changes between [C/α] with [X/H] that are expected to be at the level of a factor 2-3.

5.
The 155 H I-selected SLFSs ( N 14.6 log 16.2 H I <  ) probe a wide range of metallicities from [X/H] < − 4 to [X/H] ; − 0.2. The metallicity posterior pdf is negatively skewed, with a prominent tail extending well below [X/H] < − 3. It has a main peak around the median value [X/H] ; − 2.4 and another smaller peak around −0.6 dex solar (the dip around −1.1 is observed in both binned and unbinned data). 6. The 24 H I-selected pLLSs ( N 16.2 log 17.2 H I <  ) probe a range of metallicity from [X/H] < − 4.2 to about [X/H] ; − 1. The sample is still too small to robustly characterize the distribution, and the lack of pLLSs with [X/H] > − 1 could be mainly due to small number statistics. However, both the median (−2.1 dex) and IQR metallicities imply an overall increase in metallicities compared to the SLFSs despite the lack of [X/H] > − 1 pLLSs. 7. Our sample of H I-selected LLSs ( N 17.2 log 19 H I <  ) consists of 16 absorbers that we combine with the HD-LLS survey to reach a size sample of 62 LLSs. The full range and median of the LLS metallicity pdf are quite similar to those of the pLLSs. Combining the pLLSs +LLSs, the median metallicity is −2.2 and the IQR is only 1 dex, compared to −2.42 and 2 dex, respectively, for the SLFSs. The pLLSs+LLSs are more frequent in the metallicity range − 3.2  [X/H]  − 1.2 but less frequent at very low metallicity [X/H]  − 3.5 than the SLFSs, showing that there is a shift in the metallicity enrichment properties of these absorbers below and above N log 16.2 H I  (at 2.2 z 3.6, this corresponds to overdensities of δ ; 140-50). , the overall metallicity trend observed with N H I continues in these regimes: the median/mean metallicity increases and IQR decreases with increasing N H I . For the SLLSs, the median and IQR metallicities are −1.9 and 1.1 dex, while these are −1.4 and 0.65 dex for the DLAs. DLAs have therefore rarely VMP gas (i.e., [X/H] < − 2.4 gas), while at N log 20 H I  , VMP absorbers are not rare at 2.2 z 3.6. 9. The fractions of extremely metal-poor systems with [X/H] < − 3.5 for the SLFSs, pLLSs, and LLSs are about the same, around 3%-10% (90% confidence interval). For the SLLSs, it is somewhat smaller with a fraction in the range 0. absorbers, the latter being similar to the amount of pristine gas found in the diffuse IGM at similar redshifts. On the other hand, supersolar absorbers at any N H I at 2.2  z  3.6 are very uncommon (<1%). 11. Using paired absorbers with velocity separations of Δv  500 km s −1 along the same QSO sight lines, we find that there is a large scatter in the metallicity from about 0.2 to >2 dex. For about half of the paired absorbers, there is evidence for metallicity variations over Δv  500 km s −1 of a factor 2-3, while for the other half of a factor >140. The larger variations (on average a factor 50 difference in metallicities) are, however, observed more often in paired absorbers with Δv > 150 km s −1 than in smaller velocity separated absorbers (where the metallicity variation is on average a factor 5). It is therefore plausible that the smaller variations correspond to the change within a single galaxy halo, while the large variations correspond to differences between galaxy halos. Both the large metallicity range and metallicity variation between paired absorbers imply that the transport of metals from their formation sites (galaxies) into the CGM and the IGM is very inhomogeneous. 12. The photoionized gas associated with pLLSs, LLSs, and SLLSs contributes to the cosmic baryon budget Ω g /Ω b ; 7%, and the SLFSs contribute another 7%. From the first KODIAQ survey (Lehner et al. 2014), the O VI-bearing hot collisionally ionized gas associated with LLSs and SLLSs is likely to contribute about the same level. These are the second-largest contributors to the cosmic baryon budget at high redshift behind the LYAF. 13. We show that about 18%, 6%, and 4% of the metals ever produced at 2.2  z  3.6 are in photoionized gas associated with SLFSs, pLLSs, and LLSs, respectively. The SLLSs account for another 18%, and therefore about 45% of the metals at 2.2  z  3.6 in photoionized absorbers with N 14.6 log 20.3 H I <  . Combining the SLFSs, pLLSs, LLSs, and SLLSs, their comoving metal mass density of the photoionized gas probed by these absorbers is ρ m = 5.2 × 10 5 M e cMpc −3 . Another possible  5% of the cosmic metal budget may also be in the form of highly ionized metals in collisionally ionized gas with N log 17.8 H I  (Lehner et al. 2014). Therefore over 50% of the metals produced at 2.2  z  3.6 are in the diffuse ionized gas.
To study the cosmic evolution of these absorbers, we combine our results with the CCC survey that explores the properties of similar absorbers at z < 1 (Lehner et al. 2018Wotta et al. 2019). The cosmic evolution of absorbers from 2.2  z  3.6 to z  1 with N log 15 H I  is summarized as follows.
14. We find that from 2.2 z 3.6 to z < 1, there is an overall increase of the metallicity of the gas probed by SLFSs, pLLSs, LLSs, SLLSs, and DLAs by a factor ∼8. 15. While the metallicity threshold for the VMP absorbers increases by about 1 dex from 2.2 z 3.6 to z < 1, a similar fraction of about 50% of absorbers with N log 18 H I  are VMP at low and high z. 16. Although there is plenty of primitive gas at z  1, i.e., gas that has largely not been polluted, the fraction of pristine SLFSs, pLLSs, and LLSs with [X/H] < − 3 at z  1 is < 1% at the 90% confidence level over the range N 15 log 19 H I < < . In contrast, 10%-25% of similar absorbers at 2.2  z  3.6 have [X/H] < − 3, implying that although the transport of metals outside galaxies is still very inhomogeneous at z < 1, all regions with N log 15 H I  have been polluted to some level by z < 1. 17. The hydrogen column density (N H ) is a factor 10-15 smaller at z < 1 than at 2.2  z  3.6. The contributions to baryonic budget from the SLFSs, pLLSs, and LLSs are about 5 times smaller at z < 1 than at 2.2  z  3.6. 18. The photoionized metals associated with the SLFSs, pLLSs, and LLSs take about 6% of the cosmic metal budget at low redshift. This is about a factor 4.4 decrease compared to the contribution of the same absorbers at 2.2  z  3.6, i.e., at low redshift most of the metals are found in or near galaxies and in the intracluster medium (see Péroux & Howk 2020). This is a major change in the distribution of metals in the universe from z ∼ 2.8 to z ∼ 0.5.
We set our empirical results in the context of the cosmological zoom-in simulations using simulations from the FOGGIE project, and the main conclusions from that study are as follows.
19. Contrary to cruder resolution simulations (especially in the CGM), a striking feature of the FOGGIE cosmological zoom-in simulations is that the behavior of the metallicity as a function of N H I is broadly similar to the observed empirical relationship: as N H I increases, the overall metallicity increases and the dispersion of the metallicity decreases. However, as in other cosmological simulations, FOGGIE appears to have too many metals at any N H I (e.g., supersolar-metallicity gas is not uncommon in these simulations, but observationally it is), implying that metals are more homogeneously distributed than observed. 20. In the FOGGIE simulations, outflowing absorbers with the highest metallicities tend to be at or within the virial radius of these central halos and probe active or recent galaxy outflows, while outflowing absorbers with the lowest metallicity are found well outside the virial radius of the central galaxies and therefore the remnants of ejected gas. On the other hand, inflowing absorbers can be found at all radii, with their metallicity generally increasing the closer they are to their central galaxy. VMP absorbers with [X/H] < − 2.4 are excellent probes of inflowing gas in these simulations. We finally note that recent high-resolution simulations of the IGM show that some of the observed pristine absorbers could also directly originate in the IGM, i.e., in regions well outside the virial radii of galaxies.
We thank the referee for constructive comments that improved our manuscript. We are grateful to Anna Wright for providing the satellite catalogs used in Figure 31 based on the methods of Pontzen & Tremmel (2018). We are grateful to Saloni Deepak for identifying a few typographical errors in the cosmic baryon and metal budget sections prior to publication. The main support for this research was made by NASA through the Astrophysics Data Analysis Program (ADAP) grant NNX16AF52G. Additional support was provided by NSF grant award No. 1516777. Support for the development of the CLOUDY ionization models was provided by NASA through grants HST-AR-12854 and HST-AR-15634 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. C.C. and B.W.O. acknowledge support by NSF grants no. AST-1517908 and AST-1908109 and NASA ATP grants NNX15AP39G and 80NSSC18K1105. This project has received funding from the