# KODIAQ-Z: Metals and Baryons in the Cool Intergalactic and Circumgalactic
Gas at $2.2\lesssim z\lesssim 3.6$

Nicolas Lehner Department of Physics, University of Notre Dame, Notre Dame, IN
46556  Claire Kopenhafer Department of Physics and Astronomy, Department of
Computational Mathematics, Science, and Engineering, Michigan State
University, East Lansing, MI 48824  John M. O’Meara W.M. Keck Observatory
65-1120 Mamalahoa Highway Kamuela, HI 96743  J. Christopher Howk Department of
Physics, University of Notre Dame, Notre Dame, IN 46556  Michele Fumagalli
Dipartimento di Fisica G. Occhialini, Università degli Studi di Milano
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INAF - Osservatorio
Astronomico di Trieste, via G. B. Tiepolo 11, 34143 Trieste, Italy  J. Xavier
Prochaska UCO/Lick Observatory, Department of Astronomy & Astrophysics,
University of California Santa Cruz, 1156 High Street, Santa Cruz, CA 95064
Kavli Institute for the Physics and Mathematics of the Universe (WIP), 5-1-5
Kashiwanoha, Kashiwa, 277-8583, Japan  Ayan Acharyya Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218  Department of Physics &
Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD
21218  Research School of Astronomy and Astrophysics, Australian National
University, Weston Creek, ACT 2611, Australia  Brian W. O’Shea Department of
Computational Mathematics, Science, and Engineering, Department of Physics and
Astronomy, National Superconducting Cyclotron Laboratory, Michigan State
University, East Lansing, MI 48824  Molly S. Peeples Space Telescope Science
Institute, 3700 San Martin Drive, Baltimore, MD 21218  Department of Physics &
Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD
21218  Jason Tumlinson Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218  Department of Physics & Astronomy, Johns Hopkins
University, 3400 N. Charles Street, Baltimore, MD 21218  Cameron B. Hummels
Department of Astronomy, California Institute of Technology, Pasadena, CA
91125

###### Abstract

We present the KODIAQ-Z survey aimed to characterize the cool, photoionized
gas at $2.2\lesssim z\lesssim 3.6$ in 202 H I-selected absorbers with
$14.6\leq\log N_{\rm H\,I}<20$, i.e., the gaseous interface between galaxies
and the intergalactic medium (IGM). We find that the $14.6\leq\log N_{\rm
H\,I}<20$ gas at $2.2\lesssim z\lesssim 3.6$ can be metal-rich gas
($-1.6\lesssim{\rm[X/H]}\lesssim-0.2$) as seen in damped Ly$\alpha$ absorbers
(DLAs); it can also be very metal-poor (${\rm[X/H]}<-2.4$) or even pristine
gas (${\rm[X/H]}<-3.8$) not observed in DLAs, but commonly observed in the
IGM. For $16<\log N_{\rm H\,I}<20$ absorbers, the frequency of pristine
absorbers is about 1%–10%, while for $14.6\leq\log N_{\rm H\,I}\leq 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 $\Delta
v<500$ ${\rm km\,s}^{-1}$ imply that the metals are poorly mixed in
$14.6\leq\log N_{\rm H\,I}<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\lesssim z\lesssim 3.6$. We find the mean metallicity increases with
$N_{\rm 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 $\sim$8 from
$2.2\lesssim z\lesssim 3.6$ to $z<1$, but the contribution of the
$14.6\leq\log N_{\rm H\,I}<19$ absorbers to the total metal budget of the
universe at $z<1$ is half that at $2.2\lesssim z\lesssim 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 ${\rm[X/H]}$ with $N_{\rm H\,I}$ similar to
our observational results. Very metal-poor absorbers with ${\rm[X/H]}<-2.4$ at
$z\sim 2$–3 in these simulations are excellent tracers of inflows, while
higher metallicity absorbers are a mixture of inflows and outflows.

Circumgalactic medium (1879) — Damped Lyman-alpha systems — Intergalactic
medium (813) — Lyman limit systems (981) — Metallicity (1031) — Quasar
absorption line spectroscopy (1317)

††software: Astropy (Astropy Collaboration et al., 2018), emcee (Foreman-
Mackey et al., 2013), Matplotlib (Hunter, 2007), PyIGM (Prochaska et al.,
2017a)††facilities: Keck(HIRES)

##  1 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\sim 2$) are therefore of prime importance to
understanding their role in the evolution of galaxies. 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 unprecedented high spatial resolution,
even in the more diffuse regions of the simulation boxes (e.g., van de Voort
et al. 2019; Hummels 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
inter-relationships with galaxies.

Over the last few years, our group has engaged in several H I-selected surveys
at both low and high redshifts to target absorbers in the H I column density
range $15\lesssim\log N_{\rm H\,I}\lesssim 19$. This column density interval
spans a range of physical environments from the diffuse IGM (Ly$\alpha$
forest, hereafter LYAF) to the denser portions of the CGM and the edges of
galaxy disks. Following Lehner et al. (2018), we describe absorbers as strong
Ly$\alpha$ forest systems (SLFSs, $\log N_{\rm H\,I}=14.5$–16.2), as partial
Lyman systems (pLLSs, $\log N_{\rm H\,I}=16.2$–17.2), or as Lyman limit
systems (LLSs, $\log N_{\rm H\,I}=17.2$–19). Below this range is the LYAF
($\log N_{\rm H\,I}\lesssim 14.5$) and above it are the super-LLSs (SLLSs,
$19.0-20.3$, a.k.a. the sub-damped Ly$\alpha$ absorbers) and the damped
Ly$\alpha$ absorbers (DLAs, $\geq 20.3$). With overdensities between the
diffuse IGM (probed by the LYAF) and galaxies (probed by DLAs), the SLFSs,
pLLSs, and LLSs spanning $15\lesssim\log N_{\rm H\,I}\lesssim 19$ 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 $14.6\lesssim\log N_{\rm
H\,I}\lesssim 20$ over cosmic time (e.g., the COS CGM compedium, CCC, at low
redshift—Lehner et al. 2018, 2019; Wotta et al. 2019, and see also Lehner et
al. 2013; Wotta et al. 2016, the HD-LLS survey—Prochaska et al. 2015;
Fumagalli et al. 2016b, the KODIAQ-Z survey—Lehner 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 and on the level of metal mixing in the
diffuse regions of the universe. The H I 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. 2009, 2013; Prochaska et al. 2017b; M.
Berg et al. 2021, in prep.), 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 ${\rm[X/H]}\lesssim-1.4$ at $z\lesssim 1$, M. Berg et
al. 2021, in prep., and ${\rm[X/H]}\lesssim-3$ at $2\lesssim z\lesssim 3.5$,
Fumagalli et al. 2011a, 2016a; Crighton et al. 2016; Lofthouse et al.
2020)111For the metallicity, we use conventional squared-bracket notation
${\rm[X/H]}\equiv\log N_{\rm X}/N_{\rm H}-\log({\rm X/H})_{\sun}$, where X is
an $\alpha$-element, unless otherwise stated.. There are, however, also
counter-examples at both low and high redshifts where very low metallicity
absorbers are found in vicinity of galaxies (e.g., Ribaudo et al. 2011;
Fumagalli et al. 2016a). Fully characterizing the metallicity distributions of
the SLFSs, pLLSs, and LLSs at low and high $z$ directly probes the metallicity
enrichment (or lack of it) 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).

In previous papers, our group presented a sample of 261 $z<1$ absorbers with
$15\lesssim\log N_{\rm H\,I}\lesssim 19$ (Lehner et al. 2018, 2019; Wotta et
al. 2019, hereafter, ; ; ). This sample includes an unexpectedly large
fraction of metal poor absorbers, many of which have metallicities
${\rm[X/H]}<-2$ implying little chemical enrichment since $z\sim 2$–3.

The finding that metal-poor absorbers with ${\rm[X/H]}\leq-1.4$ represents
about half of the population of $z\lesssim 1$ absorbers with $15\lesssim\log
N_{\rm H\,I}\lesssim 19$ was unanticipated. This is because neither the LYAF
(owing, 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, and see also Tripp
et al. 2005 for a rare example of a primitive SLLS at $z\sim 0.06$). CCC also
reveals that at $z<1$ there is no evidence of pristine gas, implying $>99\%$
gas with $15\lesssim\log N_{\rm H\,I}\lesssim 19$ has been polluted by
$z\lesssim 1$, even if it is only at very low level. Low-redshift photoionized
gas with $15\lesssim\log N_{\rm H\,I}\lesssim 19$ has a wide range of
metallicities observed ($-3\lesssim{\rm[X/H]}\lesssim+0.4$) and can show large
metallicity variations (factor up to 25) over small redshift/velocity
separation ($\Delta v<300$ ${\rm km\,s}^{-1}$), implying these absorbers
sample a highly inhomogeneous medium probed by these absorbers. The CCC
results beg the questions: are these properties unique to $z\lesssim 1$? how
do they evolve from earlier times, across the peak of cosmic star formation at
$z\sim 2.5$.

To address these questions, we present here the KODIAQ-Z survey (Keck Database
of Ionized Absorbers toward Quasars), a survey designed to identify H I
absorbers with $\log N_{\rm H\,I}\gtrsim 14.5$ at $z>2$ in Keck HIRES spectra
available in the KODIAQ database (O’Meara et al., 2015, 2017) and determine
their metallicities. Owing to the higher sensitivity of the Keck HIRES data
than the Hubble Space Telescope (HST) COS Origins Spectrograph data, we can
probe low $N_{\rm H\,I}$ absorbers at high redshift and still be sensitive to
low metallicity 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 $\log N_{\rm
H\,I}\lesssim 14.5$ (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\lesssim z\lesssim 3.6$, totaling 202 absorbers, 196 of them with
$14.5\lesssim\log N_{\rm H\,I}\lesssim 19$. 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 (Lehner et al. 2016).
However, in this paper, we undertake a systematic search and characterization
of H I-selected SLFSs.222Simcoe et al. (2004) surveyed the metallicity of the
IGM at $2.2\lesssim z\lesssim 2.8$. While they focused on absorbers with $\log
N_{\rm H\,I}<14.5$, they also included about 30–40 SLFSs as part of their
survey. 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 $14.5\lesssim\log
N_{\rm H\,I}\lesssim 20$. The HD-LLS survey is directly complementary to
KODIAQ-Z since most of the HD-LLS absorbers have $\log N_{\rm H\,I}\gtrsim
17.5$ while KODIAQ-Z largely samples gas with lower $N_{\rm 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_{\rm 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\gtrsim 2.2$–3.6. Our paper is organized as
follows. In §2, we describe the survey design of KODIAQ-Z, the search of the H
I-selected absorbers, and the sample of SLFSs, pLLSs, and LLSs. In §3, we
present the methods to derive the column densities of the metals and H I for
each absorber, while in §4 we present basic statistical properties of our
sample ($z$, $N_{\rm H\,I}$, and the Doppler parameter $b$). In §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 I-selected
absorbers are presented in §6. In §7, we derive the physical properties of the
absorbers, including their density ($n_{\rm H}$), total column density
($N_{\rm H}$), temperature, and linear scale ($l\equiv N_{\rm H}/n_{\rm H}$),
and we model the neutral fraction as a function of $N_{\rm H\,I}$ in the range
$14.6\leq\log N_{\rm H\,I}\leq 20$. In §8, we estimate the cosmic metal and
baryon budgets of these absorbers. In §9, we compare the KODIAQ-Z results to
the properties of the gas in cosmological zoom simulation 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 on the connection of the
$14.5\lesssim\log N_{\rm H\,I}\lesssim 20$ gas depending on its metallicity
with galaxies. In §10, we discuss the results, and in §11, we summarize our
main conclusions.

##  2 Survey Design, Database, and Observations

###  2.1 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\lesssim z_{\rm em}\ \lesssim 4.5$, limiting the
sample size to 235 QSOs. We also searched in the spectra of higher redshift
($z_{\rm 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\gtrsim 2$ is required to have wavelength coverage
of at least Ly$\alpha$ and Ly$\beta$.333The Keck HIRES spectra are normalized
and therefore we cannot rely on the break at the Lyman limit to estimate
$N_{\rm H\,I}$ for the pLLSs and LLSs, relying solely on the Lyman series, and
hence limiting the sample we can estimate accurately $N_{\rm H\,I}$. In each
of the QSO spectra, we search for absorbers with H I column density $\log
N_{\rm H\,I}\gtrsim 14.5$. This contrasts with our KODIAQ-Z I survey where our
search of H I absorbers was limited to absorbers with $\log N_{\rm
H\,I}\gtrsim 16.2$ (to match our initial lower $z$ survey, L13) and was by no
means a systematic and complete examination of the KODIAQ DR1 database
(O’Meara et al., 2015).

A major difference (besides the redshift probed) with the low redshift survey
CCC (CCC I) is that all the HIRES spectra are flux-normalized for the
coaddition of the individual exposures and echelle orders (O’Meara et al.,
2015, 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; specifically, the break at the Lyman limit for
absorbers with $\log N_{\rm H\,I}\gtrsim 16.5$ cannot be used. This limits the
sample size of absorbers with $\log N_{\rm H\,I}\gtrsim 17.2$ because it
requires enough weak Lyman transitions that are not too contaminated to derive
$N_{\rm H\,I}$ (see §3.2). For absorbers with $\log N_{\rm H\,I}\gtrsim 18$,
damping wings start to be observable in Ly$\alpha$ absorption, and therefore
this transition can be used to derive $N_{\rm 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 $\log N_{\rm H\,I}\gtrsim 15.3$ owing to the need to
be sensitive to metallicities ${\rm[X/H]}\lesssim-1$, the typically higher
signal-to-noise ratios (SNRs) and higher spectral resolution of the Keck HIRES
spectra allow us to probe lower $N_{\rm H\,I}$, down to the LYAF threshold of
about $\log N_{\rm H\,I}\simeq 14.5$, and still be sensitive to metallicity
with ${\rm[X/H]}\lesssim-2$. Although with very high SNR Keck HIRES spectra, a
strong limit (${\rm[X/H]}\lesssim-2$) on the metallicities can be placed down
to $\log N_{\rm H\,I}\simeq 13.6$ (e.g., Simcoe et al. 2004), our lower limit
on $N_{\rm H\,I}$ was selected to correspond to the limit between the diffuse
IGM and denser IGM/diffuse CGM. Following Schaye (2001a), a 14.5 dex limit on
the H I column density corresponds to an overdensity of $\delta\gtrsim 5$ at
$z\sim 2.8$ (about the average redshift of the absorbers in our sample, see
§4.1).

###  2.2 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$\alpha$) down to
917 Å. In Fig. 1, the colored profiles show single-component Voigt models for
absorbers with a Doppler parameter $b=25$ ${\rm km\,s}^{-1}$ and $\log N_{\rm
H\,I}=15,16,17,18$. That $b$-value is selected based on previous profile
fitting of strong H I absorbers at similar redshifts (e.g., KODIAQ-Z I); 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$\alpha$ must have flux $F_{\lambda}\leq 0$ $\gtrsim
0.4$ Å in contiguous pixels; 2) at the average velocity/redshift ($v_{c}$) of
the identified Ly$\alpha$, the flux in the Ly$\beta$ absorption must also
reach $F_{\lambda}\leq 0$; 3) at $v_{c}$, we also require there is substantial
absorption in at least Ly$\gamma$ (and Ly$\delta$ if there is wavelength
coverage of that transition). Additional conditions were added to ensure that
when a DLA or a 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$\alpha$, Ly$\beta$, and Ly$\gamma$. However, since the damping wings of
Ly$\alpha$ can be used to determine $N_{\rm H\,I}$ if $\log N_{\rm
H\,I}\gtrsim 18.5$, we also undertook a separate search that relies solely on
Ly$\alpha$ where the flux $f_{\lambda}\simeq 0$ over a wavelength width of
$\gtrsim 2$ Å. A posteriori, only a few of these very strong H I absorbers
ended up in our sample owing to the fact that the KODIAQ spectra are
normalized prior to coaddition, and the normalization was often not reliable
enough to model the damping wings for these absorbers, especially the strong
SLLSs and DLAs (DLAs being not an issue since they are not part of our
survey).

Figure 1: Example of an absorber identified by our automatic searching tool.
For each identified candidate-absorber, a similar figure is created that shows
the normalized spectra (in black) of the H I transitions from the 1215 to 917
Å transitions as function of the velocity in the redshift rest-frame. The
colored profiles are Voigt profiles for an absorber with $b=25$ ${\rm
km\,s}^{-1}$ and $\log N_{\rm H\,I}=15,16,17,18$ and are used to help assess
the level of contamination in each transition. Ultimately for this candidate
two absorbers were identified at $z=2.55125$ and $2.55216$ (corresponding to
$v\simeq-50$ and $+17$ ${\rm km\,s}^{-1}$ in this figure).

For each candidate absorber, our search tool created a stack of H I
transitions as shown on Fig. 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 to 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, for which we show their redshift
distribution in Fig. 2. In this figure, we also show the redshift distribution
of the QSOs. About 83% of the candidate absorbers are at $2.2\leq z_{\rm
abs}\leq 3.6$; candidate absorbers at $z_{\rm abs}\leq 2$ are detected solely
based on Ly$\alpha$.

Figure 2: Distribution of the redshifts for the QSOs and the candidate
absorbers from the automatic search in KODIAQ DR2.

To assemble our final sample for metallicity determinations, we considered
both the ability to estimate $N_{\rm 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 set of stack of velocity profiles with Ly$\alpha$,
Ly$\beta$, Ly$\delta$, 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 if there is enough metal-ion coverage
to determine the metallicity. After all the measurements are done (detailed in
§3), the final sample size sample consists of 202 absorbers, i.e., $\approx
7\%$ of the candidates from the initial automatic search were absorbers where
we could place a robust estimate on $N_{\rm H\,I}$ and metal-line column
densities (including strong upper limits). The final KODIAQ-Z sample includes
155 SLFSs, 24 pLLSs, 16 LLSs, and 7 SLLSs.

###  2.3 Comparison Samples

One of the goals of our survey is to determine how the metallicity changes
with $N_{\rm 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 $\log N_{\rm H\,I}\geq 19$ (see ; ;  and references therein).
At $z\gtrsim 2$, we use two surveys: 1) the HD-LLS survey (Prochaska et al.,
2015; Fumagalli et al., 2016b), a survey of absorbers with $17.3\leq\log
N_{\rm H\,I}<20.3$ at $z=2$–4.39, $\sim$95% of the absorbers having H I column
densities in the range $18\leq\log N_{\rm H\,I}<20.3$ (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 R12-DLA survey. 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, mostly dictated
our use of these surveys. 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\lesssim z\lesssim 3.6$ (see below). The recent
XQ-100 survey of SLLSs by Berg et al. (2021) shows also the same metallicity
distributions as in HD-LLS.

For the R12-DLA, we adopt the metallicities from the $\alpha$-elements (or Zn
if not available). For the HD-LLS, we re-estimate all the metallicities using
the same EUVB adopted in CCC in order to avoid mismatch in the comparisons of
the metallicities although as we will show in §5.3, the systematic error
arising from the EUVB in the photoionization modeling of the absorbers is much
smaller (in fact negligible) at $z>2$ than at lower redshift, especially for
the strong LLSs and SLLSs that are part of the HD-LLS survey.

##  3 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 profile
fitting–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 the account
that the contamination, blending, and saturation effects are more severe than
for the metal lines (see Fig. 1).

###  3.1 Metal Column Densities

Although we only estimate the column densities of the metal species once the H
I transitions were modeled (and hence the candidate absorber became 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 unit velocity,
$\tau_{a}(v)=\ln[F_{\rm c}(v)/F_{\rm obs}(v)]$, where $F_{c}(v)$ and $F_{\rm
obs}(v)$ are the modeled continuum and observed fluxes as a function of
velocity. The AOD, $\tau_{a}(v)$, is related to the apparent column density
per unit velocity, $N_{a}(v)$, through the relation $N_{a}(v)=3.768\times
10^{14}\tau_{a}(v)/(f\lambda(\mbox{\AA})$) ${\rm cm}^{-2}\,({\rm
km\,s^{-1}})^{-1}$, where $f$ is the oscillator strength of the transition and
$\lambda$ is the wavelength in Å. The total column density is obtained by
integrating the profile over the pre-defined velocity interval,
$N=\int_{v_{\rm min}}^{v_{\rm max}}N_{a}(v){\rm d}v$, where $[v_{\rm
min},v_{\rm max}]$ are the boundaries of the absorption. We computed the
average line centroids through the first moment of the AOD $v_{a}=\int
v\tau_{a}(v){\rm d}v/\int\tau_{a}(v){\rm d}v$ ${\rm 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_{\rm
min},v_{\rm max}]$ corresponds to a single absorbing H I complex as shown,
e.g., in Fig. 3, but this does not mean necessarily 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 difference between central velocities between two H I
components are at least separated by $|\Delta v|\gtrsim 40$ ${\rm km\,s}^{-1}$
and $\log N_{\rm H\,I}\lesssim 18$. As we will see below the mean Doppler
parameter of the H I components is around 27 ${\rm km\,s}^{-1}$, corresponding
to a temperature of $T<4\times 10^{4}$ K, typical for photoionized gas. This
means 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 $\sim 50$ ${\rm
km\,s}^{-1}$. In that case, we integrate the metal ions over a similar
velocity width as the H I when that information is available from weak H I
transitions; if the latter is not available (e.g., SLLS where weak H I
transitions are not covered, which is the case for $\lesssim 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 this 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 Figs. 1 and 3, even deep in the LYAF, the original continuum
placement (see O’Meara et al. 2015, 2017) is quite good and can be directly
used in many cases to estimate the physical parameters of the absorption lines
accurately. We therefore took full advantage that the Keck HIRES spectra are
already normalized to use the AOD method in a fully automated way. For each
identified absorber, we made typically two iterations to estimate the physical
parameters of the lines ($v_{a}$, $N_{a}$, and equivalent width
$W_{\lambda}$). In the first iteration, the velocity range is estimated based
on the stack of velocities profiles used to preview the data. Our program
created figures that show for each transition the normalized and apparent
column density profiles (such as Fig. 3 in CCC I), allowing us to refine the
integration velocity range and determine if the continuum needed any local
refinement. In the second iteration, the refined velocity range is applied to
the original or re-normalized spectra. If we determine that the continuum
model needed to be re-estimated, we either use the automated Legendre
polynomial fitting described in CCC I or selected by hand the continuum region
to be fitted by low-order Legendre polynomials.

We integrate the equivalent widths using the same $[v_{\rm min},v_{\rm max}]$
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 if the
absorption is detected at the $\geq 2\sigma$ level. If it is not, we quote a
2$\sigma$ upper limit on the column density, which is simply defined as twice
the 1$\sigma$ error derived for the column density assuming the absorption
line lies on the linear part of the curve of growth. The 1$\sigma$ error is
determined by integrating the spectrum over a similar velocity interval to
that of a detected ion or very weak H I transition or otherwise $\pm 20$ ${\rm
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 $\geq 2\sigma$ significance, the
next step is to assess if 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 Fig. 3 and Appendix). Other ions are also considered such as
Fe II, Fe III, Al II, Al III or atoms such as O I, but in most cases these are
not detected (Fig. 3 and Appendix References). For C IV $\lambda\lambda$1548,
1550, Si IV $\lambda\lambda$1393, 1402, and O VI $\lambda\lambda$1031, 1037,
Si II $\lambda\lambda$1193, 1260, 1304, 1526, these are either doublets or
ions with several transitions, and both contamination and saturation can be
readily checked by comparing their $N_{a}(v)$ profiles. Ions with
$\lambda\lesssim 1215$ Å are more often likely to be contaminated than ions
with transitions at longer wavelength owing to the dense LYAF at these
redshifts. For example, for the absorber shown in Fig. 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$ ${\rm km\,s}^{-1}$ seen in Fig. 3 (that is
not analyzed as part of this work because its H I column density is $\log
N_{\rm H\,I}<14.5$) 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 $\lambda$1334 is detected; more often
than not it is contaminated by the LYAF.

On the other hand, ions such as C III $\lambda$977 and Si III $\lambda$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,
comparing both the central velocities of the peak optical depth and shape of
the velocity structure of the absorption profiles across transitions from the
same or similar ions.444Of course if the investigated ion is not in the same
gas-phase than an ion that is uncontaminated (e.g., O VI in some cases), these
criteria cannot be used. 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 than 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 changed slightly 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 modeling555This upper limit can be differentiated from a non-
detection upper limit as the error bars of the former are not the standard
non-detection upper limit errors $+0.18,-0.30$ dex..

For most of the metal ions in our survey, saturation is not a major issue
owing to the lower column H I densities than those of DLAs and even SLLSs
(there are only 7 SLLSs in our sample and no DLA), lower metallicities on
average than $z<1$, and to the high spectral resolution of the Keck data. To
determine if there is saturation, we proceed in a similar way that 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 $\Delta\log(\lambda f)\simeq 0.3$
and the difference in apparent column densities of the weak and strong
transitions is $\Delta(\log N_{a})\equiv\log N^{\rm weak}_{a}-\log N^{\rm
strong}_{a}\leq 0.13$ dex (see W16, ), we are able to correct for the mild
saturation using the method described in Savage & Sembach (1991) and W16. 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.666Contamination and saturation can
be a priori both present at the same time and deceitful. However, it is
unlikely to occur frequently over the same pixels, and therefore each can be
disentangled reliably in most cases.

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, averages velocities, and
redshifts for the metal ions and H I (see §3.2) are summarized in Table 2. We
also provide all the results in a machine readable format (see Appendix).

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 range of the profiles are shown in red).
Another absorber separated by 65 ${\rm 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 Fig. 4)

###  3.2 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 non-detections). 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, to the entire
Lyman series from 1215 Å all the way down to 913 Å (depending on wavelength
coverage, redshift of the absorber, presence of a SLLS/DLA along the
sightline, and strength of the absorption), spanning a difference in
$f\lambda$ of over 2,700. However, the entire Lyman series is in the LYAF (and
Ly$\beta$ forest), and therefore many transitions can be and are often heavily
contaminated, requiring a close look at all the H I absorption features to
assess the level of contamination. Therefore our first task was to
systematically look simultaneously at several H I profiles (from 1215 to 915
Å) to assess if there is enough information (i.e., portion of the profiles of
H I transitions that are uncontaminated) to derive $N_{\rm H\,I}$. If it was
judged that the contamination could be dealt with, we proceeded 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).

Owing to the LYAF contamination, each transition was manually handled, and we
often also reassessed the continuum that can be more uncertain than in the
redder part by modeling it with low-order Legendre polynomials. As mentioned
above, we treated as much as possible an absorber as a single H I component.
For the example shown in Fig. 3, there are two components separated by 65
${\rm km\,s}^{-1}$, which correspond to two absorbers at $z_{\rm abs}=2.55125$
and 2.55216. We therefore integrated the AOD profiles over each component (for
the absorber at $z_{\rm abs}=2.55125$, that corresponds to the velocity
interval of about $[-30,40]$ ${\rm km\,s}^{-1}$). For absorbers with $\log
N_{\rm H\,I}\gtrsim 18$, this was not always possible especially if there is
not coverage of the weaker Lyman series transitions (for these absorbers, the
damping wings of Ly$\alpha$ and Ly$\beta$ can be used to determine $N_{\rm
H\,I}$, see below). However, these absorbers are rare in our sample, $\lesssim
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_{\rm 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 $\chi^{2}$
minimization and $\chi^{2}$ error derivation approach outlined by Sembach &
Savage (1992). A single component COG is assumed. Our procedure solves for
$\log N$ and $b$ independently in estimating the errors. We started with all
the H I transitions for which we estimated the equivalent widths and did not
appear a priori to be contaminated based on the AOD analysis. We checked the
results from the COG model and then removed only the transitions clearly
departing from the model. We ran our COG program again with the new set of
transitions to obtain the final estimate of $N_{\rm H\,I}$ with the COG
method.

Figure 4: Example of the same absorber shown in Fig. 3 where we show the Voigt
profile fit (in red the composite profile fit, and in blue each individual
component) to the individual transitions of H I (black spectra). Despite some
contamination in each transition, there is enough uncontaminated regions in
each profile to estimate robustly the H I column density in each component.

Using the normalized and masked H I spectra, we also Voigt-profile fit the H I
transitions using an in-house PF program described in Lehner et al. (2011,
2014, 2018), which was updated and refined from Fitzpatrick & Spitzer (1997).
A major difference from the previous methods is that the model profiles were
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 more)—column density ($N_{i}$), Doppler
parameter ($b_{i}$), and central velocity ($v_{i}$)—are input as initial
guesses and were subsequently varied to minimize $\chi^{2}$ of the fit to all
the fitted H I transitions. As for the COG, we started with all the H I
transitions that did not appear a priori contaminated and iterated, removing
any transitions that a posteriori appear contaminated (or sometimes revisiting
the masked pixels). In Fig. 4, we show an example of profile fitting that
consists of a two-component fit corresponding to the absorbers at $z_{\rm
abs}=2.55125$ and 2.55216. As it can be seen in each panel of Fig. 4, there is
some evidence of contamination especially in the component at 65 ${\rm
km\,s}^{-1}$ ($z_{\rm abs}=2.55216$), but even in that component there is
enough information to derive robustly $N_{\rm H\,I}$ when using the different
H I transitions simultaneously.

In general, there is a good agreement between these different methods to
estimate $N_{\rm H\,I}$ for absorbers with $14.5\lesssim\log N_{\rm
H\,I}\lesssim 17.2$. The three methods explore different parameters and
different transitions, allowing us to robustly estimate $N_{\rm H\,I}$. When
multiple approaches are valid, we simply average the values from the different
methods and propagate the errors accordingly. For systems with $\log N_{\rm
H\,I}\gtrsim 17.2$, the PF results are most of the times solely adopted as a
consequence of the need to either fit the Ly$\alpha$ damping wings for
absorbers with $\log N_{\rm H\,I}\gtrsim 19$ or combination of weaker and
stronger H I transitions in order to correctly model the absorption that
depends more strongly on both $N$ and and $b$ in this regime. The eight
systems with $17.5\leq\log N_{\rm H\,I}\leq 18.65$ 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$\alpha$. These were all
analyzed as part of Lehner et al. (2014) and KODIAQ-Z I, and we refer the
reader to these papers for a full description, noting, very conservative
errors were adopted with a minimum error on the $N_{\rm H\,I}$ for any LLS of
$\sigma_{\log N}=0.15$ using this methodology. However, the metal lines
associated with these absorbers were all reanalyzed as part of this work.

The adopted $N_{\rm 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, and 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.

##  4 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_{\rm 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_{\rm H\,I}$ but at lower
redshift). As we show below the $N_{\rm H\,I}$ and Doppler parameter
distribution also motivate the use of photoionization in the ionization
modeling to determine the metallicity of the absorbers.

###  4.1 H I Column Density and Redshift Distributions

Figure 5: Distribution of the redshifts and H I column densities for the
absorbers in the KODIAQ-Z sample. The top and right panels show the redshift
and $\log N_{\rm H\,I}$ distributions.

In Fig. 5, we show the distributions of the absorber redshifts and H I column
densities as well as the scatter plot between these two quantities. There are
202 absorbers in KODIAQ-Z and their full $N_{\rm H\,I}$ and $z$ intervals are
$14.43\leq\log N_{\rm H\,I}\leq 19.87$ and $2.23\leq z_{\rm abs}\leq 4.35$.
However, 99% of the absorbers have $14.57\leq\log N_{\rm H\,I}\leq 19.87$ and
$2.23\leq z_{\rm abs}\leq 3.57$ with means and standard deviations
$\langle\log N_{\rm H\,I}\rangle=15.8\pm 1.1$ (the median being $\log N_{\rm
H\,I}=15.5$) and $\langle z_{\rm abs}\rangle=2.79\pm 0.31$ (the median being
$z_{\rm abs}=2.75$), respectively. There is an evident peak in the redshift
distribution at $2.4\lesssim z_{\rm abs}\lesssim 2.6$ where there are about 65
absorbers; otherwise between $z_{\rm abs}\simeq 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_{\rm 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;
as $z$ increases above 3.6, the level of contamination from the LYAF increases
(and indeed at $3.5\lesssim z_{\rm abs}\lesssim 4.4$, only for 5 absorbers
were we able to robustly derive $N_{\rm 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\leq z\lesssim 3.6$
when comparing the absorbers as a function of $N_{\rm H\,I}$ over the same
redshift interval.

Owing to the selection of the absorbers, the $N_{\rm H\,I}$ distribution
should follow to some extent the column density distribution of H I (Prochaska
et al., 2014), which is why the $N_{\rm H\,I}$ distribution is not flat. From
approximately $10^{15}$ to $10^{16.5}$ cm-2, we see in Fig. 5 the expected
drop in the number of absorbers as $N_{\rm H\,I}$ increases. At $\log N_{\rm
H\,I}<15$, our sample is not complete as there is a sharp drop in the number
of absorbers. This is because as $N_{\rm 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. For absorbers with $16.5\lesssim\log N_{\rm
H\,I}\lesssim 17.5$, we have to depend more on weak Lyman series transitions
to estimate $N_{\rm H\,I}$, limiting the number of absorbers especially as $z$
increases beyond $z\gtrsim 3.2$. As described in §3.2, absorbers with
$17.5<\log N_{\rm H\,I}\lesssim 18.65$ are rare because they are essentially
on the flat part of the COG, dramatically restricting the number of absorbers
where we can estimate $N_{\rm H\,I}$. And indeed although absorbers with
$18.7\lesssim\log N_{\rm H\,I}\lesssim 20$ should be less frequent than the
previous category, there are more of these strong LLSs/SLLSs in our sample
owing to that we can use the damping of Ly$\alpha$ and Ly$\beta$ to derive
$N_{\rm H\,I}$ (see also §3.2).

###  4.2 H I Column Density and Doppler Parameter Distributions

Figure 6: Distribution of the Doppler parameters ($b$) and H I column
densities in the KODIAQ-Z sample obtained from the profile fitting of
individual components. The right and top panels show the $b$ and $\log N_{\rm
H\,I}$ distributions. The dotted line shows the lowest fitted $b$-value
(corresponding $T\lesssim 10^{4}$ K); 90% of the components have $13.3\leq
b\leq 40$ ${\rm km\,s}^{-1}$. Note that for the Keck HIRES spectra, the
instrumental $b$-value is about 4–5 ${\rm km\,s}^{-1}$, well below the minimum
$b$-value, i.e., narrower H I absorbers than $b<12$ ${\rm km\,s}^{-1}$ could
have been detected if present.

In Fig. 6, we compare the Doppler parameter, $b$, and $N_{\rm H\,I}$ derived
from the profile fitting in the individual components of H I. There is no
strong trend between $b$ and $N_{\rm H\,I}$, although broad components with
$b>40$ ${\rm km\,s}^{-1}$ are more frequent at $\log N_{\rm H\,I}<16$ ($b=40$
${\rm 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).777The Doppler parameter for H I is related to the gas temperature via
$b_{\rm H\,I}=(2kT/m_{\rm H}+b^{2}_{\rm nt})^{0.5}$, where $k$ is the
Boltzmann constant, $T$ the temperature of the gas, $m_{\rm H}$ is the mass of
hydrogen, and $b_{\rm nt}$ is the (unknown) non-thermal component to the
broadening. We, however, note that the majority of these broad components are
found in more complex velocity profiles (2 or more components) and/or lower
SNRs spectra, where a broad component reduces the $\chi^{2}$.

Considering the entire sample, the mean and standard deviation of $b$ is
$\langle b\rangle=30\pm 11$ ${\rm km\,s}^{-1}$ and median is 28 ${\rm
km\,s}^{-1}$. About 90% of the components have $13.3\leq b\leq 40$ ${\rm
km\,s}^{-1}$ where $\langle b\rangle=27\pm 6$ ${\rm km\,s}^{-1}$ and the
median being 27 ${\rm km\,s}^{-1}$. That latter value corresponds to a
temperature of the gas of $T<4\times 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$
${\rm km\,s}^{-1}$ and a mean of $28\pm 8$ ${\rm km\,s}^{-1}$ (CCC I and Fig.
7 therein). This contrasts from the weaker LYAF absorbers with $\log N_{\rm
H\,I}\lesssim 14$ where the median and mean $b$-values are systematically
larger by 15%–30% at $z\lesssim 0.5$ than at $1.5\lesssim z\lesssim 3.6$,
implying a larger fraction of the low-$z$ IGM is hotter and/or more
kinematically disturbed than the high-$z$ IGM (Lehner et al., 2007).

##  5 Determining the Metallicity of the Absorbers

###  5.1 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; Prochaska et al. 2017b; ; ; ; ;
Crighton et al. 2013; Fumagalli et al. 2016b; Cooper et al. 2021; Zahedy et
al. 2021). In §4.2, we show the KODIAQ-Z absorbers have temperatures on
average $T\lesssim 4\times 10^{4}$ K, consistent with the gas being
photoionized (the same result was found for the low redshift absorbers with
similar $N_{\rm H\,I}$, CCC I). As illustrated in the velocity profiles
provided in the appendix, 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 (although we emphasize again that sometimes metal ions owing to their
smaller $b$ can have more than 1 component within the breadth of the H I
profiles even though the H I absorption can be fitted with a single
component—that also implies that the thermal broadening is a large contributor
to the broadening).

In contrast to the low redshift, for the ionizing radiation field and typical
lower metallicities at $z\sim 2$–3 (about 0.1% solar or ${\rm[X/H]}=-2$, see
below, ; ), even strong transitions like C II $\lambda$1334 and Si II
$\lambda$1260 are often not detected, so we have to rely more on intermediate
(Si III, C III, Si IV) and high (C IV) ions to determine the metallicity.
Because the current survey probes lower $N_{\rm H\,I}$ absorbers than our
previous surveys at high redshift (Lehner et al. 2014; KODIAQ-Z I), 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 from the
very broad and strong O VI absorption frequently found in strong LLSs, SLLSs,
and DLAs (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.

Figure 7: Comparison of the several spectral energy distributions for the
ionizing background radiation field. The HM05 is adopted in this work,
although we systematically compare it with the results from HM12. We also show
the recent EUVB models from Khaire & Srianand (2019) for comparison. The key
ions used to constrain the metallicities are shown where the dotted vertical
bars show their minimum ionization energy while their range of ionization
energy is shown with the horizontal bar.

###  5.2 Photoionization Modeling

To determine the metallicity for each absorber, we follow the same methodology
and make similar assumptions described and discussed in CCC II and CCC III. As
laid out in the previous section, owing 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, 2001,
2012). 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
Fig. 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 3 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., 2021). And indeed we show a posteriori
that the systematic errors arising from using a different EUVB are not as
large as observed at $z\lesssim 1$. We note that for stronger absorbers ($\log
N_{\rm H\,I}\gtrsim 18$) at $z\sim 3$, 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 summarized in Table 3 in CCC II.

The main variables for the photoionization models are the ionization
parameter—$U\equiv n_{\gamma}/n_{\rm H}={\rm H}$, the ionizing photon
density/total hydrogen number density (neutral + ionized)—and the metallicity,
${\rm[X/H]}$. We assume solar relative— heavy element abundances from Asplund
et al. (2009) but allow for possible variation between between C and
$\alpha$-elements (where $\alpha$ 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 detail). ${\rm[C/\alpha]}$ is therefore not necessarily
solar in the gas probed by absorbers studied here (see Aguirre et al. 2004,
and also Lehner et al. 2013; CCC II 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 )
and CCC II.888The code presented in Fumagalli et al. (2016b) and CCC II has
been incorporated into the PyIGM Python package and is available at
https://github.com/pyigm/pyigm (Prochaska et al., 2017a). We perform the MCMC
analysis on an absorber-by-absorber basis. For each absorber we apply Gaussian
priors on $N_{\rm H\,I}$ and redshift999For the metal ions, the lower (upper)
limits, the likelihoods are based on a rescaled Q-function (cumulative
distribution function).. 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 PDFs), 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 ${\rm[C/\alpha]}$ (i.e., for which we could estimate
${\rm[C/\alpha]}$), which absorbers require ${\rm[C/\alpha]}=0$ (e.g., those
with no detections for all the ions or no carbon ions), and which absorbers
may require a Gaussian prior on ${\rm[C/\alpha]}$ (absorbers with poorly
constraints on carbon ions and/or $\alpha$-elements).

Figure 8: The $\log U$ median values as a function of $N_{\rm H\,I}$ for 126
absorbers in the KODIAQ-Z sample. The solid line shows a linear fit to the
data with $14.6\leq\log N_{\rm H\,I}\leq 17.2$ (with a slope $-0.46\pm 0.10$
and intercept $+5.56$). The crosses show the mean values of the median values
in each interval $N_{\rm H\,I}$ interval indicated by the horizontal bar.

When the constraints from the metal ions are not optimal, we adopt a Gaussian
prior on $\log U$ (CCC II). These concern only SLFSs and pLLSs; the 16 LLSs
and 7 SLLSs have enough metal-ion constraints for the ionization models to
converge. Fifty eight percent of the SLFSs and pLLSs do not require a Gaussian
prior on $\log U$. In Fig. 8, we show the $\log U$ median values as a function
of $N_{\rm H\,I}$ for the absorbers that can be modeled assuming a flat prior
on $\log U$. As observed at low $z$, (; ), a strong anti-correlation is
observed between the ionization parameter and $N_{\rm H\,I}$, with a slope
steeper by 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_{\rm H\,I}$ range into two bins, $14.4\leq\log N_{\rm H\,I}<16$ and
$16.2\leq\log N_{\rm H\,I}<17.2$. 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 each PDF is well-fit with a normal distribution with
$\langle\log U\rangle=-1.43\pm 0.59$ and $-2.06\pm 0.53$. We adopt these
results for our Gaussian prior on $\log U$ in the $14.4\leq\log N_{\rm
H\,I}<16$ and $16.2\leq\log N_{\rm H\,I}<17.2$ intervals, respectively.

For the Gaussian prior on ${\rm[C/\alpha]}$, we similarly use absorbers with
relatively good constraints from carbon and $\alpha$-element ions to determine
a mean and standard deviation of
$\left\langle{\rm[C/\alpha]}\right\rangle=-0.31\pm 0.42$. We apply these for
the Gaussian prior on the ${\rm[C/\alpha]}$ ratio for absorbers with poorly-
constrained ${\rm[C/\alpha]}$. We note that for a few (13) absorbers,
${\rm[C/\alpha]}$ is clustered around $-1\lesssim{\rm[C/\alpha]}\lesssim-0.9$.
Not including those in the mean would of course increase the mean
${\rm[C/\alpha]}$ value, but it would not change the posterior PDF of
${\rm[C/\alpha]}$ for these absorbers. We discuss in §§5.6, 6.5 in more detail
the ${\rm[C/\alpha]}$ 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 sightline
name, the absorber redshift ($z_{\rm abs}$), $\log N_{\rm H\,I}$, metallicity,
ionization parameter ($\log U$), and ${\rm[C/\alpha]}$ 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 PDFs provided by the MCMC sampling of the Cloudy models. While those
PDFs can be well-behaved and nearly Gaussian, some are not as well behaved
(see corner plots in the appendix).

For the HD-LLS absorbers used to compare with KODIAQ-Z, we apply the same
overall methodology and same HM05 EUVB in the photoionization models in order
to reduce any systematic errors between the two surveys. We, however, note
that the absorption in absorber in the HD-LLS and KODIAQ-Z is not treated in
the same way. In the HD-LLS survey, a velocity interval of $\pm 500$ ${\rm
km\,s}^{-1}$ defines an absorber, which is driven by the methods to derive
$N_{\rm H\,I}$ that use either the break at the Lyman limit or the damped
wings from Ly$\alpha$. For the metal lines, they sum the absorption components
identified within that interval even though it is rarely larger than 200 ${\rm
km\,s}^{-1}$. As detailed in §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 $\pm 200$–500 ${\rm km\,s}^{-1}$, KODIAQ-Z can have several
absorbers. However, as $N_{\rm H\,I}$ increases and overlaps with the HD-LLS
survey, i.e., for $\log N_{\rm H\,I}\gtrsim 18$, we similarly lose our ability
to derive the column densities in individual components, especially when we
have information only from Ly$\alpha$ and Ly$\beta$. 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 (§3). However, a posteriori we
will show that this different treatment does not lead to large different
results.

###  5.3 Uncertainties on the Metallicities

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.

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.

###  5.4 Ionizing Background

The first one is caused by the uncertainty in the ionizing EUVB as alluded to
previously. At low redshift, CCC II shows that on average the metallicity is
systematically shifted on average by about $+0.4$ dex for the SLFSs/pLLS and
$+0.2$ dex for the LLSs when the ionizing EUVB changes from HM05 to HM12. We
have done a similar experiment with 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 Fig. 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
$\log U$ (that 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 find $\langle{\rm[X/H]}_{\rm
HM12}-{\rm[X/H]}_{\rm HM05}\rangle=+0.10\pm 0.17$. The difference is also not
as systematic as low $z$ since there are several cases where ${\rm[X/H]}_{\rm
HM12}\leq{\rm[X/H]}_{\rm HM05}$.

Figure 10: Difference between the median metallicities for the KODIAQ-Z and
HD-LLS absorbers derived using the HM12 and HM05 EUVBs against the H I column
densities of the absorbers.

As illustrated in Fig. 10 where the differences of the median metallicities
between HM12 and HM05 vs. $N_{\rm H\,I}$ are shown for both the KODIAQ-Z and
HD-LLS absorbers, that difference decreases with increasing $N_{\rm H\,I}$ as
observed for the low redshift absorbers. In fact, for $\log N_{\rm
H\,I}\gtrsim 17.2$, $\langle{\rm[X/H]}_{\rm HM12}-{\rm[X/H]}_{\rm
HM05}\rangle\simeq 0.0\pm 0.1$, i.e., the HM12 and HM05 EUVBs essentially give
the same metallicities and any small differences are essentially random. For
the SLFSs and pLLSs, $\langle{\rm[X/H]}_{\rm HM12}-{\rm[X/H]}_{\rm
HM05}\rangle\simeq+0.12\pm 0.15$. These results are not entirely surprising
when considering the different EUVBs in Fig. 7 with their difference being
smaller than at low $z$ over similar ionizing energies (Gibson et al., 2021).
Based on Fig. 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.

###  5.5 Multiple Ionization Phases

Throughout we also assume a single ionization phase. In cases where a higher
ion (O VI in our survey) is clearly under produced by the photoionizing model,
this ion is not included in the final model. In that 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. However, as we explore lower
$N_{\rm H\,I}$ absorbers in this present survey, we find in many cases that a
single ionization gas-phase can 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 gas-phase 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
structures for the absorbers in our sample would drastically change the
metallicities.

Figure 11: The ratio of ${\rm[C/\alpha]}$ against the H I column density
(left) and neutral fraction (right) from the KODIAQ-Z ($2.2\lesssim z\lesssim
3.6$), HD-LLS ($2.2\lesssim z\lesssim 3.6$), and CCC ($z\lesssim 1$) surveys.
The median values of the ${\rm[C/\alpha]}$ posterior PDFs are adopted as the
central values with 68% CI.

###  5.6 Non-solar Relative Abundances

Other uncertainties may arise from dust or non-solar relative abundances.
Since carbon, silicon, and oxygen ions are typically used to constrain the
ionization modeling at high redshift, we explicitly take account for possible
non-solar ${\rm[C/\alpha]}$ by letting this ratio vary in the range
$-1\leq{\rm[C/\alpha]}\leq+1$. In Fig. 11, we show ${\rm[C/\alpha]}$ as a
function of $N_{\rm H\,I}$ (left panel) and the neutral fraction $f_{\rm
H\,I}\equiv N_{\rm H\,I}/(N_{\rm H\,I}+N_{\rm H\,II})$ (right panel) for the
KODIAQ-Z and HD-LLS absorbers. We also show in this figure the results from
CCC at $z\lesssim 1$ (CCC III). For ${\rm[C/\alpha]}$-$N_{\rm H\,I}$, a
Spearman rank-order test shows no strong evidence for a correlation between
${\rm[C/\alpha]}$ and $N_{\rm H\,I}$ with a correlation coefficient $r_{\rm
S}=0.14$ and a $p$-value $=0.04$. On the other hand, the same statistical test
shows a positive monotonic correlation between ${\rm[C/\alpha]}$ and $f_{\rm
H\,I}$ with a correlation coefficient $r_{\rm S}=0.44$ and a $p$-value $\ll
0.05\%$. However, removing the 13 KODIAQ-Z absorbers with low
${\rm[C/\alpha]}\lesssim-0.9$ eliminates any evidence of correlation between
${\rm[C/\alpha]}$ and $f_{\rm H\,I}$ with a correlation coefficient $r_{\rm
S}=-0.05$ and a $p$-value $=0.88$.

So far the CCC, HD-LLS, and KODIAQ-Z were treated as a single sample.
Considering the samples individually would not change the conclusions for the
lack of correlation between ${\rm[C/\alpha]}$ and $N_{\rm H\,I}$. However,
while neither CCC nor HD-LLS shows any correlation between ${\rm[C/\alpha]}$
and $f_{\rm H\,I}$, for KODIAQ-Z, the Spearman rank-order test shows a strong
positive monotonic correlation between ${\rm[C/\alpha]}$ and $f_{\rm H\,I}$
with a correlation coefficient $r_{\rm S}=0.72$ and a $p$-value $\ll 0.05\%$.
In that case, removing the 13 absorbers with low ${\rm[C/\alpha]}\lesssim-0.9$
still yields $r_{\rm S}=0.55$ and a $p$-value $\ll 0.05\%$. Therefore in the
KODIAQ-Z sample, ${\rm[C/\alpha]}$ is more affected by ionization correction
than in the other samples, most likely because often only C IV is detected
(the limit on C II is often not very constraining and C III is often
contaminated) and larger ionization correction are applied (for CCC, C II and
C III are used to constrain ${\rm[C/\alpha]}$, while for HD-LLS the ionization
corrections are typically smaller). Based on CCC and HD-LLS where we always
find ${\rm[C/\alpha]}\geq-0.6$, ${\rm[C/\alpha]}$-values with
${\rm[C/\alpha]}<-0.6$ are most likely predominantly caused by ionization
rather than nucleosynthesis effects. Excluding the 13 absorbers with low
${\rm[C/\alpha]}\lesssim-0.9$, we find that
$\langle{\rm[C/\alpha]}\rangle=-0.20\pm 0.31$ in the KODIAQ-Z sample compared
to $\langle{\rm[C/\alpha]}\rangle=-0.05\pm 0.30$ in the CCC survey and
$\langle{\rm[C/\alpha]}\rangle=-0.05\pm 0.24$ in the HD-LLS survey. The
similarity between CCC and HD-LLS is quite remarkable since different $z$,
$N_{\rm H\,I}$, and ${\rm[X/H]}$ are probed. The $1\sigma$ 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 ${\rm[C/\alpha]}$ most likely caused by
ionization. In conclusion, the variable C/$\alpha$ 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/$\alpha$ (see also §6.5).

Finally, we note that for stronger LLSs and SLLSs 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
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).

##  6 Metallicity of the IGM/CGM Absorbers at $2.2\lesssim z\lesssim 3.6$

###  6.1 Metallicity Distributions

Figure 12: Posterior metallicity PDF of the SLFSs in the KODIAQ-Z sample. The
light-colored regions indicate the contribution from the upper limits.  Figure
13: Posterior metallicity PDFs of the pLLSs in the KODIAQ-Z sample. The light-
colored regions indicate the contribution from the upper limits.  Figure 14:
Comparison of the posterior metallicity PDFs of the LLSs in the KODIAQ-Z and
HD-LLS sample. The light-colored regions indicate the contribution from the
upper limits.

In Figs. 12, 13, and 14, we show the posterior probability distribution
functions (PDFs) of the metallicities of the SLFSs, pLLSs, and LLSs for
KODIAQ-Z. In Fig. 14, we also overplot the LLSs from HD-LLS survey, showing
similar PDFs between the two surveys, and hence implying we can combine both
samples to improve the statistics for the LLSs. These posterior PDFs are
constructed by combining the normalized metallicity PDFs of all of the
absorbers within a given H I column density regime. In Table 4, 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
${\rm[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_{\rm H\,I}$ increases, the overall median metallicity increases and the IQR
decreases. For pLLSs+LLSs, the median metallicity is $-2.2$ and IQR only
$\sim$1 dex from $-2.8$ to $-1.7$ dex. These features are distinctly
illustrated in Fig. 15 where we compare the PDFs of the SLFSs and pLLSs+LLSs:
there are more frequently pLLSs+LLSs in the metallicity range
$-3.2\lesssim{\rm[X/H]}\lesssim-1.2$ than SLFSs, but SLFSs become more
frequent at extremely low metallicity ${\rm[X/H]}\lesssim-3.5$. Compared to
the SLFS and pLLS PDFs, the LLS PDF is not as negatively skewed with about the
same values for the mean and median.

Figure 15: Comparison of the posterior metallicity PDFs of the SLFSs
(KODIAQ-Z) and pLLSs+LLSs (KODIAD-Z+HD-LLS). The light-colored regions
indicate the contribution from the upper limits.

For completeness and if one ignores the H I column density dependence of the
metallicity PDFs between the SLFSs, pLLSs, and LLSs, we show in Fig. 16 the
resulting metallicity PDF of the absorbers with $14.5<\log N_{\rm H\,I}<19$ at
$2.2\lesssim z\lesssim 3.6$. Since the SLFSs are about 2/3 of the sample, the
metallicity PDF of the $14.5<\log N_{\rm H\,I}<19$ absorbers is not too
different to that of the SLFSs with a similarly skewed distribution toward low
metallicities, but also a second weaker peak around $-0.6$ dex.

Figure 16: Posterior metallicity PDF of the absorbers with $14.5\lesssim\log
N_{\rm H\,I}<19$ in the combined sample of KODIAQ-Z and HD-LLS. The light-
colored regions indicate the contribution from the upper limits.

###  6.2 Metallicity vs. $N_{\rm H\,I}$

As shown in the previous section, there is already evidence for an overall
increase of the metallicities with $N_{\rm H\,I}$ in the range $14.5\leq\log
N_{\rm H\,I}\leq 19$. We now explore further how the metallicity varies with
$N_{\rm H\,I}$ considering the higher $N_{\rm H\,I}$ absorbers (the SLLSs and
DLAs). From Table 4, the trend described above with increasing median/mean
metallicity with increasing $N_{\rm H\,I}$ continues to apply in the SLLS and
DLA regimes. The IQR and standard deviation for the SLLSs are about similar to
that 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 PDFs of the other weaker absorbers, is fully
consistent with a Gaussian distribution according to the Kolmogorov–Smirnov
(KS) test with a $p$-value $\ll 0.1\%$ (see also Rafelski et al. 2012).

These conclusions are strengthened considering the cumulative distribution
functions (CDFs) of the metallicity PDFs (see CCC II for the description of
the estimation of the CDFs including upper and lower limits). We show in Fig.
17 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_{\rm H\,I}$
increases: very metal poor (VMP) absorbers with ${\rm[X/H]}\leq-2.4$ are more
likely to be found in the SLFS, pLLS, LLS regimes than in the SLLS or DLA
regime (${\rm[X/H]}_{\rm VMP}=-2.4$ is the $2\sigma$ lower bound of the DLA
metallicities). Notably the median metallicity of the SLFS is right at that
threshold of ${\rm[X/H]}_{\rm VMP}=-2.4$, implying that 50% of these absorbers
are VMP while only $<5\%$ of the DLAs may be VMP by definition.

Figure 17: Cumulative probabilities of the SLFS, pLLS, LLS, SLLS, and DLA
metallicity PDFs. The SLFSs and pLLSs are from KODIAQ-Z, the LLSs and SLLSs
from KODIAQ-Z and HD-LLS, and the DLAs from R12-DLA.

Thus far, we have separately analyzed the absorbers in different $N_{\rm
H\,I}$ categories. However, we can also simply plot the metallicities of the
absorbers as a function of $N_{\rm H\,I}$, which is shown in Fig. 18. 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 PDFs, and the error bars
represent the 68% CI. For the upper and lower limits, the down and upward
triangles give the $90^{\rm th}$ and $10^{\rm th}$ quantiles while the
vertical bar gives the 80% CI. For the DLAs, the best estimates with their
$1\sigma$ error bars are shown. The horizontal dashed line at $[{\rm X/H}]=0$
represents solar metallicity. The horizontal dotted line at ${\rm[X/H]}=-2.4$
represents the VMP limit as defined above. Note that the lack of data between
$\log N_{\rm H\,I}=20$ and 20.3 is purely artificial, resulting from our grid
of photoionization models stopping at $\log N_{\rm H\,I}=20$.

Fig. 18 further demonstrates the striking metallicity changes as a function of
$N_{\rm H\,I}$. Absorbers with ${\rm[X/H]}\gtrsim-2.4$ are observed at any
$N_{\rm H\,I}$, but absorbers with ${\rm[X/H]}<-2.4$ are observed only at
$\log N_{\rm H\,I}\lesssim 19.8$. While the VMP definition implies this is VMP
gas for DLAs, for low $N_{\rm H\,I}$ absorbers the VMP definition is the
median value for SLFSs and not even $1\sigma$ below the median of pLLSs and
LLSs. As already pointed out above, the frequency of VMP absorbers also
increases as $N_{\rm H\,I}$ decreases, which is quantitatively summarized in
Table 6. We therefore also define extremely metal-poor absorbers as absorbers
with ${\rm[X/H]}<-3$, which are mostly observed at some level in absorbers
with $\log N_{\rm H\,I}\lesssim 18$, see Table 6.101010We emphasize that VMP
DLAs with ${\rm[X/H]}<-2.4$ exist at $z\sim 2$–4, but are found only in
targeted searches (e.g., Cooke et al., 2011, 2016)—they are uncommon and they
do not show in sample as small as presented here since they are only $<5\%$ of
the population. Extremely metal-poor absorbers have metallicities that are so
low that only Pop III stars may have contaminated them if they have some
metals (see, e.g., 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 non-binned distribution
shown in Fig. 18. In that figure, there is a lack of data points at
$-1.1\lesssim{\rm[X/H]}\lesssim-1$ for H I column density absorbers with
$14.5\lesssim\log N_{\rm H\,I}\lesssim 16.4$. This is reminiscent of the
metallicity dip at a similar value of ${\rm[X/H]}\simeq-1$ observed in the
pLLS+LLS PDF at $0.45\lesssim z\lesssim 1$ (CCC III, and see §6.3 for the
cosmic evolution of the metallicity).

Figure 18: Metallicities of the SLFSs ($14.5<\log N_{\rm H\,I}<16.2$), pLLSs
($16.2\leq\log N_{\rm H\,I}<17.2$), LLSs ($17.2\leq\log N_{\rm H\,I}<19$),
SLLSs ($19\leq\log N_{\rm H\,I}<20.3$), and DLAs ($\log N_{\rm H\,I}\geq
20.3$) at $2.2.\lesssim z\lesssim 3.6$ as a function of $N_{\rm H\,I}$. The
blue data from the KODIAQ-Z survey (this paper). The orange data are from the
HD-LLS survey (Fumagalli et al. 2016b) but we re-estimated 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).

###  6.3 Cosmic Evolution of the Metallicity

Figure 19: Comparison of the posterior metallicity PDFs of the absorbers with
$15\leq\log N_{\rm H\,I}<19$ at high and low redshifts. The light-colored
regions indicate the contribution from the upper limits. For CCC, The light
colored histogram indicates the contributions from upper and lower limits.

In Lehner et al. (2016), for the first time, we studied the redshift evolution
of the metallicity absorbers with $16.2\lesssim\log N_{\rm H\,I}<19$, 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_{\rm H\,I}$ regimes below $\log N_{\rm H\,I}\lesssim 16.2$. In Fig. 19, we
show the metallicity PDFs at low ($z<1$) and high ($2.2\leq z\leq 3.6$)
redshift for H I-selected absorbers with $15\lesssim N_{\rm H\,I}<19$. There
is some overlap between the two distributions between
$-3\lesssim{\rm[X/H]}\lesssim-1.4$, but evidently there are more absorbers
with ${\rm[X/H]}\lesssim-2$ at high $z$ while many more with ${\rm[X/H]}>-1.2$
at $z<1$. In fact the latter value corresponds to the median metallicity at
$z<1$ for $16.2\lesssim\log N_{\rm H\,I}<19$ absorbers, and only 10% of the
high redshift absorbers have metallicities ${\rm[X/H]}\geq-1.2$. The medians
and means at $z<1$ are about $+1$ dex higher than at $2.2\leq z\leq 3.6$. Both
distributions are negatively skewed with similar IQRs$\,\simeq 1.5$ at low and
high $z$.

Figure 20: The metallicity as a function of $N_{\rm H\,I}$ at low $z$ (CCC)
and high $z$ (KODIAQ-Z, HD-LLS, R12-DLA). For the $z<1$ absorbers, we subtract
the metallicity of each absorber by $0.9$ dex (see text for more details). For
the black circles, lighter symbols are upper limits on the metallicity. For
the blue crosses, lighter symbols are upper limits if ${\rm[X/H]}\leq-1$ and
lower limits if ${\rm[X/H]}>-1$. For clarity, we did not plot the error bars
(which are available in Fig. 18 and Fig. 9 of CCC III for the low redshift
sample). Median metallicities are plotted, but that we treat upper and lower
limits as Fig. 18, i.e., the lower and upper values represent the 10th and
90th percentiles, respectively.

In Fig. 20, we compare the metallicities of the absorbers (following the same
methodology used to produce Fig. 18) at $z<1$ and $2.2\leq z\leq 3.6$ over the
H I column density range $14.5\lesssim N_{\rm H\,I}<22$. For the low redshift
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, and $-1$ dex for lower $N_{\rm 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 dispersion at low and high $z$ show striking
overlap at any given $N_{\rm H\,I}$ (where there is a good sampling in both
datasets). 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\leq z\leq 3.6$ and $z<1$ absorbers is an increase of their
metallicity by a factor $\sim$8.

Figure 21: The metallicity of different absorbers as a function of the
redshift. For HD-LLS and KODIAQ-Z surveys, lighter symbols are upper limits on
the metallicity. For CCC, lighter symbols are upper limits if
${\rm[X/H]}\leq-1$ and lower limits if ${\rm[X/H]}>-1$. Median metallicities
are plotted for the CCC, KODIAQ-Z, and HD-LLS surveys, but we treat upper and
lower limits as Fig. 18, i.e., the lower and upper values represent the 10th
and 90th percentiles, respectively.

In Fig. 21, 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 $\log N_{\rm H\,I}<19$ absorbers
with $-0.15\lesssim{\rm[X/H]}\lesssim+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 ${\rm[X/H]}\lesssim-3$, while there is at
$2.2\leq z\leq 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\leq z\leq
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 $\log N_{\rm H\,I}\lesssim 18$ VMP
absorbers at low and high $z$ lie in the respective VMP regime.

In Fig. 21, there is also some evidence that the metallicity distribution
somewhat changes around $z\sim 2.85$. Considering the KODIAQ-Z and HD-LLS
sample with $\log N_{\rm H\,I}<19$ absorbers at $z>2.85$, the fraction of
absorbers with ${\rm[X/H]}\geq-1$ is very small ($4/112\simeq 3.6\%$) while at
$z\leq 2.85$ it is a factor $\sim$4 larger ($17/130\simeq 13.1\%$). To explore
that potential of the change of the metallicity with $z$, in Fig. 22, we show
the comparison between the metallicity CDFs of the absorbers with
$14.5\leq\log N_{\rm H\,I}<19$ from KODIAQ-Z and HD-LLS above and below
several redshift thresholds, $z_{\rm th}=2.8$, 2.9, 3.0. The large number of
absorbers with ${\rm[X/H]}\gtrsim-1$ at $z\lesssim 2.9$ compared to higher $z$
is further demonstrated in this figure. However, this figure also shows that
at $z\gtrsim 2.9$, and in particular at $z\geq 3.0$, there is a much larger
fraction of metal-poor absorbers with ${\rm[X/H]}-2.8$ than at lower
redshifts, demonstrating that within the redshift interval $2.2\lesssim
z\lesssim 3.6$, there is a change in the metallicity PDFs with an overall
increase of the metallicity at $2.2\lesssim z<2.8$ compared to $2.8\lesssim
z\lesssim 3.6$.

Figure 22: Comparison of the metallicity CDFs of the absorbers in KODIAQ-Z and
HD-LLS with $14.5\leq\log N_{\rm H\,I}<19$ above and below a given redshift
threshold, $z_{\rm th}=2.8$, 2.9, 3.0. The numbers between parentheses
indicate the number of absorbers in each redshift interval.

###  6.4 Metallicity Variation over Small Velocity Scales

Owing to the treatment of the absorbers where we separate individual H I
components as much as possible (see §3), our KODIAQ-Z survey provides also
information on the metallicity variation over small redshift/velocity
separation between absorbers. Following CCC III, we define here paired
absorbers as absorbers that are closely separated in redshift space so that
$\Delta v\equiv|(z^{2}_{\rm abs}-z^{1}_{\rm abs})/(1+z^{1}_{\rm abs})\,c|<500$
${\rm km\,s}^{-1}$. We have a sample of 37 such paired absorbers.

In Fig. 23, we show in the left panel the absolute metallicity difference
($\Delta{\rm[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 $\Delta{\rm[X/H]}$. As observed at low $z$
(CCC III), we find no obvious relation between $\Delta v$ and the level of the
metallicity variation. At any velocity separation, there is a large scatter
from about 0.2 dex to $>2$ dex. Using the survival method to account for the
lower limits (Feigelson & Nelson, 1985), we estimate a mean
$\langle\Delta{\rm[X/H]}\rangle=1.20\pm 0.16$ (error on the mean) and the IQR
in $\Delta{\rm[X/H]}$ is $[0.38,1.81]$. There appear to be two main clusters
of data separated at $\Delta{\rm[X/H]}\simeq 1$ with about 50% of the sample
in each bin and with a mean $\langle\Delta{\rm[X/H]}\rangle=0.39\pm 0.05$ and
$2.15\pm 0.13$ below and above this limit, respectively. Therefore for the
about half of the paired absorbers, there is evidence for metallicity
variations over $\Delta v\lesssim 500$ ${\rm km\,s}^{-1}$ of a factor 2–3 and
for the other half of a factor $>140$.

Owing to the large metallicity variation over $\Delta v\lesssim 500$ ${\rm
km\,s}^{-1}$, a posteriori it is justified to treat these paired absorbers as
individual absorbers. The other direct consequence is, 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, 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 (; ; see also Kacprzak
et al. 2019; Zahedy et al. 2021). We discuss the implication of these findings
in §10.4 in more detail.

Figure 23: Left: scatter plot of the difference between the higher and lower
median metallicities for the closely-redshift separated absorbers (a.k.a.
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.

###  6.5 The ${\rm[C/\alpha]}$ ratio as a Function of the Metallicity

Figure 24: The ${\rm[C/\alpha]}$ ratio against the metallicity of the
absorbers from the KODIAQ-Z ($2.2\lesssim z\lesssim 3.6$), HD-LLS
($2.2\lesssim z\lesssim 3.6$), and CCC ($z\lesssim 1$) surveys. The median
values of the ${\rm[C/\alpha]}$ and ${\rm[X/H]}$ posterior PDFs are adopted
with 68% CI. Note that ${\rm[C/\alpha]}$-values with ${\rm[C/\alpha]}<-0.6$
are most likely predominantly caused by ionization corrections rather than
nucleosynthesis effects (see text for more details).

In the photoionization modeling of the absorbers to determine their
metallicities, the ${\rm[C/\alpha]}$ ratio is allowed to vary in the range
$-1\lesssim{\rm[C/\alpha]}\lesssim+1$ to take into account possible
nucleosynthesis effects that likely arise from the time-lag in the production
between $\alpha$-elements and carbon (e.g., Cescutti et al. 2009; Mattsson
2010). The general trend between ${\rm[C/\alpha]}$ and ${\rm[\alpha/H]}$ seen
in stars (e.g., Akerman et al., 2004; Fabbian et al., 2010) and SLLSs/DLAs
(e.g., Pettini et al., 2008; Penprase et al., 2010; Cooke et al., 2011) is a
decrease of ${\rm[C/\alpha]}$ from about ${\rm[C/\alpha]}\simeq 0$ (up to
$+0.2$ dex) at ${\rm[\alpha/H]}\gtrsim 0$ to ${\rm[C/\alpha]}\simeq-0.6$ at
${\rm[\alpha/H]}\simeq-0.5$ where ${\rm[C/\alpha]}$ plateaus around that
metallicity value up to a metallicity of ${\rm[\alpha/H]}\simeq-2$ where an
upturn is observed in ${\rm[C/\alpha]}$ with an increase of ${\rm[C/\alpha]}$
to ${\rm[C/\alpha]}\gtrsim+0.2$ at ${\rm[\alpha/H]}\leq-2.8$.

As discussed in CCC III at low redshift and Lehner et al. (2016) at high
redshift, the SLFSs, pLLSs, and LLSs do not quite follow this overall trend.
In fact, little trend is observed between ${\rm[C/\alpha]}$ and
${\rm[\alpha/H]}$ for these absorbers with a large scatter in
${\rm[C/\alpha]}$ at any metallicity. While no obvious trend is observed
between ${\rm[C/\alpha]}$ and ${\rm[\alpha/H]}$ for these absorbers, CCC III
note three characteristics: 1) a lower floor level for the distribution at
about ${\rm[C/\alpha]}\simeq-0.6$, similar to that observed for the stars, H
II regions, and SLLSs/DLAs; 2) an upper ceiling of the distribution at
${\rm[C/\alpha]}\simeq+0.5$, but with a concentration of data around
$0<{\rm[C/\alpha]}\sim+0.3$, similar again to the highest values observed in
stars, H II regions, and SLLSs/DLAs; and 3) a tentative upturn in
${\rm[C/\alpha]}$ at ${\rm[X/H]}\lesssim-2$, which would be similar to the one
observed in stars and SLLSs/DLAs.

With KODIAQ-Z and our re-modeling 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 Fig. 24, we show ${\rm[C/\alpha]}$ as a
function of ${\rm[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
${\rm[C/\alpha]}$, see Fig. 11, which we did not highlight in Fig. 24).
Excluding the 13 absorbers with ${\rm[C/\alpha]}\lesssim-0.9$ (because those
are more largely affected by ionization, see 5.3 and below), the addition of
the KODIAQ-Z and HD-LLS data confirm the first two characteristics described
above and summarized in CCC III. However, at ${\rm[X/H]}\lesssim-2$, there is
no clear upturn in ${\rm[C/\alpha]}$ for the absorbers with $\log N_{\rm
H\,I}<19$, although the scatter in the ${\rm[C/\alpha]}$ distribution becomes
smaller. On the other hand, considering the SLLSs from HD-LLS (diamond symbols
in Fig. 24), the upturn in ${\rm[C/\alpha]}$ is apparent when the metallicity
becomes smaller than ${\rm[X/H]}\lesssim-2.2$.

Since the derived ratio of ${\rm[C/\alpha]}$ is controlled by both ionization
corrections and nucleosynthesis effects, it is plausible that the ionization
corrections are washing away some of the trends that would appear if only
nucleosynthesis effects were at play (as for example for the DLAs). In §5.6,
we show that ${\rm[C/\alpha]}$ in the KODIAQ-Z sample is more affected by
ionization than in the two other samples in view of the correlation between
${\rm[C/\alpha]}$ and $f_{\rm H\,I}$. However, no such correlation is observed
in CCC and HD-LLS, and the overall dispersion is similar in these three
surveys. The mean values for ${\rm[C/\alpha]}$ in CCC and HD-LLS are also the
same (see §5.6), despite widely different ionization correction levels.
Therefore the ${\rm[C/\alpha]}$ ratio with ${\rm[C/\alpha]}>-0.6$ determined
from the ionization modeling does not seem to have a strong dependence with
the level of the ionization correction. Nevertheless, the ionization
correction essentially adds noise in the ${\rm[C/\alpha]}$ distribution at the
level of $\pm 0.3$ dex at any $N_{\rm H\,I}$ and metallicity, which is large
enough to hide subtle changes of ${\rm[C/\alpha]}$ with ${\rm[X/H]}$ at the
expected level of 0.2–0.5 dex. Hence, the scatter observed in Fig. 24 must be
largely caused by ionization while the range
$-0.6\lesssim{\rm[C/\alpha]}\lesssim+0.5$ where most of the absorbers lie may
be controlled by nucleosynthesis effects based on its similarity observed in
DLAs, stars, and H II regions.

##  7 Physical Properties of the SLFSs, pLLSs, LLSs, and SLLSs at $2.2\lesssim
z\lesssim 3.6$

The photoionization models discussed in §5 also predict the physical
properties of the absorbers. In particular, these include the neutral fraction
of the gas ($f_{\rm HI}$), the density of the gas ($n_{\rm H}$), the total H
column density ($N_{\rm H}\equiv N_{\rm H\,I}+N_{\rm H\,II}$), the length-
scale of the cloud ($l\equiv n_{\rm H}/N_{\rm 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; CCC III). A 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.

###  7.1 Densities, H Column Densities, Temperatures, and Length Scales

Figure 25: Posterior PDFs of the total H column density (left), H density
(middle), and path-length (right) for SLLSs, LLSs, pLLSs, and SLFSs for the
entire sample from the KODIAQ-Z survey and HD-LLSs.

We summarize the physical quantities ($n_{\rm H}$, $N_{\rm H}$, $l$, and $T$)
derived from our modeling for the KODIAQ-Z sample in Table 7. Values of $\log
n_{\rm H}$ with a trailing colon represents absorbers for which we adopted a
Gaussian prior on $\log U$, and therefore on $\log n_{\rm H}$. In Table 8, we
summarize the statistical properties of the physical properties derived for
the entire sample (median, mean, and dispersion of $N_{\rm H\,I}$, $n_{\rm
H}$, $N_{\rm H}$, $l$, and $T$) and a restricted sample where only a flat
prior on $\log U$ 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, SLLSs did
not require a Gaussian prior on $\log U$. However, for about 45% of the pLLSs,
we had to use a Gaussian prior on $\log U$, and the results are essentially
the same in both sample. Hereafter, we therefore consider the results from the
entire sample. In Fig. 25, we show the posterior PDFs of $N_{\rm H}$, and
$n_{\rm 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 PDFs contrast remarkably from the metallicity PDFs 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 small dispersion, and there is a clear
trend of decreasing temperatures as $N_{\rm H\,I}$ increases of about 0.1 dex
decreases between each absorber category from the SLFSs, pLLSs, LLSs, to SLLSs
(see Table 8). For the SLFSs, the average temperature is around $3\times
10^{4}$ K, while it is a factor 2 less on average for the SLLS. In §4.2, we
show that for about 90% of the H I components $\langle b\rangle=27\pm 6$ ${\rm
km\,s}^{-1}$. This places an upper limit on the temperature from the
observations of $T<4\times 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 non-thermal component, while for the
other absorbers, non-thermal motions take over the thermal ones. Compared to
the CCC III, the predicted temperatures of the $z<1$ SLFSs, pLLSs, and LLSs
are about a factor 1.4 times smaller than at $2.2\lesssim z\lesssim 3.6$.

Considering the SLFSs, pLLSs, and LLSs, the total H column densities decrease
with decreasing $N_{\rm H\,I}$, but their distributions overlap within
$1\sigma$ dispersion (that are large for these quantities, see Table 8 and
Fig. 25). The differences between the mean $\log N_{\rm H\,I}$ for SLFSs
compared with the pLLSs and the pLLSs to the LLSs are a factor 4 and 28,
respectively. In contrast, the differences in total H are significantly
smaller with differences of a factor $\sim 4$ and 6.5 between the SLFSs and
pLLSs and pLLSs and LLSs, respectively. Thus, despite their lower $N_{\rm
H\,I}$ values, the contribution 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
and see §8). The $N_{\rm H}$ column densities for SLLSs are similar to those
for LLSs, $N_{\rm H}$ is about the same even though $N_{\rm 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 CCC III, $N_{\rm H}$ at the $z<1$ SLFSs, pLLSs,
and LLSs are about a factor $\sim 10$–15 smaller than at $2.2\lesssim
z\lesssim 3.6$.

The estimated hydrogen densities increase as $N_{\rm H\,I}$ increases. The
SLLSs has the largest dispersion on $n_{\rm H}$, which is caused by a long
tail at low and high $n_{\rm H}$ values (see Fig. 25). While there is an
increase in $n_{\rm H}$ with increasing $N_{\rm H\,I}$, the bulk of the
hydrogen density for the pLLSs, LLS, and LLSs are within about $\log n_{\rm
H}\simeq-3.4$ and $-1.8$, while for SLFSs it is between $\log n_{\rm
H}\simeq-4.2$ and $-2.5$. At $z\sim 2.8$ (the mean redshift of KODIAQ-Z), the
mean cosmic density is about $\bar{n}_{\rm H}\simeq 10^{-5}$ cm-3 (Schaye,
2001b). The pLLSs, LLSs, and SLLSs therefore probe over-densities of
$\delta\equiv n_{\rm H}/\bar{n}_{\rm H}-1\sim 40$–1600 while for the lower
column SLFSs $\delta\sim 5$–320. A large part of the observed pLLSs, LLSs, and
SLLSs largely probe virialized (mini)halos ($\delta>100$, Schaye 2001b). For
SLFSs, there is more a mixture of systems probing over-densities below and
above 60, at the transition between the IGM and virialized halos (Schaye,
2001b). Compared to the CCC III, $n_{\rm H}$ at the $z<1$ SLFSs, pLLSs, LLSs,
SLLSs are about a factor $\sim 1.5$–2.5 smaller than at $2.2\lesssim z\lesssim
3.6$. However, owing to the redshift dependence of the mean cosmic density
that is $\propto((1+z)/4)^{3}$ (Schaye, 2001b), the density contrast at $z<1$
for the SLFSs and other stronger H I absorbers is always $\delta>100$, even
for the lower densities probed by SLFSs (although we note that at $z<1$, only
$\log N_{\rm H\,I}>15$ SLFSs could be probed).

The length-scales, which are directly related to $n_{\rm H}$ and $N_{\rm H}$,
have the largest dispersion among all these physical quantities. For the SLFSs
and pLLSs, the interquartile range (IQR) of length-scale of the absorbers is
about $4\lesssim l\lesssim 150$ kpc. For the LLSs, IQR for $l$ is $18\lesssim
l\lesssim 288$ kpc and for the SLLSs it is $8\lesssim l\lesssim 80$ kpc. In
contrast at $z<1$, the IQR of length-scale of these absorbers is $0.3\lesssim
l\lesssim 3.2$ kpc (CCC III). At $2.2\lesssim z\lesssim 3.6$, the SLFSs,
pLLSs, LLSs, and even the SLLSs may probe a variety of structures from gas
embedded in virialized halos to large-scale 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 intra-
group gas (e.g., Fumagalli et al., 2016a, 2017; Lofthouse et al., 2020).

###  7.2 Neutral Fractions

The behaviour of the neutral fraction with $N_{\rm H\,I}$ can be used to
constrain the cosmic baryon and metal budgets (see §8). In Fig. 26, we show
the neutral fraction of the absorbers estimated from the Cloudy models as a
function of $N_{\rm H\,I}$ from $\log N_{\rm H\,I}\ \simeq 14.6$ to 20 from
KODIAQ-Z and HD-LLS (see §5.2, and see also the recent work for the SLLSs from
Berg et al. 2021). There is a clear correlation between $f_{\rm H\,I}$ and
$N_{\rm H\,I}$. This is confirmed by the Spearman rank-order test that shows a
strong positive monotonic correlation between $f_{\rm H\,I}$ and $N_{\rm
H\,I}$ with a correlation coefficient $r_{\rm S}=0.9$ and a $p$-value $\ll
0.05\%$. However, the slope of the correlation is shallower between the SLFS,
pLLS, LLS regimes compared to the LLS, SLLS regimes. The dashed red line in
Fig. 26 is a linear fit to the data with $14.6\leq\log N_{\rm H\,I}<19$,
yielding $\log f_{\rm H\,I}=(0.51\pm 0.04)\log N_{\rm H\,I}-12.00$. The slope
is quite similar to that found at $z<1$ for absorbers with $15.1\leq\log
N_{\rm H\,I}<19$ (CCC III). The dash-dotted magenta line is a linear fit to
the data with $17.2<\log N_{\rm H\,I}\leq 20$, $\log f_{\rm H\,I}=(1.00\pm
0.08)\log N_{\rm H\,I}-20.76$, which is consistent with the linear fit derived
in the HD-LLS survey for the LLSs+SLLSs (Fumagalli et al., 2016b).

To model this overall trend between $f_{\rm H\,I}$ and $N_{\rm 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 about this method, 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 Regression package within
scikit-learn (Pedregosa et al., 2011; Buitinck et al., 2013) to model the data
with a squared-exponential 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 single
data point, and we use the $1\sigma$ 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 Fig. 26 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 black curve is well-fitted by a third order
polynomial, $\log f_{\rm H\,I}=-0.0104x^{3}+0.6268x^{2}-11.62x+64.23$ (where
$x=\log N_{\rm H\,I}$).

Figure 26: The neutral fraction of H I versus $N_{\rm H\,I}$ using the
KODIAQ-Z and HD-LLS surveys. The median values of the posterior PDFs are
adopted with 68% CI. Light blue data were estimated using a Gaussian prior on
$\log U$, while the other data have a flat prior on $\log U$. The black curve
is a GP-derived model to the entire $N_{\rm H\,I}$ range shown in the figure
(the dark area corresponds to the dispersion derived from the GP model). The
dashed red line is a linear fit to the data with $14.6\leq\log N_{\rm
H\,I}<19$ (i.e., the SLFSs, pLLSs, and LLSs). The dash-dotted magenta line is
a linear fit to the data with $17.2<\log N_{\rm H\,I}\leq 20$ (i.e., the LLSs
and SLLSs). For all the models, we consider only data with a flat prior on
$\log U$ (although the results are similar if the entire dataset was
considered instead).

##  8 Cosmological Budgets

###  8.1 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 photo- and collisionally-ionized gas found in $\tau_{\rm LL}>1$
systems is very likely the second largest contributor to the baryon budget at
$z\sim 2$–3. However, this result could be affected by small-number statistics
as Lehner et al. sample included only 7 LLSs and 8 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\lesssim z\lesssim 3.6$.

The mass gas density relative to the critical density, $\Omega_{\rm g}$, can
be estimated via (e.g., Tytler 1987; Lehner et al. 2007; O’Meara et al. 2007;
Prochaska et al. 2010; CCC III):

$\Omega_{\rm g}=\frac{\mu_{\rm H}m_{\rm H}H_{0}}{\rho_{c}\,c}\int\frac{N_{\rm
H\,I}}{f_{\rm H\,I}}\mathcal{F}(N_{\rm H\,I})dN_{\rm H\,I}\,,\\\ $ (1)

where $m_{\rm H}=1.673\times 10^{-24}$ g is the mass of hydrogen, $\mu_{\rm
H}=1.3$ corrects for the presence of helium, $c=2.9979\times 10^{10}$ cm s-1
is the speed of light, $\rho_{c}=3H^{2}_{0}/(8\pi G)=8.62\times 10^{-30}$ g
cm-3 is the $z=0$ critical density, $H_{0}=2.20\times 10^{-18}$ s-1 is the
Hubble constant,111111We adopt here the Planck 2015 results (Planck
Collaboration et al. 2016). and $\mathcal{F}(N_{\rm HI})$ is the column
density distribution function. For $\mathcal{F}(N_{\rm HI})$, we adopt the
spline function model from Prochaska et al. (2014) at $z\approx 2.5$ (see
their Table 2 and Fig. 7). For the neutral fraction, $f_{\rm H\,I}$, we adopt
the GP model from §7.2. Integrating the above equation over the different
$N_{\rm H\,I}$ regimes of the SLFSs, pLLSs, LLSs, and SLLSs, we find that
$\Omega_{\rm g}\simeq(2.7,1.1,1.4,1.1)\times 10^{-3}$ for the SLFSs, pLLSs,
LLSs, and SLLSs, respectively. The contributions of $\Omega_{\rm g}$ to the
total baryon density are $\Omega_{\rm g}/\Omega_{\rm b}\simeq
5.5\%,2.2\%,2.8\%,2.3\%$ for the SLFSs, pLLSs, LLSs, and SLLSs, respectively
($\Omega_{\rm 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).121212Although we remind the reader that
here the SLLSs include only absorbers with $19\leq\log N_{\rm H\,I}\leq 20$
owing to the grid of ionization models stopping at $\log N_{\rm H\,I}=20$.

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 to
0.4–10% to baryonic budget. As noted in §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_{\rm H\,I}$ decreases to
$\log N_{\rm H\,I}<17$ compared to stronger H I absorbers. The observed column
densities for these narrow O VI components are also well-reproduced by single-
phase photoionization models as those described in §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
contribute $\Omega_{\rm g}/\Omega_{\rm b}\simeq 7\%$ and the SLFSs contribute
to another 5.5% to the total baryon budget.131313We separate the SLFSs from
the other absorbers because, e.g., Kim et al. (2001) included those in the
LYAF contribute where they estimate that 90% of all baryons reside in the LYAF
at $1.5<z<4$. In contrast, the neutral gas that is mostly found in DLAs yield
only $\Omega_{\rm g}/\Omega_{\rm b}\simeq 2\%$ at $z\sim 2$–5. (e.g.,
Prochaska et al., 2005; Wolfe et al., 2005; Lehner et al., 2014; Péroux &
Howk, 2020). At $z<1$, CCC III find fractional contributions to baryonic
budget for the SLFSs, pLLSs, and LLSs that are about 10 times smaller. At
$z<1$, the SLFSs, pLLSs, and LLSs are probing much larger overdensities, and
at these redshifts only $\sim 5\%$ of the baryons are thought to be in the CGM
of galaxies (Shull et al., 2012).

###  8.2 Cosmological Metal Budget

The metal density of the absorbers can be estimated with the baryon density
via $\Omega_{\rm m}=\Omega_{\rm g}Z$, where $Z=10^{{\rm[X/H]}}Z_{\sun}$ is the
metallicity of the gas in mass units and $Z_{\sun}=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 about 0.5 dex bin of $N_{\rm H\,I}$ from
$\log N_{\rm H\,I}=14.5$ to $\log N_{\rm H\,I}=20$. We note that the linear
mean metallicities for SLFSs, pLLSs, and LLSs are $0.035,0.032,0.028Z_{\sun}$,
respectively. We find $\Omega_{\rm m}=(1.4,4.7,5.9)\times 10^{-7}$ for the
SLFSs, pLLSs, and LLSs, respectively. The latter value is consistent with that
derived by Fumagalli et al. (2016b), $\Omega_{\rm m}=5.1\times 10^{-7}$. It is
not surprising that our value is somewhat larger since KODIAQ-Z added a few
LLSs on the high metallicity side, including one with ${\rm[X/H]}\simeq-0.4$
(see, e.g, Fig. 18) and since the estimate of $\Omega_{\rm m}$ that depends on
the mean linear metallicity is much more sensitive to the positive tail end of
the logarithmic metallicity PDF. Since we do not have SLLSs with $\log N_{\rm
H\,I}>20$ in our sample, we use the result from the HD-LLS survey where
Fumagalli et al. (2016b) derived $\Omega_{\rm m}=1.6\times 10^{-6}$ for the
SLLSs. For reasons (sample variance, using the arithmetic mean) discussed in
detail in Fumagalli et al., these values are uncertain by a factor $\sim 3$
(see also discussion in Berg et al. 2021). We note that the comoving metal-
mass density can be estimated via $\rho_{\rm m}=\Omega_{\rm m}/\rho_{c}$
(e.g., Péroux & Howk 2020). Therefore, for the combined photoionized SLFSs,
pLLSs, and LLSs, we find $\rho_{\rm m}=1.5\times 10^{5}$ M⊙ cMpc-3, while
combined with the SLLSs, this increases to $3.5\times 10^{5}$ M⊙ cMpc-3.

The total expected cosmic metal density from the output of stars at $z\sim
2.8$ is about $\Omega^{\rm exp}_{\rm m}\simeq 9\times 10^{-6}$ or expressed in
terms of cosmic metal mass density $\rho^{\rm exp}_{\rm m}\simeq 10^{6}$ M⊙
cMpc-3 (Péroux & Howk, 2020). Péroux & Howk (2020) also recently estimated for
the DLAs, i.e., for the neutral gas, that $\Omega_{\rm m}\simeq 3.5\times
10^{-6}$ and for the dust $\Omega_{\rm m}\simeq 1\times 10^{-6}$ at
$2.2\lesssim z\lesssim 3.6$, , i.e., about 40% of the metals are in the
neutral gas and about 10% are in form of dust at $2.2\lesssim z\lesssim 3.6$.
The SLLSs account for 17%, while the SLFSs, pLLSs, and LLSs account for 1.5%,
5.1%, 6.4%, and all combined they account for about 30% of the metals at
$2.2\lesssim z\lesssim 3.6$. We note that the recent estimates from Péroux &
Howk (2020) for the SLLSs and LLSs are about within $1\sigma$ 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 $\gtrsim 7\%$ (Lehner et al., 2014).

While most of the baryons at $2.2\lesssim z\lesssim 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 90% 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 $\log N_{\rm
H\,I}\gtrsim 14.5$ absorbers. Using the recent results on $\Omega_{\rm C\,IV}$
from D’Odorico et al. (2013, 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 $\Omega_{\rm m}\gtrsim
3.4\times 10^{-7}$ at $z\sim 2.8$ for the “IGM” or about $\gtrsim 4\%$ of the
total cosmic metal budget (although we note this may include contributions
from higher H I column density absorbers than LYAF absorbers). Overall the
ionized universe probed by absorbers $\log N_{\rm H\,I}<20$ contains about 40%
of the metals at $2.2\lesssim z\lesssim 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\sim 2.8$.

At $z<1$, $\Omega_{\rm m}=(4.6,6.5,25.9)\times 10^{-7}$ for the SLFSs, pLLSs,
and LLSs, respectively, (CCC III).141414In CCC III, we inadvertently used the
mean logarithmic metallicity, instead of the linear mean, which does not
change the value for the LLSs, but does for the SLFSs and pLLSs, which are,
respectively, $0.177Z_{\sun}$ and $0.197Z_{\sun}$. This increases $\Omega_{\rm
m}$ by a factor 5 for these absorbers. As the LLSs dominate the metal budget,
the effect on the total $\Omega_{\rm m}$ for the SLFSs, pLLSs, and LLSs is
only an increase by a factor 1.3, i.e., $\Omega_{\rm m}\simeq 3.7\times
10^{-6}$ for these absorbers at $z<1$.  In terms of comoving mass density, the
combined SLFSs, pLLSs, and LLSs at $z<1$ have $\rho_{\rm m}=4.7\times 10^{5}$
M⊙ cMpc-3. Therefore the amount of metals in SLFSs, pLLSs, and LLSs has
increased by a factor 3 from $z\simeq 2.2$–3.6 to $z\lesssim 1$. However, the
total expected cosmic metal density at $z\sim 0.5$ is about $\Omega^{\rm
exp}_{\rm m}\simeq 6\times 10^{-5}$ or cosmic metal mass density $\rho^{\rm
exp}_{\rm m}\simeq 8\times 10^{6}$ M⊙ cMpc-3 (Péroux & Howk, 2020), a factor 7
larger than at $z\sim 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 2 times lower than the contribution of
the same absorbers at $2.2\lesssim z\lesssim 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 intra-cluster medium,
demonstrating a major shift in the distribution of the metals from 2.3 Gyr to
$\gtrsim 9$ Gyr post-Big Bang.

##  9 Comparison to the FOGGIE Simulations

###  9.1 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 comprised of
cosmological zoom 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 §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
(van de Voort et al., 2019; 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$,
owing to the spatial resolution being at a fixed comoving size). This spatial
resolution profile enables both the development of small relatively-dense and
cool cloudlets and 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 $\sim
10^{12}$ M⊙ at $z=0$ and have no mergers of mass ratios more than 1:10 after
$z=2$ (see also Zheng et al. 2020 and Lochhaas 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 sightline 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 cutoff151515Adjusting 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 a 11% increase.. An interval is
defined encompassing the sightline’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
${\rm km\,s}^{-1}$ of each other. These steps are repeated until the total
column density of the remaining, unmasked cells falls below $\log N_{\rm
H\,I}=12$. SALSA then calculates the total column density and column-density
weighted 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 sightlines 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, Smith, &
Silvia, 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 sightline, within 100 pc of each other, and
have line-of-sight velocities within 10 ${\rm 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.

###  9.2 Comparison of the Metallicities

Figure 27: 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,
R12-DLA surveys (where we limit these surveys to $2\leq z\leq 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 absobers
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
$\log N_{\rm X}>12$. 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.  Figure 28: Comparison of the
metallicity CDFs for the different categories of absorbers in the KODIAQ-Z+HD-
LLS+R12 observations and FOGGIE simulations at $2\leq z\leq 3$ (the sample of
pLLSs in this redshift interval for the observation sample is too small to
produce a reliable CDF and therefore is not shown here).

The column density and average metallicity of our extracted FOGGIE H I
absorbers are shown in Fig. 27. The plot bounds and dotted and dashed lines
are the same as in Fig. 18. The left panel shows all H I absorbers in cyan,
with the observational samples overplotted. Gray-scale 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 selection—is
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
${\rm[X/H]}\gtrsim-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 a 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_{\rm H\,I}$
decreases seen in the middle panel is caused by the column density requirement
of $\log N_{\rm X}>12$ 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 agreeement 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_{\rm 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_{\rm
H\,I}$, with the spread tending to decrease as $N_{\rm H\,I}$ increases. This
is remarkable because neither EAGLE nor FIRE zoom simulations could reproduce
the similar trend seen at low redshift (; ). In Fig. 28, we compare the
metallicity CDFs for the various absorbers at $2\lesssim z\lesssim 3$ in the
FOGGIE simulations and observations, further demonstrating the similar trends
of increasing metallicity with increasing $N_{\rm H\,I}$.161616The FOGGIE H I
absorbers can have metallicity as low as $\log(Z/Z_{\odot})=-8$ (the FOGGIE
metallicity floor). There are absorbers with $\log(Z/Z_{\odot})<-5$ classified
as SLFSs and pLLSs. These absorbers are not excluded from the CDFs in Fig. 28,
so these CDFs do not start at zero.  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.

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 substantial higher metallicity (the
FOGGIE median metallicity is 0.5 dex higher or a factor 3 higher than in the
R12-DLA sample).