The nature of damped HI absorbers probed by cosmological simulations: satellite accretion and outflows

TL;DR

This study uses advanced cosmological simulations to analyze the distribution and origins of high column density neutral hydrogen absorbers around galaxies, revealing the roles of satellites, outflows, and galaxy components across redshifts.

Contribution

It demonstrates that high-resolution simulations can accurately reproduce observed absorber distributions and clarifies the contributions of galaxy discs and halos to damped extsc{H}{i} absorption.

Findings

Simulations match observed absorber distributions across redshifts.

DLA and sub-DLA absorbers originate from both galaxy discs and halos.

High-resolution is essential to link absorbers at large distances to dwarf satellites.

Abstract

We use state-of-the-art cosmological zoom simulations to explore the distribution of neutral gas in and around galaxies that gives rise to high column density \ion{H}{i} \mbox{Ly-$\alpha$} absorption (formally, sub-DLAs and DLAs) in the spectra of background quasars. Previous cosmological hydrodynamic simulations under-predict the mean projected separations $(b)$ of these absorbers relative to the host, and invoke selection effects to bridge the gap with observations. On the other hand, single lines of sight (LOS) in absorption cannot uniquely constrain the galactic origin. Our simulations match all observational data, with DLA and sub-DLA LOS existing over the entire probed parameter space ($-4\lesssim $[M/H]$\lesssim 0.5$, $b<50$ kpc) at all redshifts ($z\sim 0.4 - 3.0$). We demonstrate how the existence of DLA LOS at $b\gtrsim 20-30$ kpc from a massive host galaxy require high numerical resolution, and that these LOS are associated with dwarf satellites in the main halo, stripped metal-rich gas and outflows. Separating the galaxy into interstellar ("\ion{H}{i} disc") and circumgalactic ("halo") components, we find that both components significantly contribute to damped \ion{H}{i} absorption LOS. Above the sub-DLA (DLA) limits, the disc and halo contribute with $\sim 60 (80)$ and $\sim 40 (20)$ per cent, respectively. Our simulations confirm analytical model-predictions of the DLA-distribution at $z\lesssim 1$. At high redshift ($z\sim 2-3$) sub-DLA and DLAs occupy similar spatial scales, but on average separate by a factor of two by $z\sim 0.5$. On whether sub-DLA and DLA LOS sample different stellar-mass galaxies, such a correlation can be driven by a differential covering-fraction of sub-DLA to DLA LOS with stellar mass. This preferentially selects sub-DLA LOS in more massive galaxies in the low-$z$ universe.

Tables

Table 1. Table 1: Simulation-based detection fractions of DLA-LOS enclosed within the F08+K17 model 1 ​ σ 1 𝜎 1\sigma -, 2 ​ σ 2 𝜎 2\sigma -, and 3 ​ σ 3 𝜎 3\sigma -confidence regions.

Redshift range	$1\sigma$	$2\sigma$	$3\sigma$
	$\left(0.68\right)$	$\left(0.95\right)$	$\left(0.99\right)$
$3.0-2.0$	$0.41$	$0.74$	$0.77$
$2.0-1.4$	$0.41$	$0.81$	$0.85$
$1.4-0.8$	$0.58$	$0.89$	$0.91$
$0.8-0.4$	$0.78$	$0.97$	$0.98$

Table 2. Table 2: Same as Table 1 , but including sub-DLA LOS.

Redshift range	$1\sigma$	$2\sigma$	$3\sigma$
	$\left(0.68\right)$	$\left(0.95\right)$	$\left(0.99\right)$
$3.0-2.0$	$0.41$	$0.70$	$0.74$
$2.0-1.4$	$0.39$	$0.76$	$0.81$
$1.4-0.8$	$0.42$	$0.77$	$0.81$
$0.8-0.4$	$0.54$	$0.84$	$0.87$

Equations

$${\dot{\rho }}_{\star }={ϵ}_{\mathrm{ff}}\frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}},$$

$$M_Z=2.09M_{\mathrm{O}}+1.06M_{\mathrm{Fe}}$$

$$N_{\text{textsc}Hi}=\sum _i{n}_{\text{textsc}Hi,i}\Delta x_i$$

$$[M/H]=\frac{{\sum }_i{n}_{\text{textsc}Hi,i}{[M/H]}_i}{{\sum }_i{n}_{\text{textsc}Hi,i}},$$

Full text

††pubyear: 2019††pagerange: The nature of damped HI absorbers probed by cosmological simulations: satellite accretion and outflows–The nature of damped HI absorbers probed by cosmological simulations: satellite accretion and outflows

1 Introduction

The chance alignment of neutral gas clouds intervening the line of sight (LOS) towards background quasi-stellar objects (quasars) imprints characteristic absorption features in the quasar power-law continua. In neutral hydrogen (Hi), the most Hi-rich absorbers are the Damped Lyman-$\alpha$ Absorbers (DLAs; $\log_{10}[\mathrm{N}_{\mathrm{\textsc{Hi}}}~{}(\mathrm{cm}^{-2})]\geq 20.3$, Wolfe et al. (1986)) and the sub-DLAs ($19.0\leq\log_{10}[\mathrm{N}_{\mathrm{\textsc{Hi}}}~{}(\mathrm{cm}^{-2})]<20.3$, e.g., Péroux et al. (2003); Zafar et al. (2013)). Unless otherwise specified, we refer to sub-DLAs and DLAs uniformly as damped H i absorbers. Both classes of absorbers are always accompanied by low-ionization metal line complexes (Prochaska et al., 2003; Noterdaeme et al., 2012; Rafelski et al., 2014) which suggests an association with a gaseous medium affected by chemical enrichment. Yet, the origin and nature of damped H i absorbers and their relation to galactic environments remains highly debated.

Absorption velocities and asymmetric line-profiles with leading edges (Prochaska & Wolfe, 1997), as well as models of disc formation (Mo et al., 1998), support the idea that DLAs preferentially probe rapidly rotating galaxy discs (Wolfe et al., 1986). However, using cosmological hydrodynamic simulations, Haehnelt et al. (1998) demonstrated that the absorption kinematics could equally well be explained by irregular proto-galactic clumps in dark matter halos subject to a combination of rotation, random motions, gas flows, and mergers. A picture is emerging in which damped H i absorbers trace neutral gas on scales of tens of kpc in galaxies (Møller & Warren, 1998; Christensen et al., 2007; Pontzen et al., 2008; Monier et al., 2009; Fynbo et al., 2010, 2011; Meiring et al., 2011; Rao et al., 2011; Krogager et al., 2012; Péroux et al., 2012; Fynbo et al., 2013b; Krogager et al., 2013; Rahmati & Schaye, 2014; Rahmani et al., 2016; Rhodin et al., 2018).

A number of analytical models have been successful in reproducing observed bulk properties of DLAs (Fynbo et al., 2008; Krogager et al., 2017, from hereon the “F08+K17 model") (also Freudling et al., 2019, in prep.). By construction, these models depend on mean scaling relations of galaxies, and cannot capture process related to galactic sub-structure or gas structures in the intergalactic medium (IGM). Large-scale cosmological hydrodynamic simulations have the potential to mitigate these shortcomings. However, such simulations struggle with reproducing the observed absorption kinematics (Bird et al., 2014; Bird et al., 2015); lack resolution to capture cold gas in halos (McCourt et al., 2018; Hummels et al., 2018); and require very strong stellar feedback in order to produce cold gas at large impact parameters, destroying the star forming gas discs in the process (Liang et al., 2015). Furthermore, work using cosmological simulations have called into question the detection of high impact parameter $(b~{}\mathrm{[{\,\rm kpc}]}$) damped H i absorbers, suggesting instead that observations suffer from a selection bias towards the most luminous (but unrelated) galaxy in the projected sky vicinity (Rahmati & Schaye, 2014). We are thus left wondering what LOS through galactic environments that give rise to damped H i absorbers.

In this work, we employ high resolution cosmological zoom simulations to build probability-functions of the distribution of damped H i absorption in and around galaxies. These probability-functions are matched to the most recent compilation of spectroscopically confirmed absorber-galaxy pairs (Møller & Christensen, 2019; Christensen et al., 2019, from hereon, the “MC19" compilation) and the analytical F08+K17 model, and are used to explore the physical origin of the damped H i absorption. Throughout the analysis, we confine the comparison to a circular beam of $50~{}\mathrm{kpc}$ radius, centered on the galaxy. This value is observationally motivated, and selected to enclose the current compilation of detections. The paper is organized as follows: Section 2 describes our simulations, with particular emphasis given to the feedback recipes and resolution, Section 3 presents our results, and Section 4 summarizes our conclusions.

2 Simulations

2.1 Simulation setup

We carry out a cosmological hydrodynamic+$N$-body zoom-in simulation of a Milky Way mass galaxy using the adaptive mesh refinement (AMR) code RAMSES (Teyssier, 2002), assuming a flat $\Lambda$-cold dark matter cosmology with $H_{0}=70.2{\,\rm{km\,s^{-1}}}{\,\rm Mpc}^{-1}$, $\Omega_{\rm m}=0.272$, $\Omega_{\Lambda}=0.728$, and $\Omega_{\rm b}=0.045$. From a dark matter only simulation, with a simulation cube of size $L_{\rm box}=85{\,\rm Mpc}$ at $z=0$, a halo of $R_{200,{\rm m}}=334{\,\rm kpc}$ (radius within which the mass density is 200 times the mean matter density) and $M_{200,{\rm m}}=1.3\times 10^{12}{\,\rm M}_{\odot}$ was selected for re-simulation at high resolution. Particles within $3R_{200,{\rm m}}$ at $z=0$ were traced back to $z=100$, and the Lagrangian region they defined was regenerated at high resolution, still embedded within the full lower-resolution volume, using the MUSIC code (Hahn & Abel, 2011). The simulation was then run to $z=0$, with outputs every $\Delta a=0.01$. These are the same initial conditions as the “m12i” halo from Hopkins et al. (2014) and Wetzel et al. (2016), drawn from the volume used in the AGORA galaxy formation comparison project (Kim et al., 2014; Kim et al., 2016).

The dark matter particle mass in the high resolution region is $m_{\rm dm}=3.5\times 10^{4}{\,\rm M}_{\odot}$, with a gas mass resolution of $7070{\,\rm M}_{\odot}$. Star formation is sampled using $10^{4}{\,\rm M}_{\odot}$ particles, with stellar evolution reducing this by up to $50\%$, see below. The adaptive mesh is allowed to refine if a cell contains more than eight dark matter particles. This allows the local force softening to closely match the local mean inter-particle separation, which suppresses discreteness effects (e.g., Romeo et al., 2008). A similar criterion is employed for the baryonic component, where the maximum refinement level is set to allow for a mean constant physical resolution of $\sim 20{\,\rm pc}$ in the dense interstellar medium. As such, this simulation is a significant improvement over previous cosmological zoom simulation used to study cold gas in and around galaxies (e.g. Liang et al., 2015).

The adopted star formation and feedback physics is presented in Agertz et al. (2013) and Agertz & Kravtsov (2015, 2016). Briefly, star formation is treated as a Poisson process occurring on a cell-by-cell basis according to the star formation law,

[TABLE]

where $\dot{\rho}_{\star}$ is the star formation rate density, $\rho_{\rm g}$ the gas density, $t_{\rm ff}=\sqrt{3\pi/32G\rho_{\rm g}}$ is the local free-fall time and $\epsilon_{\rm ff}$ is the local star formation efficiency per free-fall time of gas in the cell. The efficiency is computed following the relation from Padoan et al. (2012), derived from simulations of star formation in magnetized supersonic turbulence111$\epsilon_{\rm ff}=0.5\exp(-1.6t_{\rm ff}/t_{\rm dyn})$, where the dynamical time is $t_{\rm dyn}=L/2\sigma$, and $\sigma$ is the local velocity dispersion compute using neighbouring gas cells over a region of size $L=$ 3 grid cells per spatial dimension..

Each formed star particle is treated as a single-age stellar population with a Chabrier (2003) initial mass function. We account for injection of energy, momentum, mass and heavy elements over time from Type II and Type Ia supernovae (SNe), stellar winds and radiation pressure (allowing for both single scattering and multiple scattering events on dust) on the surrounding gas. Each mechanism depends on the stellar age, mass and gas/stellar metallicity (through the metallicity dependent age-mass relation of Raiteri et al., 1996), calibrated on the stellar evolution code STARBURST99 (Leitherer et al., 1999).

Furthermore, to accurately account for SN feedback we adopt the SN momentum injection model recently suggested by Kim & Ostriker (2015, see also ). A SN explosion is considered resolved when its cooling radius222The cooling radius in gas of density $n$ and metallicity $Z$ scales as $r\approx 30(n/1{\,\rm{cm^{3}}})^{-0.43}(Z/{\,\rm Z}_{\odot}+0.01)^{-0.18}{\,\rm pc}$ for a supernova explosion with energy $E_{\rm SN}=10^{51}$ erg (e.g. Cioffi et al., 1988; Thornton et al., 1998). is captured by at least 6 grid cells. In this case the explosion is initialized in the energy conserving phase by injecting the relevant energy ($10^{51}{\,\rm erg}$ per SN) into the nearest grid cell. If this criterion is not fulfilled, the SN is initialized in its momentum conserving phase, i.e. the total momentum generated during the energy conserving Sedov-Taylor phase is injected into the cells surrounding a star particle. It can be shown (e.g. Blondin et al., 1998; Kim & Ostriker, 2015) that at this time, the momentum of the expanding shell is approximately $4\times 10^{5}(E_{\rm SN}/10^{51}{\,\rm erg})^{16/17}(n/1~{}{\rm cm}^{-3})^{-2/17}{{\,\rm M}_{\odot}{\,\rm{km\,s^{-1}}}}$.

We track iron (Fe) and oxygen (O) abundances separately, and advect them as passive scalars. When computing the gas cooling rate, which is a function of total metallicity, we construct a total metal mass as

[TABLE]

according to the mixture of alpha (C, N, O, Ne, Mg, Si, S) and iron (Fe, Ni) group elements relevant for the sun (Asplund et al., 2009). The code accounts for metallicity dependent cooling by using tabulated cooling functions of Sutherland & Dopita (1993) for gas temperatures of $10^{4-8.5}$ K, and rates from Rosen & Bregman (1995) for cooling down to lower temperatures. Heating from the ultraviolet background radiation is accounted for by using the model of Haardt & Madau (1996), assuming a reionization redshift of $z=8.5$. Self-shielding is modelled following the model of Aubert & Teyssier (2010). Finally, we follow Agertz et al. (2009) and adopt an initial metallicity of $Z=10^{-3}{\,\rm Z}_{\odot}$ in the high-resolution zoom-in region in order to account for enrichment from unresolved population III star formation (e.g. Wise et al., 2012).

2.2 Model verification and analysis

By accounting for the above stellar feedback budget, Agertz & Kravtsov (2015, 2016) simulated the cosmological assembly of Milky Way-mass galaxies, demonstrating that it led to a realistic late-type galaxy matching observed properties such as disc size, the presence of a thin and thick stellar disc, stellar and gas surface density profiles, the Kennicutt-Schmidt relation, the stellar mass-gas metallicity relation (and its evolution), and a specific angular momentum typical of spiral galaxies of the Milky Way mass (stellar mass of $\approx 5\times 10^{10}{\,\rm M}_{\odot}$).

The simulation considered in this work reaches a significantly higher mass and spatial resolution, and matches the above mentioned observational relations (to be presented in Agertz et al. in prep). Furthermore, we find that Hi column density distributions, $f(N_{\textsc{Hi}})$, at all redshifts feature a turnover at column densities $N_{\textsc{Hi}}\sim 10^{21}~{}{\rm cm^{-2}}$, as observed both at high and low redshifts (Zwaan et al., 2005; Noterdaeme et al., 2009), and in the ISM of individual local galaxies (see analysis by Erkal et al., 2012). Our simulation does therefore not suffer from an excess of high column density Hi ($N_{\textsc{Hi}}\gg 10^{21}~{}{\rm cm^{-2}}$) in the inner few kpc of galaxies, as found by Erkal et al. (2012) in their high resolution simulations with inefficient stellar feedback333Using simulations at lower numerical resolution, Altay et al. (2011) found that the $\textsc{Hi}-{\rm H}_{2}$ transition was sufficient to explain the observed turnover, in contrast to Erkal et al. (2012) (see also Altay et al., 2013) who argued that establishing a realistic star formation-feedback cycle is crucial for explaining the column density distribution of DLAs.. As such, our zoom simulation is a relevant platform for studying the nature of DLA and sub-DLA LOS in and around massive disc galaxies.

Using the dark matter halo finder HOP (Eisenstein & Hut, 1998), we identify all dark matter halos444Substructures are not robustly identified using HOP, but as we are mainly interested in central galaxies, this is no issue. within the simulated high resolution volume, in all simulation snapshots (see Section 2.1). For each halo, we generate two-dimensional maps of $N_{\textsc{Hi}}$ and metallicity ([M/H]). We obtain Hi from the neutral gas density in each cell by correcting for its molecular hydrogen content using the model by Krumholz et al. (2008, 2009), as implemented in Agertz & Kravtsov (2015).

All maps cover an area of $200\times 200{\,\rm kpc}^{2}$, and are centered on the dark matter and stars using a shrinking sphere algorithm (Power et al., 2003). The pixels have the same sizes as the finest cell sizes available in the simulation at the given snapshot, typically $\Delta x\sim 20{\,\rm pc}$. In each pixel, the quantities have been computed as

[TABLE]

and

[TABLE]

where the sum runs over all cells along a given pixel LOS along the simulation volume (with a 200 kpc depth), $n_{\textsc{Hi},i}$ is the local cell density of neutral hydrogen, [M/H]i the local cell gas metallicity, and $\Delta x_{i}$ the size of a cell. For each snapshot we compute maps from the three Cartesian directions. In subsequent sections we will present results both for single snapshots and directions to illustrate our findings, as well as results averaged over all analysed directions and different redshift ranges. Finally, quoted impact parameters $b$ are defined as the projected distance from the center of the main galaxy to a point of interest.

Numerical studies of DLAs have generally focused on large simulations volumes in order to capture a statistical sample of galaxies, at the sacrifice of numerical resolution (see Section 1). Indeed, modern cosmological simulations of galaxy formation in domains with sizes of over tens of Mpc, aimed at reproducing the stellar luminosity function, e.g. EAGLE (Schaye et al., 2015) and Illustris (Vogelsberger et al., 2014), have been restricted to spatial resolutions of $\sim 0.5-1{\,\rm kpc}$. In contrast, we are here interested in a specific set of observations relating DLAs and sub-DLAs, at impact parameters up to $b\sim 50{\,\rm kpc}$, to confirmed host galaxies, here predominantly massive and star forming (see Section 3.1), making a comparison to high resolution zoom simulations more appropriate.

The importance of high numerical resolution is illustrated in Figure 1. The top row shows $N_{\textsc{Hi}}$ and [M/H] maps, as well as $b$ vs. [M/H], at $z=1$ for the most massive galaxy in our zoom region ($M_{\star}\sim 5\times 10^{10}{\,\rm M}_{\odot}$). The bottom row shows the same quantities, but for the simulation run at 32 times lower spatial resolution ($\Delta x\sim 640{\,\rm pc}$), and dark matter particle masses of $1.4\times 10^{6}{\,\rm M}_{\odot}$, compatible with the aforementioned large volume simulations. While the low resolution simulation also leads to an extended disc, with a similar stellar mass at $z=0$ (made possible by the stellar feedback driven galactic outflows), the lack of detail in the circumgalactic medium (CGM) compared to the fiducial high resolution model is striking. The fiducial model resolves a greater number of low mass satellites, all found to contribute to DLA and sub-DLA detections (see also Section 3.2). In addition, galaxy interactions, inflows and feedback driven outflows create complex gas structures with large covering fractions (see Sections 3.2 and 3.4). A significant fraction of these gas structures show high levels of metal-enrichment for which feedback and satellite accretion is key; most prominently observed as potential damped H i detections at $b>20{\,\rm kpc}$ in the $b-$[M/H] plane (Section 3.1). All of the above are crucial in order to interpret the observations presented in the next following analysis.

3 Results

3.1 Comparing simulations to analytical models

By combining scaling relations, the simple analytical F08+K17 model was constructed, which successfully describes DLAs as resulting from random LOS through inclined H i discs of the same population of galaxies traced by the Lyman-break technique (LBGs) (Fynbo et al., 2008; Krogager et al., 2017). This model combines a $z\sim 3$ UV luminosity function (Reddy et al., 2008); a $z\sim 3$ metallicity-luminosity relation (Pettini et al., 2001); a Holmberg relation to describe the size of the H i-disc as a function of the galaxy luminosity (following Wolfe et al., 1986); and a prescription for the metallicity gradient. In light of its success, with our simulations we now ask if the assumed scaling relations are reproduced in the simulations and if so at which redshifts.

To answer this, we use HOP to identify the most massive dark matter halo, and follow its evolution in redshift from $z=3$ to $z=0.4$. This halo grows in stellar mass from $9.6\leq\log_{10}[\mathrm{M}_{\star}~{}(\mathrm{M}_{\odot})]\leq 10.8$, which matches the LBG population stellar mass distribution $\log_{10}[\mathrm{M}_{\star}~{}(\mathrm{M}_{\odot})]=9.87\pm 0.53$ (Reddy & Steidel, 2009) from whose luminosity function the F08+K17 model was constructed. We compute $b~{}(\mathrm{kpc})$, $[\mathrm{M/H}]$, and $\log_{10}[N_{\textsc{Hi}}~{}(\mathrm{cm^{-2})}]$ at each redshift and along each side of the simulation-box using the method described in Section 2.2. Binning in redshift intervals, we calculate the angle-averaged detection-fractions ($f_{\mathrm{detections}}$) at each locus. The results are presented in Fig 2 and consists of probability-functions in the $b~{}(\mathrm{kpc})-[\mathrm{M/H}]$ -and the $b~{}(\mathrm{kpc})-\log_{10}[N_{\textsc{Hi}}~{}(\mathrm{cm^{-2})}]$ parameter space, with $f_{\mathrm{detections}}$ normalized for each panel. We show these results together with predictions from the F08+K17 model (with the mean relation in black, and $1\sigma$-; $2\sigma$- and $3\sigma$ intrinsic dispersions (grey) contours), and the MC19 compilation for a sequence of redshifts (columns). Row one shows a simulation snapshot centered on the galaxy; rows two- and three show the space spanned by $b~{}(\mathrm{kpc})-[\mathrm{M/H}]$ with a cut on the column-density at the sub-DLA and at the DLA lower limits, respectively; and row four spans the $b~{}(\mathrm{kpc})-\log_{10}[N_{\textsc{Hi}}~{}(\mathrm{cm^{-2})}]$ parameter space.

Figure 2 reveals significant detection probabilities of LOS above the DLA column density threshold extending to impact parameters of $50~{}\mathrm{kpc}$ from the host center, across all redshifts. However, the distribution clearly changes with redshift. Qualitatively, this can be understood as follows; due to the cosmic density evolution with redshift, $\rho(z)\propto(1+z)^{3}$, at higher redshifts the higher density environments gives rise to more interactions and more tidal and/or turbulent disturbances. At $z>2$, this inhibits gas from settling into coherent and long-lived structures such as extended H i discs. We therefore see $f_{\mathrm{detections}}$ across a broad range of grid-points in the $b~{}(\mathrm{kpc})-[\mathrm{M/H}]$ -and the $b~{}(\mathrm{kpc})-\log_{10}[N_{\textsc{Hi}}~{}(\mathrm{cm^{-2})}]$ parameter spaces alike. This is particularly clear in the highest redshift panel ($z=2-3$, Fig 2 bottom left), for which our simulations generate detections across the whole parameter-space.

Towards progressively lower redshift, merger-rates and densities drop, and the H i gas is allowed to settle. This includes the formation of an H i disc which encodes information on a metallicity-gradient, most prominently observed as the diagonal stripe of high detection-probability (Fig. 2 third row, rightmost panel, but see also Section 3.2 for a direct observation of the disc in this parameter-space). In the same row, we also find a region of diffuse (lower detection probability per unit grid-point) $f_{\mathrm{detections}}$ but extending to large impact parameters. This is indicative of a significant contribution from the circum-galactic medium and the extended halo to the DLA cross-section (see Sections 3.2 and 3.3).

In Table 1 we quantify the agreement between the F08+K17 model and our simulation by calculating the summed detection-fraction enclosed within the $1\sigma$- $2\sigma$- and $3\sigma$ limits. The DLA LOS detection fractions in simulations converge to model predictions towards progressively lower redshifts, reaching a good agreement at $z\sim 0.8$.

Relaxing the threshold on the column-density to include detection-fractions from sub-DLA LOS in $b~{}(\mathrm{kpc})-[\mathrm{M/H}]$ (Figure 2, second row), we find that $f_{\mathrm{detections}}$ becomes less concentrated and covers a larger region in the parameter-space at all redshifts. Phenomenologically, this suggests that sub-DLA LOS fill the void between islands of DLA-$f_{\mathrm{detections}}$, and could be related to a density gradient such that sub-DLAs on average trace lower-density gas. We pursue this by separating the galaxy into a halo and a disc (Section 3.2) and quantify statistical differences between the average DLA and sub-DLA LOS in Sections 3.3 and 3.4.

In Table 2 we demonstrate how the inclusion of sub-DLA column densities affects the match between simulation $f_{\mathrm{detections}}$ and the analytical disc model predictions. We find that such an inclusion results in a worse match at all redshifts. This most probably reflects the underlying Holmberg relation in the construction of the model, which works well down to the lower limit of DLA column densities. Indeed, this statistical discrepancy also matches the observational results of Rhodin et al. (2018), who found an excess of sub-DLAs beyond the predicted $1\sigma$ model region. We therefore conclude that in order to match simulations and observational data below the formal DLA limit, the analytical disc-model must be amended or expanded.

3.2 The fraction of DLA/sub-DLA systems in discs and halos

The relative contributions of interstellar (“disc”) and circum-galactic (“halo”) gas to DLA and sub-DLA detections are not yet known. Despite employing an infinitely thin circular disc to define a geometry on which to distribute the DLA cross-section, the F08+K17 model does not claim to describe galactic H i discs. In simulations, such a separation is complicated as these gaseous components are ill-defined, especially for the earliest stages of disc formation. This is especially true for gas in the interface between the disc and the halo, related to accretion and outflows. The results below should hence be treated as indicative.

We adopt the following simple approach for separating neutral disc gas from halo gas, and apply it to our most massive galaxy in the zoom-in region. For each simulation snapshot, the angular momentum vector is computed for all baryons residing within $3r_{1/2}$ of the central galaxy, where $r_{1/2}$ is the (baryonic) half mass radius. This defines a disc plane, for which we compute the density profile $n(r)$ of neutral gas, where $r$ is the spherical radius. Such profiles are nearly exponential, with a fall-off below $n(r)\sim 10^{-2}\,{\rm cm}^{-3}$ due to the inability of gas to self-shield against the background UV field (also found in the radiative transfer calculations by Aubert & Teyssier, 2010). We define the size of the disc to be the radius $r_{\rm d}$ where $n(r)$ falls below $n=10^{-3}\,{\rm cm}^{-3}$, but note that our conclusions are not very sensitive to cut-off densities in the range to $n=10^{-4}-10^{-2}\,{\rm cm}^{-3}$. At $z<1$, $r_{\rm d}$ fall in the range $\sim 15-30$ kpc. For each simulation snapshot, all neutral gas inside of the cylindrical radius $r_{\rm d}$ and below an altitude above the disc plane of $h=8$ kpc is defined as ’disc’, and the gas outside of this domain is defined as ’halo’.

At high redshifts ($z\gtrsim 1-2$), mergers and and galaxy interactions yield mostly transients discs - if at all. To allow for a cleaner disc-halo separation, we restrict our analysis to $z<1$, and in order to convey that a robust (as opposed to a transient) disc is identified, we stack the results between $z=0.4-1$. In this way, we ensure that our methodology is not biased to the conditions of a single snapshot. The results of the component separation are shown in in Figure 3, for which the top row displays an H i column density map of the galaxy (left), the H i disc (center) and halo (right) at $z=0.6$. Row two and three display the associated detection fractions in the $b-$[M/H] plane, with lower column-density thresholds set to sub-DLA and DLA limits, respectively. Our simple approach allows a clear disc component to be extracted, with the residual halo gas containing satellites, stripped filamentary gas, and gas at smaller impact parameters residing just outside, but soon to be accreted onto, the disc.

The disc component is visually distinct in the $b-$[M/H] relation, both for DLA LOS alone, and when sub-DLA LOS are included, and extends all the way to the $3\sigma$ confidence relation of the analytical F08+K17 model. In fact, $99.6\%$($88\%$) of all simulated DLA sightlines lie within the $2$(1)$\sigma$ relation, indicating that the F08+K17 model does not necessarily represent the physical nature of turbulent discs, but illustrates why this model has been highly successful in capturing observational properties of DLAs (Krogager et al., 2017). For the above disc definition, $\sim 40\%$ ($60\%$) of LOS with $\log_{10}[\mathrm{N}_{\textsc{Hi}}~{}(\mathrm{cm}^{2})]>19$ intersect ’halo’ (’disc’) gas, and $\sim 20\%$ ($80\%$) for $\log_{10}[\mathrm{N}_{\textsc{Hi}}~{}(\mathrm{cm}^{2})]>20.3$. By allowing for a disc definition with $h=4$ kpc, rather than 8 kpc, we find identical results, indicating that our results are robust to reasonable parameter changes. In summary, our experiments illustrate that for the extended disc galaxy of the type analysed here, compatible with $L_{\star}$ galaxies ($M_{\star}\sim{\rm few}\times 10^{10}{\,\rm M}_{\odot}$), DLAs predominantly originate in extended H i discs, but with halo gas playing a role at all impact parameters - especially at high impact parameters and in the sub-DLA regime.

Both simulations and observations indicate a significant amount of metal rich DLA/sub-DLA LOS at high impact parameters ($b\gtrsim 25$ kpc). From the simulated component separation, it is clear that most of the gas contributing to these LOS is of halo origin. This can result from a high metal column, or a low H i column. We interpret this as a clear sign of metal rich gas stripped from dwarf satellites, as the shallow potential wells in these galaxies allow for efficient unbinding of enriched gas via outflows, ram pressure and tidal stripping. This scenario is supported by Figure 1 which shows that the metallicity map does not necessarily follow the H i map but instead shows higher metal content extending from the satellites; and visually by the presence of satellites and stripped tidal tails in the halo component of Figure 3.

3.3 A redshift evolution in the mean impact parameters

The results of Section 3.1 suggest that there is a redshift evolution in the internal distribution of sub-DLA and DLA LOS through a galaxy, and that this evolution may be different for sub-DLAs and DLAs. We know from observations that damped H i absorbers follow a column-density dependent distribution of impact parameters (Christensen et al., 2007; Monier et al., 2009; Fynbo et al., 2010, 2011; Meiring et al., 2011; Rao et al., 2011; Rahmani et al., 2016). Marginalizing over redshifts, these observations state that DLA LOS are found systematically closer to their host galaxy than sub-DLA LOS by a factor of $\sim 2$ (Rao et al., 2011; Rhodin et al., 2018). Likewise, investigations of the $\Delta V_{90}-[\mathrm{M/H}]_{\mathrm{abs}}$ relation (e.g. Ledoux et al., 2006; Møller et al., 2013; Neeleman et al., 2013) have shown that DLAs have smaller mean velocity widths and lower metallicities than sub-DLAs (e.g. Som et al., 2015). Taken together, these observations suggest that the two classes of absorbers on average trace different parts of their hosts.

We now explore whether we can quantify such average differences. The simplest non-parametric statistic to quantify the difference is the mean impact parameters of sub-DLA and DLA LOS. From here onwards, these will be referred to as $\langle b_{\mathrm{sub-DLA}}\rangle$ and $\langle b_{\mathrm{DLA}}\rangle$ respectively. These quantities are particularly interesting from an observational point of view, as impact parameters are amongst the simplest observables and readily available for 43 absorber-galaxy pairs in the MC19 compilation. In addition, the analysis presented in Section 3.1 suggest that DLA LOS can be found at large impact parameters with significant probabilities.

We first calculate $\langle b_{\mathrm{sub-DLA}}\rangle$ and $\langle b_{\mathrm{DLA}}\rangle$. To ensure that the galaxy analysed in previous sections is not an outlier, we allow HOP (see Section 2.2) to identify the seven most massive halos. This allows us to cover diverse galaxy environments, a stellar mass range of $6.7<\log_{10}[\mathrm{M}_{\star}(\mathrm{M}_{\odot})]<10.8$ which covers most observations, and to span a redshift range $0.4<z<3.0$. Using the information on $\langle b_{\mathrm{sub-DLA}}\rangle$ and $\langle b_{\mathrm{DLA}}\rangle$ together with the associated redshift and stellar mass from the three Cartesian directions and for each of the halos, we present the results in Figure 4. The top panel displays the individual data of $\langle b_{\mathrm{sub-DLA}}\rangle$ (orange) and $\langle b_{\mathrm{DLA}}\rangle$ (grey) together with their binned-mean relations in the respective colours, as a function of redshift. The results indicate a co-evolution in mean DLA -and sub-DLA impact parameters towards progressively smaller values with time. However, perhaps the most striking feature is the large scatter in the mean values at any redshift. The scatter is set from the range of galaxy-sizes spanned by the seven halos, and the effect of a fixed $50~{}{\,\rm kpc}$ beam. To cancel the size-dependence and enhance any differential change in the mean impact parameters with redshift, in the bottom panel, we therefore proceed to plot the ratio of the means; $\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle$. A horizontal line implies a scale-free fraction, which at unity translates to identical mean impact parameter sub-DLA and DLA LOS. For reference, we display such a $1:1$ relation as a red dashed line.

To ensure that any correlation in $\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle$ with redshift is indeed real, and not a consequence of an underlying correlation with stellar mass, we calculate the Pearson (PC) coefficient, the Spearman rank-order correlation (SC) coefficient, and their associated null-hypothesis’s p-values. Expressed $\mathrm{Test}~{}(\mathrm{parameter1},~{}\mathrm{parameter2})=[\mathrm{coefficient},~{}\mathrm{p-value}]$ we find $\mathrm{PC}(\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle,~{}\log_{10}[\mathrm{M}_{\star}(\mathrm{M}_{\odot})])=[0.073,\sim 10^{-2}]$ and $\mathrm{SC}(\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle,~{}\log_{10}[\mathrm{M}_{\star}(\mathrm{M}_{\odot})])=[0.144,\sim 10^{-5}]$, whereas significantly higher correlation coefficients are identified for $\mathrm{PC}(\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle,~{}z)=[0.396,\sim 10^{-31}]$, $\mathrm{SC}(\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle,~{}z)=[0.399,\sim 10^{-32}]$. With consistent results from a Pearson and a Spearman test, both strongly suggesting that redshift is indeed the primary driver in the impact parameter evolution, we proceed to ask whether it is driven by a particular halo.

We find consistent correlation coefficients of the same magnitude using both the Pearson and the Spearman test (six positive and one negative correlation) for the individual halos. The magnitude of these coefficients for individual halos are similar to that of the population as a whole, but at elevated p-values in a broad range $\sim[10^{-33}:10^{-1}]$. We therefore proceed to explicitly calculate linear fits to individual halos (Figure 4, grey lines). These fits yield similar gradients that systematically favor an evolution in the mean impact parameter ratio with time. With a fit to all data, we find a mean gradient of $0.16\pm 0.05$ per unit redshift, i.e. at a formal $3.2\sigma$ significance (Figure 4, black solid line).

This mean evolution in $\langle b_{\mathrm{DLA}}\rangle/\langle b_{\mathrm{sub-DLA}}\rangle$-ratio with redshift is remarkably consistent with the analysis presented in Section 3.1. At $z\sim 3$ the ratio intercepts the 1:1 correlation. This is in part caused by the presence of multiple galaxies within the $50~{}\mathrm{kpc}$ projected radius beam - each one generating its own internal DLA and sub-DLA LOSs - which on average will tend to average out any difference between the mean impact parameters discovered in the beam. In part, it is caused by the interaction between galaxies in the beam, which causes irregular gas flows that inhibit the formation of coherent structures, rendering a stochastic distribution. Moving down in redshift, the distribution of sub-DLA and DLA LOS diverge such that, on average, they are separated by a factor of two below redshifts $z\sim 0.5$.

It is striking that such a factor of $\sim 2$ difference in mean projected separation between sub-DLA and DLA LOS is retrieved in observations at redshifts $z<1.0$ (Rao et al., 2011; Rhodin et al., 2018). Indeed, even the absolute mean impact parameters derived in these studies ($\langle b\rangle_{\mathrm{DLA}}^{\mathrm{empirical}}\sim 15~{}{\,\rm kpc}$, $\langle b\rangle_{\mathrm{sub-DLA}}^{\mathrm{empirical}}\sim 30~{}{\,\rm kpc}$) are marginally consistent with some of the high means retrieved for individual galaxies in our simulations at comparable redshifts (see Figure 4, top panel). As per Section 3.1 and Figure 2 bottom row, however, we find significant detection fractions of both sub-DLA and DLA LOS to projected separations $\sim 50~{}{\,\rm kpc}$ at all redshifts, and marginal inconsistencies in the absolute mean values between simulations and observations may reflect (observational) low number statistics and that cosmological zoom simulations have not converged on the spatial scales required to resolve cold gas in the CGM (McCourt et al., 2018; Hummels et al., 2018). We also note that Rahmati & Schaye (2014, R14) identify an anti-correlation between median impact parameter and $N_{\textsc{Hi}}$, for which they report $\langle b\rangle_{\mathrm{sub-DLA}}^{\mathrm{R14}}\sim 7~{}\mathrm{pkpc}$ and $\langle b\rangle_{\mathrm{DLA}}^{\mathrm{R14}}\sim 2~{}\mathrm{pkpc}$. Albeit a large intrinsic scatter in the data, with 15-85 per cent percentiles which we estimate to extend from $2-29~{}\mathrm{kpc}$ and $0.5-4~{}\mathrm{kpc}$ for sub-DLAs and DLAs, respectively (Rahmati & Schaye, 2014, private communication), their median values are significantly lower than those identified in observations and in our simulations. Additionally, despite a consistent factor of two difference in the means, their analysis is based on redshift $z=3$, whereas our analysis of the redshift evolution suggests that such a clear separation is only in place at $z\lesssim 0.7$.

3.4 Do sub-DLA and DLA sight-lines probe galaxies of different mass?

An alternative explanation to the idea that sub-DLA and DLA LOS trace different relations within a galaxy, is that they on average probe intrinsically different galaxy populations; sub-DLA LOS on average probing more massive galaxies than do DLA LOS (e.g. Khare et al., 2007; Kulkarni et al., 2010). Similar conclusions were also reached with deep ground-based imaging of sub-DLA hosts, which suggest that these systems preferentially sample more luminous - and therefore more massive galaxies than DLAs (Meiring et al., 2011). By measuring the stellar mass of the absorbing galaxies from matching spectral energy distribution galaxy templates to photometry, Augustin et al. (2018) and Rhodin et al. (2018) directly confirmed the apparent anti-correlation between the host’s stellar mass and the H i column-density of the absorber.

If real, this correlation is remarkable as we expect both high- and low column-density LOS to pass through a given galaxy. In addition, the shape of the column-density frequency distribution function ($f(N_{\textsc{Hi}},X)$) forces the frequency of sub-DLAs to outnumber that of DLAs. From Prochaska et al. (2014, P14) and with integration limits $[19.00,20.30]_{\mathrm{sub-DLA}}$ and $[20.30,25.00]_{\mathrm{DLA}}$, one finds $N_{\mathrm{DLA}}/N_{\mathrm{sub-DLA}}\sim 0.38$, which together with the apparent correlation above implies an excess of massive gas-rich galaxies relative to low-mass galaxies, at odds with naive expectations from galaxy luminosity functions. A viable alternative is that the correlation is driven by sample selection effects.

Indeed Dessauges-Zavadsky et al. (2009) argue for a metallicity bias at low-redshifts, while at $z>1.7$ DLA and sub-DLA selections show statistically consistent metallicity distributions. Recent work on dust attenuation and quasar surveys confirms that we expect a dust-bias to preferentially act on high $\mathrm{N}_{\textsc{Hi}}$, massive, metal-rich, dusty galaxies selected against quasar LOS (Vladilo & Péroux, 2005; Zafar & Watson, 2013; Fynbo et al., 2013a; Noterdaeme et al., 2015; Krogager et al., 2015; Zafar et al., 2015; Krogager et al., 2016; Fynbo et al., 2017). In Rhodin et al. (2018) we argue that the detection void in the low mass, low H i column-density parameter space could be explained by a functional dependence of the ratio in covering fraction of sub-DLA to DLA LOS with stellar mass, such that sub-DLAs are preferentially detected in more massive galaxies.

To test this prediction, we calculate the covering fraction ($f_{\mathrm{cov}}$) of sub-DLA and DLA LOSs in each halo, and repeat the statistical analysis described in Section 3.3. We present the angle-averaged results in Figure 5. In the top panel, we plot the data for individual galaxies $f_{\mathrm{cov}}$ in DLA (grey) and sub-DLA (orange) LOS, together with binned-mean relations in the respective colours, as a function of stellar mass. $f_{\mathrm{cov}}$ for sub-DLA and DLA LOS correlate with stellar mass such that larger stellar masses, on average, are associated with larger $f_{\mathrm{cov}}$. At stellar mass $\log_{10}(\mathrm{M}_{\star}~{}[\mathrm{M}_{\odot}])\sim 10$ the entire central region sees an increase in neutral gas, which leads to a rapid increase in the final covering fractions. For the most massive disc galaxy in our sample, this increase correlates with rapid disc growth at $z\sim 1.5$. This increase in gas disc sizes of galaxies forming in Miky Way-mass halos, around this redshift and stellar mass, is found in larger samples of zoom simulations (Kretschmer et al. in prep), as well as observationally in terms of the stellar size-mass relation (see compilation by Agertz & Kravtsov, 2015).

Calculating the Pearson correlation coefficient; the Spearman rank-order correlation coefficient, and their associated null-hypothesis’s p-value, we find $\mathrm{PC}(z,f_{\mathrm{DLA}}/f_{\mathrm{sub-DLA}})=[-0.162,\sim 10^{-6}]$ and $\mathrm{SC}(z,f_{\mathrm{DLA}}/f_{\mathrm{sub-DLA}})=[-0.164,\sim 10^{-6}]$, whereas $\mathrm{PC}(\log_{10}(\mathrm{M}_{\star}~{}[\mathrm{M}_{\odot}]),f_{\mathrm{DLA}}/f_{\mathrm{sub-DLA}})=[-0.273,\sim 10^{-15}]$ and $\mathrm{SC}(\log_{10}(\mathrm{M}_{\star}~{}[\mathrm{M}_{\odot}]),f_{\mathrm{DLA}}/f_{\mathrm{sub-DLA}})=[-0.244,\sim 10^{-12}]$. With consistent results, we conclude that a correlation in differential covering fraction with stellar mass indeed dominates over a correlation with redshift by a factor of $\sim 2$. In Figure 5 lower panel, we show the ratio of DLA to sub-DLA covering fractions to envisage the differential scaling with stellar mass, and colour-code the data by redshift.

Having established the subordinate role of redshift to stellar mass in the covering fraction ratio, we proceed to ask whether the correlation is driven by a particular halo. Individual halos display similar correlation coefficients to the that of the population as a whole, but at elevated p-values in the range $\sim[10^{-3}:10^{-1}]$, suggesting that each halo contributes to the final p-value and correlation-strength of the population. We therefore proceed to calculate linear fits to the data of individual halos (overplotted as grey lines). With consistent trends, we find a mean differential covering fraction ratio per unit logarithmic stellar mass of $-0.14\pm 0.04$. This corresponds to a weak functional dependence at a $3.5\sigma$ significance which we over-plot as a black solid line. Qualitatively, this correlation with stellar mass is consistent with the results of Faucher-Giguère et al. (2015) who report that the covering fraction of gas at different column-density thresholds is sensitive to individual dark matter halo masses for LBGs at $z\sim 2$.

Despite weak, if real and acting across the three decades in stellar mass $\sim[8.0<\mathrm{M}_{\star}~{}[\mathrm{M}_{\odot}]<11.]$ identified in observations, we expect this gradient to account for $3\times 0.14\sim 0.4$, i.e. a $40\%$ measured effect. We consider a more realistic value to be that produced within the $1\sigma$ dispersion in the stellar mass distribution of confirmed hosts. For the sample used in Rhodin et al. (2018), with mean and dispersion in stellar mass of $9.87\pm 0.66$, we expect a measured effect of $\sim 20\%$. Whereas we expect this effect to contribute to the observed anti-correlation, a direct comparison to the current set of data is beyond the scope of this work, as such a comparison has to account for sample selection functions and observational strategy bias.

As a reference-point, we instead calculate the predicted ratio of the number of DLAs to sub-DLAs, based on $f(N_{\textsc{Hi}},X)$. With the P14 parametrization, $f(N_{\textsc{Hi}},X)_{\mathrm{P14}}$, we expect a ratio $N_{\mathrm{DLA}}/N_{\mathrm{sub-DLA}}=0.38$. We over-plot this estimator as a dotted purple line. Using the $f(N)$ statistic as an unbiased proxy for the cosmic average covering fraction ratio in DLAs to sub-DLAs, we find an over-density of DLA (or equivalently, an under-density of sub-DLA) LOS within a 50 kpc projected radius beam from the mass center of the hosting halo. From this, we conclude that the majority of sub-DLA are located at separations $b>50~{}\mathrm{kpc}$.

4 Conclusions

In this work, we have used a cosmological zoom-in hydrodynamics simulation of the assembly of a Milky Way mass galaxy to study beams of $50~{}\mathrm{kpc}$ radius centered on galaxies to study the environments capable of producing sub-DLA and DLA LOS. Building probability functions of observables and logging their redshift evolution, we have matched these to the most up-to-date compilation of spectroscopically confirmed absorbing galaxies, and to analytic models. The main results can be summarised as follows:

• •

In agreement with recent results published by other groups, we find that increased numerical resolution leads to an excess of H i coverage in galaxy halos. This increases the mean impact parameters of damped H i absorbers and produces detection probabilities at all observed values.

• •

High numerical resolution and an effective feedback prescription allow us to reproduce the column density distribution function of blind quasar surveys and for individual galaxies. This validates the application of these simulations to study sub-DLA and DLAs in a galactic context.

• •

At redshift $z\lesssim 1$ we statistically match the distribution of impact parameters with metallicity based on the analytical F08+K17 model for DLAs. For higher redshifts, and when including the contribution from sub-DLAs, we find a greater statistical miss-match and excess detection probability beyond the envelope of the analytical model.

• •

To understand the physical origin of the detection-fractions, we perform a disc-halo separation using the baryonic angular momentum vector and the density profile of neutral gas. We find that $\sim 80\%$ of DLAs originate in extended H i discs, with the remaining $20\%$ originating in the halo. Including sub-DLAs increases the relative contribution of the halo to $40\%$. We note that this separation relies on the existence of a clearly defined H i disc.

• •

We find a $3.2\sigma$ significance redshift evolution in the mean impact parameters of sub-DLA and DLA LOS. At high redshifts the mean impact parameters overlap, but separate at successively lower redshifts with sub-DLA LOS on average identified at a factor two larger projected separations than DLAs at $z\sim 0.4$. At all redshifts, our simulations can match the observational data.

• •

We find a $3.5\sigma$ significance anti-correlation in differential covering-fraction of DLA to sub-DLA LOSs with stellar mass. This suggests an observational selection effect that causes the preferential detection of sub-DLAs in more massive galaxies in the low-redshift universe.

Despite the remarkable agreement with observations - which may suggest that damped H i absorption predominantly form in typical gas-rich disc galaxies - we caution against over-interpreting the results. Whereas our simulation qualitatively capture and explain the current set of observations, the results should be anchored in the next generation of zoom-in simulations performed across a large grid to account for environmental dependence and sample the underlying dark matter halo distribution.

Acknowledgements

This work was supported by grant ID DFF-4090-00079. The cosmic Dawn center is funded by the DNRF. This work used the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. OA acknowledges support from the Swedish Research Council (grant 2014- 5791). OA and FR acknowledge support from the Knut and Alice Wallenberg Foundation.