Low-z Lyman Limit Systems in the FIRE Simulations

Low-Redshift Lyman Limit Systems as Diagnostics of Cosmological Inflows and Outflows


We use cosmological hydrodynamic simulations with stellar feedback from the FIRE project to study the physical nature of Lyman limit systems (LLSs) at . At these low redshifts, LLSs are closely associated with dense gas structures surrounding galaxies, such as galactic winds, dwarf satellites, and cool inflows from the intergalactic medium. Our analysis is based on 14 zoom-in simulations covering the halo mass range M at , which we convolve with the dark matter halo mass function to produce cosmological statistics. We find that the majority of cosmologically-selected LLSs are associated with halos in the mass range M. The incidence and HI column density distribution of simulated absorbers with columns cm are consistent with observations. High-velocity outflows (with radial velocity exceeding the halo circular velocity by a factor ) tend to have higher metallicities () while very low metallicity () LLSs are typically associated with gas infalling from the intergalactic medium. However, most LLSs occupy an intermediate region in metallicity-radial velocity space, for which there is no clear trend between metallicity and radial kinematics. The overall simulated LLS metallicity distribution has a mean (standard deviation) [X/H] (0.4) and does not show significant evidence for bimodality, in contrast to recent observational studies but consistent with LLSs arising from halos with a broad range of masses and metallicities.

galaxies: formation — galaxies: evolution — galaxies: haloes — quasars: absorption lines — intergalactic medium — cosmology: theory

1 Introduction

The circum-galactic medium (CGM) plays a vital role in determining the evolution of galaxies. Accretion from the intergalactic medium (IGM) passing through the CGM provides the fuel to sustain star formation in galaxies over a Hubble time (e.g. Kereš et al. 2005; Prochaska & Wolfe 2009; Bauermeister et al. 2010). At the same time, galactic feedback processes driven by stars and active galactic nuclei are observed to expel gas out of galaxies and into the CGM as outflows (e.g., Heckman et al. 2000; Shapley et al. 2003; Martin 2005; Weiner et al. 2009; Steidel et al. 2010; Greene et al. 2011; Jones et al. 2012; Rubin et al. 2014; Cicone et al. 2014). Cosmological models of galaxy formation furthermore show that strong outflows are necessary to regulate galaxy growth (e.g., Kereš et al. 2009; Anglés-Alcázar et al. 2014; Hopkins et al. 2014; Somerville & Davé 2015) and explain the metal enrichment of the IGM (Aguirre et al. 2001; Oppenheimer & Davé 2006; Booth et al. 2012). At low redshift, Lyman limit systems (LLSs; loosely defined as systems optically thick at the Lyman limit) trace overdense structures in the halos of galaxies and are thus useful observational probes of inflows and outflows (e.g., Ribaudo et al. 2011; Tripp et al. 2011).

In the past few years, the observational census of LLSs has been greatly improved by new analyses of Hubble Space Telescope (HST) observations. Ribaudo et al. (2011) measured the incidence of LLS and column density distribution of HI absorbers, using archival data from the Faint Object Spectrograph (FOS) and Space Telescope Imaging Spectrograph (STIS) instruments. Thanks to the installation of the Cosmic Origins Spectrograph (COS) on HST, the increased sample of low-redshift HI-selected LLSs makes it now possible to measure the unbiased metallicity distribution of LLSs. Lehner et al. (2013) measured this metallicity distribution and found evidence for a metallicity bimodality. This measurement included all systems with neutral hydrogen column density . Here and in the rest of this paper, logarithms are to the base 10, is in units of cm, and is in units of . Lehner et al. (2013) tentatively identified the metal-rich branch as likely tracing galactic winds, recycled outflows, and gas tidally-stripped from galaxies. They also noted that the metal-poor branch has properties consistent with the cold accretion streams predicted by cosmological simulations (e.g., Kereš et al. 2005, 2009; Dekel & Birnboim 2006; Dekel et al. 2009; Brooks et al. 2009; van de Voort et al. 2011; Faucher-Giguère et al. 2011). In an analysis that doubles the Lehner et al. absorber sample, Wotta et al. (2016) find broadly consistent evidence for a metallicity bimodality, but note that the bimodality primarily arises from the (more numerous) partial LLSs with column in the range . Quiret et al. (2016) reported hints of a bimodality in the metallicity distribution of super LLSs (SLLSs) with (which they refer to as sub damped Ly absorbers) at , but with lower statistical significance and with no evidence of a metallicity bimodality at higher redshift. Fox et al. (2013) found evidence that metal-rich low- LLSs tend to show higher O VI columns and broader O VI profiles than metal-poor LLSs.

At the same time, the incidence and metallicity properties of high-redshift LLSs have become increasingly well constrained by ground-based observations (Fumagalli et al. 2011, 2013, 2016; Lehner et al. 2014; Cooper et al. 2015; Prochaska et al. 2015). Despite this tremendous observational progress, forward modeling of IGM absorbers using cosmological simulations remains critical to understanding the physical nature of LLSs, since spectroscopy only provides direct information on the line-of-sight structure of absorbers.

Overall, cosmological simulations have proved very successful at reproducing the observed column density distribution of HI absorbers (Hernquist et al. 1996; McQuinn et al. 2011; Altay et al. 2011; Altay et al. 2013; Rahmati et al. 2013; Gurvich et al. 2016), demonstrating how this technique can be used to inform the physical nature of absorption systems. In general, the agreement between cosmological simulations and observations is best for the Lyman- forest (e.g., Bird et al. 2013), which is relatively little affected by feedback processes from galaxies (e.g., Kollmeier et al. 2006). On the other hand, simulations predict that a large fraction of circum-galactic LLSs trace galactic winds (Fumagalli et al. 2014; Rahmati et al. 2015; Faucher-Giguère et al. 2015; Liang et al. 2016; Faucher-Giguère et al. 2016) and the distribution of LLSs is therefore sensitive to feedback processes. Most of the other circum-galactic LLSs arise in streams of infalling gas (e.g., Faucher-Giguère & Kereš 2011; Fumagalli et al. 2011; van de Voort et al. 2012; Goerdt et al. 2012; Faucher-Giguère et al. 2015) or on the outskirts of galaxies. To date, however, most simulation analyses have focused on the properties of LLSs at high redshift (see also Kohler & Gnedin 2007).

Our goal in this paper is to study the properties of simulated LLSs at , with an eye toward interpreting recent HST measurements and developing inflow-outflow diagnostics. We use a set of cosmological zoom-in simulations from the FIRE project (“Feedback In Realistic Environments”)1, which implement models for stellar feedback on the scale of star-forming regions. These simulations have been shown to correctly reproduce the star formation histories of galaxies below (Hopkins et al. 2014), mass-metallicity relations (Ma et al. 2016) at all redshifts where observations are available, and LLS covering fractions in the halos of galaxies (Faucher-Giguère et al. 2015; Faucher-Giguère et al. 2016). In these simulations, galaxy-scale outflows are self-consistently generated from the local injection of feedback momentum and energy on small scales (Muratov et al. 2015, 2016). These successes of the FIRE simulations make them particularly well suited to address the physical nature of LLSs.

For this paper, we define partial LLSs as systems with , LLSs as systems with , and SLLSs as systems with . Since partial LLSs and LLSs trace similar structures in the CGM, we will often group them together and refer to them simply as LLSs. Throughout, we assume a standard flat CDM cosmology with , , and , consistent with the latest observational constraints (e.g., Planck Collaboration et al. 2015).

The structure of this paper is as follows. In §2, we describe in more detail the simulations used in our analysis and our methodology to extract observational statistics. In section 3, we convolve our zoom-in simulations results with the dark matter halo mass function to predict cosmological statistics, including the incidence of LLSs and the HI column density distribution of CGM absorbers. We quantify relationships between LLS column density, metallicity, and radial velocity relative to central galaxies in §4. In that section, we also compute the overall metallicity distribution of low-redshift LLSs predicted by our simulations. We discuss our results and conclude in §5. The Appendix summarizes convergence tests done to validate our analysis.

2 Methodology

2.1 Simulation sample

Our analysis is based on cosmological zoom-in (e.g., Porter 1985; Katz & White 1993) simulations from the FIRE project. As in previous FIRE papers, all SPH simulations presented in this paper were run with exactly the same code as the original FIRE simulations presented in Hopkins et al. (2014). Therefore, we only summarize the key points here. Our CGM analysis methods are similar to the ones used in Faucher-Giguère et al. (2015) and more details can be found in that paper.

Our simulations were run with the P-SPH (Hopkins 2013) hydrodynamics solver implemented in the GIZMO code (Hopkins 2015). P-SPH uses a pressure-entropy formulation of the smooth particle hydrodynamics (SPH) equations that eliminates the artificial surface tension at contact discontinuities found in traditional density-based SPH formulations (e.g., Agertz et al. 2007) and resolves the major historical differences between SPH and grid-based codes. The gravity solver in GIZMO is a modified version of the GADGET-3 gravity algorithm (Springel 2005) which implements the adaptive softening method of Price & Monaghan (2007) (which we use for the gas) and a modified softening kernel that represents the exact solution for the potential of the SPH smoothing kernel following Barnes (2012).

In our simulations, gas is allowed to cool to K via atomic and molecular line emission, in addition to the standard processes described by Katz et al. (1996). The FIRE simulations are designed to resolve the Toomre mass/length of gas within galaxies, and therefore the most massive giant molecular clouds (GMCs) in which most star formation occurs in typical galaxies (e.g., Murray 2011). In the simulations, star formation occurs only in dense, locally self-gravitating, and molecular/self-shielding gas. Where all these criteria are met, the gas is converted into stars on a local free-fall time. The density threshold for star formation is set cm in most of our simulations, though in practice the self-gravity criterion is generally the most stringent. Stellar feedback is modeled by implementing energy, momentum, mass, and metal return from radiation pressure, photoionization, photoelectric heating, supernovae, and stellar winds following the STARBURST99 stellar population synthesis model (Leitherer et al. 1999). We explicitly follow the chemical abundances of nine metal species (C, N, O, Ne, Mg, Si, S, Ca, and Fe). During the hydrodynamic simulations, ionization balance of all elements is computed using the cosmic ultraviolet (UV) background model of Faucher-Giguère et al. (2009) and we apply an on-the-fly correction for dense, self-shielded gas based on a local Jeans-length approximation.

(M) (M) (kpc) (M) (pc) (M) (pc)
m11.4a 2.6e11 6.2e9 180 3.3e4 9 1.7e5 140
m11.9a 8.4e11 3.0e10 250 3.4e4 9 1.7e5 140
MFz0_A2 1.0e13 2.7e11 630 3.0e5 9 1.4e6 140

and are the dark matter and stellar masses inside the present-day virial radius, , defined as in Bryan & Norman (1998). and are the baryonic and dark matter particle masses. The simulations use adaptive gravitational softening lengths for the gas but fixed softening lengths for the dark matter. is the minimum force softening length for the gas and is Plummer-equivalent gravitational softening length for the dark matter.

Table 1: Parameters of the New Simulations

We analyze a sample of 14 simulations run to and spanning the halo mass range . From previous work, our sample includes the seven main simulations in Hopkins et al. (2014) (m09, m10, m11, m12v, m12q, m12i, and m13) and four simulations from Chan et al. (2015) (m10h1297, m10h1146, m10h573, and m11h383). We ran three new simulations (m11.4a, m11.9a, MFz0_A2) with to better populate that halo mass range for our analysis. Simulation MFz0_A2 follows the same main halo as the MFz2_A2 MassiveFIRE simulation (Feldmann et al. 2016), but include a larger high-resolution region and is evolved to . The simulation parameters for the three new runs are detailed in Table 1. With the exception of m13 and MFz0_A2, all simulations in our sample have a gas particle mass  M, minimum adaptive baryonic smoothing/force softening length , dark matter particle mass  M, and dark matter force softening length . For dwarf galaxies in the sample, the resolution parameters are significantly better. For example, the m09 and m10 runs have a gas particle mass M and correspondingly smaller minimum softening lengths.

2.2 Extraction of absorber statistics

For each simulation, we center our analysis on the most massive main halo within the zoom-in region, and include 20 snapshots evenly spaced in redshift over and 11 evenly spaced snapshots over . We use Amiga Halo Finder (AHF) (Knollmann & Knebe 2009) and the virial overdensity definition from Bryan & Norman (1998) to identify halos. We denote the virial radius . Figure 1 shows HI and metallicity maps for six of our simulated halos at . From the panels on the right, which show velocity vectors overlaid on top of the metallicity maps, we anticipate a complex relationship between instantaneous velocity and metallicity in the CGM of our simulated halos at this redshift, a point to which we will return to quantitatively in §4.

Figure 1: Left: HI column density maps at for simulations m10, m10h573, m11h383, m11.9a, m12i, and m13 in reading order from top left to bottom right. Black contours indicate cm (partial Lyman limit systems and higher). Right: Same as left, but metallicity in solar units for the same halos as on the left. The arrows indicate projected velocity vectors in units of the halo circular velocity. The HI column density and metallicity shown in these maps are for the strongest absorber along each sight line, calculated as described in §2.2. Overall, the circum-galactic gas dynamics and their relationship to gas metallicity are complex.

To extract the properties of individual CGM absorbers, we first grid the SPH particle data onto cubic 3D meshes, each with cells. The full SPH kernel is used to interpolate particle data onto the grid. In cases where the grid does not resolve individual particles, the cloud-in-cell method is used to distribute the particle mass to the nearest 8 grid points in a mass-conserving manner. We verified the convergence of our results with respect to grid resolution (see Appendix A). At , our simulations indicate that the vast majority of LLSs occur within the virial radius of halos. To capture this halo cross section, we use grids of side length for each main halo. In any given snapshot, we find that up to of LLSs in our grids are found away from the central galaxy, but almost all are associated with nearby galaxies that have their own main halos. Thus, such LLSs are accounted for by our convolutions over the dark matter halo mass function described in §3. The HI density for each cell is computed using a self-shielding-corrected photoionization rate (Rahmati et al. 2013), which is a function of the total hydrogen density and the ultraviolet background. We do this by post-processing our simulated grids, but have verified that computing the HI density on the particles and then depositing that density into cells gives similar results.

Figure 2: Example sight line profiles, as a function of line-of-sight velocity, used to identify and quantify the properties of CGM absorbers. This is a random sight line through the main halo of our m11.9a simulation, with zero LOS velocity corresponding to the galaxy velocity. Top: of HI column density. The horizontal solid black line indicates the average value of over the sight line, and the dashed lines show intervals. Middle: Metallicity profile in solar units. Bottom: Gas mass-weighted radial velocity relative to the central galaxy. For each sight line, we select the strongest HI absorber (highlighted in red) using the shape of the profile and extract the corresponding metallicity and radial velocity.

For each line of sight (LOS), we extract the properties of the strongest HI absorber identified in velocity space. To select the strongest HI component along each LOS, we use the LOS profile binned as a function of LOS velocity. The strongest component is identified by the peak along the LOS, and its width is set by the points at which drops below its mean value along the LOS. We experimented with different quantitative criteria for automatically identifying strongest absorbers along each LOS but found that there is in general no ambiguity in our analysis. This is because there is typically a single strongest absorber along each LOS, well separated in velocity space from other column density peaks. This approach is a good proxy for how LLSs are identified in observations (e.g., Lehner et al. 2013) and is more accurate than, for example, estimating absorber metallicities by projecting the total gas and metal masses across the simulation volume. The latter approach can be significantly biased when the metallicity is not uniform along the LOS, which is often the case in the multiphase CGM. We do not analyze the properties of weaker absorbers. Figure 2 illustrates this procedure for a random LOS through the m11.9a halo at .

For each absorber identified in this way, we evaluate the HI column density , metallicity [X/H], and the gas mass-weighted radial velocity of the component relative to the central galaxy. We define the metallicity [X/H] of an absorber as , where is the metal mass fraction of the absorber gas and is the solar metal mass fraction (Asplund et al. 2010).

3 Cosmological H I Statistics

Our zoom-in simulations provide us with predictions for the distribution of absorbers within the virial radius of galaxies as a function of halo mass and redshift. Observations that randomly select quasars on the sky, on the other hand, probe a cosmological distribution of absorbers that receives contributions from halos of all masses. In this section, we convolve our zoom-in results with the dark matter halo mass function to derive predictions for the cosmological incidence and column density distribution of HI CGM absorbers.

3.1 Weighting halos as a function of mass and redshift

Figure 3: Top: Incidence of dark matter halos (virial cross sections per sight line) per unit redshift and halo mass bin. Middle: Covering fraction of LLSs () within a projected virial radius as a function of halo mass, as predicted by our zoom-in simulations. A redshift-dependent fit to the covering fractions is shown at and in solid blue and red respectively, with the shaded regions showing the standard deviation of the fit. Bottom: Incidence of LLS per unit redshift per unit . The majority of randomly selected LLS are associated with halos in the mass range . In each panel, the solid dots correspond to the covering fractions averaged over three orthogonal sky projections and the error bars (usually smaller than the solid dot) show the sample standard deviation among the different projections.

Consider a prescribed range. We define to be the mean number of absorbers in this range per sight line per unit redshift. For a given redshift , we quantify the contribution of halos of different mass bins to the total using the weights


We evaluate as follows.

For any main halo, we define the halo cross section as . The standard dark matter halo mass function provides the mean number density of halos per halo mass bin as a function of redshift, . The three-dimensional halo number density can be converted into a surface density per unit redshift using the cosmological line element, . We can then simply evaluate the mean number of halos intersected per unit redshift per halo mass bin for random sight lines:


The top panel of Figure 3 shows this quantity at and 1. For the halo mass function, we used the CPMSO+AHF fit to -body cosmological simulations from Watson et al. (2013). As expected, the curve increases with decreasing redshift (as halos grow) and decreases with increasing halo mass (as halos become rarer).

To predict the contribution of different halos to the cosmological incidence of absorbers, we need to know the fraction of the halo cross section covered by the absorbers of interest. We define this quantity as the covering fraction and use our zoom-in simulations to predict its values. For each simulation snapshot, we evaluate the covering fraction by extracting absorbers along one LOS for each grid pixel with impact parameter . To capture viewing angle variance, we compute covering fractions for three orthogonal sky projections for each snapshot. For each halo, we plot in the middle panel of Figure 3 the average covering fraction and the estimated standard deviation among the projections using an error bar. For this figure we include all systems with to match the observational analysis of Lehner et al. (2013). This panel shows that the LLS covering fraction versus halo mass at peaks broadly for . A similar (but quantitatively different) trend for LLS covering fractions peaking in relatively massive halos is seen in the FIRE simulations at (Faucher-Giguère et al. 2015; Faucher-Giguère et al. 2016).

In all snapshots, the variance in (the error bars in Figure 3) owing to different viewing angles is small compared to the variance arising from time variability and halo-to-halo scatter. To demonstrate the strong time variability of the covering fractions in individual halos, we plot in Figure 4 as a function of redshift for the six halos shown in Figure 1. The curves in Figure 4 are based on the 31 simulation snapshots between and analyzed in this paper for each simulation. The time-dependence of the covering fractions arises from the time-dependence of cosmological inflows inherited from the evolution of large-scale cosmic structure and from the bursty nature of star formation and galactic winds predicted by our simulations (e.g., Muratov et al. 2015; Faucher-Giguère et al. 2015; Sparre et al. 2015). In addition to the time variability in individual halos, there could be some differences in time-averaged covering fractions between different halos, for example connected to the environments of halos. In §4.3, we perform a bootstrapping analysis to quantify the uncertainty in the metallicity distribution predicted using our set of zoom-in simulations arising from halo-to-halo scatter.

We fit our simulated covering fraction data to a function of halo mass with a redshift-dependent double power law of the form


Performing a least squares fit to all the data in the redshift interval (for ) results in the best-fit parameters = . The root mean square (r.m.s.) error for the fit is 0.34. The fit (plus and minimum the r.m.s. error) is plotted in the middle panel of Figure 3 for and . To match the LLS definition used in previous high-redshift FIRE papers (Faucher-Giguère et al. 2015; Faucher-Giguère et al. 2016), we also fitted the covering fractions for . The shape and redshift evolution of the best fit are consistent with the one above, so we held the , , , and parameters fixed and solved for the best-fit normalization for the different cut. The best-fit normalization for is , with an r.m.s. error of of 0.30. This is lower than for because higher-column systems are generally rarer.

The contribution of different halos to the cosmological incidence of absorbers can then be expressed in terms of functions describing the covering fraction as a function of redshift and halo mass, and the halo mass function:


For LLSs, is equivalent to . We show predicted by our simulations at and 1 in the bottom panel of Figure 3. The bottom panel shows that, at , the majority of randomly-selected LLSs arise from halos in the mass range , where covering fractions peak.

Figure 4: Covering fractions for systems with as a function of redshift for the six simulated halos shown in Figure 1. The burstiness of star formation and galactic winds in our simulations drives the time variability of the covering fractions.

3.2 Cosmological incidence

Figure 5: Solid curves: The cosmological incidence, , of LLSs predicted by our simulations for different HI column ranges ( in green, in orange and in blue). The error bars show observational measurements of from the LLS census of Ribaudo et al. (2011).

The total cosmological incidence of absorbers is obtained by integrating over halo mass:


We show the cosmological incidence of absorbers as a function of redshift for three different ranges in Figure 5. The first range, , is chosen to match the selection from Lehner et al. (2013). The two other ranges, and , match the observational measurements reported in Ribaudo et al. (2011). The comparison with the observational compilation of Ribaudo et al. (2011) shows that convolving the covering fractions predicted by our zoom-in simulations with the halo mass function yields cosmological incidences of LLSs that are consistent with observations at . Overall, the simulated increases slightly from to ; the significant redshift fluctuations apparent in Figure 5 reflect the strong time variability of covering fractions in individual simulations in our zoom-in sample (Fig. 4).

3.3 Column density distribution

Figure 6: H I column density distribution, , predicted using our simulations at , 0.5, and 1. The black dash-dotted lines show observational constraints from Ribaudo et al. (2011). In the LLS and SLLS regimes, our simulations agree well with the observational constraints. For lower column systems, our simulations lie below the observational constraints by up to dex but this is because many such low-column absorbers arise outside the virial radius of galaxy halos and therefore are not accurately captured by our analysis, which focuses on absorbers.

Our zoom-in simulations can also be used to predict the column density distribution function of CGM absorbers, , defined as the incidence of absorption systems per unit per unit absorption length, i.e.


The absorption length, , is defined such that the incidence of absorbers per absorption length, , is constant if the comoving surface area of the absorbers is constant across cosmic time (Bahcall & Peebles 1969):




Using the procedure of the previous section to compute for small column density intervals, we can compute :


where is the contribution to the incidence of absorbers for systems with column between and at redshift .

Figure 6 shows the column density distribution predicted by our simulations for cm. These predictions are compared to the observational constraints on the column density distribution from Ribaudo et al. (2011). Overall, the predicted column density distribution agrees well with observations for LLSs and SLLSs but the simulations under-predict the column density distribution of lower-column systems by dex at cm. Such a discrepancy at lower columns is expected since our analysis considers only absorbers that arise within the projected virial radius of main halos, i.e. it ignores absorbers that arise in the IGM outside dark matter halos.

4 Relationships between HI column, metallicity and radial kinematics

We now quantify the relationships between the HI column density of absorbers, their metallicity, and their radial kinematics relative to central galaxies.

For this analysis, we define “cosmological histograms” that enable us to study correlations between absorbers properties properly weighted by their cosmological incidence. For any prescribed interval, consider a random absorber (along the sight line to a random background quasar) with properties (e.g., HI column, metallicity, and radial velocity) described by the vector . Following the formalism of §3, the probability density function is


Here is the conditional probability density for given a redshift and halo mass . In our case, is the joint probability density function of HI column density, metallicity, and radial velocity. We use three orthogonal sky projections for each snapshot in evaluating , although the distribution varies relatively little for different viewing angles. Since is in general multi-variate, we can consider different 2D projections to quantify correlations between absorber attributes.

As a complement to the 2D cosmological histograms, we compute corresponding 2D histograms that quantify the mean halo mass contributing to different bins in parameter space:


For each bin of , this histogram gives the mean of halos contributing to the bin weighted by the cosmological incidence of absorbers.

For any given column density interval (such as LLSs), the weighting by in equation (10) is derived such that the resulting 2D histogram represents the bivariate distribution of physical properties for absorbers in that column density interval along random sight lines in the Universe. Below, we are also interested in the metallicity distribution as a function of HI column. For this, we replace the weighting factor by simply in equations (10) and (11). This is because in this case is not restricted to a prescribed interval, so we need not include the term in .

All histograms are evaluated using all simulated data. Quantities are interpolated linearly with respect to and to fill gaps between snapshots. Given the time variability of the CGM in individual halos (e.g., Fig. 4) and our modest-sized sample of simulated halos, averaging over this full redshift interval is necessary to mitigate stochastic fluctuations. Furthermore, LLS measurements at are presently limited to a few dozen systems so their distributions of physical parameters are obtained by averaging over a similar redshift interval (e.g., Lehner et al. 2013; Wotta et al. 2016).

4.1 Metallicity vs. HI column

Figure 7: Left: Cosmological 2D histogram for metallicity vs. HI column density from our simulations (representative of randomly selected LLSs). Simulated LLSs typically have a metallicity within 1 dex of the LLS mean metallicity . The metallicity distribution for lower column systems, cm, extends to lower metallicities. Right: 2D histogram showing the mean contributing to different bins on the left hand side. Overall, the contribution of more massive halos is weighted toward higher metallicities. The color bars are logarithmically scaled. We include systems with cm to provide information about their distribution integrated over galaxy halos, but caution that this distribution does not include a contribution from absorbers arising outside the virial radius of halos.

Figure 7 shows the cosmological histogram for metallicity vs. HI column and the corresponding mean halo mass histogram. As the panel on the left reveals, simulated LLSs typically have a metallicity within 1 dex of the median metallicity () of all LLSs. Lower column systems, cm on the other hand, have a metallicity distribution that extends well below , consistent with many of these absorbers representing a tail of IGM gas that has not yet been significantly enriched by galaxies. In the next section, we show that most absorbers have infall velocities and thus predominantly trace accretion of metal-poor intergalactic gas.

Overall, the mean mass histogram on the right hand side of Figure 7 reveals that more massive halos have a contribution weighted toward higher absorber metallicities. This overall trend is a CGM analog of the well-known mass-metallicity relation for galaxies (e.g., Tremonti et al. 2004; Zahid et al. 2011; Peeples et al. 2014) and which the FIRE simulations broadly reproduce (Ma et al. 2016). The tail of the metallicity distribution is weighted toward M halos.

4.2 Metallicity vs. radial kinematics

A primary motivation for our analysis is to understand how LLS metallicity can be used as a diagnostic of inflows and outflows. Figure 8 illustrates the relationship between absorber metallicity and its radial velocity relative to the central galaxy, , along with the corresponding mean mass histogram, for absorbers with (i.e., including partial LLSs). Positive radial velocities correspond to gas that is outflowing in an instantaneous sense, while negative radial velocities correspond to gas that is inflowing. In order to compare radial velocities for absorbers arising in a wide range of halo mass, velocities are normalized by halo circular velocity, , which corresponds to for  M.

Figure 8 shows that for cosmologically-selected LLSs in the redshift interval in our simulations:

  1. Very high velocity () outflows tend to have high-metallicity .

  2. Very low metallicity () LLSs tend to have radial velocities corresponding to inflows ().

  3. Most LLSs have modest absolute velocity . Overall, these LLSs have a broad metallicity distribution ranging from to , with no clear trend with radial kinematics.

Thus, our simulations suggest that very low metallicity LLSs with generally trace infalling IGM gas (consistent with observational interpretations of such low-metallicity absorbers, e.g. Ribaudo et al. 2011; Fumagalli et al. 2011) but that metallicity alone cannot robustly distinguish between inflows and outflows for the majority of LLSs with higher metallicity.

The concentration of LLSs with is easy to understand since is both the characteristic velocity of gas that is accelerated as it falls into halos from the IGM and the characteristic velocity of galactics winds in the FIRE simulations (Muratov et al. 2015). The difficulty of associating metal-enriched LLSs with inflows vs. outflows is due in large part to the importance of wind recycling in our simulations (for a previous analysis of how wind recycling shapes the galaxy stellar mass function, see Oppenheimer et al. 2010). When galactic winds recycle efficiently, metal-rich wind ejecta that initially have positive radial velocity later fall back onto the source galaxy with . Muratov et al. (2015) showed that the powerful galactic winds driven by high-redshift galaxies transform into galactic fountains by (or earlier) in massive galaxies in the FIRE simulations. Ma et al. (2016) furthermore demonstrated that the FIRE galaxies, except dwarfs with stellar mass M, retain most of the metals they have produced in their halos. Anglés-Alcázar et al. (2016) confirm the importance of wind recycling in the FIRE simulations using a particle tracking analysis that traces the full history of baryons.

As a check of our above result based on instantaneous radial kinematics that absorbers with can generally be associated with the accretion of fresh gas from the IGM, we have used the particle tracking pipeline from Anglés-Alcázar et al. (2016) and found the following: in the main halo of our m12i simulation at nearly all the gas associated with LLSs and classified as fresh accretion based on its full time history has metallicity , with more than half having metallicity . On the other hand, nearly all the wind LLSs in this halo have metallicity . We plan to more systematically extend this particle tracking analysis in future work. We note that previous authors found using particle tracking that wind recycling can constitute a significant (and possibly dominant) fraction of the instantaneously infalling CGM gas at late times in other simulations as well (Ford et al. 2014; Christensen et al. 2016). However, the magnitude of this effect depends on the properties of the (uncertain) wind models so previous results do not necessarily quantitatively apply to the FIRE simulations.

The ambiguities in using LLS metallicity to distinguish between inflows and outflows in the FIRE simulations at low redshift are in contrast to some earlier simulation analyses at high redshift (), in which maps of the CGM appear to show a much stronger correlation between instantaneous radial kinematics and metallicity than the halos shown in Figure 1 or than quantified in this section (e.g., Shen et al. 2012). Quantitatively, the analysis of a single main halo in the ErisMC simulation by Shen et al. (2012) revealed a factor of difference in mean metallicity between instantaneously inflowing and outflowing halo gas at . The difference between our results and the high redshift results of Shen et al. (2012) is likely at least partially attributable to the fact that both cosmological inflows and galactic winds are much stronger at high redshift. Furthermore, high redshift accretion of cool intergalactic gas dense enough to give rise to LLS absorption tends to be collimated in streams (e.g., Kereš et al. 2005; Dekel et al. 2009; Fumagalli et al. 2011; Faucher-Giguère & Kereš 2011) between which outflowing gas can propagate relatively unimpeded. By , Figure 1 shows that well-defined cool streams of infalling gas have largely disappeared in the simulations. We therefore caution that our present conclusions may not apply at .

Figure 8: Left: Cosmological 2D histogram for LLS metallicity vs. radial velocity relative to the central galaxy (representative of randomly selected LLSs). High velocity () outflows tend to have high metallicities and very low metallicity () LLSs tend to have infall radial velocities. However, most LLSs occupy an intermediate central region in metallicity-radial velocity space. For these more typical LLSs, there is no clear trend between metallicity and radial kinematics. Metal-enriched inflows arise in the FIRE simulations as a result of galactic winds that efficiently recycle at low redshift. Right: Mean mass histogram showing the mean contributing to different bins on the left hand side. The color bars are logarithmically scaled.
Figure 9: Top: Overall metallicity distribution for LLSs, averaging over halos of all masses (solid curve). The shaded region shows the 95% confidence interval estimated by sampling from the full simulation set with replacement. The multiple peaks apparent in the fiducial distribution are not significant. Bottom: Same LLS metallicity distribution as in the top panel (black) and divided into inflow (blue) and outflow (red) components based on instantaneous radial velocity. Outflowing LLSs have on average slightly higher metallicity than inflowing ones, but the two distributions overlap strongly, with no clean separation in metallicity.

4.3 The metallicity distribution of Lyman limit systems

To directly compare with the LLS metallicity distribution recently measured by Lehner et al. (2013), we project the 2D histogram in Figure 8 onto the metallicity axis. The result is shown in Figure 9. The simulated metallicity distribution has a mean (standard deviation) (0.4).

To test the statistical significance of features in the metallicity distribution, we performed a bootstrapping analysis. We recalculated the metallicity histogram 200 times, each time using a new sample of simulations drawn with replacement from our original set of simulations. The shaded region in the top panel of Figure 9 indicates the 95% confidence interval for each metallicity bin. From this analysis, we conclude that apparent multiple peaks in the metallicity distribution computed using our entire simulation sample are consistent with statistical noise inherent to the finite number of simulations in our sample. Thus, we do not find significant evidence for a multimodal metallicity distribution in our simulations. It is possible, though, that the simulated distribution has weak multimodal features within the noise of our analysis.

The lack of a significant metallicity bimodality in our simulations is in contrast with the possible metallicity bimodality observed by Lehner et al. (2013) and in the expanded sample of Wotta et al. (2016), which suggests a well-defined dip at . In our simulations, certain individual snapshots show rather well defined bimodal metallicity distributions (for example, the m11.9a snapshot shown in Fig. 1) with the low- and high-metallicity peaks roughly associated with inflows and galactic ejecta. However, inflows and outflows are highly time variable in our simulations and such features are transients. The main result of Figure 9 is that bimodal features in the metallicity distribution do not survive cosmological averaging over the broad range of halos in our simulations. This result is consistent with the broad range of halo mass contributing to the cosmological incidence of LLSs. Over the halo mass range M that contributes most to the cosmological LLS incidence, the interstellar gas phase and stellar metallicities span a total range of dex in our simulations (taking into account both the systematic increase with stellar mass and scatter around the mean mass-metallicity relation). These simulations are in good agreement with observational measurements of the mass-metallicity relations (Ma et al. 2016). Since the metallicity of galactic winds is similar to the metallicity of the interstellar medium (Muratov et al. 2016), the circum-galactic LLS metallicity distribution predicted by our simulations should therefore have a total width comparable to this. As a result, any clean bimodal “dip” in the metallicity distribution at a given halo mass (if present) may be expected to be substantially washed out in the cosmological average.

Accretion of fresh gas from the IGM could in principle give rise to a distinct very low metallicity peak in the LLS metallicity distribution, but such a distinct component is not present in our simulated cosmological distribution. This is likely because even gas that is first accreting onto a central galaxy from the IGM tends to be pre-enriched by ejecta from surrounding lower-mass galaxies. Smooth IGM accretion generally originates from similar large-scale structures as infalling satellite galaxies, and we have previously noted the effects of galactic winds from satellites on the observational properties of infalling gas at high redshift (Faucher-Giguère et al. 2015; Faucher-Giguère et al. 2016).

In the updated metallicity analysis of Wotta et al. (2016), the statistical evidence for a bimodality is strongest for the subsample of partial LLSs, with . However, 15 out of 44 absorbers in that sample only have lower or upper limits on their metallicity. Allowing the lower limit absorbers to have arbitrarily high metallicity and the upper limit absorbers to have arbitrarily low metallicity could significantly flatten the distribution and potentially reduce the statistical significance of the dip at . It will be important to measure more precisely the metallicities of absorbers with only limits currently available (or to rigorously account for the limits in the statistical analysis) to more conclusively determine statistical significance of the observed dip. At , the observational analyses of Fumagalli et al. (2016) and Lehner et al. (2016) reveal a broad unimodal LLS metallicity distribution.

Another discrepancy between our simulations and the metallicity distributions measured by Lehner et al (2013) and Wotta et al. (2016), which cannot be explained by the presence of limits in the observational sample, is that our simulations contain significantly fewer very low metallicity LLSs than the observations. For example, our simulated metallicity distribution contains very few () LLSs with but about half of Wotta et al.’s low-metallicity branch have such low metallicity. Wotta et al. explore systematic uncertainties in metallicity estimates from the assumed UV background model and find that using a UV background model with a harder spectrum due to a larger contribution from quasars would increase their inferred metallicities by up to dex. This effect could explain part of the discrepancy at low metallicity if pushed to its extreme, but would introduce some tension at the high-metallicity end.

Differences in metal mixing may also contribute to the low-metallicity discrepancy. In our SPH simulations, metals are advected by particles but do not diffuse between particles, so that metal mixing occurring below the resolution limit is not captured. Underestimating subgrid metal mixing in the simulations should (if anything) overestimate the number of very low metallicity LLSs, rather than underestimate it. We have verified that the “missing” low-metallicity systems cannot be explained by the absence of subgrid metal mixing in our simulations by analyzing a re-simulation of our m12i halo that includes a subgrid turbulent metal diffusion model similar to Shen et al. (2010). This re-simulation did not produce more absorbers. On the other hand, if metals in the real CGM are poorly mixed, then cosmological simulations (both SPH and grid-based) could overestimate their mixing by forcing new metals to be distributed to at least one resolution element. In our simulations, metals produced by stellar evolution are distributed to the gas particles in the SPH smoothing kernel, corresponding to a gas mass M for the m11.4a and m11.9a simulations whose parameters are listed in Table 1. In an observational analysis of CIV absorbers, Schaye et al. (2007) concluded that those intergalactic absorbers likely arise from a large population of compact metal clumps with masses M, which we are unable to resolve in our simulations. If the metals in LLSs are similarly poorly mixed, then our simulations could overestimate the metallicity of many LLS sight lines within which metals in reality would not have had enough time to mix. This effect could potentially “hide” a population of low-metallicity LLSs.

To further investigate the relationship between metallicity and inflows vs. outflows, we divide the total metallicity distribution into instantaneously inflowing () and outflowing () components in the bottom panel of Figure 9. The instantaneously inflowing LLSs have a mean (standard deviation) (0.4), while the instantaneously outflowing LLSs have a mean (standard deviation) (0.4). Consistent with our results from the previous section, infalling LLSs are on average of slightly lower metallicity than outflowing ones, but the distributions are not cleanly separated in metallicity. This suggests that metallicity alone cannot be used as a diagnostic to distinguish instantaneous inflows from instantaneous outflows for LLSs with [X/H] , i.e. the majority of LLSs. However, as discussed in §4.2, the lowest metallicity systems with [X/H] do tend to be inflows.

5 Conclusions

We have used cosmological hydrodynamic zoom-in simulations from the FIRE project to investigate the physical properties of circum-galactic absorption at , with particular emphasis on Lyman limit systems and their use as diagnostics of cosmological inflows and galactic winds. The FIRE simulations self-consistently generate galactic winds from energy and momentum injection on the scale of star-forming regions by stellar feedback. We analyzed 14 simulations covering the halo mass range  M at , which we convolved with the dark matter halo mass function to obtain cosmological statistics representative of absorbers randomly selected along unbiased sight lines to background quasars. Our main conclusions are as follows:

  1. When convolved with the dark matter halo mass function, the FIRE simulations are consistent with the observed cosmological incidence and HI column density distribution of Lyman limit systems at .

  2. The majority of HI-selected LLSs are associated with relatively massive halos in the mass range M. Dwarf halos with M have extremely small () LLS covering fractions within a projected virial radius.

  3. The LLS covering fractions of individual halos are strongly time variable. The strong variability results from a combination of time-variable inflows (including accreting satellites) and bursty outflows.

  4. Simulated LLSs typically have a metallicity within 1 dex of the mean metallicity . The metallicity distribution of lower column systems, cm, extends well below . This is consistent with many of these low-column absorbers representing a tail of IGM gas that has not yet been significantly enriched by galaxies.

  5. LLSs with large radial outflow velocity () tend to have high metallicities , while very low metallicity () LLSs tend to have infall radial velocities.

  6. Most LLSs have moderate radial velocities (). For these more common LLSs, there is no strong trend between metallicity and instantaneous radial kinematics.

  7. When separating LLSs in groups with (inflows) and (outflows), the inflowing LLSs have a slightly lower mean metallicity () than the outflowing LLSs (). However, both distributions have a standard deviation of [X/H] , causing them to overlap strongly.

  8. Overall, we find no significant evidence for a bi-modality of LLS metallicity in our simulations. This result is in tension with observations that suggest two metallicity branches clearly separated at .

  9. The simulated metallicity distribution also lacks a population of low-metallicity LLSs () that is prominent in observations, with only of all simulated LLSs having that low metallicity. The existence of such low-metallicity LLSs may indicate that metals are poorly mixed in the observed CGM on scales below the resolution of our simulations.

Overall, our simulations indicate that very low metallicity LLSs are predominantly associated with cosmological inflows at , but that higher metallicity LLSs trace gas with roughly equal probability of having instantaneous inflow or outflow kinematics. This result is largely a consequence of the prevalence of gas recycling in the FIRE simulations, which causes a large fraction of metal-rich galactic wind ejecta to later fall back onto galaxies (Anglés-Alcázar et al. 2016). Thus, metallicity is a powerful diagnostic of pristine intergalactic inflows but in general cannot robustly distinguish between recent outflows and inflows of recycled gas. It will be interesting to sharpen this result by using particle tracking to more accurately identify the physical nature and history of gas elements in our simulations. Such an analysis may reveal that metallicity is a more powerful diagnostic of gas that has been previously processed by galaxies vs. not than of instantaneous inflows vs. instantaneous outflows. It will also be important to firm up the statistical significance of our analysis using a larger sample of simulated halos.

We also plan to extend our analysis of LLSs to the high-redshift Universe, where observational measurements of LLSs and their metallicities are also available (Lehner et al. 2014; Fumagalli et al. 2016; Lehner et al. 2016). At high redshift, many LLSs arise outside the virial radius of galaxy halos (e.g. Cooper et al. 2015) and the present zoom-in approach will not be adequate to study cosmological statistics. We plan to pursue full-volume cosmological simulations with the FIRE resolution and physics to address such IGM questions. Another promising observational diagnostic of inflows and outflows that will be worthwhile to investigate using simulations is the distribution in azimuthal angle relative to the galactic disk plane of strong MgII absorbers (e.g., Bordoloi et al. 2011; Bouché et al. 2012; Kacprzak et al. 2012; Lan et al. 2014).


The authors are grateful to Nicolas Lehner, Chris Wotta, Chris Howk, J. X. Prochaska, Cameron Liang, Joop Schaye, and Andrey Kravtsov for discussions regarding the observed Lyman limit system metallicity bimodality. We are also grateful to Benedikt Diemer for providing a set of Python cosmology modules used in this work. ZH, CAFG, and DAA were supported by NSF through grants AST-1412836, AST-1517491 and DGE-0948017, by NASA through grant NNX15AB22G, and by STScI through grants HST-AR-14293.001-A and HST-GO-14268.022-A. DK and TKC were supported in part by NSF grant AST-1412153 and Cottrell Scholar Award from the Research Corporation for Science Advancement. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NASA ATP grant NNX14AH35G, and NSF grants AST-1411920 and AST-1455342. EQ was supported by NASA ATP grant 12-ATP-120183, a Simons Investigator award from the Simons Foundation, and the David and Lucile Packard Foundation. The simulations analyzed in this paper were run on XSEDE computational resources (allocations TG-AST120025, TG-AST130039, and TG-AST140023), on the NASA Pleiades cluster (allocation SMD-14-5189), on the Northwestern Quest computer cluster, and on the Caltech Zwicky computer cluster.

Appendix A Convergence Properties

Figure 10: Grid resolution convergence for the LLS () covering fractions within a projected virial radius in the m12i simulation. We tested grids equivalent to our fiducial resolution and one higher resolution level ( and grid points on grids with side lengths of ), but spanning a single quadrant of the virial radius. The covering fractions are well converged for the fiducial grid resolution used in our analysis.
Figure 11: Same as Figure 10 but for the LLS metallicity distribution in the m12i simulation.

To test our results for convergence with respect to grid resolution, we regridded a single quadrant of m12i at both fiducial resolution and at one higher resolution level over the full redshift range . If we define the origin as the center of the m12i galaxy, our new cubic meshes have and grid points over a volume with and . This produces grids with and equivalent resolution over the full halo. The results of our grid resolution convergence tests are shown in Figure 10 for LLS covering fractions and in Figure 11 for the metallicity distribution. Both the covering fractions and the metallicity distribution are well converged with respect to grid resolution. This is expected since the SPH smoothing lengths of LLS gas are comparable to the grid cell size at the fiducial grid resolution used for our analysis. We note that the LLS covering fractions plotted in Figure 10 for the m12i simulation are somewhat more time variable than in Figure 4 in the main text because here we average over only a single quadrant of the simulation, whereas Figure 4 averages over the entire projected virial area.

Figure 12: Overall cosmological LLS metallicity distribution recomputed by artificially increasing (red) and decreasing (blue) by a factor of 5 the weights of the most massive halos (m13 and MFz0_A2), which are simulated with coarser mass resolution. The good agreement with the fiducial distribution (black) indicates that these contribute relatively little to the cosmological LLS metallicity distribution. The lower resolution of these simulations is thus unlikely to significantly affect our results.

The most massive simulated halos in our sample, m13 and MFz0_A2, were run at one resolution level lower (gas particles masses M) than the rest of our simulation sample. To test whether our results are sensitive to possible resolution effects in these simulations, we artificially varied their LLS covering fractions by a factor of 5 up and down and re-computed the cosmological metallicity distribution. Figure 12 shows the results. Even if these massive halos were subject to significant resolution limitations, our main results would be affected only very slightly. This is because halos in this mass range are rare and contribute negligibly to global LLS statistics.


  1. Project website: http://fire.northwestern.edu


  1. Agertz O., Moore B., Stadel J., Potter D., Miniati F., Read J., Mayer L., Gawryszczak A., Kravtsov A., Nordlund Å., Pearce F., Quilis V., Rudd D., Springel V., Stone J., Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MNRAS, 380, 963
  2. Aguirre A., Hernquist L., Schaye J., Weinberg D. H., Katz N., Gardner J., 2001, ApJ, 560, 599
  3. Altay G., Theuns T., Schaye J., Booth C. M., Vecchia C. D., 2013, MNRAS, 436, 2689
  4. Altay G., Theuns T., Schaye J., Crighton N. H. M., Dalla Vecchia C., 2011, ApJL, 737, L37
  5. Anglés-Alcázar D., Davé R., Özel F., Oppenheimer B. D., 2014, ApJ, 782, 84
  6. Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2016, ArXiv e-prints
  7. Asplund M., Grevesse N., Sauval A. J., Scott P., 2010, Astrophys. Space Sci., 328, 179
  8. Bahcall J. N., Peebles P. J. E., 1969, ApJ, 156, L7
  9. Barnes J. E., 2012, MNRAS, 425, 1104
  10. Bauermeister A., Blitz L., Ma C.-P., 2010, ApJ, 717, 323
  11. Bird S., Vogelsberger M., Sijacki D., Zaldarriaga M., Springel V., Hernquist L., 2013, MNRAS, 429, 3341
  12. Booth C. M., Schaye J., Delgado J. D., Dalla Vecchia C., 2012, MNRAS, 420, 1053
  13. Bordoloi et al. 2011, ApJ, 743, 10
  14. Bouché N., Hohensee W., Vargas R., Kacprzak G. G., Martin C. L., Cooke J., Churchill C. W., 2012, MNRAS, 426, 801
  15. Brooks A. M., Governato F., Quinn T., Brook C. B., Wadsley J., 2009, ApJ, 694, 396
  16. Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  17. Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov a. L., Faucher-Giguère C. a., Quataert E., 2015, MNRAS, 454, 2981
  18. Christensen C. R., Davé R., Governato F., Pontzen A., Brooks A., Munshi F., Quinn T., Wadsley J., 2016, ApJ, 824, 57
  19. Cicone et al. 2014, A&A, 562, A21
  20. Cooper T. J., Simcoe R. A., Cooksey K. L., O’Meara J. M., Torrey P., 2015, ApJ, 812, 58
  21. Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
  22. Dekel A., Birnboim Y., Engel G., Freundlich J., Goerdt T., Mumcuoglu M., Neistein E., Pichon C., Teyssier R., Zinger E., 2009, Nature, 457, 451
  23. Faucher-Giguère C.-A., Feldmann R., Quataert E., Kereš D., Hopkins P. F., Murray N., 2016, MNRAS, 461, L32
  24. Faucher-Giguère C.-A., Hopkins P. F., Kere D., Muratov A. L., Quataert E., Murray N., 2015, MNRAS, 449, 987
  25. Faucher-Giguère C. A., Kereš D., 2011, MNRAS, 412, 118
  26. Faucher-Giguère C.-A., Kereš D., Ma C.-P., 2011, MNRAS, 417, 2982
  27. Faucher-Giguère C. a., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
  28. Feldmann R., Quataert E., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2016, ArXiv e-prints
  29. Ford A. B., Davé R., Oppenheimer B. D., Katz N., Kollmeier J. A., Thompson R., Weinberg D. H., 2014, MNRAS, 444, 1260
  30. Fox A. J., Lehner N., Tumlinson J., Howk J. C., Tripp T. M., Prochaska J. X., O’Meara J. M., Werk J. K., Bordoloi R., Katz N., Oppenheimer B. D., Davé R., 2013, ApJ, 778, 187
  31. Fumagalli M., Hennawi J. F., Prochaska J. X., Kasen D., Dekel A., Ceverino D., Primack J., 2014, ApJ, 780, 74
  32. Fumagalli M., O’Meara J. M., Prochaska J. X., 2011, Science, 334, 1245
  33. Fumagalli M., O’Meara J. M., Prochaska J. X., 2016, MNRAS, 455, 4100
  34. Fumagalli M., O’Meara J. M., Prochaska J. X., Worseck G., 2013, ApJ, 775, 78
  35. Fumagalli M., Prochaska J. X., Kasen D., Dekel A., Ceverino D., Primack J. R., 2011, MNRAS, 418, 1796
  36. Goerdt T., Dekel A., Sternberg A., Gnat O., Ceverino D., 2012, MNRAS, 424, 2292
  37. Greene J. E., Zakamska N. L., Ho L. C., Barth A. J., 2011, ApJ, 732, 9
  38. Gurvich A., Burkhart B., Bird S., 2016, arXiv:1608.03293
  39. Heckman T. M., Lehnert M. D., Strickland D. K., Armus L., 2000, ApJS, 129, 493
  40. Hernquist L., Katz N., Weinberg D. H., Miralda-Escudé J., 1996, ApJL, 457, L51
  41. Hopkins P. F., 2013, MNRAS, 428, 2840
  42. Hopkins P. F., 2015, MNRAS, 450, 53
  43. Hopkins P. F., Keres D., Onorbe J., Faucher-Giguere C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581
  44. Jones T., Stark D. P., Ellis R. S., 2012, ApJ, 751, 51
  45. Kacprzak G. G., Churchill C. W., Nielsen N. M., 2012, ApJL, 760, L7
  46. Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
  47. Katz N., White S. D. M., 1993, ApJ, 412, 455
  48. Kereš D., Katz N., Fardal M., Davé R., Weinberg D. H., 2009, MNRAS, 395, 160
  49. Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  50. Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
  51. Kohler K., Gnedin N. Y., 2007, ApJ, 655, 685
  52. Kollmeier J. A., Miralda-Escudé J., Cen R., Ostriker J. P., 2006, ApJ, 638, 52
  53. Lan T.-W., Ménard B., Zhu G., 2014, ApJ, 795, 31
  54. Lehner N., Howk J. C., Tripp T. M., Tumlinson J., Prochaska J. X., O’Meara J. M., Thom C., Werk J. K., Fox a. J., Ribaudo J., 2013, ApJ, 770, 138
  55. Lehner N., O’Meara J., Fox A., 2014, ApJ, 788, 119
  56. Lehner N., O’Meara J. M., Howk J. C., Prochaska J. X., Fumagalli M., 2016, arXiv:1608.02588
  57. Leitherer C., Schaerer D., Goldader J. D., Delgado R. M. G., Robert C., Kune D. F., de Mello D. F., Devost D., Heckman T. M., 1999, ApJS, 123, 3
  58. Liang C. J., Kravtsov A. V., Agertz O., 2016, MNRAS, 458, 1164
  59. Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, MNRAS, 456, 2140
  60. McQuinn M., Oh S. P., Faucher-Giguère C.-A., 2011, ApJ, 743, 82
  61. Martin C. L., 2005, ApJ, 621, 227
  62. Muratov A. L., Keres D., Faucher-Giguere C.-A., Hopkins P. F., Ma X., Angles-Alcazar D., Chan T. K., Torrey P., Hafen Z. H., Quataert E., Murray N., 2016, arXiv:1606.09252
  63. Muratov A. L., Keres D., Faucher-Giguere C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691
  64. Murray N., 2011, ApJ, 729, 133
  65. Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265
  66. Oppenheimer B. D., Davé R., Kereš D., Fardal M., Katz N., Kollmeier J. A., Weinberg D. H., 2010, MNRAS, 406, 2325
  67. Peeples M. S., Werk J. K., Tumlinson J., Oppenheimer B. D., Prochaska J. X., Katz N., Weinberg D. H., 2014, ApJ, 786, 54
  68. Planck Collaboration Ade P. A. R., Aghanim N., Arnaud M., Ashdown M., Aumont J., Baccigalupi C., Banday A. J., Barreiro R. B., Bartlett J. G., et al. 2015, arXiv:1502.01589
  69. Porter D. H., 1985, PhD thesis, California Univ., Berkeley.
  70. Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  71. Prochaska J. X., Wolfe A. M., 2009, ApJ, 696, 1543
  72. Prochaska et al. 2015, ApJS, 221, 2
  73. Quiret S., Péroux C., Zafar T., Kulkarni V. P., Jenkins E. B., Milliard B., Rahmani H., Popping A., Rao S. M., Turnshek D. A., Monier E. M., 2016, MNRAS, 458, 4074
  74. Rahmati A., Pawlik A. H., Raičevic̀ M., Schaye J., 2013, MNRAS, 430, 2427
  75. Rahmati A., Schaye J., Bower R. G., Crain R. A., Furlong M., Schaller M., Theuns T., 2015, MNRAS, 452, 2034
  76. Ribaudo J., Lehner N., Howk J. C., 2011, ApJ, 736, 42
  77. Ribaudo J., Lehner N., Howk J. C., Werk J. K., Tripp T. M., Prochaska J. X., Meiring J. D., Tumlinson J., 2011, ApJ, 743, 207
  78. Rubin K. H. R., Prochaska J. X., Koo D. C., Phillips A. C., Martin C. L., Winstrom L. O., 2014, ApJ, 794, 156
  79. Schaye J., Carswell R. F., Kim T.-S., 2007, MNRAS, 379, 1169
  80. Shapley A. E., Steidel C. C., Pettini M., Adelberger K. L., 2003, ApJ, 588, 65
  81. Shen S., Madau P., Aguirre A., Guedes J., Mayer L., Wadsley J., 2012, ApJ, 760, 50
  82. Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581
  83. Somerville R. S., Davé R., 2015, Annu. Rev. Astron. Astrophys., 53, 51
  84. Sparre M., Hayward C. C., Feldmann R., Faucher-Giguère C.-A., Muratov A. L., Kereš D., Hopkins P. F., 2015, arXiv:1510.03869
  85. Springel V., 2005, MNRAS, 364, 1105
  86. Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289
  87. Tremonti C. A., Heckman T. M., Kauffmann G., Brinchmann J., Charlot S., White S. D. M., Seibert M., Peng E. W., Schlegel D. J., Uomoto A., Fukugita M., Brinkmann J., 2004, ApJ, 613, 898
  88. Tripp T. M., Meiring J. D., Prochaska J. X., Willmer C. N. A., Howk J. C., Werk J. K., Jenkins E. B., Bowen D. V., Lehner N., Sembach K. R., Thom C., Tumlinson J., 2011, Science, 334, 952
  89. van de Voort F., Schaye J., Altay G., Theuns T., 2012, MNRAS, 421, 2809
  90. van de Voort F., Schaye J., Booth C. M., Haas M. R., Dalla Vecchia C., 2011, MNRAS, 414, 2458
  91. Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G., 2013, MNRAS, 433, 1230
  92. Weiner B. J., Coil A. L., Prochaska J. X., Newman J. A., Cooper M. C., Bundy K., Conselice C. J., Dutton A. A., Faber S. M., Koo D. C., Lotz J. M., Rieke G. H., Rubin K. H. R., 2009, ApJ, 692, 187
  93. Wotta C. B., Lehner N., Howk J. C., O’Meara J. M., Prochaska J. X., 2016, ApJ, 831, 95
  94. Zahid H. J., Kewley L. J., Bresolin F., 2011, ApJ, 730, 137
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description