THE INCIDENCE OF LOW-METALLICITY LYMAN-LIMIT SYSTEMS AT z ∼ 3.5: IMPLICATIONS FOR THE COLD-FLOW HYPOTHESIS OF BARYONIC ACCRETION^*

We obtained higher resolution spectra along 15 quasar sightlines (Table 2) selected as described above using the Magellan Echellete Spectrograph (MagE, Marshall et al. 2008) on the 6.5 m Magellan Clay telescope. MagE covers optical wavelengths from 3100 ${\rm{\mathring{\rm A} }}$ to 1 μm. At ${z}_{\mathrm{LLS}}$ = 3.3, the 912 ${\rm{\mathring{\rm A} }}$ Lyman break is redshifted to 3900 ${\rm{\mathring{\rm A} }}$. With a 085 slit, MagE has a resolution ${\mathcal{R}}=4950$ (or FWHM = 60.7 $\mathrm{km}\;{{\rm{s}}}^{-1}$). Observations were made on UT 2013 March 17/19 and UT 2013 May 05 with typical seeing of 06–08. The typical total integration time along a sightline is 3600 s, although there is some variation due to quasar magnitude and sightline observability.

Data were reduced using the MASE pipeline (Bochanski et al. 2009), using GJ 620.1 B/HIP 80300 as a standard for calibration. MASE is an IDL software package designed for reducing MagE data and performs the full extraction and calibration process. We manually construct a cubic-spline fit to the continuum of each spectrum.

Figure 1 shows portions of the SDSS and MagE spectra of J083832 around several of the LLS metal lines for comparison. In this example there is no statistically significant metal absorption from the LLS in the SDSS spectrum (FWHM ≈ 150 $\mathrm{km}\;{{\rm{s}}}^{-1}$) but the MagE spectrum has clear absorption lines.

Zoom In Zoom Out Reset image size

Several of the LLSs have no metal absorption in their corresponding MagE spectra. We observed one such object, J124957, with the Magellan Inamori Kyocera Echellette spectrograph (MIKE, Bernstein et al. 2003) at the same telescope on UT 2013 March 20 for 2400 s using a 1'' slit. MIKE is a double echelle spectrograph. The blue arm covers 3200–5100 ${\rm{\mathring{\rm A} }}$, and the red arm covers 4900 ${\rm{\mathring{\rm A} }}$ to 1 μm. With a 1'' slit, the blue/red arms have resolutions of 28,000/22,000 (FWHM = 10.7 $\mathrm{km}\;{{\rm{s}}}^{-1}$/13.6 $\mathrm{km}\;{{\rm{s}}}^{-1}$). All metal lines used in the current survey are in the red portion of the spectrum. MIKE data were reduced using MIKE Redux,⁵ a series of IDL tools that encompass all calibrations and extractions.

We also make use of an 1800 s medium-resolution infrared spectrum of the same object taken with the Folded-port Infrared Echellette (FIRE, Simcoe et al. 2013), on the 6.5 m Magellan Baade telescope, observed as part of a different survey. FIRE has a bandpass covering 0.8–2.5 μm, at a resolution of ${\mathcal{R}}=6000$ (FWHM = 50 $\mathrm{km}\;{{\rm{s}}}^{-1}$). For data acquisition and reduction details see Matejek & Simcoe (2012).

The primary motivation for the MIKE observation was the possibility of identifying deuterium absorption associated with the H i, which can indicate that gas is unprocessed (e.g., Crighton et al. 2013). Unfortunately, the H i absorption from the LLS along this sightline proved too broad to distinguish deuterium absorption. However, the higher resolution MIKE spectrum enables more sensitive column-density measurements, and the FIRE spectrum allows us to measure several ions not covered by the optical instruments.

We are interested in how our sample of metal-poor LLSs compares to the global LLS population. As a comparison sample, we studied spectra from the blind LLS survey of Fumagalli et al. (2013), also conducted with MagE. Since their sample is at slightly lower redshifts than ours, we selected the ten highest-redshift absorbers from their survey (excluding several DLAs and absorbers close to the QSO redshift) to achieve the best redshift overlap with our metal-poor sample. This method of choosing objects also avoids introducing any selection bias with regard to metallicity. The ten LLSs we examined have a median redshift of ${z}_{\mathrm{LLS}}$ = 3.04. For the remainder of the paper, we refer to this dataset as the "metal-blind sample," and to our observations as the "metal-poor sample." We applied identical analysis techniques to reduced spectra in both samples.

${N}_{{\rm{H}}\;{\rm{I}}}$ is notoriously difficult to measure in the LLS regime; the $\mathrm{Ly}\alpha$ curve-of-growth is flat and higher order transitions, including the Lyman break, are saturated by construction. We found that for the purpose of calculating metallicities, ionization modeling techniques can marginalize over a wide range of H i column densities, at least within the LLS range, as we discuss in Section 3.3.4.

Rather than attempting to find H i explicitly using Voigt profiles, we determined a range of viable column densities for each system, listed in Table 2. We used modified versions of IDL software from the XIDL ⁶ library in conjunction with the Voigt profile fitting packages VPFIT and RDGEN.⁷

For each system, the plausible range of ${N}_{{\rm{H}}\;{\rm{I}}}$ is determined by fitting various different aspects of the H i absorption signature (Figure 2). We estimated the Doppler b parameter, a measure of the width of the Voigt profile, using higher order Lyman series transitions. For weaker systems where the Lyman limit is not fully saturated, we were able to constrain ${N}_{{\rm{H}}\;{\rm{I}}}$ using the flux at the Lyman limit and/or fitting higher order Lyman series lines. For the strongest absorption systems, we fit the weak damping wings on the Lyα profile. For systems of middling strength, the saturation of the Lyman limit coupled with the non-existence of Lyα damping wings restricts the range of possible ${N}_{{\rm{H}}\;{\rm{I}}}$.

Zoom In Zoom Out Reset image size

The typical ${N}_{{\rm{H}}\;{\rm{I}}}$ uncertainty (defined as max(${N}_{{\rm{H}}\;{\rm{I}}}$) – min(${N}_{{\rm{H}}\;{\rm{I}}}$) for a given system) for both samples is 0.7 dex. This median total deviation is akin to an error bar of ±0.35 dex. The best constrained system had a total deviation of 0.3 dex, while the least constrained had a deviation of 1.7 dex, although the maximum is an outlier, with the next largest being 1.3 dex. These errors are incorporated into our metallicity uncertainty as bounds on a flat prior distribution of allowed ${N}_{{\rm{H}}\;{\rm{I}}}$, defining the range where we explore possible values for our solutions.

We measure column densities for ionic species using the apparent optical depth method (AODM, Savage & Sembach 1991). For each absorber, we integrated over a fixed velocity width in order to maintain consistency between different species/lines (rounded to the nearest pixel). For ions with multiple available lines, we performed an error-weighted average of the detections. The absorption profiles of the LLSs comprising our metal-poor and metal-blind samples are shown in Figures 3 and 4, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each line where there was a non-detection, we calculated a 3σ upper limit to the column density corresponding to the error in the spectrum over ±1 resolution element using 100,000 Monte Carlo realizations. For each iteration, we added to the flux a value drawn from a Gaussian distribution with width equal to the error at that pixel, then measured the column density of the resulting mock profile, discarding realizations where the column density was negative. If an equal amount of flux were scattered above and below the continuum level, this would result in a positive column density because the relationship between flux and apparent optical depth is exponential (see Fox et al. 2005). This produced a distribution of the largest possible column densities consistent with the observed lack of absorption. We adopted the column density larger than 99.7% of other trials as the 3σ upper limit. To mitigate the effects of poor continuum fits resembling absorption with low column density, we set the flux to unity and repeated the process if less than 50% of the trials produced a positive column density.

Similar to the H i analysis, we dismiss all metal absorption lines contaminated by interloping absorbers as well as lines obscured by large amounts of noise. The Si ii λ1260 transition, which is a powerful diagnostic for low-metallicity absorbers due to its large oscillator strength, was unavailable for ten of the 17 systems due to confusion with the Lyα forest.

All AOD column density measurements assumed that the absorption profiles reside on the linear portion of the curve-of-growth and are unsaturated. While this assumption is expected to hold for all of the lines in our metal-poor sample, there is no absorption strength cut for the metal-blind sample, so we need to test for saturation. In low- and medium-resolution spectra, absorption profiles can be saturated without clearly exhibiting zero flux, since the instrument blurs the absorption profile. In Table 1 we list the measured column densities, and below we discuss the identification of saturated lines.

Table 1.  Details on Specific Absorbers

Ion	λrest	$\mathrm{log}{N}_{\mathrm{AODM}}$	$\mathrm{log}{N}_{\mathrm{adopted}}$
Metal-poor Sample
J080853–070940 z = 3.545
C iv	1548	13.38 ± 0.13	13.39 ± 0.11
C iv	1550	13.42 ± 0.21	⋯
O i	1302	<13.65	<13.65
Si ii	1304	<13.31	<13.30
Si ii	1526	<13.45	⋯
Si iv	1393	12.93 ± 0.12	12.93 ± 0.12
Si iv	1402	<13.36	⋯
J083832+200142 z = 3.47595
Al ii	1670	12.73 ± 0.04	12.73 ± 0.04
C ii	1334	14.29 ± 0.02	>14.29
C iv	1548	13.65 ± 0.05	13.70 ± 0.04
C iv	1550	13.85 ± 0.07	⋯
Si ii	1526	13.49 ± 0.10	13.49 ± 0.10
Si iv	1393	13.52 ± 0.03	13.52 ± 0.03

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

3.2.1. Testing for Saturation

The primary interest of our study is the metallicity, $[{\rm{M}}/{\rm{H}}]$, of the absorbers. However, this requires knowledge of elemental abundance ratios, rather than the information on specific ions that we measure. Rather than assuming ionization conditions to convert between ionic and atomic column densities, we used the software package Cloudy (version 13.02, last described by Ferland et al. 2013) to solve for the ionization and metallicity simultaneously. With Cloudy, we modeled the ionization conditions of the absorbers over a range of metallicities, obtained ionic column densities for each model, and determined which models best matched the observed column densities using Monte Carlo simulations. Ionization conditions are typically described by the ionization parameter U, a proxy for hydrogen density ${n}_{{\rm{H}}}$, defined as

$$\begin{eqnarray}&&U=\displaystyle \frac{{\rm{\Phi }}}{{n}_{{\rm{H}}}c},\end{eqnarray} \tag{ 1 }$$

where c is the speed of light and ${n}_{{\rm{H}}}$ = n_{H i} + n_{H ii} + ${n}_{{{\rm{H}}}_{2}}.$ The flux of H i-ionizing photons, Φ, is given by

$$\begin{eqnarray}&&{\rm{\Phi }}=4\pi {\displaystyle \int }_{{\nu }_{\mathrm{LL}}}^{\infty }\displaystyle \frac{{J}_{\nu }}{h\nu }d\nu ,\end{eqnarray} \tag{ 2 }$$

where J_ν is the specific intensity of the incident radiation and ν_LL is the frequency corresponding to 1 Ryd.

Specifically, we used Cloudy to calculate the column densities of different species as a function of $\mathrm{log}U$ over a grid of H i column densities, redshifts, and metallicities. The models were generated assuming a plane-parallel geometry of a uniform and isothermal layer of photoionized gas, with the shape of the ionizing radiation spectrum derived from a combination of the cosmic microwave background (CMB) and the cosmic ultraviolet background (UVB) spectrum from Haardt & Madau (2012, CUBA software). The CMB is much weaker than the UVB at all relevant wavelengths, and omitting it caused no appreciable change. We adopted a solar relative abundance pattern for the Cloudy models.

The UVB spectrum includes emission from galaxies and QSOs, as well as a sawtooth absorption pattern due to the He ii Lyman series (Madau & Haardt 2009). Although observations of the hydrogen ionizing flux disagree with the normalization of this spectrum (Faucher-Giguère et al. 2008; Becker & Bolton 2013), Cloudy models adjust the normalization according to the input value of U, so this is inconsequential. Efforts to adjust the shape of the spectrum have found that best-fit models typically require fairly small modifications that ultimately translate to a difference in metallicity of ≲0.2 dex (Crighton et al. 2015), although there are outliers that are best fit with larger modifications to the shape. Since our modeling is based on a small number of ions, we keep the number of fit parameters to a minimum to avoid overfitting and do not let the shape of the spectrum vary in our models. Additionally, this ionizing background (without adjustment) was used in most of the works we compare to in Section 5.3, so any uncertainties in the shape of the spectrum should minimally affect comparison with other observations.

For each LLS, we interpolated on this grid to the absorber redshift, producing for each ionic metal species the column density as a function of $[{\rm{M}}/{\rm{H}}]$, $\mathrm{log}U,$ and ${N}_{{\rm{H}}\;{\rm{I}}}$. In Figures 5 and 6 we sketch the ionization modeling technique for LLSs with and without metal-line absorption, respectively. These examples assume a value of ${N}_{{\rm{H}}\;{\rm{I}}}$ intermediate to the allowed range, so one dimension is missing from this schematic representation of the modeling process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.1. Dependence of Results on N_{H i}

Although many systems require us to marginalize over a fairly broad range in ${N}_{{\rm{H}}\;{\rm{I}}}$, ionization modeling within the bounds we consider is surprisingly insensitive to the particular H i column density. Just as Figures 5 and 6 are projections onto an assumed value of ${N}_{{\rm{H}}\;{\rm{I}}}$, in Figure 7 we project along two different values of $\mathrm{log}U$ (at fixed $[{\rm{M}}/{\rm{H}}]$ = −2.5), to clarify how several properties vary with ${N}_{{\rm{H}}\;{\rm{I}}}$.

Zoom In Zoom Out Reset image size

In the top panel, we plot the model column density as a function of ${N}_{{\rm{H}}\;{\rm{I}}}$ for several different ions, as well as the total hydrogen column density ${N}_{{\rm{H}}}$, scaled to fit on the same plot. For a given value of $\mathrm{log}U,$ ${N}_{{\rm{H}}}$ does not scale very rapidly with ${N}_{{\rm{H}}\;{\rm{I}}}$; over the two orders of magnitude of ${N}_{{\rm{H}}\;{\rm{I}}}$ shown, ${N}_{{\rm{H}}}$ only changes by about 0.3 dex. The column density curves for C iv and Si iv are comparably flat. The variation in the low-ionization metals with ${N}_{{\rm{H}}\;{\rm{I}}}$ is more appreciable, but is still one order of magnitude smaller than the variation in ${N}_{{\rm{H}}\;{\rm{I}}}$. We do not show the scaling with $[{\rm{M}}/{\rm{H}}]$, although it is as expected—the metal column densities increase by an order of magnitude when $[{\rm{M}}/{\rm{H}}]$ is increased by one. Since ${N}_{{\rm{H}}}$ and most of the metal column density curves are relatively flat functions of ${N}_{{\rm{H}}\;{\rm{I}}}$, the model parameters U and $[{\rm{M}}/{\rm{H}}]$ corresponding to the model that matches measured ionic column densities from an LLS do not vary strongly with ${N}_{{\rm{H}}\;{\rm{I}}}$.

Comparing the model column density curves for the two different ionization parameters plotted, we see that the low ions are not strongly dependent on $\mathrm{log}U.$ C iv and Si iv column densities depend much more strongly on how ionized the gas is, as can also be seen in the right-hand plots of Figures 5 and 6.

We estimate from the Cloudy output how a metallicity measurement (at fixed $\mathrm{log}U$) based on a single ion varies with ${N}_{{\rm{H}}\;{\rm{I}}}$. For ion x corresponding to atom X, with measured column density N_x and model ionization fraction f_x = N_x/N_X, we can write the metallicity as

$$\begin{eqnarray*}[{\rm{X}}/{\rm{H}}] & = & \mathrm{log}{N}_{{\rm{X}}}/{N}_{{\rm{H}}}-\mathrm{log}{({N}_{{\rm{X}}}/{N}_{{\rm{H}}})}_{\odot }\\ & = & \mathrm{log}\displaystyle \frac{{f}_{{\rm{H}}\;{\rm{I}}}{N}_{{\rm{x}}}}{{f}_{{\rm{x}}}{N}_{{\rm{H}}\;{\rm{I}}}}-\mathrm{log}{({N}_{{\rm{X}}}/{N}_{{\rm{H}}})}_{\odot }\end{eqnarray*}$$

where f_{H i} is the H i ionization fraction in the corresponding model. Noting that the measured column density is a constant, we can differentiate with respect to $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}:$

$$\begin{eqnarray*}&&\displaystyle \frac{\partial [{\rm{M}}/{\rm{H}}]}{\partial \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}}=\displaystyle \frac{\partial \mathrm{log}{f}_{{\rm{H}}\;{\rm{I}}}}{\partial \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}}-\displaystyle \frac{\partial \mathrm{log}{f}_{{\rm{x}}}}{\partial \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}}-1.\end{eqnarray*}$$

The derivative depends only on ionization model output, and is shown in the bottom panel of Figure 7 for several ions.

We first note that all of the slopes are negative, indicating that metallicity decreases with ${N}_{{\rm{H}}\;{\rm{I}}}$ as expected. The Si iv and C iv slopes approach zero above ${N}_{{\rm{H}}\;{\rm{I}}}$ ∼ 18, indicating that modeled metallicities should be very robust in this range of column densities—given N_{C iv} or N_{Si iv} and $\mathrm{log}U,$ the metallicity is independent of ${N}_{{\rm{H}}\;{\rm{I}}}$. Si ii and C ii have larger derivatives with ${N}_{{\rm{H}}\;{\rm{I}}}$, as expected from the column density curves, but weaker dependence on $\mathrm{log}U.$

This suggests that uncertainty of ${N}_{{\rm{H}}\;{\rm{I}}}$ is manageable. Treating the metallicity slope as a proxy for the model uncertainty in the resulting metallicity and integrating the C ii derivative from ${N}_{{\rm{H}}\;{\rm{I}}}$ = 17.5 to 18.5, where the model uncertainty is worst, results in only a ∼ 0.5 dex uncertainty in metallicity. In practice, the uncertainties are significantly less since models are fit using multiple ions.

3.3.2. Matching Ionization Models to LLSs

To compare how well different ionization models fit the data, we define the likelihood function

$$\begin{eqnarray*}&&{\mathcal{L}}([{\rm{M}}/{\rm{H}}],{N}_{{\rm{H}}\;{\rm{I}}},\mathrm{log}U)={{\rm{\Pi }}}_{i}{{\ell }}_{i}({N}_{i})\end{eqnarray*}$$

where ℓ_i(N_i) is the likelihood for each ion measured given N_i, the model column density obtained by interpolating the grid to the corresponding values of $[{\rm{M}}/{\rm{H}}]$, ${N}_{{\rm{H}}\;{\rm{I}}}$, and $\mathrm{log}U.$ For detections, ℓ_i is taken to be a Gaussian with a mean and standard error given by the AODM measurements. For upper and lower limits, we let ℓ_i be unity if the model column density is within the range allowed by the limits, and ℓ_i decays as a Gaussian with σ = 0.05 beyond the allowed range. In practice, we used $\mathrm{log}{\mathcal{L}}$ to avoid computational instabilities resulting from small likelihoods.

Motivated by the implementation of Crighton et al. (2015), we used the Python module emcee (Foreman-Mackey et al. 2013) to perform a Markov Chain Monte Carlo (MCMC) sampling of the parameter space. This approach allows for calculation of the posterior probabilities of the metallicity and ionization parameter while marginalizing over the possible values of ${N}_{{\rm{H}}\;{\rm{I}}}$. We used flat priors with $\mathrm{log}U,$ $[{\rm{M}}/{\rm{H}}]\in [-4,-1],$ and ${N}_{{\rm{H}}\;{\rm{I}}}$ within the viable range determined for each LLS. We ran 1000 iterations of 100 walkers sampling the parameter space, discarding the first 100 iterations as "burn-in" to allow the walkers to explore the full posterior distribution and to remove the signature of the walkers' initial conditions. We constructed the posterior distribution from the remaining 900 iterations. An example posterior distribution (for J123525) is shown in the top left of Figure 8.

Zoom In Zoom Out Reset image size

Although many systems required marginalization over a wide range in ${N}_{{\rm{H}}\;{\rm{I}}}$, often larger than an order of magnitude, we found that the ionic column densities and posterior probabilities are not strongly dependent on ${N}_{{\rm{H}}\;{\rm{I}}}$, consistent with expectations from the discussion in Section 3.3.1.

3.3.3. Metallicity Upper Limits

3.3.4. Applicability of Single-cloud Models

Larger values of the ionization parameter U generally correspond to lower derived abundances for the ions we considered, so systematically overestimated ionization parameters would drive our results to artificially low metallicities. Fumagalli et al. (2011a) compiled results from published high-redshift LLSs (z > 1.5, their Figure S6) and determine that all systems have $-3\lt \mathrm{log}U\lt -1,$ with most systems in the range of $-3\lt \mathrm{log}U\lt -2.5.$ These systems generally have smaller redshifts than our sample. For all of the absorbers with detected heavy-element absorption in our metal-poor sample, we find $-2.7\lt \mathrm{log}U\lt -1.9,$ with most clustered around $\mathrm{log}U=-2.3$ (Figure 11).

Zoom In Zoom Out Reset image size

There are several plausible explanations for the moderately larger ionization parameters we measured. The ionization parameters in Fumagalli et al. (2011a) were compiled from studies using varying modeling techniques and ionizing radiation spectra. In addition, the complete literature sample comprises a relatively small sample, roughly the same size as our survey, often with redshifts and H i column densities different than our sample, suggesting that they may not form a uniformly selected comparison group for our results. Indeed, Fumagalli et al. (2011a) used the values from the literature only to show that $\mathrm{log}U=-3$ is a justifiable lower limit. The absorbers in Crighton et al. (2015) are at z = 2.5 and have sub-LLSs H i column densities.

Lehner et al. (2013) studied absorbers with $16.2\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 18.5$ at z ≲ 1 and found $\mathrm{log}U=-3.3\pm 0.6.$ Werk et al. (2014) looked at LLSs located near L* galaxies at z ≲ 0.2 (overlapping the sample from Lehner et al. 2013) and measured a mean ionization parameter $\mathrm{log}U=-2.8.$ A trend of LLS ionization parameter moderately increasing with redshift is supported by the studies previously mentioned and by our sample, although the total sample size is small.

The hydrogen density, ${n}_{{\rm{H}}}$, for each LLS follows from the ionization parameter U from Equation (1). To calculate Φ, we used Equation (2) using the shape of the UVB spectrum from Haardt & Madau (2012), renormalized to match the observed H i photoionization rates from Becker & Bolton (2013). Under the uniformity assumption, the sizes of the gas clouds can be estimated via ℓ = ${N}_{{\rm{H}}\;{\rm{I}}}$/(χ_{H i}${n}_{{\rm{H}}}$), where the neutral hydrogen ionization fraction χ_{H i} is output by the ionization models. In Table 2 we list the range of densities corresponding to the ionization parameter posterior distribution, and the range of sizes that follow assuming the central value of $\mathrm{log}U$ over the range of viable ${N}_{{\rm{H}}\;{\rm{I}}}$. Typical densities are of the order of 10⁻³ ${\mathrm{cm}}^{-3}$, and typical lengths are a few tens to hundreds of kiloparsecs, with a median of 160 kpc.

A uniform cloud with the median estimated size of 160 kpc at z = 3.5 would have a velocity gradient of ∼60 $\mathrm{km}\;{{\rm{s}}}^{-1}$ due to Hubble expansion, significantly larger than the typical Lyα forest Doppler parameter of b ∼ 30 $\mathrm{km}\;{{\rm{s}}}^{-1}$ (Hu et al. 1995). However, this is also roughly the size of a MagE resolution element, and we routinely see absorption that is ∼3 resolution elements wide (e.g., Figure 5), suggesting that the LLSs we measure often have multiple absorbing components (as is sometimes evident from metal absorption profiles).

Since we make no attempt to measure the multiple component structure of the LLSs using the medium-resolution data, our estimated cloud sizes represent an aggregate measurement of all components contributing to each LLS. Allowing for variation in the ionization state/density, individual components may have sizes that are of the order of 1–10 times smaller than the aggregate cloud sizes; although the absorbers' aggregate sizes are large, their individual substructures will be smaller, and could plausibly have sizes consistent with expectations for accretion flows.

In Figure 12, we plot $[{\rm{M}}/{\rm{H}}]$ over the overdensities $\delta ={n}_{{\rm{H}}}/{\bar{n}}_{{\rm{H}}}$ of the absorbers, where ${\bar{n}}_{{\rm{H}}}$ is the cosmic mean baryon density at the redshift of the absorber: ${\bar{n}}_{{\rm{H}}}$ = Ω_bρ_c(1 + z)³/m_p, where Ω_b = 0.04 is the cosmic baryon density relative to ρ_c, the critical density of the universe, and m_p is the proton mass. The systems with detections have typical overdensities ranging from δ ∼ 10–100. This figure also portrays how conservative the limits are: the Type 1 upper limits correspond to limiting overdensities at least two-fold greater than the detections, implying that the metallicities may be significantly lower than the limits we adopted.

Zoom In Zoom Out Reset image size

The high-resolution MIKE spectrum of J124957 confirmed that the lack of metal lines in the MagE spectrum of the same object was not a resolution effect, and the IR FIRE spectrum likewise shows no absorption at expected locations. In Figure 13 we compare the SDSS, MagE, and MIKE spectra for J124957, as well as showing cutouts of several portions of the FIRE spectrum where absorption is expected. It is evident that the four-fold increase in resolution provided by MIKE does not reveal any weak lines in this particular case.

Zoom In Zoom Out Reset image size

We obtained stricter upper limits on the column density with the high-resolution spectrum, allowing us to measure a metallicity upper limit of $[{\rm{M}}/{\rm{H}}]$ < −2.90, 0.2 dex less than the upper limit measured with MagE. The limits we measure from the FIRE spectrum are not as strict, since it has lower resolution. When we use column density limits from both MIKE and FIRE, the ions measured with FIRE do not influence the result. Using only ions from the FIRE spectrum, we find an upper limit of $[{\rm{M}}/{\rm{H}}]$ < −2.60.

It remains to be seen whether high-resolution observations could reveal weaker lines in other examples, but in light of current sensitivities and to maintain consistency among the sample, we use the metallicity limit obtained from the MagE spectrum in subsequent analysis.

A complete description of the distribution of metallicities needs to incorporate information provided by both detections and upper limits. To that end, we employ survival analysis methods developed to deal with censored datasets containing a mixture of measurements and limits.

For univariate data, the Kaplan–Meier estimator provides a general, non-parametric maximum likelihood estimate of the population from which a censored sample was drawn. For details on the application of the Kaplan–Meier method to similar datasets, see Simcoe et al. (2004) and references therein.

Briefly, the Kaplan–Meier method constructs a cumulative distribution function (CDF) for the sample, handling ambiguities introduced by upper limits by including them in probability calculations only when they are guaranteed to be unambiguous. For example, an upper limit of $[{\rm{M}}/{\rm{H}}]$ < −2.5 is guaranteed to be less than $[{\rm{M}}/{\rm{H}}]$ = −2.4, but may not be less than $[{\rm{M}}/{\rm{H}}]$ = −2.6 and will not be treated as such. The resulting cumulative distribution is a piecewise function that remains constant at upper limits and jumps at detections.

The Kaplan–Meier method requires the sample to satisfy two criteria. First, the upper limits must be independent. This is clearly true for our sample, where each measurement is drawn from a different absorber and most are from different sightlines. Second, the probability of a measurement being censored must be uncorrelated with the value of the measurement itself. While our sample likely does not strictly meet this criterion, since lower metallicity systems are more likely to result in non-detections, there is a characteristic to our survey that preserves the randomness of the censoring: all targets were selected using the same criterion, so the priors on metallicity are uniform across the sample. Hence, the selection method should adequately randomize the censoring. Also, since an LLS observation resulting in a metallicity limit partially depends on the signal-to-noise ratio of the spectrum, the brightness of the quasar and observing conditions also serve as randomizing factors. From the discussion in Section 3.3.1, metal column densities generally do not depend strongly on ${N}_{{\rm{H}}\;{\rm{I}}}$, so ${N}_{{\rm{H}}\;{\rm{I}}}$ should not bias whether or not a system gives a metallicity upper limit. Since the metallicity measurements and limits we find are not segregated such that all limits fall below detections, the Kaplan–Meier method is applicable.

We calculate the Kaplan–Meier distributions for our samples using ASURV Rev 1.2 (Isobe & Feigelson 1990; Lavalley et al. 1992), which implements the methods discussed in Feigelson & Nelson (1985). The resulting cumulative distributions are shown in Figure 14.

Zoom In Zoom Out Reset image size

We also extrapolate our results for each sample independently to estimate the entire LLS population of z ∼ 3.73, the mean redshift of the metal-poor sample. Since the metal-blind sample is at lower redshift and metallicities evolve with redshift, we can apply a shift to the entire CDF to bring it to the same redshift as the metal-poor sample. Taking the mean slope of the IGM and DLA metallicity with redshift, $[{\rm{M}}/{\rm{H}}]$ ∝ 0.28z (see Section 5.1), we shift the metal-blind CDF by the difference between the mean redshifts of the samples.

Alternatively, we may estimate the full LLS CDF from the metal-poor sample if we assume that all systems with metal lines detected in SDSS have higher metallicity than systems without detected metal lines. This assertion is demonstrably false for some individual cases, but lacking the (forthcoming) fully unbiased set of H i-selected LLS metallicities, it can represent a first attempt at generalization. In this case we may construct the CDF for the entire SDSS sample as

$$\begin{eqnarray*}{P}_{\mathrm{SDSS}} & = & {{xP}}_{\mathrm{MagE}}+(1-x)\\ x & = & \displaystyle \frac{{\rm{Number}}\;{\rm{of}}\;{\rm{SDSS}}\;{\rm{LLSs}}\;{\rm{meeting}}\;{\rm{metal-blind}}\;{\rm{criteria}}}{{\rm{Number}}\;{\rm{of}}\;{\rm{SDSS}}\;{\rm{LLSs}}}.\end{eqnarray*}$$

In other words, we scale the CDF by the fraction of LLSs matching our selection criterion (x), then add to it the fraction of LLSs that do not. Since, under our assumption, all LLSs not meeting our criteria have metallicities larger than those that do, this approximates the CDF for metallicities in the range probed by our metal-poor sample. In our scan of the SDSS spectra, we found x = 0.48. We stress that variations in signal-to-noise ratio and ionization fractions may result in some SDSS LLSs that failed to meet our initial selection criteria having lower metallicity than others that did.

In Figure 14 we compare these extrapolations to the metal-poor CDF and to each other. They diverge at high metallicity, because the assumptions on the extrapolation from the metal-poor sample place an unrealistic floor on the corresponding CDF, which is most relevant at higher metallicity. They are in fairly good agreement for metallicities below $[{\rm{M}}/{\rm{H}}]$ = −2.5; both extrapolations suggest ∼20% of LLS at z ∼ 3.73 have $[{\rm{M}}/{\rm{H}}]$ < −2.5, a value roughly 1σ above the measured IGM abundance.

In order to assess whether or not low-metallicity LLSs are consistent with being cold flows, we compare the metallicity CDFs of both the metal-poor and metal-blind samples with a toy model parent CDF. Motivated by the bimodal metallicity distribution found at low redshift (Lehner et al. 2013), we assume a mixed Gaussian model where the absorbers are drawn from a combination of two different parent populations, one being the IGM (representing potential accretion flows) and the other having more highly enriched gas that has been polluted by a local host galaxy. We refer to absorbers drawn from the IGM distribution as cold-flow candidates (CFCs).

We assume that the parent distribution in $[{\rm{M}}/{\rm{H}}]$ for CFCs is the same as the IGM's, which we interpolate from the measurements of Simcoe (2011) to the mean redshift of the sample. Note that this study used the same ionizing background spectrum to measure metallicities, so any systematics from uncertainty in a particular realization of the UVB spectrum are common to both studies.

The parent $[{\rm{M}}/{\rm{H}}]$ distribution of enriched CGM gas is less clear, but for this study we associate this phase with DLA abundances. The exact physical structures giving rise to DLAs are neither fully understood nor expected to be uniform, but they are thought to be locally enriched. Since DLAs have systematically lower abundances than H ii regions measured in emission (e.g., Sanders et al. 2015), this is a conservative choice to represent the non-CFC branch (using larger metallicities for this branch would give a larger fraction of CFCs). We model enriched gas using a lognormal distribution with parameters given by DLA measurements, since DLAs are thought to originate from gas in galaxies (e.g., Rafelski et al. 2011). We use DLA metallicities from Rafelski et al. (2012, and references therein). Compared to LLSs, DLA metallicities have been analyzed more extensively since DLAs are predominantly neutral and tend to have small ionization corrections that do not require modeling. In addition, the Lyα damping wings allow for accurate H i column densities in moderate-resolution spectra.

To compare with the metal-poor (blind) sample we average over all DLAs between z = 3.26–4.37 (2.90–3.25). The model metallicity probability distribution is then

$$\begin{eqnarray*}&&p([{\rm{M}}/{\rm{H}}])={f}_{\mathrm{IGM}}\;{p}_{\mathrm{IGM}}([{\rm{M}}/{\rm{H}}])+(1-{f}_{\mathrm{IGM}}){p}_{\mathrm{DLA}}([{\rm{M}}/{\rm{H}}]),\end{eqnarray*}$$

where p_IGM and p_DLA are Gaussian metallicity distributions with parameters summarized in Table 3, and ${f}_{\mathrm{IGM}}$ is the fraction of the model distribution drawn from the IGM that we are estimating (i.e., the fraction of CFCs).

Table 3.  Metallicity Distribution Parameters

	z	$[{\rm{M}}/{\rm{H}}]$IGM	σIGM	$[{\rm{M}}/{\rm{H}}]$DLA	σDLA	${f}_{\mathrm{IGM}}$	68% c.i.	95% c.i.
Metal-poor	3.73	−3.36	0.8	−1.69	0.48	0.71	0.60–0.84	0.44–0.96
Metal-blind	3.04	−3.12	0.8	−1.49	0.52	0.48	0.34–0.62	0.18–0.79

Download table as:  ASCIITypeset image

We performed Kolmogorov–Smirnov tests to find which IGM fractions are allowed by the two measured distributions (Figure 15). We list the values of ${f}_{\mathrm{IGM}}$ allowed within 68% and 95% confidence intervals (c.i.) in Table 3. Performing a least-squares fit, we find a best fit to the CDFs with ${f}_{\mathrm{IGM}}$ = 0.71 and 0.48 for the metal-poor and metal-blind samples, respectively. We adopt these best-fit ${f}_{\mathrm{IGM}}$ values, with errors given by the 68% confidence intervals, for the remainder of the paper. We caution that the small sample sizes enable intrinsic variation within the LLS population to appreciably influence the results.

Zoom In Zoom Out Reset image size

In addition to fitting the CDF, we also tested several fits to the measured metallicities directly and found approximately the same maximum likelihood values and ranges for ${f}_{\mathrm{IGM}}$. We compared mock histograms of metallicities drawn from mixed Gaussian models to the observed histogram, and performed an MCMC analysis using a likelihood function similar to that used for ionization modeling, with Gaussians for measured metallicities and one-sided Gaussians for limits.

Since 48% of the SDSS sample meets our metal-poor criteria, assuming only systems passing our initial cuts can be CFCs implies that the range of acceptable values for ${f}_{\mathrm{IGM}}$ for the entire z ∼ 3.7 LLS population is 0.34 ± 0.06. This is somewhat less than ${f}_{\mathrm{IGM}}$ for the metal-blind sample, suggesting that the assumption regarding the SDSS sample may be questionable, as discussed in Section 4.3, although sample variance may also account for the disagreement.

In Figure 16, we show the best-fit CDF and 68% confidence intervals for both samples, as well as the best-fit probability distribution function with the relative contributions from the enriched and unenriched parent populations.

Zoom In Zoom Out Reset image size

Of particular interest is Ω_CFC: the mass fraction of the universe (relative to the critical density) at z ∼ 4 contained in LLSs with IGM metallicities (i.e., CFCs). This quantity can be compared to simulations to test whether the global mass contained in cold flows agrees. Using our measurements and ionization models we make an order-of-magnitude estimate of this quantity for comparison with simulations.

The ratio of the CFC mass density to the cosmological critical density (at z = 0), ρ_c, is given by

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{\mathrm{CFC}}=\displaystyle \frac{{H}_{0}}{c}\displaystyle \frac{\mu {m}_{{\rm{H}}}}{{\rho }_{{\rm{c}}}}{\displaystyle \int }_{0}^{\infty }{N}_{{\rm{H}}}({N}_{{\rm{H}}\;{\rm{I}}}){f}_{\mathrm{CFC}}({N}_{{\rm{H}}\;{\rm{I}}}){{dN}}_{{\rm{H}}\;{\rm{I}}}\end{eqnarray*}$$

where ${N}_{{\rm{H}}}$(${N}_{{\rm{H}}\;{\rm{I}}}$) is the total hydrogen column density of an absorber, f_CFC(${N}_{{\rm{H}}\;{\rm{I}}}$) is the frequency distribution of CFCs as a function of neutral hydrogen column density, μ is the reduced mass of the gas, and m_H is the mass of the hydrogen atom.

To compute the integral, we need to assume a form for f_CFC(${N}_{{\rm{H}}\;{\rm{I}}}$). POW10 find f(${N}_{{\rm{H}}\;{\rm{I}}}$) for z = 3.7 LLSs can be fit by

$$\begin{eqnarray*}{f}_{\mathrm{LLS}}({N}_{{\rm{H}}\;{\rm{I}}})=\left\{\begin{array}{ll}{10}^{-4.85}{N}_{{\rm{H}}\;{\rm{I}}}^{-0.8} & 17.5\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\leqslant 19\\ {10}^{2.75}{N}_{{\rm{H}}\;{\rm{I}}}^{-1.2} & 19\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 20.3.\end{array}\right.\end{eqnarray*}$$

To estimate f_CFC(${N}_{{\rm{H}}\;{\rm{I}}}$), we multiply f_LLS by the fraction of LLSs that are CFCs, ${f}_{\mathrm{IGM}}$. (Note ${f}_{\mathrm{IGM}}$ is the fraction of LLS that are CFCs and f_CFC is the frequency distribution of these LLSs.) We take ${f}_{\mathrm{IGM}}$ = 0.34 ± 0.06 from Section 5.1.

We assume a reduced mass μ = 1.3, appropriate for absorbers with 75% H and 25% He by mass. The total hydrogen column density is readily found from the H i column density N_H = (n_H/n_{H i})${N}_{{\rm{H}}\;{\rm{I}}}$, with the ratio of ionized-to-total hydrogen coming from the ionization model. We use the median value derived from ionization solutions in both of our samples: $\langle {n}_{{\rm{H}}\;{\rm{I}}}/{n}_{{\rm{H}}}\rangle =0.0063.$ The median value for each sample separately is similar. Due to our small sample size, we cannot adequately measure and include in the calculation any variation in f_CFC and $\langle {n}_{{\rm{H}}}/{n}_{{\rm{H}}\;{\rm{I}}}\rangle$ with H i column density.

We also need to set bounds on the integral. For the lower bound, we use $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=17.5,$ the stated sensitivity limit of the LLS survey in POW10. Although lower column density absorbers are more numerous, larger column density absorbers dominate the mass density, so the result is largely insensitive to the lower bound; this may not be the case in calculations based on larger datasets, where the H i neutral fraction may be found to vary with ${N}_{{\rm{H}}\;{\rm{I}}}$. As an upper bound we take $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=19.5,$ a typical value where systems transition from being considered LLSs to sub-DLAs.

With these parameters, we obtain Ω_CFC = 0.0017. Comparing this to the cosmic baryon density, Ω_b = 0.04, we find roughly 5% of baryons at z ∼ 3.7 are contained in CFC LLSs. Note that this calculation is sensitive to both the maximum column density set for the integral and the H i ionization fraction; increasing (decreasing) the maximum ${N}_{{\rm{H}}\;{\rm{I}}}$ by 0.2 dex increases (decreases) the result by a factor of ∼1.5, and the H i ionization fractions for individual systems can vary from the median by a factor of ∼5. As we discuss in Section 5.4, our result is fairly consistent with simulated results.

In Figure 17, we compare the metallicities of our samples with DLA metallicities from Rafelski et al. (2012, and references therein). DLAs generally have higher metallicities than both our metal-poor and metal-blind samples, and there are suggestions that DLAs have a metallicity floor, which many LLSs are below. It is clear from this comparison that LLSs and DLAs at high redshift differ significantly in their metallicity distributions.

Zoom In Zoom Out Reset image size

Also shown are low-redshift LLSs from Lehner et al. (2013, and references therein). They categorized their H i systems as LLS ($16.2\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 19$) or super LLS (SLLS, $19\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 20.3$). Their LLS sample is mostly composed of systems with $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 17.5,$ the cutoff for our survey, so some differences are expected.

The low-redshift LLS population clearly exhibits a metallicity bimodality, with most of the lower-metallicity branch below most DLAs, although the difference between the LLS and DLA populations is not as emphasized as at higher redshift. While a bimodal model fits our high-redshift sample well, the two populations blend together more smoothly than they do at lower redshift, where Lehner et al. see very few systems at intermediate abundances. This may be a result of many of the lower abundances at high redshift producing upper limits rather than measurements.

Several LLSs drawn from the literature (Table 4) that exhibit low metallicities and are claimed as potential evidence of cold flows are also included in Figure 17. Our work corroborates the finding of low-metallicity LLSs and provides some statistical context for the population from which they are drawn. The two high-redshift metallicity upper limits are from Fumagalli et al. (2013). Using high-resolution spectra, they were able to model and subtract the Lyα forest to obtain column-density upper limits for C iii λ977 and Si iii λ1206, which provide tighter constraints than the ions available in our medium-resolution spectra (see their Figure S5).

Table 4.  Low-metallicity LLSs in the Literature

QSO	${z}_{\mathrm{LLS}}$	$[{\rm{M}}/{\rm{H}}]$	Source
PG1630+377	0.274	−1.71 ± 0.06	Lehner et al. (2013)
J144535+291905	2.44	−2.0 ± 0.17	Crighton et al. (2013)
HE 0940–1050	2.917	−2.93 ± 0.13	Levshakov et al. (2003)
J113418+574204	3.411	< −4.2	Fumagalli et al. (2013)
Q0956+122	3.096	< −3.8	Fumagalli et al. (2013)

Download table as:  ASCIITypeset image

By comparing our sample with structures having analogous properties in simulations, we explore the agreement between simulations and observations and gain insight into the nature of metal-poor LLSs. Fumagalli et al. (2011b) simulated absorption profiles produced by cold flows at z ∼ 2.3 and found that much of the gas is ionized by the UV background, appearing mostly as LLSs. They determined that DLAs have higher metallicities than LLSs and SLLSs, with DLA metallicities fairly consistent with observed systems, and the authors suggested metal-poor LLSs may therefore be an observational signature of cold-flow accretion. Their simulations predict that most of the cold-flow observational signatures are LLSs with $17\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 18,$ with a peak metallicity for cold-stream LLSs of one-hundredth solar.

While our measured metallicities tend to be somewhat lower, enrichment of the IGM between z = 3.7 and z = 2.3 may account for some of the difference. An additional difference is that our metal-poor sample is mostly composed of absorbers with $18\lt \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 19,$ the higher end of the LLS column-density distribution.

A limitation of the simulations employed in Fumagalli et al. (2011b) and similar studies is the simulated volume. In Fumagalli et al. (2011b), zoom-in simulations of seven halos were considered. This restricts analysis of the covering fraction of LLSs to roughly a few times the virial radius of each galaxy. Although these simulations allowed for a descriptive picture of the neutral-gas content immediately around the seven relatively massive galaxies presented, it remains unclear if a random quasar sightline is more likely to intersect a LLS in the immediate vicinity of one of these systems or elsewhere in the cosmic web. Such questions can only be addressed with larger simulation volumes. Full-volume cosmological simulations can complement the analysis of Fumagalli et al. (2011b) by probing the full range of LLSs that are probed through quasar-selected samples.

To demonstrate this point, we briefly consider the global distribution of neutral hydrogen in the full-volume cosmological simulations presented in Bird et al. (2014). These simulations were run using the cosmological hydrodynamical simulation code AREPO (Springel 2010) in a periodic box of size 10 h⁻¹ Mpc. The simulations contain 512³ dark matter particles and a similar number of baryon resolution elements yielding a mass resolution of 1.4 × 10⁵ M_⊙—about an order of magnitude larger than that presented in Fumagalli et al. (2011b). The simulation physics is the same as in the Illustris Simulation (Vogelsberger et al. 2014), which importantly includes star formation-driven winds at a level that allows for appropriate evolution of the galaxy stellar mass function (Genel et al. 2014; Torrey et al. 2014). Neutral hydrogen fractions are obtained assuming a uniform UV background, with self-shielding corrections (Rahmati et al. 2013). Complete simulation and post-processing details are presented in Bird et al. (2014).

We examine the simulations by applying a similar technique to that used in both Fumagalli et al. (2011b) and Bird et al. (2014). Specifically, we project the neutral hydrogen and mass-weighted average metallicity onto a two-dimensional grid. We do this only along one projected direction but do not expect our results to be modified if we considered other projections. We employ a grid of 16,000 by 16,000 cells, which results in converged neutral-hydrogen column-density distribution functions. Projecting the full 10 Mpc box onto a single grid could boost the LLS number count by adding multiple, well-separated systems of low column density. To minimize this effect, we use ten slices each with a thickness of 1 Mpc. We treat each pixel as an independent line of sight.

In Figure 18, we show a map of the neutral-hydrogen column density through one slice of the simulated box at z = 3.5, truncated at $\mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}=17.5$ (to emphasize LLSs), and the corresponding map of metallicities. Lower column density material appears as semi-transparent blue in the column-density map. There is LLS-level column density material tracing the cosmic web extending beyond the virial radii of galaxies in the simulated volume.

Zoom In Zoom Out Reset image size

Since we do not have explicit projected offsets for the observed quasar sightlines, it is hard to know which part of the IGM/CGM is being probed. It is possible that some fraction of our observed CFCs intersect material still in the IGM, well outside the halo virial radius. In those cases, it is not immediately clear over what timescale the observed neutral gas will fall through the virial radius nor whether it will stay neutral (or be shock heated) as it is accreted. Additionally, most of the simulated LLSs within a halo virial radius are enriched, with the metal-poor LLSs tending to trace the cosmic web. This suggests that either the feedback prescription in the simulation artificially contaminates inflowing material or much of our observed sample is inter-halo material, as opposed to intra-halo.

We can directly address global statistics of the simulated LLS population. Adopting standard column density limits ($17.5\leqslant \mathrm{log}{N}_{{\rm{H}}\;{\rm{I}}}\lt 19.5$), we find roughly 7% of the simulation baryons reside in LLSs, and it has been previously shown that the column density distribution function for neutral hydrogen in these models is reasonably consistent with observations (Bird et al. 2014). The simulated LLS CDF at redshift z = 3.5 is shown in Figure 19. Our metal-blind sample CDF (when shifted to z = 3.5)⁸ is in fairly good agreement with the simulated LLS CDF with a peak metallicity for LLSs around one-hundredth solar. Although we are sampling a significantly larger volume, this result is similar to that presented in Fumagalli et al. (2011b). An extrapolation of the full LLS population from the metal-blind sample is in agreement at low metallicities, but diverges at higher metallicities where the simple extrapolation is inadequate.

Zoom In Zoom Out Reset image size

We consider all LLSs with $[{\rm{M}}/{\rm{H}}]$ < −2.5 to be CFCs. The fraction of LLSs constituting CFCs is approximately the ratio of the LLS and CFC mass densities, Ω_CFC/Ω_LLS = 0.312. This is in accordance with the observational result for ${f}_{\mathrm{IGM}}$ extrapolated to the entire LLS population.

However, we find the derived Ω_CFC to be smaller for the simulations by a factor of ∼2. Given that the simulations have both a similar ${N}_{{\rm{H}}\;{\rm{I}}}$ distribution function (Bird et al. 2014) and a metal-poor fraction when compared with observations, it is likely this offset is driven by the applied ionization corrections. In our Ω_CFC calculation based on observations, we assumed that the hydrogen neutral fraction, n_{H i}/${n}_{{\rm{H}}}$, is independent of ${N}_{{\rm{H}}\;{\rm{I}}}$ due to an insufficient amount of data to measure any such relationship. Since the integrals are heavily weighted toward larger ${N}_{{\rm{H}}\;{\rm{I}}}$ systems, unaccounted-for variations in neutral fraction could systematically and significantly influence the calculation. Hence, the disagreement between simulated and observed values of Ω_CFC and Ω_LLS may be resolved by improved observational statistics and does not necessarily undermine the agreement on the fraction of LLSs that are metal-poor.

Additionally, we suspect the prevalence of metal-poor LLSs in the simulations can be influenced by (i) the specifics of the adopted feedback model and (ii) the mass and spatial resolution. Further investigations on both of these fronts are warranted, alongside developing an understanding of the dependence of the hydrogen neutral fraction on ${N}_{{\rm{H}}\;{\rm{I}}}$. Full-volume cosmological simulations that are able to simultaneously reproduce the ${N}_{{\rm{H}}\;{\rm{I}}}$ distribution function as well as the low-metallicity tail of the LLS distribution presented in this paper will be helpful in identifying the fraction of low-metallicity LLSs residing within halo virial radii and, thus, the true mass density of cold flows.