Damped Lyman Alpha Systems

Damped Lyman Alpha Systems in Galaxy Formation Simulations

Andrew Pontzen, Fabio Governato, Max Pettini, C.M. Booth, Greg Stinson, James Wadsley, Alyson Brooks, Thomas Quinn, Martin Haehnelt
Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK
Astronomy Department, Box 351580, University of Washington, Seattle, WA 98195, USA
Department of Physics, Institute for Computational Cosmology, University of Durham, South Road, Durham, UK
Sterrewacht Leiden, University of Leiden, P.O. Box 9513, 2300 RA Leiden, The Netherlands
Department of Physics and Astronomy, McMaster University, Hamilton, ON L8S 4M1, Canada
Email: apontzen@ast.cam.ac.uk
Accepted 2008 July 30. Received 2008 July 30; in original form 2008 April 19

We investigate the population of damped Lyman alpha systems (DLAs) in a recent series of high resolution galaxy formation simulations. The simulations are of interest because they form at some of the most realistic disk galaxies to date. No free parameters are available in our study: the simulation parameters have been fixed by physical and observational constraints, and thus our work provides a genuine consistency test. The precise role of DLAs in galaxy formation remains in debate, but they provide a number of strong constraints on the nature of our simulated bound systems at because of their coupled information on neutral  H i densities, kinematics, metallicity and estimates of star formation activity.

Our results, without any parameter-tuning, closely match the observed incidence rate and column density distributions of DLAs. Our simulations are the first to reproduce the distribution of metallicities (with a median of ) without invoking observationally unsupported mechanisms such as significant dust biasing. This is especially encouraging given that these simulations have previously been shown to have a realistic stellar mass-metallicity relation. Additionally, we see a strong positive correlation between sightline metallicity and low-ion velocity width, the normalization and slope of which comes close to matching recent observational results. However, we somewhat underestimate the number of observed high velocity width systems; the severity of this disagreement is comparable to other recent DLA-focused studies.

DLAs in our simulations are predominantly associated with dark matter haloes with virial masses in the range . We are able to probe DLAs at high resolution, irrespective of their masses, by using a range of simulations of differing volumes. The fully constrained feedback prescription in use causes the majority of DLA haloes to form stars at a very low rate, accounting for the low metallicities. It is also responsible for the mass-metallicity relation which appears essential for reproducing the velocity-metallicity correlation. By the majority of the neutral gas forming the DLAs has been converted into stars, in agreement with rough physical expectations.

quasars: absorption lines – galaxies: formation – methods: numerical
pubyear: 2008

1 Introduction

One of the most difficult and most important questions to ask of any cosmological simulation is whether, even if it matches some known properties of the observed Universe, the route via which it obtained those results is physically meaningful. It is tempting to argue that, with the degree of parameter-tuning available to the modern simulator (stemming from our inability to maintain a sufficient dynamic range, uncertainty in gas physics and in particular star formation and feedback prescriptions), attempts to match a small number of observed properties can succeed without representing a physical route to that success. A sensible test of any suite of simulations, therefore, is to scrutinize its predictions for observed relations which were not considered in the process of planning those simulations. Success or failure of the simulation to match such relations cannot be equated to success or failure of the simulation and its predictions as a whole, but it can lend weight in either direction.

In this paper we will apply such an approach to the series of galaxy formation simulations most recently described in Governato et al. (2007), Brooks et al. (2007) and Governato et al. (2008) (henceforth G07, B07 & G08 respectively). The outputs of these simulations contain more realistic disc galaxies than have previously been achieved. In particular, the simulated galaxies form rotationally supported disks falling on the Tully-Fisher and baryonic Tully-Fisher relations, have a distribution of satellites compatible with local observations (G07) and have reasonable stellar mass-metallicity relations (B07).

It should be noted, however, that in common with other galaxy formation simulations, the mass in the bulge component of the G07 galaxies is overestimated (see also Eke et al., 2001). The problem is essentially one of angular momentum loss; it seems that to prevent this, both high numerical resolution (Kaufmann et al., 2007) and better models of feedback from supernova explosions are important elements (G08). The exaggerated bulges cause rotation curves to decline unrealistically over a few disk scale lengths to of their peak value. However, these problems appear to be shrinking in magnitude as resolution increases – and, regardless, the simulations under consideration form galaxies at which are as realistic as current computational and modeling power will allow, so that a detailed study of their properties during formation is justified.

In this paper, we will investigate at the predominantly neutral gas which gives rise to damped Lyman alpha systems (DLAs). These are systems with column densities of  H i  in excess of , seen in absorption against more distant luminous sources (generally quasars); for a recent review see Wolfe et al. (2005). The particular limit is historical, corresponding to the column densities expected if the Milky Way were to be viewed face on (Wolfe et al., 1986), but simple physical arguments suggest it makes a convenient distinction between trace  H i in the ionised intergalactic medium (IGM, below the limit) and clouds which are predominantly composed of  H i (above the limit, see Wolfe et al., 2005). The latter clouds must absorb, in an outer layer, the majority of incident photons capable of ionising hydrogen () and are therefore termed “self-shielding”.

The existence of an intergalactic ultraviolet (UV) field arising from the cumulative effect of external galaxies and quasars (e.g. Haardt & Madau, 1996) plays many roles in understanding the state of these clouds – not only does it affect the ionisation levels in the optically thin transition regions, but it also contributes substantially to the heating budget via photo-ejection of electrons from atoms. Consequently, gas cannot cool to form neutral clouds without first collapsing in the presence of a gravitational mass with virial velocity of tens of (Rees, 1986; Quinn et al., 1996), suggesting such clouds are associated with dark matter haloes and hence protogalaxies. This rough physical argument is verified by previous simulations (many of which are listed in Table 1), although in the simulations of Razoumov et al. (2006) and Razoumov et al. (2008) DLAs often bridge multiple haloes, extending into the intervening IGM. (We discuss this issue further in Section 5.5.)

Reference(s) Type SF Ionization/RT Max Vol Gas Res
Katz et al. (1996b) SPH None Plane Correction
Gardner et al. (1997a) SPH None Plane Correction
Gardner et al. (1997b)
Haehnelt et al. (1998) SPH None Den. Cut N/A
Gardner et al. (2001) SPH Yes, weak FB Plane Correction
Cen et al. (2003) Eulerian Yes, with FB Hybrid
Nagamine et al. (2004a) SPH Multiphase/GW Eq. Thin/MP
Nagamine et al. (2004b)
Razoumov et al. (2006) Adpt Eulerian None Non-Eq. Live RT/post-processor
Nagamine et al. (2007) SPH Multiphase/GW Eq. Thin/MP
Razoumov et al. (2008) Adpt Eulerian Basic Non-Eq. Thin/post-processor
This work SPH Yes, with FB Eq. Thin/RT post-processor
Table 1: Selected previous simulations of DLAs. Notes on methods: The largest volume simulated for the study, in comoving units. The best gas resolution achieved in the study, which may not have been achieved in the largest volume. For SPH (Lagrangian) simulations we give the smallest particle mass; for Eulerian simulations we give the finest grid resolution (in physical units at ). UV background in optically thin limit; sightlines post-processed using plane parallel radiative transfer and ionisation equilibrium. UVB optically thin, but in post-processing all gas particles assumed fully neutral for number densities . This study used a re-sampling technique to study the high resolution dynamics of a limited number of haloes and did not collect cosmological statistics.FB = Feedback, i.e. the deposition of energy into the ISM due to supernova explosions. By “Weak” FB is meant thermal injection only, which is generally recognized to be insufficient (see Section 5.1) The optical depth for each cell was calculated and used to approximate a local attenuation to the UVB. The multiphase method keeps track of the fraction of a gas particle in a cold cloud phase in pressure equilibrium with the warm medium; the cold clouds are assumed to be fully self shielded while the ambient medium is regarded as optically thin. Phenomenological galactic winds (GW) are added, causing bulk outflows from all haloes. Adaptive Eulerian: the grids which keep track of gas properties are automatically refined as regions collapse on small scales. A comprehensive treatment of radiative transfer was used by this paper, using eight angular elements in the live simulation and 192 in a post-processing stage. See Sections 2 and 3.1.

Irrespective of their physical nature, it is simple to show directly that DLAs contain the majority of  H i over all redshifts (e.g. Tytler, 1987), which suggests they have an important role to play in the global star formation history. Curiously, however, a wide range of diagnostics provide compelling evidence that the star formation rates in typical DLAs are small (). These include the low characteristic metallicities and dust depletions (for a review see Pettini, 2006), the extreme rarity of detectable optical counterparts and the faintness of Ly emission (a large number of studies are summarised in Wolfe et al., 2005). In DLAs with detectable molecular hydrogen (H), the local UV can be estimated from pumping into high energy rotational levels (e.g. Srianand et al., 2005), results suggesting star formation rates comparable to the present day Milky Way (). However, since few DLAs are associated with detectable H absorption, these measurements are not necessarily representative of the wider population. The estimated cooling rates from the recently developed C ii] technique (Wolfe et al., 2003) lead to star formation rate density estimates of , although the exact interpretation of these results is complicated by various assumptions in the method and the unknown area of a typical DLA system.

A further constraint on the hosts of DLAs is given by kinematic information encoded in unsaturated metal absorption lines. Both the neutral gas (traced by low-ion transitions such as Si ii) and any surrounding ionised gas (traced by high-ion transitions such as C iv) may be probed. The exact relationship between the high- and low-ion regions is not entirely certain (e.g. Fox et al., 2007, and references therein). Emphasis in our work, and most other studies, is placed on the low-ion profiles since these presumably reflect the kinematics of the gas giving rise to the DLA itself.

The earliest systematical survey of DLA low-ion velocity profiles was conducted by Prochaska & Wolfe (1997), who suggested that the observed kinematics could arise from a population of thick cold rotating disks with a distribution of rotational velocities similar to that observed in local disk galaxies. This is at odds with physical intuition and the prevailing Cold Dark Matter (CDM) hierarchical cosmogony in the sense that it requires the halo mass function, and hence power spectrum of fluctuations, to remain almost unchanged over Gyr between and . A view more compatible with the standard model seems desirable. Haehnelt et al. (1998) were quick to apply numerical simulations of structure formation to show that a fiducial population of CDM haloes were not incapable of producing velocity profiles similar to those observed – however, some details have proved more problematic as we describe below.

One way of quantifying the simplest property of the kinematics is to assign to each DLA a “velocity width” , roughly measuring the Doppler broadening of any unsaturated low-ion transition (see Section 3.2 for a more precise definition). The velocity width may be presumed to give an indication of the underlying virial velocity of the system responsible, although there is no guarantee that a particular sightline will sample the entire range of the velocity dispersion within the system; the simulations of Haehnelt et al. (1998) suggested that with a large scatter.

Simulations such as those by Katz et al. (1996b); Gardner et al. (1997a, 2001); Nagamine et al. (2004a) focused instead on the cross-sectional size of haloes as DLAs. Taken with a halo mass function (which throughout this work, we will produce for a CDM “concordance” scenario) such trends can be used to predict the overall rate of occurrence of DLA systems – a prediction which, in most simulations, agrees roughly with observations. Including the results of Haehnelt et al. (1998), one can further attempt to predict the relative proportions of systems of differing velocity width. In general, simulations underestimate the incidence of high velocity widths; this can be traced to the majority of the cross-section being assembled from low mass () haloes. Similar difficulties are also encountered in varying degrees by semi-analytic models of DLAs. These are, of course, not suited to producing the details of the line widths but they can help to pin down which physical considerations affect them (see e.g. Kauffmann, 1996; Maller et al., 2001; Johansson & Efstathiou, 2006).

Overall, the extent of the implications of the velocity mismatch is somewhat unclear. Some authors have argued that a fundamental difficulty with the CDM scenario has been uncovered (Prochaska & Wolfe, 2001). However, numerical modeling of the DLA population is intrinsically troublesome. The dynamic range of processes involved in making the final population is tremendous, with the relevant scales ranging from cosmological to stellar. Therefore we suggest that failure to match – or, indeed, success in matching – specific observations should be interpreted conservatively. We now describe some specific uncertainties in understanding the simulated DLA population, and outline how these are dealt with in our work. For comparison, Table 1 gives details of previous simulations and their approach to these problems.

  1. Star Formation. Although Gardner et al. (1997a) argue that star formation at has little impact on the column densities of individual systems, it is not a priori clear how the kinematics and cross-section are affected by the supernovae feedback. Investigations in this direction have been made by Nagamine et al. (2004a, b), using a phenomenological galactic wind prescription. In the present paper, our star formation is fully prescribed by physical models and observations, leaving no parametric freedom; see Section 2.

  2. Self-Shielding. DLAs are known to contain gas which is largely self-shielding (see above); thus the fiducial “uniform UV background” which is used in many cosmological simulations proves inadequate. We use a simple radiative transfer post-processor to correct for this (Section 3.1). We assess both the algorithm’s reliability and the severity of neglecting radiative transfer in the live simulation in Section 5.4.

  3. Cosmological Sampling. The total DLA cross-section is thought to be evenly spread through many orders of magnitude of parent halo mass (Gardner et al., 1997a, b); observations of low redshift DLAs appear to confirm such a view (Zwaan et al., 2005). Presumably one must allow a sufficient volume in a simulation to study a number of protocluster regions as well as maintain sufficient resolution to resolve low mass dwarf galaxies. Techniques involving extrapolation into unresolved regimes (such as fitting functional forms to the DLA cross-section as a function of halo mass: Gardner et al., 1997a, b) provide useful signposts but naturally must be regarded with caution. Our simulations, and those by Razoumov et al. (2006) and Nagamine et al. (2004a) resolve all haloes of relevance only by separately simulating a number of boxes of varying size; thus a way of combining results from the various boxes must be found. In our case, this involves using the halo mass function to re-weight sightlines in calculating cosmological properties (Section 3.3).

The remainder of this paper is structured as follows. In Sections 2 and 3 we describe respectively the series of simulations in use for our study and our methods for extracting DLA sightlines and producing quantities representative of a cosmological ensemble. We give results for these quantities and their underlying relations in Section 4. Section 5 describes a number of consistency checks and runs with altered parameters which shed further light on the origin of some of our relations. Finally, we summarise and discuss future work in Section 6. Two appendices contain technical detail for completeness.

We adopt the following conventions. Except where specified, all quoted measurements are given in physical units; where cosmological parameters enter calculations, these are based on the standard cosmology used for the simulations: , , , , , . We briefly investigated the effect of incorporating the favoured parameters from the fifth year WMAP results (Dunkley et al., 2008), but the overall differences are expected to be minor (Section 5.3).

2 Simulations

Tag Usable Vol (comoving)
Table 2: The simulations used in this work. The first column is the tag which we use to refer to each simulation. For all except “Cosmo”, a subsample of the full box is simulated in high resolution; no results are taken from outside this region. The second and third columns refer respectively to the mean gas and dark matter particles within the region, the fourth to the gravitational softening length (in physical units) and the final column gives the comoving volume of the region. The separate boxes are generated from entirely different sets of initial conditions; the “Dwf” and “MW” simulations are designed to form respectively a dwarf and Milky Way type galaxy at while the “Large” and “Cosmo” boxes are more statistically representative.

The simulations are successors to the galaxy formation simulations described in G07, with higher resolution and a more physical star formation feedback prescription as described by Stinson et al. (2006) (S06, see also below). They are computed using the SPH code Gasoline (Wadsley et al., 2004), and include gas cooling and the effects of a uniform ultraviolet (UV) background (following Haardt & Madau, 1996) in an optically thin approximation. We later perform post-processing to account for self-shielding effects (Section 3.1). We also ran test simulations which included an approximate treatment of self-shielding within the live simulation (Section 5.4).

The SPH smoothing lengths are defined adaptively to use the 32 nearest neighbours in all averaging calculations, except that a minimum smoothing length of times the gravitational softening () is imposed. The star formation and feedback recipe (S06) is based on the algorithms originally described by Katz (1992). In short, gas particles can only form stars if , and the local hydrodynamic flow is converging. The rate of star formation in such regions is assumed to follow the Schmidt (1959) law () with index . This is a rather natural choice of , since it implies the cold gas is turned into stars over some multiple of the local dynamical timescale:


Following the evolution of each star particle consistently with the Kroupa et al. (1993) IMF, a fixed fraction of the supernova energy created at each timestep is deposited thermally into the surrounding gas particles. To emulate the physics of the processes responsible for distributing this energy to the gas, radiative cooling is disabled in particles within a radius which must be determined. Traditionally this can involve further parameters, but the S06 method obviates the need for these by modeling the physical processes according to a prescription based on blast wave models (e.g. Chevalier, 1974; Ostriker & McKee, 1988): this sets the radius of the local ISM affected as a function of the local density and temperature. The free parameters, consisting of the constant of proportionality in the Schmidt law () and the efficiency of supernova energy deposition into the ISM () are tuned such that isolated galaxy models match the Kennicutt (1998b) law and cold gas fractions observed in present day disk galaxies. This leaves no free parameters in the context of this study (setting and – see Stinson et al. 2006) but produces galaxies which satisfy a wide range of observational constraints (see Introduction). As part of the supernova algorithm, metals are also deposited into the ISM (see Section 4.2.3 for more details), but note that our simulations do not yet contain any contribution to radiative cooling from metals (which is expected to be physically dominant below temperatures ).

As in G07 & G08, the prescriptions are implemented in a fully cosmological context, based on a fiducial model with parameters given at the end of Section 1. The majority of our results are derived from four simulations (Table 2). Three of these take advantage of the “volume renormalization” technique, in which progressively higher resolution is employed towards the central galaxy of interest (Katz & White, 1993). The “Dwf” and “MW” simulations simulate the same region as the DWF1 and MW1 runs in G07 and G08, forming at a dwarf and Milky Way type galaxy respectively. The third simulation, “Large”, is a compromise between high resolution and volume, while the fourth simulation, “Cosmo”, is of a full box populated with gas at uniform resolution. The smoothing length is constant in physical units for , and evolves comovingly for ; its final fixed value for each simulation is given in Table 2. Otherwise, the physics in each of these runs is identical.

3 Processing Pipeline

Figure 1: A two-dimensional representation of our three-dimensional radiative transfer scheme, which involves building a nested grid and integrating optical depths from a cell along that grid to the edges of the box. A combination of local fine spatial resolution and speed is achieved by traversing the grid at the lowest available level until this level is closed. At such a point, the algorithm jumps to a higher level in the grid. Subsequently, when lower levels become available they are ignored since these high density regions are not directly connected to the original high density region. This is primarily a computational simplification, but it also (in an admittedly ad hoc fashion) prevents long distance unphysical shadowing which would otherwise arise from our limited angular resolution – see text for details.

The processing pipeline is based on SimAn111www.ast.cam.ac.uk/~app26/siman/ , an object-oriented C++/Python-based environment for analyzing simulations of arbitrary file format with support for OpenGL real-time visualization. Such a framework allows us to understand our data rapidly, especially since it includes support for stereoscopic glasses – giving a true 3D rendering of the data. It is written in such a way as to allow analysis of data from an interactive python session, hiding a large number of details from the scientific user.

3.1 Self-Shielding

We employ a basic radiative transfer post-processor to assess the effect of self-shielding. The simulation is divided into nested grids so that the lowest level cells contain close to 32 particles, to correspond with the SPH scheme of G07. We use a cosmic ultraviolet background (UVB) following a revised version of Haardt & Madau (1996) (Haardt 2006, private communication); this is initially assumed to be present throughout the simulation volume.

In each cell, using the radiation field described above, the ionization state of the hydrogen and helium gas is determined using an equilibrium solver algorithm based on Katz et al. (1996a). This requires the temperature and density of each SPH particle, derived directly from the simulation output. We keep the temperature, not the internal energy, fixed during this process. While not ideal, it is necessary to make some such assumption to determine the thermal state in post-processing. As a test we verified that keeping instead the internal energy fixed makes little noticeable difference to our results. On the other hand, shielding can significantly drop the heating rate, so that an assessment of the dynamical effect is essential; we find it to be minor in terms of our final results, although uncertainties remain (see Sections 5.2 and 5.5).

With the ionisation state determined, the attenuation of the UV field due to the  H i, He i and He ii ions is then calculated. This translates into an optical depth (as a function of frequency) for the cell. In subsequent iterations, the attenuation of the incoming UVB radiation is calculated by integrating the cell optical depths along the six directions defined by the orientation of the grid to the edge of the box. Because the grid is refined adaptively, this process is somewhat complicated by the need to “level-jump”, i.e. move to higher levels of the nested structure when the lower levels run out. Note that the grid walk stays at the lowest level in any directly connected region of the origin cell (including when crossing boundaries of higher levels). As it moves out of the high resolution region, higher levels are selected as appropriate. The walk does not re-descend, even if lower levels again become available, instead using the averaged properties of regions defined by the current level. This is primarily a matter of computational speed; however, it has the side-effect of preventing long distance shadowing which can arise in low angular-resolution radiative transfer mechanisms. Admittedly we have not formulated this in terms of any limiting procedure yielding a well-defined physical calculation, but heuristically the averaging over larger solid angles as the integration proceeds away from the target cell is quite correct. We have illustrated this situation in Figure 1.

This entire process is iterated over the full simulation, convergence being assessed by changes in the UVB field and ionization state between steps. We found that the system defined above had some oscillatory behaviour on its approach to convergence, which led us to introduce a damping term which in effect averages the optical depth between iterations. This results in faster convergence, but does not essentially alter the scheme described.

Our self-shielding process is crude when compared to recent radiative transfer codes (for a comparison see Iliev et al., 2006); however, we do not believe this makes our results obsolete: see Section 5.4. A further complication must be considered, which is the failure of the scheme as described to account for any dynamic effects of the changing ionisation fraction and heating rates. We have assessed the severity of this problem using a local attenuation approximation: see Section 5.4.

Figure 2 shows a map of the neutral hydrogen in a (physical) cube centred on the major progenitor of a Milky Way type galaxy in our “MW” box. It is coloured such that DLAs appear in red and Lyman limit systems appear in green/yellow. The locations and virial radii (see Section 3.2, below) of dark matter haloes exceeding in virial mass are overplotted.

3.2 Calculation of Halo and Sightline Properties

We start by locating individual dark matter haloes in the simulation using the grid-based code AHF222http://www.aip.de/People/aknebe/AMIGA (Knebe et al., 2001; Gill et al., 2004). We then mark a spherical region of radius (the radius within which the density exceeds the mean density of the simulation at the given redshift by a factor 178, by analogy with spherical top-hat collapse models) as belonging to the halo of mass . We also record the masses in gas, stars and  H i  within the virial radius. The virial velocity is defined as usual, .

In the resampled simulations, we immediately discard any haloes which have been contaminated by low resolution particles from the outer regions. We also discard any haloes with fewer than 200 gas particles or 1 000 dark matter particles. We determined this limit empirically by examining the point at which halo properties diverged in simulations at different resolutions – for more details, see Section 5.1.

Figure 2: The neutral column density of  H i  in a cube centred on the major progenitor to a Milky Way type galaxy (box MW). The colours are such that DLAs () appear in dark red and Lyman limit systems () appear in green and yellow. The circles indicate the projected positions and virial radii of all dark matter haloes with . All units are physical. A stereoscopic version and an animated version of this plot are available at www.ast.cam.ac.uk/~app26/.
Figure 3: Example low-ion (Si ii) absorption profiles from random sightlines through selected haloes: in reading order, these have virial velocities of , , , (making the first halo extremely rare in cosmological sightline samples, see Section 4.1). The first is taken from our “Cosmo” box, the second is the “MW” major progenitor and the remaining two are smaller haloes from the “MW” simulation. Note that the velocity axes are scaled differently in the respective plots. The zero velocity offset corresponds to the motion of the centre of mass of the halo (unobservable in practice), illustrating a qualitative shift in behaviour from multiple clumps (higher virial velocities) to single clumps (lower virial velocities). The vertical dotted lines indicate the velocity offsets where the cumulative optical depth reaches and of its maximum value; is given by the difference in their position in velocity space. The pixel size is approximately , although internally we use a higher resolution. For the purposes of this plot, the profiles are normalized to correspond to a definite transition, although this is not necessary for our computations. This chosen transition, along with the sightline  H i  column density and metallicity, is indicated in each panel. Noise is added to simulate S:N ; again this is only for illustrative purposes and is not part of our pipeline.

Line-of-sight properties can be calculated from particle-based simulations by projecting quantities onto a grid. For instance, projecting all gas particles in a simulation onto the plane allows the column density along the direction to be estimated by summing the mass in a grid square and dividing by its area. This is the approach taken in previous SPH simulations of DLA properties. However, in initial numerical experiments we found that sightline properties in our simulations were not robust to changes in the somewhat arbitrary grid resolution.

We have instead calculated all quantities using a true SPH approximation; for details, see Appendix A. This results in a varying spatial resolution of sightlines which is automatically consistent with the simulation data. The minimum smoothing length allowed is times the softening (given in Table 2), meaning we can resolve spatial gradients over as little as in our highest resolution simulation, although a more typical effective resolution is .

For our main results, sight-lines are projected through the simulation in random orientations and at random sky-projected offsets from the centre of a halo up to its virial radius. We verified that extending this search area to twice the virial radius had no impact on our results (this can also be seen directly from Figure 2). This confinement of our DLA cross sections is discussed in Section 5.5.

For each sightline, we measure directly the column density in neutral hydrogen. If this column density exceeds the DLA threshold (), it is added to our catalogue. If not, it is immediately discarded; but we keep track of the numbers of all sightlines taken so that we may calculate


where is the search area, is the total number of random sightlines calculated and is the number of such sightlines which exceed the DLA threshold. In this way, we obtain a representative DLA cross-section for each halo without assuming any particular projection.

We also produce an absorption line profile for a low-ion transition such as those of Si ii. We assumed that the relative abundances of heavy elements were solar and that Si ii  was perfectly coupled to  H i, so that for solar metallicity and (Lodders, 2003). In other words, given the metallicity of each gas particle, we take . Although an approximation, we found that the effect of relaxing the assumption of the Si ii  –  H i  coupling was minor; see Section 5.2.

Example profiles are shown in Figure 3, in which we have chosen four haloes and displayed five random sightlines from each. For the purposes of this plot, we choose one of SiII1808, 1526, 1304 or 1260 according to which transition has maximum optical depth closest to unity. The plots are centred such that corresponds to the motion of the centre of mass of the parent dark matter halo. We also added, for display purposes, gaussian noise such that the signal-to-noise ratio (S/N) is . The range contributing to (see Section 4.2.2) is shown in Figure 3 by vertical lines at either end. Absorption arising in haloes with virial velocities is often composed of multiple clumps whereas for smaller haloes there tends to be one main, central  H i  clump, moving with the centre of mass of the halo. Visual inspection of the different systems suggests that the former systems are less dynamically relaxed, presumably because they have formed more recently and have longer dynamical timescales; however, we did not verify this systematically.

3.3 Cosmological Sampling

We aim to make statements about the agreement or otherwise of our simulations with cosmological observations of DLAs and therefore need to construct a representative global sample of absorbers. Our approach is related to that of Gardner et al. (1997a), in that we correct our limited sample by reweighting in accordance with the halo mass function. However, we emphasize that in our case this is to make allowance for our limited statistics and combine results from separate boxes, whereas in Gardner et al. (1997a) this approach was used to extrapolate behaviour into an unresolved low-mass regime.

Consider any measurable property of DLA absorbers, . The aim of the process discussed below is to construct the distribution function, . (We will adopt the custom of using the absorption distance , where : populations with constant physical cross sections and comoving number densities maintain constant line densities as they passively evolve with redshift.) Given the monotonic invertible relations and , one has


where is the cross section of a DLA of mass and is the halo mass function333We adopt the numerically calibrated version of the Sheth & Tormen (1999) halo mass function given by Reed et al. (2006)., so that the physical number density of haloes with mass in the range is , and . However, for an ensemble of haloes of the same mass, will be scattered – although it may be possible to define a mean value . Similarly, will in fact be scattered around some suitably chosen average for a given halo mass; furthermore even for a single halo most properties will vary depending on the particular sightline taken through the halo to the distant quasar444We regard both of these effects as stochastic scatter, although presumably a complete theory would account for the exact value of in terms of a sufficient number of parameters pertaining to the halo and sightline..

The easiest approach is to make the replacements , in equation (3). However, even if significant deviation of from is rare, one may be interested in variations of over several orders of magnitude – these rare fluctuations can therefore contribute significantly to the “tail” of the observed values.

Accordingly, the method by which we perform the reweighting is non-parametric. Since at first this can seem a little obscure, we offer two descriptions. Below, we have described the method in a heuristic manner. In Appendix B we have outlined how this method can be derived by discretizing a well defined integral, which allows for an exact interpretation of our results should one be necessary.

First, we bin haloes (in logarithmic bins) by their virial mass. Within each bin , ranging over virial mass , one may calculate a mean DLA cross-section and a total physical number density :


The contribution of DLA systems per unit physical length from the bin is , and consequently one may calculate:


To investigate the distribution of properties observed we have two approaches. The simplest is to extract a representative set by discarding selected sightlines. This is ideal for comparing correlations between sightline parameters (e.g. velocity & metallicity, Section 4.2.3), where the observed data sets consist of only separate sightlines. For each halo , we need to select a number of sightlines proportional to (cf equation (33)) where is the cross-section of the specific halo, represents the mass bin to which halo belongs and is the total number of simulated haloes in the mass bin . The actual sightlines chosen from each halo are determined by a pseudo-random but deterministic approach, for stability of results.

Figure 4: The DLA cross-section of haloes which meet our resolution criteria in the Dwf (plus symbols), MW (dots), Large (cross symbols) and Cosmo (tripod symbols) boxes plotted against their virial mass. There is a (resolution-independent) sharp cut-off at below which the cross-section for DLA absorption is negligible. Haloes with no DLA cross-section are shown artificially at . The fit to the equivalent results for two models in Nagamine et al. (2004a) is given by the dotted and dash-dotted lines (their models P3 and Q5 respectively). The major progenitors to the “MW” (Milky Way like) and “Dwf” (Dwarf type) galaxies are indicated.

This method is not ideal, however, because it throws away potentially useful information. In particular, when constructing distribution functions for a property , we use an alternative method which uses the cosmological halo mass function to weight, rather than select, results. The set of all sightlines is binned by the property , indexed by . One then has


where the sum is over all sightlines from any halo in any box, denotes the halo associated with the sightline and is the bin size. The constants are not hard to calculate, but obscure the technique and are presented in full in Appendix B, equation (31).

Of course, the above methods require that line-of-sight effects (such as the coalignment of multiple DLAs along a single sightline) are not a dominant effect. We verified that extending all our sightlines through the entire boxes made no significant difference to any of our results. This is because the DLA cross-sections per halo are small so that even when strong clustering is taken into account, the probability of a double-intersection along a line-of-sight is very low. (Note also that velocity widths, Section 4.2.2, are always determined from unsaturated line profiles, so that trace metals in the IGM are too weak to change the measured widths.)

One should also require that environmental variations of halo properties (in effect systematic variations with parameters other than halo mass) are unimportant. This requirement is harder to verify, but we neither detect variations of DLA properties of differing regions in our “Cosmo” box, nor any strong correlations with indirect measures of environment such as the halo spin parameters. More systematic work in the form of the Gimic project, resimulating carefully chosen regions of the Millennium Simulation (Springel et al., 2005), also suggests that environmental considerations are largely controlled by the local variation of the finite volume halo mass function (Crain et al., in prep).

4 Results

4.1 Masses of Haloes Hosting DLAs

Figure 5: The data of Figure 4 multiplied by the halo mass function from Reed et al. (2006) to give a total line density for each representative system. Haloes with zero cross-section are shown artificially at .

Taking all haloes of sufficient resolution (Section 5.1), we plot their DLA cross-section against virial mass in Figure 4. Our DLAs have cross-sections of between and . In agreement with rough physical expectations, the cross-section increases as a function of mass. In their regions of overlap, the results agree between boxes: although the Cosmo box has a vastly larger volume than our other boxes, and hence exhibits more scatter by probing rarer haloes, we verified that by restricting the number of haloes at a given mass the distributions of haloes from differing simulations were in agreement in each mass bin. Two extreme outliers in the Cosmo box were individually inspected and found to be systems in the middle of major mergers.

A notable feature is that the cross-section drops off extremely quickly for (i.e. ). This behaviour is accompanied by a sharp break in the mass of neutral hydrogen in such haloes (Figure 6), and we attributed it to the UV field which prohibits cooling of gas in lower mass haloes (see Rees, 1986; Efstathiou, 1992; Quinn et al., 1996, and Section 5.5).

Figure 6: The total mass of neutral hydrogen plotted against the virial mass for our haloes. The symbols are as described in the caption of Figure 4.

Our cross-sections are overall somewhat larger than previous simulation works have suggested (e.g. Nagamine et al., 2004a; Gardner et al., 1997b) and furthermore are not compatible with a single power-law link to the halo mass. The results of two models from Nagamine et al. (2004a) are overplotted on Figures 4 and 5 for comparison (runs “P3” and “Q5” here refer to weak and strong galactic winds in the cited work). Our enhanced cross-sections for are apparently a consequence of our particular feedback implementation – see Section 5.1.

From our cross-sections, we can calculate directly the line density of DLAs via equation (5). For this, we find , in good agreement with the observationally determined value from SDSS DR5555www.ucolick.org/~xavier/SDSSDLA/DR5/ using the method described for SDSS DR3 in Prochaska et al. (2005). In Figure 5, we have shown how haloes of different masses contribute to this total line density by plotting for each halo; is simply the integral under the curve defined by the locus of these points. The major contributors to the total line density are haloes of masses . At lower and higher masses the contribution is cut off by the rapidly decreasing cross-sections or exponential roll-off in the halo mass function respectively.

As expected from our previous discussion, our results have a peak at which contrasts with the flatter results from fiducial power-law cross-sections. Consequently, at first glance, it appears that the area under our locus of points must be larger than that under the N04 curves and hence a substantial disagreement in line density is inevitable; however the cut-off for N04 is at rather lower masses (), and this brings the total line density in N04 considerably closer to the observed value.

Plotting the total  H i mass against the virial mass (Figure 6) and DLA size against the total  H i mass of a halo (Figure 7) gives an alternative view of our cross-sections. A striking feature of the latter plot is a bifurcation, particularly notable in the “Cosmo” box but also traced by the “Large” box, in which haloes of a fixed  H i mass can have different cross-sections. The primary physical distinction between the upper and lower branches is in halo mass: the former trace  H i-rich haloes with and the latter trace a population of  H i-poor haloes with . This is reminiscent of recent claims of bimodality in observed DLAs (Wolfe et al., 2008). However we limit ourselves, for the moment, to more general considerations, noting that the high mass, low branch DLAs are extremely rare in our simulations (making up less than 2% of the total cross-section). Our method for generating cosmological samples (Section 3.3) will propagate through any such bimodalities into our final results without difficulty, assuming the Cosmo box has a representative selection of haloes. See Section 6.5 for a further discussion of bimodality.

Figure 7: As Figure 4, except the cross-sections are now plotted against the total mass of neutral hydrogen in each halo. Haloes with no DLA cross-section are shown artificially at . The dotted lines are of timescales for  H i  depletion through star formation; from top to bottom , and yr. These illustrate constraints on our locus of points by considering both and the lifetime of a typical DLA (see text for details).

A guide relationship can be estimated by fixing a particular star formation rate. As explained in Section 2, in our simulations the instantaneous star formation rate is given by the Schmidt law. From this may be estimated a timescale for conversion of neutral gas into stars,


With the assumption that , one has the approximate relation


Our locus lies roughly along the line years; this is not unexpected, especially if our simulations are to match observations that .666In other words, the total mass of stars at in a given comoving volume is roughly equal to the total mass of neutral hydrogen in that same volume at ; see equation (9) and the ensuing discussion. This suggests that a large fraction of the DLA cross-section should be converted to stars by , giving an upper limit of to the star formation timescale (unless star formation proceeds in rapid discrete bursts, which is not the case in our simulations). Assuming DLAs are not short-lived objects, or achieved by very fine balancing of rapid gas cooling and star formation, one would also expect . The constraint is obeyed, suggesting that this stable model is reasonable.

4.2 Column Densities, Velocity Widths and Metallicities

4.2.1 Column Density Distribution

Figure 8: The simulations’ DLA column density distribution (solid line) compared to the observed values from SDSS DR5 (points with error bars, based on Prochaska et al. 2005; see main text for explanation). The dashed, dashed-dotted and dotted lines show the contribution from , and haloes respectively (these are not directly observable distributions, but give guidance as to how our cross-section is composed).

One of the best constrained quantities, observationally, is the neutral hydrogen column density distribution . This is defined such that gives the number of absorbers with column densities in the range and absorption distance . Applying the reweighting method described in Section 3.3 to our sample yields an estimate for the cosmological column density distribution, shown by the solid line in Figure 8. This can be compared directly to the observed distribution given by the points with error bars, which are derived from SDSS DR5 (see previous section for an explanation). The matching of the normalization and approximate slope of the observed column density distribution can be seen as a genuine success of the simulations: we emphasize that no fine tuning has been applied to achieve this result. Furthermore, our results appear to have converged at the resolution of the simulations used (see Section 5.1).

From the column density distribution, one may express the total neutral gas mass in DLAs in terms of the fiducial definition


where is the proton mass, is the critical density today, gives the fraction of the gas in elements heavier than hydrogen and is an upper limit for the integration, which is discussed in the next paragraph. gives the fraction of the redshift zero critical density provided by the comoving density of DLA associated gas measured at redshift . (This is different from the more natural definition of time-dependent s which express a density at any given redshift in terms of the critical density at that redshift. Only in the Einstein-deSitter universe will these definitions coincide.)

Although the calculation should take , this is not possible for the observational sample owing to the rapidly decreasing number of systems at the high column density limit. Prochaska et al. (2005) discussed how different assumptions for the functional form of the column density distribution can lead to different values of . The discrepancies are small for the two best functional fits to the observational data (a double power law or a Schechter function with exponential roll-off at high column densities). However, these extrapolations are actually only constrained by a few points at high column densities; a more robust approach – albeit less physically transparent – is to calculate directly from summing the total neutral hydrogen in the observed sample of DLAs, which for the SDSS DR5 sample is roughly equivalent to using the upper limit .

Using this limit, we obtain , which can be compared with the result from SDSS DR5 in the combined bin , . As expected from Figure 8, the results are in fair agreement; the slight mismatch is driven by the overestimation of our high column density points ().

Our weighting approach predicts the form of the distribution for much rarer, higher column density systems than the sample limited observations allow. Significant contributions to are made by these rare systems unless . In fact, directly measuring the slope for our simulations shows that it slowly decreases from for to a constant value of for . Thus a correction to is expected if we allow to extend to arbitrarily high values. Performing the calculation with gives . This value is not directly comparable to observational estimates, but shows that an observer living in our simulations would underestimate by about due to missing contributions from the rare high column density systems. A further discussion of this issue is given in Section 6.4.

4.2.2 Velocity Width Distribution

Figure 9: The simulations’ DLA velocity distribution (solid line) compared to the observed values based on the sample described by Prochaska et al. (2003) (see text for details; shown by points with error bars). The dashed, dashed-dotted and dotted lines are as described in the caption of Figure 8.

We have already discussed some qualitative features of the low-ion velocity profiles generated in our simulations (Section 3.2, Figure 3). These are important because they provide a direct observational measure of the kinematics of the DLAs, and therefore have the ability to substantially constrain the nature of the host haloes. We now turn to the comparison of our characteristic velocities with the observed quantitative distribution.

We assign a velocity width to each generated profile using the fiducial technique (Prochaska & Wolfe, 1997). This inspects the “integrated optical depth” and assigns the velocity width where and . The result is a representative velocity width for the sightline, produced without any of the difficulties associated with fitting multiple Voigt profiles.

To a good approximation, the only dependence on the particular low-ion transition chosen is in the overall normalization of the optical depths from the relative abundances and oscillator strengths. The measure of the velocity width is invariant under such rescalings, so that only for display purposes (i.e. Figure 3) do we need to choose a particular ion and transition. Observers require to choose unsaturated lines because, while our simulations calculate directly, spectra only determine which cannot be inverted (in the presence of noise) for .

Using our weighting technique, we calculate the cosmological distribution of velocity widths, . This is plotted in Figure 9 (solid line), along with observational constraints (points with error bars) calculated from data provided by Jason Prochaska (private communication, 2008) based on the compilation of high resolution ESI, HIRES and UVES777Echellette Spectrograph and Imager on the Keck Telescope; High Resolution Echelle Spectrometer on the Keck Telescope; Ultraviolet and Visual Echelle Spectrograph on the Very Large Telescope spectra described in Prochaska et al. (2003). The dataset consists of 113 observed DLAs with (we directly verified that the velocity width redshift evolution over this range is negligible, so that we use the full range of observed systems to bolster our statistics). We normalize the line density of DLAs in the observational sample to match the observed (Section 4.2.1).

Figure 10: The data of Figure 9 (thick solid line) plotted cumulatively for comparison with similar plots in the literature. The dotted and dash-dotted lines show the Razoumov et al. (2008) models N1 and H1 respectively (see text for details). The shaded band gives the approximate Poisson errors on the observational data, but we caution these are in fact strongly correlated.

Overall, the simulations reproduce the approximate pattern of observed velocity widths, with a peak in the distribution at ; the agreement is fair for velocity widths . However, it is now a familiar feature of CDM simulations that they produce too few high velocity width systems (see Introduction) and our results are no exception. We provide a discussion of this discrepancy in Section 6.3.

Recent velocity width results by Razoumov et al. (2008) and Razoumov et al. (2006) have been presented by plotting a cumulative line density, . A direct comparison between our results and those from such work can be drawn from Figure 10, in which we similarly replot our velocity width distribution cumulatively (thick solid line). As well as the observational data (thin solid line) we have overplotted results from two simulations described by Razoumov et al. (2007; R08). These two models, R08.N1 and R08.H1, differ in their size and hence resolution; the former covers a volume of approximately (comoving), resolving grid elements of minimum side length (physical) while the corresponding values for the latter are and . Comparing our simulations with those of R08, the extent of the disagreement between simulated and observed results is similar (discounting the lack of low velocity systems in R08.H1, which arises from the coarser resolution), although our own results are actually somewhat closer for low velocity widths (). This is possibly linked to our somewhat higher than normal cross-sections for haloes of virial velocities , but because much of R08’s DLA cross-section lies outside haloes it is hard to provide a concrete explanation for such disparities (see Section 5.5).

We caution that any cumulative measure has strongly correlated errors so that the plot can present a somewhat distorted picture of the discrepancy (one is really interested in its gradient in linear space). In our plot we have shown the Poisson errors on the observed data as a shaded band, but these only represent the diagonal part of the covariance.

4.2.3 Metallicity

Figure 11: The simulated distribution of metallicities (thick solid line). The points with error bars show the observed metallicities from the sample of DLAs based on Prochaska et al. (2007), normalized to the observed . The dashed, dashed-dotted and dotted lines are as described in the caption of Figure 8.

The metallicity, i.e. the ratio by mass of elements heavier than helium to hydrogen, is an important diagnostic of observed DLA sightlines. Since metals are deposited in the interstellar medium (ISM) through supernova explosions, the metallicity of a region is determined by the interplay between the integrated star formation rate and bulk motions of gas including galactic inflows and outflows. In general, one observes a positive correlation between the mass of a galaxy and its metallicity both at (e.g. Tremonti et al., 2004; Lee et al., 2006) and at higher redshifts (e.g. Savaglio et al., 2005; Erb et al., 2006). This trend is also observed in our simulations (Brooks et al., 2007, henceforth B07). Debate over a precise account of the origin of the relation continues, but the analysis of B07 shows that, for our simulations, the predominant effect is reduced star formation in low mass haloes, with the dynamics of outflows providing an important correction to the simple closed box model.

There are many uncertainties in simulating metallicities, which reflect not only the mass formed in stars but also the dynamics of the gas, its accumulation from the IGM and mixing within the halo. B07 showed that, for our simulations, high resolution is required () to attain convergence for the star formation history (and consequently the metallicities)888This requirement is weaker than the analagous result in recent work by Naab et al. (2007) who require particles for convergence; this probably relates to the neglect of feedback effects in Naab et al. (2007) – feedback in our model tends to prevent collapse to very high densities, and thus imposes an effective star formation scale independent of the resolution; see also Section 5.1. . We verified (using an independent analysis code) that our own results suffered similarly, but that our other results converge at lower resolutions (see Section 5.1). Whenever we deal with results concerning metallicities, we therefore have to discard a number of haloes (those with ) which are included elsewhere. The effect of this is to reduce the number of haloes in each mass bin, and therefore exacerbate worries that we do not have a fully representative set of haloes for all mass scales. Nonetheless, we do not find any evidence that this is causing systematic effects, as we verified that our other distributions are not significantly changed by this restriction.

Each of our sightlines is assigned a metallicity as follows. The simulations keep track of two sets of elements, the -capture and Fe-peak elements. The enrichment process follows the model of Raiteri et al. (1996), adopting yields for Type Ia and Type II supernovae from Thielemann et al. (1986) and Weaver & Woosley (1993) respectively. For our purposes the differences between the two metal groups are minor: observations show corrections are typically (e.g. Ledoux et al., 2002), and exploratory work suggested our simulations were similarly insensitive. We choose to ignore these differences for the sake of simplicity and use the total mass density in metals. For quantitative results, we assume a solar metallicity fraction of by mass (Lodders, 2003).

A cosmological distribution of DLA metallicities is generated using our fiducial technique (Section 3.3). The result is shown by the thick solid line in Figure 11 and closely matches the observational constraints, displayed as points with error bars. These are derived from the Prochaska et al. (2003) sample described in Section 4.2.2, restricted to . We use metallicities based on -capture elements (primarily Si), updated using data from Prochaska et al. (2007).

Our simulated results are in good agreement with the observed distribution. They exhibit a roughly gaussian distribution (as a function of log metallicity); the best fit parameters give a median metallicity of with a standard deviation . A similar fit to the observed data gives and . Given the uncertainties in abundances and yields, the differences are extremely minor. This success is unusual: previous simulations (e.g. Cen et al., 2003; Nagamine et al., 2004b) tended to substantially overproduce metals in DLA systems, and were only been able to approach the observed result if metal rich DLAs could be hidden using substantial dust biasing. However, such a scheme is unsupported by observational evidence, and our simulations now suggest it is unnecessary; for a further discussion see Section 6.4.

We also investigated the relationship between the mean metallicity of the gas within the virial radius of a halo and the spread of DLA metallicities derived from SPH sightlines through that halo. We found that the sightlines on average displayed slightly higher metallicities than the halo mean (by ) with a spread of . We interpret the offset by noting that DLA sightlines sample the high density gas in which star formation is occurring – even if the rate is low (Section 4.3), presumably the local star formation preferentially enriches these. Some haloes showed a detectable radius-metallicity gradient, but in general this was shallow compared with the overall scatter.

4.2.4 Correlation between Metallicity and Velocity Widths

Figure 12: The relationship between metallicity and low-ion velocity width of individual sightlines in a cosmologically weighted sample from our simulation (crosses, with linear least square bisector given by the solid line) and from the observational dataset described in the text (dots, with linear least square bisector given by the dotted line). The qualitative agreement, given the known deficiencies in the simulations and the lack of fine tuning, is remarkable. In the simulations, the origin of the correlation is an overall mass-metallicity relation; see text for further details.

Observationally, there is known to be a correlation between the low-ion velocity width and the metallicity of a DLA sightline. Ledoux et al. (2006) presented a set of observations showing a positive correlation between these parameters at the level and suggested that the relation could reflect a mass-metallicity correlation analogous to that seen in galaxy surveys; the result was confirmed, using a separate sample, by Prochaska et al. (2008).

Our matching of the metallicity relation and near-matching of the velocity width distribution gives us confidence to attempt to probe this relationship in our simulations. We employ the Monte-Carlo sample generator version of our halo mass function correction code (Section 3.3) to produce a sample of coupled velocity and metallicity measurements, matching the size of the observational comparison sample described in Section 4.2.2 restricted to as in Section 4.2.3.

The simulation results are shown as crosses in Figure 12, with the linear least-square bisector fit999The linear least-square bisector method estimates the slope of a relationship between and as the geometric mean of the slopes obtained by regression of and ; see Isobe et al. (1990). given by the solid line. The observational sample is shown by dots (the errors on each observation are small compared to the intrinsic scatter), with the linear least-square bisector fit shown as a dotted line. These fitted relationships are parametrized as and respectively.

Although the normalization of the relationship in our simulation is somewhat different from the observational sample, we emphasize that the slope and overall trend, as well as the mean metallicity of our sightlines, are correct and that qualitatively the results are in agreement. Given that the metallicity distribution (Figure 11) is in close agreement with that observed, we suggest that the quantitative disagreement arises from the previously noted underestimate of the velocity widths in our simulation (i.e. one should interpret the discrepancy in the relationships in Figure 12 as a horizontal, not vertical, displacement). It is also apparent visually that the observational results show a larger scatter about the mean relationship than does our computed sample. We tentatively suggest that this points towards an ultimate resolution of the velocity width issue which involves “boosting” the width of certain sightlines, perhaps by local effects such as outflows, while keeping the contribution by halo mass roughly as outlined in Section 4.1 (although this is not the only conceivable interpretation: see Section 6.3).

The origin of our relation between velocity widths and metallicity is the underlying mass-metallicity relation, as suggested by Ledoux et al. (2006). We verified this directly, but it can also be seen from Figures 9 and 11: the dashed, dash-dotted and dotted lines in each case show the contribution from haloes with , and respectively. Both the sightline velocities and the metallicities can be seen to be a strong function of halo mass, and it is this fact that leads to the final relationship.

4.2.5 Other correlations are weak

The correlation between column density and metallicity or velocity is weak in our simulations, in agreement with observations: contrast the strong dependencies of the velocity widths and metallicities on the underlying halo mass with the weak dependence of the column density (Figure 8). Given the historical confusion over the evidence for correlation between the column density and metallicity, we have investigated this aspect of our observational and simulated samples in a separate note (Pontzen & Pettini, in preparation).

4.3 Star formation rates

Figure 13: The star formation rates in our haloes, plotted against their virial mass. Also shown is the observationally determined characteristic LBG star formation rate from Reddy et al. (2007), and a very approximate range of star formation rates for DLAs derived using the C ii] cooling rate technique (Wolfe et al., 2003) – for more details and a list of caveats, see main text. An empirical power law fit is shown. Haloes with no star particles are shown artificially at .
Figure 14: Using our non-parametric weighting technique, the distribution of star formation rates in DLAs within our simulation is shown, with the same observed ranges indicated as in Figure 13. There is qualitative agreement between our estimates for DLA star formation rates and the observational constraints.

In the Introduction, we noted that star formation rates (SFRs) in DLAs are typically low () and yet these systems contain most of the neutral hydrogen, a necessary intermediary in the star formation process. We also see this behaviour in our simulations, and now describe how it comes about.

We associate a star formation rate with each halo by inspecting, at , the mass of star particles formed within the preceding years. For very low star formation rates (), there may have been no star particles formed in such a period; in this case, we extend the averaging period by a variable amount up to years so that we can resolve star formation rates down to .

The SFRs of our individual haloes are shown in Figure 13; the horizontal dotted line shows the star formation rate of an LBG galaxy. This is derived from Reddy et al. (2007), wherein is stated . We adopt the conversion ratio of Kennicutt (1998a), but first correct the luminosity for dust attenuation by a factor of , which is the mean correction given (Reddy et al., 2007, section 8.5). We then divide the final result by in order to consistently use the Kroupa IMF (Section 2), for which a larger number of massive, hot stars are formed relative to the Salpeter IMF assumed by Kennicutt (1998a). This gives a final characteristic LBG star formation rate of .

We also plot a shaded band indicating the approximate SFRs obtained for DLAs using the C ii] technique (Wolfe et al., 2003). This is intended to serve as a guide only – there are a number of uncertainties in the position of this band, including the intrinsic complexity of the observations and the conversion from a star formation rate per unit area (we have assumed a DLA cross-section of , which is our approximate range for the most common DLA haloes – see Figures 4 and 5).

By evaluating


(using a binning technique to discretize the integral) we derive a global star formation rate at of . The UV dust-corrected and IR estimates in Reddy et al. (2007) place this value between and , but these are again sensitive to the IMF, assuming a Salpeter form for their main results: for a Kroupa IMF one should again reduce the rate by . This caveat should be borne in mind, but overall our results are not unreasonable and a detailed understanding of remaining discrepancies is beyond the scope of the present work.

More importantly for our purposes, the star formation rate roughly scales as (see fit in Figure 13). Given the low mass end of the halo mass function is an approximate power law, , the overall dependence of the integrand of equation (10) is rather shallow and a wide range of halo masses contribute to the global star formation rate. Thus, there is no inconsistency in the view that a typical () DLA has a low star formation rate but that the DLA cross-section as a whole contributes significantly to the star formation, and is largely converted into stars by .

We generate a cosmological DLA cross-section weighted sample of these star formation rates (Section 3.3). Figure 14 shows the resulting distribution, . As expected by the range of halo masses contributing to the DLA cross-section (Figure 4), the star formation rate in most DLAs is much lower than that in visible LBGs. Qualitatively, the simulated rate of DLA star formation is consistent with or fractionally lower than the observed star formation rate estimates from the C ii technique (see above). Since we have not included upper limits in our approximate observational band, this is an acceptable result.

Figure 15: The distribution of stellar masses of our DLAs. The vertical band and vertical dotted line show respectively the range of and median value for observed LBG stellar mass from Shapley et al. (2001) (adapted for consistency with our IMF). As expected, most DLAs have formed a relatively small population of stars, consistent with their low metallicities.

Finally, we investigated the total stellar mass accrued in our DLAs (the integral of the star formation rate). This quantity can be estimated in surveys of LBG by fitting their spectral energy distribution deduced from multi-band photometry. In Figure 15 we have plotted the distribution of our DLA stellar masses and compared it with the estimated range of observed LBG’s stellar masses, using the results of Shapley et al. (2001) adapted as described above for our Kroupa IMF. Most DLAs have formed a relatively small population of stars (the distribution ranges over but peaks strongly at ), consistent with their low metallicities. There is a strong dependency on the underlying halo mass, which we have shown by plotting the contribution from different halo mass ranges. This arises not only because the baryonic mass rises roughly linearly with the virial mass, but also because the star formation in low mass haloes is substantially suppressed by feedback (see also Sections 4.2.3 and 5.1).

The simulated DLAs are clearly being consumed very slowly, and it will be a matter of considerable interest to understand the ultimate fate of the gas, especially given the apparent realism of the galaxies at . We intend to address this question fully in future work; for the present, we note that over 80% of the neutral  H i in the major progenitor at has formed stars in our Milky Way type galaxy (box MW) by . However, this accounts for only 10% of the total stars; 10% of the stars were in fact already formed at , with the remainder coming from accreted gas and stars in satellites.

5 Consistency Checks and Effects of Parameters

Brief Tag Comment
MW As Table 2
  .LR Lower Resolution
  .LR.ThF Thermal Feedback
  .LR.SS Self-Shielding
Large As Table 2
  .LR Lower Resolution
  .LR.ThF Thermal Feedback
  .SS Self-Shielding
Cosmo As Table 2
  .SS Self-Shielding
Table 3: A summary of the comparison simulations used.

In this section, we will discuss the effect of changing some details of our simulations and processing, including the resolution and star formation feedback prescription parameters. Since the effects are rather minor overall, some readers may prefer to skip directly to our conclusions (Section 6). The extra simulations used for the discussions below are summarised in Table 3.

5.1 Resolution and Star Formation Feedback

Figure 16: As Figure 4, except we now plot only two simulations (MW, Large; dots and crosses respectively) and variants thereon – see Table 3 for a description of these.

In this section, we describe a further suite of simulations which probe the dependence of our overall conclusions on resolution and feedback strength. Previous studies of DLAs have been shown to be somewhat sensitive to the resolution (e.g. Nagamine et al., 2004a; Razoumov et al., 2006) and the main difference between our simulations and previous work lies in the feedback implementation, so that both these considerations are worth exploring.

We described in Section 3.2 our resolution criteria, demanding both that and . This is quite a conservative limit, based on rerunning two of our simulations (“MW” and “Large”) at lower resolution (for details see Table 3). For haloes with coarser resolution than this stated limit, we found that the DLAs started to exhibit scatter about the locus defined by their higher resolution counterpart, and a slight tendency to underestimate the DLA cross-section (while overestimating individual column densities). These resolution effects are relatively mild compared to some previously reported effects (e.g. Gardner et al. (1997a), Nagamine et al. (2004a), Razoumov et al. (2006)). We suggest that this is because our treatment of feedback (see Section 2), which leads to efficient energy deposition from supernovae in dense areas, helps to impose an effective scale-length below which the particles composing the ISM do not collapse further. Ultimately, given the impossible task of maintaining sufficient dynamic resolution to track the actual gas components making up the ISM, this is not an unreasonable behaviour. Figure 16 shows the “MW” and “Large” simulations along with their low resolution counterparts (“MW.LR”, “Large.LR”) and demonstrates that the simulations are in good agreement when restricted according to our resolution criteria. The velocity widths from these low resolution simulations were also overall in agreement with their high resolution counterparts; however, as previously noted, the star formation histories and metallicities only converge for a more stringent resolution cut (Section 4.2.3).

We took the low resolution initial conditions and reran the simulations, replacing our normal feedback implementation with a purely thermal approach (i.e. the same energy per supernova was produced, but the cooling switch-off was not implemented). This has the well-known effect of causing the supernovae energy to be radiated away over much less than the dynamical timescale (Katz, 1992; Thacker & Couchman, 2000) and is thus the weakest conceivable mechanism. It is almost certain to be unphysical, since the fast cooling times arise from the combination of high temperature and density. This is caused by averaging properties of unresolved regions; in fact the high temperature (blast interior) and high density (undisturbed exterior) regions are presumably separated in the true parsec-scale ISM.

Our cross-sections arising from this approach lie considerably closer to those of Nagamine et al. (2004a) (Figure 16). We also find more usual results for quantities such as the metallicity; because our metallicities converge only at high resolutions (see Section 4.2.3) quantitative comparisons between our low resolution runs are a little dangerous, but we see a very much weaker mass-metallicity relation (see also Brooks et al., 2007), with metallicities all approximately . This is, as expected, in poor agreement with observations but in better agreement with older simulations. Recall that our fiducial feedback formulation causes a reduction in metallicities primarily by suppressing star formation efficiently, rather than causing any bulk outflows (Brooks et al., 2007).

We also note that the velocity widths arising from particular haloes in our thermal feedback simulations are comparable to, or somewhat () higher than, the velocity widths in our fiducial simulations. (The velocity widths agree between low resolution and standard runs.) This does not result in a closer matching of the cosmological velocity width distribution, because the proportion of intermediate and high mass haloes is reduced. We interpret this as suggesting that the origin of our velocity widths is chiefly gravitational, and that the thermal feedback runs form denser, compact  H i regions. This is no surprise: there is no physics in our simulations that can plausibly give rise to large-scale bulk galactic winds. Future improvement in DLA velocity width modeling may therefore arise from an understanding of what drives such flows, and whether cold gas can indeed exist within them. However this understanding will augment, not replace, a sufficient feedback model for the local suppression of star formation; see Section 6.3.

5.2 Si II and low-ion regions

Figure 17: The probability distribution of the ratio of the velocity width from the Cloudy-based Si ii  ionisation runs to the original velocity widths which assume . Although higher velocity widths do arise by considering these effects, the mean differences are not large enough to have a significant impact on, for example, the overall distribution of observed DLA line widths (Figure 9).

One possible explanation for the underestimated incidence of high velocity width absorbers in our simulations (Section 4.2.2) is a misidentification of the precise regions responsible for producing the low-ion profiles. Therefore, we briefly investigated the effect of relaxing our assumptions (Section 3.1) that . This assumption is based on comparing the energetics of the ion transitions; however there is no actual physical mechanism coupling the Si ii stage to  H i. This is to be contrasted with Oi, which is strongly coupled to  H i via a charge transfer mechanism (e.g. Osterbrock & Ferland, 2006). There are clear examples of DLAs in which the Oi width is significantly smaller than the Si ii width (see, for example, the lowest two panels in Figure 6 of Ledoux et al., 1998) and thus there must be contributions to the Si ii velocity structure from outside the cold  H i gas.

Because DLAs are low metallicity environments (both observationally and in our simulations), a first approximation to the silicon ionisation problem can be achieved by taking the radiation intensity from our hydrogen/helium radiative transfer code (Section 3.1) and calculating the equilibrium state of Silicon. To address this in a simple way, we calculated a grid of Cloudy101010Version 07.02.01, available from www.nublado.org and last described by Ferland et al. (1998) models which varied in the local density, temperature and incident radiation strength indexed by the single photoionisation parameter . For each particle in our simulation, we calculated the value of Siii/Si by interpolating these models. Then we recalculated our DLA sightlines (Section 3.2) using the new Siii values, rather than the old assumed values. We used exactly the same set of offsets and angles, so that the final reprocessed-catalogue can be compared on a per-sightline basis. Our sightlines extend through the entire box so that any trace gas outside our haloes could theoretically contribute to the velocity width. However, the combination of the decreasing neutral gas density and decreasing metallicities makes the intergalactic contribution to the unsaturated line widths extremely minor.

We have plotted, for all sightlines, the probability distribution of the ratio of the updated velocity widths to the original in Figure 17. A number of sightlines did, after reprocessing, display significantly higher velocity widths (as much as doubling in some cases). Thus for individual systems, these corrections can be important. However (noting that the -axis in the plot is logarithmic) it is clear that the significantly increased velocity width systems are too rare to make a significant impact on the distribution of velocity widths (Figure 9), a fact we verified explicitly.

5.3 Cosmological Model

The parameters used throughout this work were fixed when running the first of the simulations described. Since then, our knowledge of cosmological parameters has improved thanks to the expanding wealth of constraints from galaxy surveys measuring the baryon acoustic oscillations (BAO), supernovae surveys and most significantly three- and five-year WMAP results (Spergel et al., 2006; Dunkley et al., 2008). Our chosen parameters are , whereas the latest WMAP, BAO and SN results combined suggest . It is natural to question whether such differences can have a significant effect on our results.

Given our confined cross-section (Section 4.1), any effect can be split into two components: the change in the halo mass function, and the change in the individual objects. Considering the halo mass function first, since at the relative density of to CDM is suppressed by a factor of , the difference in has a negligible effect. Further, the value of has only minor consequences for the transfer function (amounting to the form of the acoustic oscillations). The most important difference, therefore, should arise from the reduction of : this will decrease the scale on which fluctuations have become non-linear, and hence reduce the number of high mass objects. However, we also need to consider the lower value of . The high mass end of the halo mass function is less affected (since it is closer to the normalization scale), but at the low mass end (corresponding to high ) there are a lower than expected number of haloes. Overall then, the effect of updating our halo mass function is to slightly reduce the line density of DLAs (to , which is marginally lower than the observed line density, but acceptable given the uncertainties in our simulations) and shift the cross-section to somewhat lower masses. However, we verified that this had negligible effect on, for example, our velocity width distribution. This is because, fortuitously, the two effects produce an almost flat value of over .

On halo scales, the shift in the value of must be the dominant effect. We do not believe this will have a large effect on our simulations, being a small correction, especially since equilibrium considerations determine the physics on these scales. A complete resolution of this question, however, awaits new simulations with the updated values.

5.4 Radiative Transfer

We have incorporated a correction for self-shielding (Section 3.1) which is coarse compared to some recent simulations (Razoumov et al., 2006, 2008). On the other hand, the complexity of the fine-structure of the ISM today (and presumably also at ) means that even the most advanced simulations cannot capture its “microscopic” (i.e. parsec scale) behaviour. Nagamine et al. (2004a) employed a sub-resolution two-phase pressure equilibrium model for the ISM (Springel & Hernquist, 2003; McKee & Ostriker, 1977), apparently obviating the need for radiative transfer since the cold clouds are assumed to be fully self-shielded and the warm ambient medium to be optically thin. This is a promising approach, but it is not entirely clear how it links to the messy observational picture in which atomic gas exists in distinct warm and cold phases, ionized gas in warm and hot phases, and photoionized regions occur within the cold phase around star formation ( H ii regions) as well as in the diffuse disk ISM.

For this work, our essential aim was to identify large scale regions in which the cosmological UV background should be significantly suppressed. We verified our code using a variety of tests comprising comparisons of equilibrium ionisation front positions with simple constant temperature Cloudy models. While our radiative transfer code performs well in simple plane-parallel setups, it resolves only 6 angular elements and its true 3D performance is therefore significantly degraded. However, we believe that for the purposes of the present study it is adequate: in both optically thin and optically thick limits, the angular resolution is of little importance; only in the transition regions will significant differences arise. We ran a final test in which we post-processed our “Large” box after rotating it through 45, thus providing a different set of angular elements to the code. None of our DLA results were significantly affected by this change.

We do not include the UV emission of star forming regions in our calculation. This would introduce a significant extra complexity to our algorithm without any real benefits for the following reason. The star formation rate of most DLAs is known observationally – and in our simulations – to be (see Section 4.3, Figure 14). For such low star formation rates the local contribution to the UVB is comparable to, or usually significantly smaller than, the field. This can be estimated, for example, using the data of Black 1987, or more usually is obtained directly from observational data – see Introduction. Although the theoretical picture of the importance of local sources is not straightforward (e.g. Miralda-Escudé, 2005; Schaye, 2006), we are confident that these direct observational upper limits justify our approximation.

There is another slight inconsistency in our approach, in that the radiative transfer post-processor identifies certain regions of gas as neutral which in the live simulations is treated as ionized; further, UV heating occurs for the uniform background even in regions which are later identified as self-shielding. To assess the magnitude of this effect, we re-simulated “MW.LR”, “Large” and “Cosmo” using a local UVB attenuation approximation (“SS” simulations in Table 3). In these simulations, for each gas particle the strength of the UVB used for ionization equilibrium and heating calculations was reduced by the mean attenuation for particles of that density in the post-processed simulation. We continue to apply our radiative transfer post-processor, but the corrections involved are smaller in the new outputs.

Figure 18: The mean (sight-line weighted) metallicity of DLAs in haloes from the “Cosmo” & “Large” boxes (tripod and cross symbols respectively) contrasted with their live self-shielding counterparts (dots and plus symbols). The most significant effect of including live self-shielding, other than a slight overall increase in DLA cross-sections, is a decrease in star formation rates in low mass haloes, causing them to display systematically lower metallicities. For clarity in this plot, we have only displayed the first 100 haloes from the “Cosmo” boxes.

In general, the differences between the “live” self-shielding and original runs were rather minor. They consisted of a slight increase in the DLA cross-section ( dex) and a reduction in the star formation rates. Consequently the cold gas metallicity of these low mass haloes was seen to drop by up to dex (Figure 18) – this was apparently the biggest change caused. In fact, this could help to resolve the slight paucity of low metallicity ([M/H]) DLAs in our simulations (Figure 11), although without a fuller set of statistics we were not able to verify this directly. Importantly, however, we reproduced all our measured distributions replacing “Cosmo” with “Cosmo.SS” and “Large” with “Large.SS”. The effects of this were dominated by an increase in overall line density to (see Section 4.2.1). This takes us further away from the observed result . However, excepting normalization, the distribution of properties was left almost unchanged. We believe these effects warrant further investigation in the future, but their direct effects seem to be fortuitously small. One should be aware that, as well as the self-shielding effect which reduces the UV heating and hence equilibrium temperature and pressure in the ISM, other inaccuracies in the detail of the ISM (such as magnetic pressure) could easily contribute opposite effects of similar magnitudes; however their study is well beyond the scope of this work.

5.5 Intergalactic DLAs?

It is clear from Figure 2 that our DLA cross-section is confined to within the virial radii of our host haloes. This is in qualitative agreement with the majority of previous DLA simulations listed in Table 1, and is an assumption of all DLA semi-analytic models that we are aware of. However, it contrasts with the recent results of Razoumov et al. (2006) and Razoumov et al. (2008) (R06 and R08 respectively) in which as much as 50% of the DLA cross-section resides outside the virial radii of haloes (see R08 Table 2 and the top right panel of R06 Figure 3).

Although rough estimates suggest that DLAs should form in dark matter haloes, it is not unreasonable for DLAs to form in overdense surrounding regions. Given the requirement for self-shielding, (Haehnelt et al., 1998), one may compare the gravitational mass required to confine such gas () to the total gas mass enclosed in the DLA (), giving the ratio


Thus, for DLAs with column densities close to the lower limit, it may be possible for the self-gravity of gas, or slight dark matter overdensities, to obviate the need for a fully collapsed dark matter halo. This is especially true if (unlike in our simulations) gas is allowed to cool efficiently to temperatures .

The bigger problem is in understanding how the gas cools to become neutral in the first place. One key difference between our simulations and those of R06, R08 is that the latter include an approximate algorithm to correct the temperature for shielding effects, whereas we keep our temperatures constant during the radiative transfer processing (Section 3.1). While our approach is clearly an approximation, the latter may overcool gas since dynamical and feedback heating effects could easily become significant as the UV field drops. The only way to correctly resolve this problem is (as in R06, but not R08) to include radiative transfer in the live simulations.

It will be interesting to see how this issue and a physical understanding of it develop with future generations of simulations and radiative transfer codes. But while it remains to be adequately resolved, we should note that (a) the effect in R06/08 is resolution dependent, and therefore the physical status is not transparent; (b) our self-shielding runs suggested the effect was rather small (Section 5.4), at least with the local approximation and (c) we have produced a DLA cross-section which matches nearly all observational constraints, which lends some indirect reinforcement to the validity of our approximation.

6 Discussion and Conclusions

6.1 Overview

We have investigated the occurrence of DLAs in a series of simulations (see Governato et al. 2007, 2008 and Brooks et al. 2007 – in this text referred to as G07, G08 and B07) which produce galaxies at with as near as currently possible to realistic physical properties (see Introduction). These high resolution simulations include a physically motivated star formation feedback prescription with only one free parameter (the efficiency of energy deposition) which is set by observations of low redshift star formation. Thus, as well as providing information on DLAs themselves, this work is an independent cross-check of the formation process of the G07/B07/G08 galaxies.

To produce cosmological statistics from a smorgasbord of boxes, we used a non-parametric weighting process, in effect correcting the halo mass function of the combined sample. This does not involve any fitting or parametrization, and hence no extrapolation: all the results presented in this paper are derived directly from resolved regions of simulations.

The picture that emerges from our simulations is of a cosmological DLA cross-section predominantly provided by intermediate mass haloes, . These ranges are somewhat higher than many early simulations suggested (e.g. Gardner et al., 1997b, 2001), and are close to the range of suggested by observations of DLA/LBG correlations (Cooke et al., 2006). Our DLAs form stars at low rates (predominantly ) and have metallicities in the range with a median of approximately , both in agreement with observations. We also investigated the distribution of total stellar masses of DLAs, finding them to be predominantly spread over , with a peak at . By the majority of the neutral gas forming the DLAs has been converted into stars, but during this time substantial merging complicates the direct identification of low redshift galaxies with high redshift absorbers. In plots showing overall halo properties we have marked the location of the major progenitors to our Milky Way like and dwarf-type galaxies (from boxes “MW” and “Dwf” respectively); they appear to have been fairly typical DLAs.

6.2 Mass-Metallicity Relationships

A tight relationship between dark matter halo mass and star formation rate (Figure 13) underlies a strong mass-metallicity relation; in Figure 12 we have shown that this is manifested by a correlation between DLA velocity widths and metallicities (for equivalent observational results see Ledoux et al., 2006; Prochaska et al., 2008). Further, the same simulations produce a realistic mass-metallicity relation for galaxies (B07), suggesting that these simulations capture the metal enrichment of collapsed systems in a meaningful way. The DLA result is not significantly affected by gradients within the disks of our forming galaxies, but the DLA sightlines do preferentially sample regions of gas with higher metallicities (by factors ) than the mass-weighted halo mean. Further work on metallicities in these simulations is ongoing, in particular regarding the importance of metal diffusion – the metals are simply advected with the SPH particles, which can lead to inaccuracies where spatial gradients are important. In connection with this work, it will be interesting to see whether the metal enrichment of the simulations’ less overdense regions also appears realistic, providing a connection to the importance (or otherwise) of outflows (Section 6.3).

Recent work has suggested that the origin of the observed relationship between kinematics as measured by Mg ii rest-frame equivalent widths () and metallicity of DLAs should not be relied upon as an indicator of an underlying variation with mass (Bouché, 2008). However, this claim is not directly relevant to fiducial velocity width determinations which are measured from unsaturated or mildly saturated transitions; MgII is a strongly saturated absorption line, and as such can easily be affected by trace gas in the outer regions of galaxies (or perhaps even the ambient IGM).

6.3 Velocity Profiles

Overall, we reproduce the spread of velocity widths seen in DLA systems (Figure 9) and approximately match the distribution of low velocity width systems (). However the long standing difficulty of underestimating the observed incidence rate of high velocity widths () persists in our simulations – although the discrepancy is rather small compared with older simulations and comparable to that seen in the recent work by Razoumov et al. (2008). It is not entirely clear why simulations in general have encountered such persistent difficulty in this respect. There are essentially three possibilities:

Firstly, it is possible that our assumption oversimplifies the identification of regions responsible for the low-ion profiles and leads to systematic underestimates of the velocity dispersion. We investigated this avenue briefly (Section 5.2), but the correction did not appear sufficient to resolve the discrepancy – although this could be sensitive to the description of the ISM on sub-resolution scales (see also below). We ran a brief test in which all gas was assumed to contribute to the velocity width – in this (completely artificial) scenario, the cosmological rate of high velocity widths is over-estimated, showing that the required motions are, at least, “available” in this sense.

Secondly, perhaps the internal gas kinematics still require some correction. Although we did not, for this work, attempt a rigorous decomposition, qualitative inspection suggested that a mixture of disk kinematics (Prochaska & Wolfe, 1998) and merging clumps (Haehnelt et al., 1998) are responsible for the final spread of velocities. We did not find that our feedback made any significant contribution via turbulent motions to the velocity widths (Section 5.1) – hardly a surprise, given that the coupling to the ISM is achieved entirely thermally. A better understanding of feedback is likely to help the kinematical situation by inducing bulk outflows. It would be interesting to see how the simulations of Nagamine et al. (2004a) perform in reproducing DLA line profiles, since they include a phenomenological model of galactic winds. Promising progress in placing such winds on a more physical footing was recently described by Ceverino & Klypin (2007), whose simulations reproduce such bulk motions starting from an explicit model for the ISM.

Thirdly, it is possible that our cross-section is associated with, on average, haloes of insufficient mass; in other words, DLAs at the high mass end could be overly compact. This in turn could be connected to the more well-known angular momentum problem in which simulated disk galaxies have underestimated scale lengths (Navarro & Steinmetz, 2000; Eke et al., 2001). Our simulations, as was noted in the Introduction, go some significant way towards a resolution of the issue for disk galaxies by combining sufficient feedback (Stinson et al., 2006) with high resolution (Kaufmann et al., 2007); our elevated DLA cross-sections for intermediate-mass haloes are connected to this same feedback (Section 5.1). As the mechanisms preventing angular momentum loss are further understood, these same mechanisms may continue to increase the DLA cross-section for high mass haloes. On the other hand, our matching of the metallicity distribution (Figure 11) provides an alternative joint constraint on the star formation history and mass of responsible haloes – and is in good agreement with observations. Any change in the mass of haloes comprising our DLA cross-section would therefore need to be compensated by a change in our virial mass – star formation relation (Figure 13).

6.4 Column Density Distribution and Dust Biasing

As we discussed in Section 4.2.1, we match the overall column density distribution well, but somewhat overestimate the number of observed systems with . This causes us to slightly overpredict value of (we obtain , whereas the SDSS data yield ). The simulated value of is obtained by imposing an upper cut-off on the column density of systems contributing to the measured because stronger systems are too rare to be observed with current observational datasets. But in our simulations, these same systems actually contribute significantly to the cosmic density of DLAs; when included in the census, they boost the total to . Thus in our picture, “invisible” systems contribute about to the  H i budget; but functional fitting of the observed data from SDSS suggests a much steeper cut-off at high column densities than our simulations, and hence that the observed value given above should be essentially correct (Prochaska et al., 2005). While the precise observational extrapolation is subjective and the fit is driven by a couple of bins with , there is a strong constraint in the lack of systems at high column density. For instance, only three systems with are observed in the SDSS DR5 survey, whereas our simulation predicts ten over the same path length (a glance at the Poisson distribution function shows that this is a significant discrepancy).

It is notable that host galaxies of gamma ray bursts (GRBs), which trace a sample weighted towards higher column densities, routinely display column densities of (Jakobsson et al., 2006) – so that such systems certainly exist but have a smaller cross-section than our simulations suggest. It is possible that we overpredict the number of high column density systems because our simulations are inadequate in the responsible regions, which are cold and dense: one should form H but the complex physics responsible is not implemented in our code. Schaye (2001) suggested that this molecular cloud formation could determine a characteristic cut-off for the column density distribution. On the other hand, observations suggest that DLAs harbour very little molecular gas, with fractions at the very most (Ledoux et al., 2003). Of course the typical DLA sightline could miss high density clumps within giant molecular cloud formations, which would mean the global fraction of H could be somewhat higher. It is therefore currently unclear how feasible this explanation could be.

It has been claimed in the past (e.g. Cen et al., 2003) that DLA dust obscuration effects should substantially reduce the number of observed high systems relative to their true abundance. However, this is hard to reconcile with the most recent observations; in particular, radio-selected samples (which are unlikely to be biased by dust considerations) show little if any evidence for increased high absorbers (Ellison et al., 2001; Jorgenson et al., 2006; Ellison et al., 2008). Further, in a companion paper (Pontzen & Pettini, 2008) we show that there is only marginal statistical evidence for dust-induced obscuration in the joint distribution of observed DLA column densities and metallicities, which are in fact nearly uncorrelated. This also arises quite naturally in our simulations, because the dependence of on the underlying halo mass is weak (Figure 8), whereas metallicities and kinematics are strongly correlated with the halo mass. Combining all the above considerations with our successful reproduction of the observed distribution of metallicities (Section 4.2.3), we believe dust-biasing is unlikely to have a major part to play in a future understanding of DLAs.

The very slight trend which is observed linking the column density to the halo mass in our simulations biases low column density observations in favour of finding low mass haloes (Figure 8). We have not explicitly studied the sub-DLA population, so that these results are not directly comparable with the work of Khare et al. (2007), in which it is suggested that sub-DLAs preferentially probe more massive haloes than DLAs. However, there is certainly a tension between these results which deserves some attention in the future.

6.5 Bimodality

Although we did not specifically set out to investigate the possibility of bimodality within the DLA population (Wolfe et al., 2008), we saw a bifurcation in our cross sections as a function of  H i mass (Figure 7). The bifurcating lower branch in Figure 7 has higher virial masses for a given  H i mass, wider velocity profiles, higher star formation rates and higher mean temperatures. It is important to emphasize that these systems, which are perhaps undergoing a starbursting phase, account for only of the DLA cross-section; our results are therefore not directly comparable to the observational claim of two populations of roughly equal importance.

We ran some careful tests to try to understand the origin of this effect, but it essentially remains work in progress: further simulations and detailed studies of the evolution of these objects will be key to a full account. The critical halo mass, , is probably too small to be directly related to the shock-stability model of Dekel & Birnboim (2006). There are no spatial correlations which would suggest the abnormal haloes are associated with protocluster regions. We verified that the results from the “Large” box were consistent with being drawn from a subvolume of the “Cosmo” box. Since the low cross-section branch objects account for such a small fraction of our DLAs, we felt justified in postponing their full study to later work.

6.6 Future Work

In general, we have commented throughout that the lack of resolution on fine ISM scales (a common feature of all current cosmological simulations, likely to persist into the foreseeable future) is a fundamental limitation. The ISM is a complex mixture of gas in different phases with variations on truly tiny scales. It is also interesting to note that observations suggest that the pressure contributed by magnetic effects in the disk of our galaxy and other local galaxies is comparable in magnitude to turbulent kinetic pressure (e.g. Heiles & Crutcher, 2005); hence magnetic effects could contribute at least as much as stellar feedback to an understanding of the formation of disk galaxies, at least at . It is not impossible to fathom that magnetic fields could have an important role to play in a better understanding of DLAs. But regardless of any particular physical processes, our feeling is that long term future focus should remain on improving our understanding of the coarse-grained behaviour of such a mixture, i.e. in developing sub-grid physics models. This conclusion was reached independently by Razoumov et al. (2008).

Even with our current set of simulations, it will be interesting to investigate in more detail the predictions of these simulations with regard to the LBG population; in Section 4.3 we showed that the total star formation rate at may be very slightly overpredicted compared to observation. But, in work not presented here, we have also investigated the shape of the luminosity function (at and higher redshifts) and found both and the faint end slope to be consistent with observations. Remaining shortcomings could easily be accounted for by details of the IMF, suggesting that the underlying populations are quite reasonable (Governato et al, in prep).

In terms of the properties of DLAs, this work opens the door to an array of further studies. In particular, we intend to address in more detail the time evolution of our DLA population and the transfer of gas through hot, warm and stellar phases in the near future. Given the realistic nature of metals in our DLAs and in our galaxies at (Brooks et al., 2007), the details of the intervening time could, amongst other things, shed light on the “missing metals” problem (Pettini, 2006). Further, it is known observationally that the metallicities, velocity widths and column density distributions evolve only slightly with decreasing redshift; will our simulations reproduce such slow evolution, which depends on finely balancing the cooling of gas, forming  H i, with the formation of stars, destroying it? If so, the manner in which this is accomplished will be of considerable interest. (For instance, do individual objects maintain the same DLA cross-section, or is the slow evolution only seen over cosmological averages?) Although flawed in many respects, simulations offer a guide to our understanding of these issues which should not be ignored – as long as it is taken with a healthy pinch of salt.


We thank Jason X. Prochaska for providing observational data, Francesco Haardt for UVB models and the referee, Peter Johansson, for a thorough and helpful report. AP is supported by a STFC (formerly PPARC) studentship and scholarship at St John’s College, Cambridge and gratefully acknowledges helpful conversations and communications with George Efstathiou, Gary Ferland, Johan Fynbo, Alexei Razoumov, Emma Ryan-Weber and Art Wolfe. FG acknowledges support from a Theodore Dunham grant, HST GO-1125, NSF ITR grant PHY-0205413 (supporting TQ), NSF grant AST-0607819 and NASA ATP NNX08AG84G. Most simulations were run at the Arctic Region Supercomputing Center. Part of the analysis was performed on Cosmos, a UK-CCC facility supported by HEFCE and STFC in cooperation with SGI/Intel utilising the Altix 3700 supercomputer.


  • Black (1987) Black J. H., 1987, in Astrophysics and Space Science Library, Vol. 134, Interstellar Processes, Hollenbach D. J., Thronson Jr. H. A., eds., pp. 731–744
  • Bouché (2008) Bouché N., 2008, pre-print (arXiv:0803.3944)
  • Brooks et al. (2007) Brooks A. M., Governato F., Booth C. M., Willman B., Gardner J. P., Wadsley J., Stinson G., Quinn T., 2007, ApJ, 655, L17 (B07)
  • Cen et al. (2003) Cen R., Ostriker J. P., Prochaska J. X., Wolfe A. M., 2003, ApJ, 598, 741
  • Ceverino & Klypin (2007) Ceverino D., Klypin A., 2007, pre-print (arXiv:0712.3285)
  • Chevalier (1974) Chevalier R. A., 1974, ApJ, 188, 501
  • Cooke et al. (2006) Cooke J., Wolfe A. M., Gawiser E., Prochaska J. X., 2006, ApJ, 652, 994
  • Dekel & Birnboim (2006) Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
  • Dunkley et al. (2008) Dunkley J. et al. 2008, ApJS, submitted (arXiv:0803.0586)
  • Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
  • Eke et al. (2001) Eke V. R., Navarro J. F., Steinmetz M., 2001, ApJ, 554, 114
  • Ellison et al. (2008) Ellison S., York B., Pettini M., Nissim K., 2008, MNRAS, Submitted
  • Ellison et al. (2001) Ellison S. L., Yan L., Hook I. M., Pettini M., Wall J. V., Shaver P., 2001, A&A, 379, 393
  • Erb et al. (2006) Erb D. K., Shapley A. E., Pettini M., Steidel C. C., Reddy N. A., Adelberger K. L., 2006, ApJ, 644, 813
  • Ferland et al. (1998) Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761
  • Fox et al. (2007) Fox A. J., Petitjean P., Ledoux C., Srianand R., 2007, A&A, 465, 171
  • Gardner et al. (1997a) Gardner J. P., Katz N., Hernquist L., Weinberg D. H., 1997a, ApJ, 484, 31
  • Gardner et al. (2001) —, 2001, ApJ, 559, 131
  • Gardner et al. (1997b) Gardner J. P., Katz N., Weinberg D. H., Hernquist L., 1997b, ApJ, 486, 42
  • Gill et al. (2004) Gill S. P. D., Knebe A., Gibson B. K., 2004, MNRAS, 351, 399
  • Governato et al. (2008) Governato F., Mayer L., Brook C., 2008, preprint (arXiv:0801.1707) (G08)
  • Governato et al. (2007) Governato F., Willman B., Mayer L., Brooks A., Stinson G., Valenzuela O., Wadsley J., Quinn T., 2007, MNRAS, 374, 1479 (G07)
  • Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
  • Haehnelt et al. (1998) Haehnelt M. G., Steinmetz M., Rauch M., 1998, ApJ, 495, 647
  • Heiles & Crutcher (2005) Heiles C., Crutcher R., 2005, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 664, Cosmic Magnetic Fields, Wielebinski R., Beck R., eds., pp. 137–+
  • Iliev et al. (2006) Iliev I. T., Ciardi B., Alvarez M. A., Maselli A., Ferrara A., Gnedin N. Y., Mellema G., Nakamoto T., Norman M. L., Razoumov A. O., Rijkhorst E.-J., Ritzerveld J., Shapiro P. R., Susa H., Umemura M., Whalen D. J., 2006, MNRAS, 371, 1057
  • Isobe et al. (1990) Isobe T., Feigelson E. D., Akritas M. G., Babu G. J., 1990, ApJ, 364, 104
  • Jakobsson et al. (2006) Jakobsson P., Fynbo J. P. U., Ledoux C., Vreeswijk, P., Kann, D. A., Hjorth, J., Priddey, R. S., Tanvir, N. R., Reichart, D., Gorosabel, J., Klose, S., Watson, D., Sollerman, J., Fruchter, A. S., de Ugarte Postigo, A., Wiersema, K., Björnsson, G., Chapman, R., Thöne, C. C., Pedersen, K., Jensen, B. L., 2006, A&A, 460, L13
  • Johansson & Efstathiou (2006) Johansson P. H., Efstathiou G., 2006, MNRAS, 371, 1519
  • Jorgenson et al. (2006) Jorgenson R. A., Wolfe A. M., Prochaska J. X., Lu L., Howk J. C., Cooke J., Gawiser E., Gelino D. M., 2006, ApJ, 646, 730
  • Katz (1992) Katz N., 1992, ApJ, 391, 502
  • Katz et al. (1996a) Katz N., Weinberg D. H., Hernquist L., 1996a, ApJS, 105, 19
  • Katz et al. (1996b) Katz N., Weinberg D. H., Hernquist L., Miralda-Escude J., 1996b, ApJ, 457, L57+
  • Katz & White (1993) Katz N., White S. D. M., 1993, ApJ, 412, 455
  • Kauffmann (1996) Kauffmann G., 1996, MNRAS, 281, 475
  • Kaufmann et al. (2007) Kaufmann T., Mayer L., Wadsley J., Stadel J., Moore B., 2007, MNRAS, 375, 53
  • Kennicutt (1998a) Kennicutt Jr. R. C., 1998a, ARA&A, 36, 189
  • Kennicutt (1998b) —, 1998b, ApJ, 498, 541
  • Khare et al. (2007) Khare P., Kulkarni V. P., Péroux C., York D. G., Lauroesch J. T., Meiring J. D., 2007, A&A, 464, 487
  • Knebe et al. (2001) Knebe A., Green A., Binney J., 2001, MNRAS, 325, 845
  • Kroupa et al. (1993) Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
  • Ledoux et al. (2002) Ledoux C., Bergeron J., Petitjean P., 2002, A&A, 385, 802
  • Ledoux et al. (1998) Ledoux C., Petitjean P., Bergeron J., Wampler E. J., Srianand R., 1998, A&A, 337, 51
  • Ledoux et al. (2006) Ledoux C., Petitjean P., Fynbo J. P. U., Møller P., Srianand R., 2006, A&A, 457, 71
  • Ledoux et al. (2003) Ledoux C., Petitjean P., Srianand R., 2003, MNRAS, 346, 209
  • Lee et al. (2006) Lee H., Skillman E. D., Cannon J. M., Jackson D. C., Gehrz R. D., Polomski E. F., Woodward C. E., 2006, ApJ, 647, 970
  • Lodders (2003) Lodders K., 2003, ApJ, 591, 1220
  • Maller et al. (2001) Maller A. H., Prochaska J. X., Somerville R. S., Primack J. R., 2001, MNRAS, 326, 1475
  • McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  • Miralda-Escudé (2005) Miralda-Escudé J., 2005, ApJ, 620, L91
  • Naab et al. (2007) Naab T., Johansson P. H., Ostriker J. P., Efstathiou G., 2007, ApJ, 658, 710
  • Nagamine et al. (2004a) Nagamine K., Springel V., Hernquist L., 2004a, MNRAS, 348, 421
  • Nagamine et al. (2004b) —, 2004b, MNRAS, 348, 435
  • Nagamine et al. (2007) Nagamine K., Wolfe A. M., Hernquist L., Springel V., 2007, ApJ, 660, 945
  • Navarro & Steinmetz (2000) Navarro J. F., Steinmetz M., 2000, ApJ, 538, 477
  • Osterbrock & Ferland (2006) Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei, 2nd. ed.  Sausalito, CA: University Science Books, 2006
  • Ostriker & McKee (1988) Ostriker J. P., McKee C. F., 1988, Reviews of Modern Physics, 60, 1
  • Pettini (2006) Pettini M., 2006, in The Fabulous Destiny of Galaxies: Bridging Past and Present, Paris:Frontier Group, 2006
  • Pontzen & Pettini (2008) Pontzen A., Pettini M., 2008, MNRAS, submitted
  • Prochaska et al. (2008) Prochaska J. X., Chen H.-W., Wolfe A. M., Dessauges-Zavadsky M., Bloom J. S., 2008, ApJ, 672, 59
  • Prochaska et al. (2003) Prochaska J. X., Gawiser E., Wolfe A. M., Castro S., Djorgovski S. G., 2003, ApJ, 595, L9
  • Prochaska et al. (2005) Prochaska J. X., Herbert-Fort S., Wolfe A. M., 2005, ApJ, 635, 123
  • Prochaska & Wolfe (1997) Prochaska J. X., Wolfe A. M., 1997, ApJ, 487, 73
  • Prochaska & Wolfe (1998) —, 1998, ApJ, 507, 113
  • Prochaska & Wolfe (2001) —, 2001, ApJ, 560, L33
  • Prochaska et al. (2007) Prochaska J. X., Wolfe A. M., Howk J. C., Gawiser E., Burles S. M., Cooke J., 2007, ApJS, 171, 29
  • Quinn et al. (1996) Quinn T., Katz N., Efstathiou G., 1996, MNRAS, 278, L49
  • Raiteri et al. (1996) Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105
  • Razoumov & Cardall (2005) Razoumov A. O., Cardall C. Y., 2005, MNRAS, 362, 1413
  • Razoumov et al. (2008) Razoumov A. O., Norman M. L., Prochaska J. X., Sommer-Larsen J., Wolfe A. M., Yang Y.-J., 2008, ApJ, 683, 149 (R08)
  • Razoumov et al. (2006) Razoumov A. O., Norman M. L., Prochaska J. X., Wolfe A. M., 2006, ApJ, 645, 55
  • Reddy et al. (2007) Reddy N. A., Steidel C. C., Pettini M., Adelberger K. L., Shapley A. E., Erb D. K., Dickinson M., preprint (arXiv:0706.4091)
  • Reed et al. (2006) Reed D., Bower R., Frenk C., Jenkins A., Theuns T., 2007, MNRAS, 374, 2
  • Rees (1986) Rees M. J., 1986, MNRAS, 218, 25P
  • Savaglio et al. (2005) Savaglio S., Glazebrook K., Le Borgne D., Juneau S., Abraham R. G., Chen H.-W., Crampton D., McCarthy P. J., Carlberg R. G., Marzke R. O., Roth K., Jørgensen I., Murowinski R., 2005, ApJ, 635, 260
  • Schaye (2001) Schaye J., 2001, ApJ, 562, L95
  • Schaye (2006) —, 2006, ApJ, 643, 59
  • Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  • Shapley et al. (2001) Shapley A. E., Steidel C. C., Adelberger K. L., Dickinson M., Giavalisco M., Pettini M., 2001, ApJ, 562, 95
  • Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
  • Spergel et al. (2006) Spergel D. N., Bean R., Doré O., Nolta M. R., Bennett C. L., Dunkley J., Hinshaw G., Jarosik N., Komatsu E., Page L., Peiris H. V., Verde L., Halpern M., Hill R. S., Kogut A., Limon M., Meyer S. S., Odegard N., Tucker G. S., Weiland J. L., Wollack E., Wright E. L., ApJS, 170, 377
  • Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  • Springel et al. (2005) Springel V., White S. D. M., Jenkins A., Frenk C. S., Yoshida N., Gao L., Navarro J., Thacker R., Croton D., Helly J., Peacock J. A., Cole S., Thomas P., Couchman H., Evrard A., Colberg J., Pearce F., 2005, Nature, 435, 629
  • Srianand et al. (2005) Srianand R., Petitjean P., Ledoux C., Ferland G., Shaw G., 2005, MNRAS, 362, 549
  • Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074 (S06)
  • Thacker & Couchman (2000) Thacker R. J., Couchman H. M. P., 2000, ApJ, 545, 728
  • Thielemann et al. (1986) Thielemann F.-K., Nomoto K., Yokoi K., 1986, A&A, 158, 17
  • Tremonti et al. (2004) 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
  • Tytler (1987) Tytler D., 1987, ApJ, 321, 49
  • Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 9, 137
  • Weaver & Woosley (1993) Weaver T. A., Woosley S. E., 1993, Phys. Rep., 227, 65
  • Wolfe et al. (2005) Wolfe A. M., Gawiser E., Prochaska J. X., 2005, ARA&A, 43, 861
  • Wolfe et al. (2003) Wolfe A. M., Prochaska J. X., Gawiser E., 2003, ApJ, 593, 215
  • Wolfe et al. (2008) Wolfe A. M., Prochaska J. X., Jorgenson R. A., Rafelski M., 2008, preprint (arXiv:0802.3914)
  • Wolfe et al. (1986) Wolfe A. M., Turnshek D. A., Smith H. E., Cohen R. D., 1986, ApJS, 61, 249
  • Zwaan et al. (2005) Zwaan M. A., van der Hulst J. M., Briggs F. H., Verheijen M. A. W., Ryan-Weber E. V., 2005, MNRAS, 364, 1467

Appendix A SPH Calculations

For convenience we will denote the SPH averaging procedure by square brackets, i.e. for any quantity