Physical drivers of galaxies' cold-gas content: exploring environmental and evolutionary effects with DARK SAGE

TL;DR

This study combines observational data and semi-analytic modeling to explore how environmental and evolutionary processes influence galaxies' cold gas content and quiescence, highlighting ram-pressure stripping's role in satellite galaxies.

Contribution

It demonstrates how the DARK SAGE model can replicate observed HI fractions and environmental effects, emphasizing the importance of ram-pressure stripping in galaxy evolution.

Findings

Ram-pressure stripping explains the anti-correlation between HI fractions and halo mass.

Model satellites are systematically gas-poor, indicating potential model limitations.

Hot gas stripping influences galaxy quiescence in dense environments.

Abstract

We combine the latest spectrally stacked data of 21-cm emission from the ALFALFA survey with an updated version of the Dark Sage semi-analytic model to investigate the relative contributions of secular and environmental astrophysical processes on shaping the HI fractions and quiescence of galaxies in the local Universe. We calibrate the model to match the observed mean HI fraction of all galaxies as a function of stellar mass. Without consideration of stellar feedback, disc instabilities, and active galactic nuclei, we show how the slope and normalisation of this relation would change significantly. We find Dark Sage can reproduce the relative impact that halo mass is observed to have on satellites' HI fractions and quiescent fraction. However, the model satellites are systematically gas-poor. We discuss how this could be affected by satellite--central cross-contamination from the group-finding algorithm applied to the observed galaxies, but that it is not the full story. From our results, we suggest the anti-correlation between satellites' HI fractions and host halo mass, seen at fixed stellar mass and fixed specific star formation rate, can be attributed almost entirely to ram-pressure stripping of cold gas. Meanwhile, stripping of hot gas from around the satellites drives the correlation of quiescent fraction with halo mass at fixed stellar mass. Further detail in the modelling of galaxy discs' centres is required to solidify this result, however. We contextualise our results with those from other semi-analytic models and hydrodynamic simulations.

Tables

Table 1. Table 1: Number of satellite galaxies in the various halo mass bins for both the observations and full model of Dark Sage (which, apart from the cut of m ∗ ≥ 10 9 ​ M ⊙ subscript 𝑚 superscript 10 9 subscript M direct-product m_{*}\!\geq\!10^{9}\,{\rm M}_{\odot} , comes from Millennium), for which results are shown in Figs.  7 & 8 . N p , max subscript 𝑁 p max N_{\rm p,max} refers to the historical maximum number of particles that each subhalo was composed of in the Millennium merger trees.

	$N_{\mathrm{sat}}$
${\log }_{10}\left(M_{\mathrm{vir}}/M_{\odot }\right)$		$\in [12,13.5)$	$\in [13.5,15)$
Observations	1 030 (11.9%)	3 486 (40.2%)	4 148 (47.9%)
Dark Sage, all	194 240 (13.3%)	704 781 (48.4%)	544 154 (37.4%)
Dark Sage, $N_{p,\max }\ge 100$	155 430 (12.9%)	585 254 (48.6%)	453 286 (37.6%)

Equations

$$R_{{\mathrm{H}}_2}(r)\equiv \frac{{\Sigma }_{{\mathrm{H}}_2}(r)}{{\Sigma }_{\mathrm{H}\text{textsc}\mathrm{i}}(r)}={\left[{\left(1-\frac{3s(r)}{4+s(r)}\right)}^{-1}-1\right]}^{-1},$$

$$s(r)\equiv \frac{\ln \left(1+0.6\chi (r)+0.01{\chi }^2(r)\right)}{0.6\tau (r)},$$

$$\chi (r)\equiv \frac{3.1}{4.1}\left[1+3.1{\left(\frac{Z(r)}{{\mathrm{Z}}_{\odot }}\right)}^{0.365}\right],$$

$$\tau (r)\equiv 0.66\frac{Z(r)}{{\mathrm{Z}}_{\odot }}\frac{c_f(r){\Sigma }_{\mathrm{gas}}(r)}{{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}},$$

$$c_f(r)\equiv \varphi {\left(\frac{Z(r)}{{\mathrm{Z}}_{\odot }}\right)}^{-\gamma },$$

$${\rho }_{\mathrm{hot},\mathrm{cen}}(R_{\mathrm{sat}})v_{\mathrm{sat}}^2\ge 2\pi G{\Sigma }_{\mathrm{gas}}(r)\left[{\Sigma }_{\mathrm{gas}}(r)+{\Sigma }_{\ast }(r)\right],$$

$${\rho }_{\mathrm{hot},\mathrm{cen}}(R_{\mathrm{sat}})v_{\mathrm{sat}}^2=\frac{G{M}_{\mathrm{sat}}(<R)m_{\mathrm{hot},\mathrm{sat}}}{8R_{\mathrm{vir},\mathrm{sat}}{R}^3},$$

Full text

††pagerange: Physical drivers of galaxies’ cold-gas content: exploring environmental and evolutionary effects with Dark Sage–LABEL:lastpage††pubyear: 2017

1 Introduction

It is widely accepted that the seeds of galaxies are born out of the gravitational fragmentation of gas within dark-matter-dominated overdensities in the Universe (á la White & Rees 1978), more commonly referred to as ‘haloes’. Our modern theoretical consensus is that galaxies then continue to accrete gas cosmologically, either directly through cold streams or via their hot circumgalactic medium, where ionised gas must cool to a predominantly neutral phase before settling into the galactic disc (see Rees & Ostriker 1977; White & Frenk 1991; Birnboim & Dekel 2003; Kereš et al. 2005; Dekel & Birnboim 2006; Benson & Bower 2011; Stevens et al. 2017). In addition, the hierarchical nature of galaxy growth means galaxies can acquire large volumes of gas at once through mergers. As the principal component of gas in galaxies, and because cold gas is the raw material for forming new stars, neutral atomic hydrogen (H i) is one of the most important ingredients in the formation and evolution of galaxies. Therefore, understanding the ensemble of internal and external astrophysical mechanisms that regulate H i content and star formation activity is critical if we are to realise a coherent picture of galaxy evolution.

Internally, beyond the direct consumption of gas in the formation of stars and the accretion of black holes, the evolution of galaxies’ gas reservoirs is consequently dictated by feedback from stellar evolution and active galactic nuclei. This feedback manifests in the form of energetic winds that eject and/or heat gas reserves, as evidenced and detailed by a variety of observations (e.g. Sharp & Bland-Hawthorn 2010; Feruglio et al. 2010; Anderson et al. 2015; Cicone et al. 2015; Nielsen et al. 2016) and simulations (e.g. Springel & Hernquist 2003; Brook et al. 2011; Costa, Sijacki & Haehnelt 2015; Taylor & Kobayashi 2015; Bower et al. 2017). Instabilities and turbulence within the interstellar medium also help regulate the distribution of gas and its ability to form stars (see, e.g., Federrath 2015; Stevens, Croton & Mutch 2016). While these secular processes all impact galaxies’ gas fractions, it has become increasingly clear that external effects, driven by environment, are of comparable importance.

The role of environment in the suppression of H i was first demonstrated by observations showing cluster galaxies to be gas-poor compared to the field (Giovanelli & Haynes 1985; Solanes et al. 2001). This was exemplified in exquisite detail by the VLA***Very Large Array Imaging of Virgo in Atomic gas survey (VIVA; Chung et al. 2009), which demonstrated the importance of the highest density environments in shutting down star formation via strong gas depletion mechanisms. Other observational (e.g. Cortese et al. 2011; Catinella et al. 2013) and theoretical (e.g. McCarthy et al. 2008; Rafieferantsoa et al. 2015; Marasco et al. 2016) efforts support this picture, with both camps generally describing external processes that are distinguishable based upon the time-scales over which they act: (i) those that act swiftly and directly upon the cold gas to remove it from the galaxy via an interaction of the interstellar and intracluster (or intragroup) media ($\sim$10s of Myr, i.e. ram-pressure stripping; Gunn & Gott 1972), or (ii) those that regulate the rate at which gas is able to accrete onto the galaxy from its dark-matter halo over more lengthy time-scales (i.e. strangulation, where galaxies consume their available gas for star formation in $\gtrsim$1 Gyr; Larson, Tinsley & Caldwell 1980). Using this distinction, Brown et al. (2017) provide strong evidence that the gas loss in massive haloes ($M_{\rm vir}>10^{13}\,{\rm M}_{\odot}$) is considerably faster than the subsequent quenching of star formation.

To zeroth order, environment can be thought of as a dichotomy between ‘central’ and ‘satellite’ galaxies. These terms are most relevant for numerical simulations, where a central galaxy belongs to the most massive subhalo of a halo (see, e.g., Springel et al. 2001; Onions et al. 2012). All remaining subhaloes host satellite galaxies. Under this definition, a central can be both isolated or have any number of satellites associated with it. By definition, only satellite galaxies are prone to gas-stripping processes. For observational data to be interpreted in this framework, group-finding algorithms are typically employed to arrange galaxies into haloes. These use a known halo mass function from simulations with abundance matching to assign halo masses, where the brightest or most massive galaxy is labelled the central and the remainder are tagged as satellites (see, e.g., Huchra & Geller 1982; Tucker et al. 2000; Yang et al. 2005; Campbell et al. 2015). The strength of environmental stripping processes depends primarily on the parent halo mass, within which the satellite resides, as evidenced by their H i fractions (Brown et al. 2017).

In recent years, single-dish surveys have begun to provide H i measurements for tens of thousands of galaxies. For example, the H i Parkes All Sky Survey (HIPASS; Meyer et al. 2004) scanned $\sim$75 per cent of the sky (30,000 deg2) and detected 21-cm emission in $\sim$5000 nearby galaxies. Its successor, the Arecibo Legacy Fast ALFA†††Arecibo L-band Feed Array survey (ALFALFA; Giovanelli et al. 2005) went even further, and the latest data release ($\alpha$.70) provides a census of H i gas for $\sim$25 000 galaxies over $\sim$7000 deg2 of sky (for a presentation of the earlier $\alpha$.40 release, see Haynes et al. 2011). In addition, the stacking of undetected H i sources based on optical position and redshift has become be a major tool for pushing 21-cm surveys beyond their nominal sensitivity limit, providing representative, statistical studies of H i as a function of galaxy properties and environment in the local Universe (Fabello et al. 2011a, b, 2012; Brown et al. 2015, 2017).

In their paper, Brown et al. (2017) compare their average gas content scaling relations with predictions from the cosmological hydrodynamic simulations of Davé et al. (2013) and the galform semi-analytic model (Gonzalez-Perez et al. 2014). They find that both sets of predictions are too gas-poor at fixed stellar mass and specific star formation rate (sSFR), although there is qualitative reproduction of ram-pressure stripping in the group and cluster regime. Marasco et al. (2016) also investigate the role of environment in dictating the H i content of galaxies within the EAGLE‡‡‡Evolution and Assembly of GaLaxies and their Environments suite of cosmological hydrodynamic simulations (Schaye et al. 2015). They use the simulations to successfully reproduce the findings of Fabello et al. (2012) and Catinella et al. (2013), demonstrating that EAGLE galaxies in more massive haloes have lower gas content at fixed stellar mass. However, rather than a continuous trend of H i depletion as a function of environment (e.g. Stark et al. 2016; Brown et al. 2017), their work finds the role of environment to be more binary, determining whether or not galaxies have any H i at all. Higher-resolution runs of these simulations have found the H i fractions to increase (Bahé et al. 2016; Marasco et al. 2016; Crain et al. 2017).

Despite significant progress, large uncertainties remain concerning where and how precisely, in terms of epoch and environment, secular and environmental processes influence galaxy gas reservoirs. We therefore use the Dark Sage semi-analytic model (Stevens et al. 2016) to investigate the contributions of various evolutionary and environmental processes in determining the gas fractions and star formation activity of both central and satellite galaxies, relating this to their stellar masses and specific star formation rates. Because semi-analytic models describe the phenomenology of hydrodynamical effects without directly modelling hydrodynamics at all, all gas-stripping processes are manually prescribed. The advantage this gives is that we can learn about the relative roles each effect has on the gas content of galaxies by turning the effects on and off. Dark Sage is especially well equipped, as the one-dimensional disc structure of galaxies is self-consistently evolved in the model, naturally leading to a spatial dependence of environmental effects within galaxies.

This paper is structured as follows. We describe the main aspects of the Dark Sage semi-analytic model in Section 2, including new features we have introduced for this work. In Section 3, we describe the observational data with which our results are closely compared. We examine the H i fractions and quiescence (lack of star formation) of galaxies within the central–satellite dichotomy in Section 4 by comparing a set of Dark Sage model variants against our observed data. The role of environment in shaping the H i fractions and quiescence of satellite galaxies is then studied in greater detail in Section 5. We discuss the context of our results with recent literature and offer concluding remarks in Section 6.

2 The semi-analytic model

The Dark Sage semi-analytic model (Stevens et al. 2016) is a heavily modified version of the models developed by Croton et al. (2006, 2016). Whereas most semi-analytic models evolve the integrated properties of singular baryonic reservoirs, Dark Sage evolves the one-dimensional structure of galactic discs in annuli of fixed specific angular momentum (similar to Stringer & Benson 2007; Dutton & van den Bosch 2009), where $j\!=\!r\,v_{\rm circ}$. Evolutionary and environmental processes thus affect each galaxy on local scales. We refer the reader to section 3 of Stevens et al. (2016) for full details, but briefly outline the key features of the model here.

• •

Cooling: Hot gas cools onto galaxies following the method introduced by White & Frenk (1991), using the cooling tables of Sutherland & Dopita (1993). The cooling gas is assumed to have an exponential surface density profile as a function of $j$ and to immediately localise itself with material in the disc of the same $j$.

• •

Star formation, evolution, and feedback: Passive star formation goes as $\Sigma_{\rm SFR}\!\propto\!\Sigma_{\rm H_{2}}$, where 43% of star-forming gas is approximated as instantly recycled (cf. Cole et al. 2000). Every star formation episode results in supernova feedback, which heats gas out of the disc. The mass of gas heated in an annulus is directly proportional to the mass of stars formed and inversely proportional to the local gas surface density. The latter approximates supernova energy as being more readily dissipated in denser gaseous regions (this model follows Fu et al. 2010, 2013). Excess energy from supernovae can eject gas out of the halo.

• •

Instabilities: A two-component Toomre $Q$ value (Toomre 1964; Romeo & Wiegert 2011) for each annulus is routinely calculated. If $Q<1$, a starburst in that annulus can occur, and any remaining unstable mass is transferred to the two adjacent annuli in proportion such that angular momentum is conserved. Sinking unstable stars in the innermost annulus are transferred to a fully dispersion-dominated bulge component, while gas can be accreted onto the central black hole.

• •

Mergers: For minor mergers, the specific angular momentum of the merging satellite is measured relative to the central galaxy to determine where in the disc the satellite’s cold gas ends up, while its stars are transferred to the central’s bulge. For major mergers, the gas disc vectors are summed to define the new plane onto which both gas discs are projected, while both stellar discs are destroyed and form a bulge. Starbursts occur in the annuli where gas originated from both systems, based on a variation of Somerville, Primack & Faber (2001). Black holes are directly fed gas during mergers phenomenologically (following Kauffmann & Haehnelt 2000).

• •

Active galactic nuclei: Radio mode accretion of hot gas onto a central black hole and subsequent feedback occurs at each time-step, whereby hot gas is prevented from cooling (Croton et al. 2016). Quasar mode feedback is triggered by mergers and instabilities. The model steps out annulus by annulus, determines whether the quasar energy is sufficient to strip the cumulative gas up to that annulus, and if so, transfers all the gas in that annulus to the hot reservoir. The remaining energy can eject gas out of the halo entirely.

By default, whenever gas is reheated by feedback, it is transferred to the hot and/or ejected reservoirs associated with the central galaxy within the halo, even if the gas originated from a satellite (but see Section 2.1.3). While the choice of whether reheated gas remains in a satellite’s hot reservoir or not can have some effect on the gas properties of model galaxies (e.g., as raised by Font et al. 2008; Lagos et al. 2014), this effect is entirely negligible for the results we present for Dark Sage.

Dark Sage has been run on the standard merger trees of the Millennium simulation (Springel et al. 2005). The simulation included $2160^{3}$ particles of mass $1.18\!\times\!10^{9}\,{\rm M}_{\odot}$ in a periodic box with a comoving length of 685 Mpc. This assumes $h\!=\!0.73$, $\Omega_{M}\!=\!0.25$, $\Omega_{\Lambda}\!=\!0.75$, $\Omega_{b}\!=\!0.045$, and $\sigma_{8}\!=\!0.9$ (Spergel et al. 2003). The merger trees carry 64 snapshots. As in sage (Croton et al. 2016), galaxies in Dark Sage are evolved on 10 sub-time-steps between each snapshot.

Galaxy catalogues built from Dark Sage are available through the Theoretical Astrophysical Observatory (Bernyk et al. 2016).§§§https://tao.asvo.org.au Here, photometric properties and light-cones for Dark Sage galaxies can be constructed. The codebase itself has now also been made publicly available.¶¶¶https://github.com/arhstevens/DarkSage

2.1 Updates to the model

For the purposes of this work, we have updated some of the default prescriptions in Dark Sage and added a few optional variants to the code as well. These were added to improve the average gas content of the entire galaxy population at low stellar masses (see Section 2.2) and to increase our ability to investigate the impact of environment on gas content (see below).

2.1.1 Atomic and molecular hydrogen

The amount of gas in the form of H i and H2 in each annulus of each Dark Sage galaxy is regularly recalculated. In Stevens et al. (2016), this was based on the mid-plane pressure of the disc (Blitz & Rosolowsky 2004), and closely followed the prescription of Fu et al. (2010). While we include that prescription in this work, we have now opted for a metallicity-based prescription as our default, where

[TABLE]

where $\Sigma(r)$ is the local surface density of the subscripted quantity [where $\Sigma_{\rm gas}(r)$ is for all gas in the disc], $Z(r)$ is the local surface density ratio of gaseous metals to all gas, and $c_{f}(r)$ is the local clumping factor. This is based on the work of McKee & Krumholz (2010) and follows the implementation of Fu et al. (2013). These equations originate from modelling spherical gas clouds subject to a photoionizing ultra-violet background, where metallicity represents the dust content of those clouds, which shield the cloud from the ionizing photons, allowing a greater fraction of the gas to be in molecular form (for further details, see Krumholz, McKee & Tumlinson 2008, 2009). Here, we assume a solar metallicity of ${\rm Z}_{\odot}\!=\!0.02$ and treat $\phi$ and $\gamma$ as free parameters, for which we find values of 3.0 and 0.3, respectively, after recalibrating the model (Section 2.2).

Note that $R_{\rm H_{2}}(r)$ can only be calculated by equation (1) provided $s(r)\!<\!2$. This requires gas to have non-zero metallicity. In Dark Sage, gas begins with zero metallicity. If passive star formation from H2 were the only means of forming stars in the model, then we would have a paradoxical chicken–egg scenario with neither chickens nor eggs, and hence star formation would never ignite with this method. This is resolved by instabilities and mergers, which both induce starbursts independent of any H2 (for full details, see Stevens et al. 2016). In other words, once a galaxy has experienced an instability or merger, which tends to happen quickly, it can form stars regularly from H2 via equation (1).

In this paper, we show results from Dark Sage using both the metallicity- and pressure-based prescriptions for determining the ratio of H2 (and H i) in gas disc annuli. We denote these as $f_{\rm H_{2}}(Z)$ and $f_{\rm H_{2}}(P)$, respectively, where $f_{\rm H_{2}}$ is formally defined as the mass ratio of H2 to all gas within an annulus, as per equation 12 of Stevens et al. (2016), and hence has a one-to-one mapping with $R_{\rm H_{2}}$.

By changing the prescription for $f_{\rm H_{2}}$, the relative contributions of each star formation channel to the final stellar populations at $z\!=\!0$ are altered for galaxies of mass $m_{*}\!\lesssim\!10^{10}\,{\rm M}_{\odot}$. This is demonstrated in Fig. 1. Compared to $f_{\rm H_{2}}(P)$, switching to $f_{\rm H_{2}}(Z)$ reduces the contribution of the passive H2 channel of star formation, and increases the contributions from both instability- and merger-driven starbursts. Regardless of the $f_{\rm H_{2}}$ prescription though, the stars in galaxies of mass $m_{*}\!\lesssim\!10^{9}\,{\rm M}_{\odot}$ are predominantly born out of the H2-dependent passive channel of star formation, whereas the stars in galaxies of mass $m_{*}\!\gtrsim\!10^{9.5}\,{\rm M}_{\odot}$ are mainly born in instability-driven starbursts. This has important consequences for some of the runs of the model we present in this paper, which we come back to in Section 4.3.

2.1.2 Ram-pressure stripping of cold gas

As was the case in Stevens et al. (2016), the cold gas in each annulus of each satellite galaxy is subject to ram-pressure stripping from the hot gas associated with the central (most massive) subhalo, based on the landmark work of Gunn & Gott (1972):

[TABLE]

where $\rho_{\rm hot,cen}(R_{\rm sat})$ is the local density of the intracluster medium (or the circumgalactic medium of the central) at the position of the satellite and $v_{\rm sat}$ is the relative speed of the satellite to the central. If this inequality is met for an annulus, all gas in that annulus is transferred to the hot reservoir of the central.

However, equation (2) does not provide a criterion for when a disc will feel ram pressure; it simply tells us if it is vulnerable to ram pressure how much gas will be stripped. In order for a disc to be vulnerable, the satellite’s hot gas on its leading side must first be stripped away. While the stripping of a satellite’s hot gas would be asymmetric in reality, Dark Sage has built-in symmetry assumptions about matter in (sub)haloes (like most semi-analytic models). Bearing all of this in mind, we have added a condition to Dark Sage that equation (2) is only applied to satellites if the total baryonic mass (cold gas + stars) of the galaxy exceeds the hot-gas mass of the subhalo. While this phenomenological condition misses a lot of the complex physics surrounding ram-pressure stripping, it provides a simple means of ensuring that some of a satellite’s protective hot gas must be lost before cold gas in the disc is susceptible to stripping. This is more important for massive satellites (see, e.g., fig. 2 of Croton et al. 2016).

2.1.3 Stripping of hot gas in subhaloes

As an alternative to hot gas being stripped in proportion to dark matter (see Croton et al. 2016), which we maintain as the default for Dark Sage, we have now included an option for a ram-pressure-stripping-like prescription of the hot gas in subhaloes (à la Font et al. 2008; McCarthy et al. 2008). Specifically, we assume equation 6 of McCarthy et al. (2008) describes the gravitational restoring force per unit area on the hot gas of the satellite, which is assumed to be distributed in a singular isothermal sphere. We then find the radius $R$ from the satellite’s centre where

[TABLE]

which equates ram pressure with the gravitational restoring force density (where subscript ‘sat’ implies the subhalo the satellite is associated with, and $M_{\rm sat}$ accounts for all matter). Any hot gas in the subhalo external to this radius is assumed to be stripped over the course of the simulation snapshot interval, meaning one tenth of that mass is stripped in a sub-time-step on which the galaxies are evolved. For the sake of simplicity, after each stripping calculation, the remaining hot gas is assumed to redistribute itself into a singular isothermal sphere again, such that equation (3) is always applicable. If left unchecked, this artificial redistribution of hot gas to larger radii would lead to continually stronger stripping on each consecutive sub-time-step. To counteract this, the stripped hot-gas mass at a given sub-time-step is not allowed to exceed that of the first sub-time-step within the same snapshot interval. Once the next snapshot is reached, the virial radius of the subhalo will have reduced, and hence so will the extent to which the hot gas is assumed to be distributed.

We also include the option of not stripping hot gas at all. With this option, we also return reheated gas from feedback to the hot component of the subhalo, rather than that of the main halo, which would otherwise be standard.

2.2 Recalibration

Before investigating how environmental and evolutionary processes can affect the gas fractions of galaxies, we first need to ensure the average gas fraction of Dark Sage galaxies is representative of the real Universe. This was always a constraint used in calibrating the model (using the data from Brown et al. 2015), but there was still an overall deficit in the gas content of galaxies at low stellar mass (see appendix A of Stevens et al. 2016). This has now been alleviated through a combination of the updates to the code above and opting for the H i fraction of galaxies to become the constraint of greatest emphasis. We further exclude any subhaloes that were never composed of at least 100 particles in their history, ensuring that we use well-resolved galaxies in calibrating the model. This is now in much finer agreement with Brown et al. (2015), as presented in Fig. 2.

The other $z\!=\!0$ constraints of the model have been maintained as well. These include the stellar, H i, and H2 mass functions, the black hole–bulge mass relation, the stellar mass–gas metallicity relation, and the Baryonic Tully–Fisher relation. See section 3 and appendix A of Stevens et al. (2016, and references therein) for further details.

We have only modified one parameter from the original version of Dark Sage; we have reduced the passive star formation efficiency from H2 by a factor of 3 to $\epsilon_{\rm SF}\!=\!1.3\!\times\!10^{-4}\,{\rm Myr}^{-1}$ (cf. table 1 of Stevens et al. 2016). This means a greater fraction of stars are formed through disc instabilities. Specifically, with the old parameter value, the dominant star formation channel transitioned from passive H2 to instabilities at a factor of $\sim$10 higher in stellar mass versus what is seen now in Fig. 1. This has a minimal effect on the stellar mass function, but it is important for when stars form in a galaxy, and for their total H i and H2 content. This parameter change has been implemented for all runs of the model presented in this paper.

2.3 Model galaxy sample

Throughout this paper, we study galaxies from Dark Sage at $z\!=\!0$, and only include those with $10^{9}\!<\!m_{*}/{\rm M}_{\odot}\!\leq\!10^{11.5}$, matching our sample of observed galaxies, described in Section 3. In principle, we could have used the $z\!=\!0.041$ snapshot from Millennium, which would have more accurately matched the redshift range for the observations. However, the evolution of galaxies is sufficiently low at $z\!\lesssim\!0.05$ such that our results would have no noticeable changes.

We note that while the H i and stellar mass measurements we present for Dark Sage galaxies are instantaneous quantities for the $z\!=\!0$ snapshot, star formation rates are average quantities over the period since the previous snapshot. The time-scale for the star formation rates is hence $\sim$260 Myr. While this is an order-of-magnitude difference in time-scale from the H$\alpha$ measurements for the observed galaxies (cf. Section 3 of this paper and section 6 of Brown et al. 2017), it is not a cause of significant concern. Benson et al. (2012) has shown that improving the temporal resolution of a simulation above Millennium’s (64 snapshots) will only result in a $\lesssim$5 per cent change in the universal average SFR from the galacticus semi-analytic model (Benson 2012) at any time. The same is true for the stellar mass function at $z\!\lesssim\!4$ (see fig. 11 of Benson et al. 2012). Stevens (2016, see appendix B) showed similar results for sage (Croton et al. 2016) using the GiggleZ simulation suite (Poole et al. 2015). Using the same simulations, we have confirmed that any sSFR-related results for Dark Sage are negligibly affected by temporal resolution, whether we increase it by a factor of 4 or decrease it by a factor of 2 (not shown here).

3 Observational data

This section describes the observational data used in this paper. Galaxies are selected according to their stellar mass ($m_{*}\!\geq\!10^{9}\,{\rm M}_{\odot}$) and redshift ($0.02\!\leq\!z\!\leq\!0.05$) from the overlap in volume between the Sloan Digital Sky Survey Data Release 7 (SDSS DR7; Abazajian et al. 2009) and the ALFALFA survey (Giovanelli et al. 2005). This yields a representative parent sample of 30,695 galaxies that is defined in Section 2 of Brown et al. (2015). We refer the reader to that paper for a more detailed description.

H i 21-cm spectra are extracted from ALFALFA data cubes using the SDSS DR7 coordinates and redshift. This provides spectra for every galaxy in the sample, regardless of whether it is a formal H i detection or not.

The MPA-JHU SDSS DR7 catalogue∥∥∥http://www.mpa-garching.mpg.de/SDSS/DR7/ and improved stellar masses from http://home.strw.leidenuniv.nl/~jarle/SDSS/ provides stellar mass and star formation rate (SFR) estimates for the entire sample. Masses are constrained following the methodology of Kauffmann et al. (2003) using stellar absorption lines and fits to the optical photometry. The catalogue follows Brinchmann et al. (2004) in deriving global SFRs by using fits to emission H$\alpha$ line fluxes and, where no or low-signal-to-noise lines are detected, the 4000-Å break strength. We note these probe completely different time-scales, but, as raised in Section 2.3, both can be safely compared with our model results.

Halo masses are assigned using the SDSS DR7 group catalogue,******http://gax.shao.ac.cn/data/Group.html detailed in Yang et al. (2007). We refer the reader to their paper for a full description. Briefly, the authors use the friends-of-friends, iterative group-finding algorithm of Yang et al. (2005) to identify galaxy groups with SDSS DR7. Groups are then assigned a halo mass via the abundance matching technique, whereby individual group luminosity or stellar mass is ranked and then matched to the halo mass function of Warren et al. (2006). For this work, we use the halo masses estimated using the stellar mass rank order.

For a full breakdown of the optical and H i data, see sections 2.1 and 2.2 of Brown et al. (2015). The sample’s stellar mass, SFR, and halo mass estimates are described in more detail in section 2 of Brown et al. (2017). The galaxy and halo properties presented in those papers assumed a Kroupa (2001) initial mass function (IMF) and a dimensionless Hubble parameter $h\!=\!0.7$. In order to be consistent with the model galaxies we compare against (see Section 2), all quantities have been converted to assume $h\!=\!0.73$ and a Chabrier (2003) IMF. We have assumed $0.66\,m_{\rm*,Chabrier}\!=\!0.61\,m_{\rm*,Kroupa}$ and $0.67\,{\rm SFR}_{\rm Chabrier}\!=\!0.63\,{\rm SFR}_{\rm Kroupa}$ (as in Madau & Dickinson 2014).

3.1 H i spectral stacking

As outlined in the Section 1, H i spectral stacking improves the statistical size of studies way beyond what is currently possible using only detections, probing the gas content of galaxies across the entire gas-rich to -poor regime. The stacking technique in this paper is based upon the technique developed by Fabello et al. (2011a) and described fully in section 3 of Brown et al. (2015). To briefly summarise, we shift the H i spectra of a given ensemble of galaxies to their rest frame velocity and co-add them regardless of their detection status, weighting each by the inverse of their corresponding stellar mass. Doing so yields a weighted average stacked profile and, therefore, the average H i fraction, $\langle m_{\mathrm{H}\,\LARGE\textsc{i}}/m_{*}\rangle$, for the galaxies in each stack.

Errors on the stacking results are estimated using a statistical delete-a-group jackknife routine whereby a random 20 per cent of spectra in a given stack are iteratively discarded without replacement, each time re-estimating the average gas fraction. This provides five separate estimates of gas fraction with each spectra discarded once. Following equation 3 in Brown et al. (2015), the error is then calculated as the standard deviation of these five average gas fraction measurements multiplied by a $\sqrt{N-1}$ factor, where $N\!=\!5$.

4 Satellites versus centrals

In this section, we break both the observed and model galaxies into centrals and satellites, and study the relative differences in their H i fractions and quiescence. For the purposes of this paper, a central is the most massive galaxy of a halo, including isolated galaxies that are the sole occupant of a halo. Satellites constitute all other galaxies. We perform multiple runs of Dark Sage with various environmental and evolutionary processes switched on and off to understand the effect each has on the gas content of galaxies. In each panel of Fig. 3, we present the mean H i fraction of the observed galaxies for fixed values of stellar mass, comparing to one of the runs of Dark Sage. Similarly, in Fig. 4, we assess H i fraction as a function of specific star formation rate (sSFR). It is important that we consider this, as Brown et al. (2015) showed that near-ultraviolet$-r$ band colour, which is a proxy for sSFR, is more tightly related to H i fraction than stellar mass is. For Dark Sage, we use bins of width 0.2 dex for Figs. 3–7. We plot the log of the mean value along each axis for each bin, for both observed and model galaxies.

In drawing conclusions from comparisons of the H i fractions of observed and model galaxies, one needs to be mindful of how the relative H i and H2 content in the model galaxies is determined. Ideally, we would simultaneously examine the results of Dark Sage with the mean H2 fractions of observed galaxies. Unfortunately, only limited samples of galaxies have their H2 content inferred from CO observations (e.g. Young et al. 1995; Leroy et al. 2009; Saintonge et al. 2011a; Boselli, Cortese & Boquien 2014a), and no H2-blind survey analogous to ALFALFA currently exists (although ongoing surveys such as xCOLDGASS will help alleviate this – Saintonge et al. in preparation). There is, however, a wealth of evidence, drawn primarily from these surveys, that closely connects the H2 in these galaxies to their star formation activity (e.g. Kennicutt 1998; Rownd & Young 1999; Bigiel et al. 2008; Saintonge et al. 2011b, 2016; Boselli et al. 2014b). We can, therefore, use the quiescent fraction of galaxies in observations and models as a first order indication of the relative levels of H2. We present the quiescent fractions of central and satellite galaxies for the various runs of Dark Sage, as a function of stellar mass, alongside the same quantity for our observed sample, in Fig. 5. For our purposes, we consider a quiescent galaxy to be one with ${\rm sSFR}\!<\!10^{-11}\,{\rm yr}^{-1}$.

4.1 Observations and the full model

As seen in all panels of Figs. 3 & 4, the difference in the mean H i fraction between the observed centrals and satellites, whether at fixed stellar mass or specific star formation rate, is $\lesssim$ 0.2 dex. Both at $m_{*}\!\lesssim\!10^{9.5}\,{\rm M}_{\odot}$ and ${\rm sSFR}\!\lesssim\!10^{-1.5}\,{\rm Gyr}^{-1}$, there is little to separate centrals from satellites. At higher stellar masses, satellites have less H i than centrals, suggesting that a significant population of satellites are subject to phenomena (stripping, strangulation, etc.) that deplete their gas reservoirs (á la Gunn & Gott 1972; Larson et al. 1980).

At high sSFR, satellites have greater H i fractions than centrals. Two effects are important for this comparison. First is the rate at which gas depletion takes place in satellites. If depletion processes are slow-acting, the sSFR and H i fraction will decrease for an individual galaxy at comparable rates. If they are of moderate speed or fast-acting, the H i fraction will be observed to decrease before sSFR (see Balogh, Navarro & Morris 2000; McCarthy et al. 2008; Vollmer et al. 2012; Brown et al. 2017). One might naïvely except then that at fixed sSFR, satellites would have lower H i fraction than centrals, which would oppose what is observed. This expectation would only be legitimate if one were comparing equivalent populations of satellites and centrals (i.e. of the same stellar mass in each bin), however. In Fig. 6, we present the mean stellar mass of satellites and centrals for equivalent bins in sSFR as used in Fig. 4. Here, we see that at fixed sSFR, centrals tend to have greater stellar mass. As such, they are biased towards lower gas fractions (cf. Fig. 3).

The complete model of Dark Sage (regardless of the prescription for the breakdown of atomic and molecular hydrogen) shows a much greater separation in the H i fractions of satellites and centrals, typically $\lesssim$ 0.65 dex (panels $a$ & $b$ of Figs. 3 & 4). Even though the model is tuned to meet the average H i fraction of all galaxies, this large separation is possible, as the satellites only make up a minority of the galaxies (see below), and because the contribution of the more-massive centrals dominates the average. As such, Dark Sage is in good agreement with observations for the H i fraction as a function of stellar mass for the centrals (only). This is largely independent of the choice of $f_{\rm H_{2}}$ prescription, although there is slightly better agreement with observations at the low-mass end when we employ the metallicity-based prescription (cf. panels $a$ and $b$ of Fig. 3). In effect, the properties of the satellites in the model are unconstrained. All results we present concerning satellites are hence predictions (or ‘postdictions’) of the model.

38.2 per cent of galaxies in the observational sample are classed as satellites, whereas the Dark Sage Millennium sample has 29.7 per cent satellites. Some difference in these values is expected, due to the nature of group finding versus subhalo finding for the observed galaxies and simulated data, respectively. In principle, subhalo finding for simulations is a more precise means of identifying satellites (although there is certainly some variation amongst codes – see Onions et al. 2012; Knebe et al. 2013; Behroozi et al. 2015). Projection effects make it possible for true centrals to be observationally classified as satellites with a group finder, but the converse is unlikely (see, e.g., Campbell et al. 2015). We return to this important point and its potential influence on our results in Section 6.

The Dark Sage galaxies also show a flattening in their mean H i fraction at ${\rm sSFR}\!\gtrsim\!10^{-1}\,{\rm Gyr}^{-1}$ (panels $a$ and $b$ of Fig. 4), whereas the H i fractions continue to increase in the observations at high sSFR, with the separation between satellites and centrals also increasing. A much smaller fraction of the model galaxies occupies these high sSFRs versus the observations, implying galaxies in the model form their first generation of stars too early, leaving them less star-forming at $z\!=\!0$ for an equivalent stellar mass. Indeed, fig. A7 of Stevens et al. (2016) demonstrates that Dark Sage galaxies are too star-forming at high redshift. Similar behaviour is exhibited in Fig. 6, where the mean stellar mass of Dark Sage galaxies as a function of sSFR begins to flatten at \sim$$10^{-1}\,{\rm Gyr}^{-1}, whereas the observed galaxies continue towards lower stellar masses. In the observations, satellites are also found to generally be of lower stellar mass for fixed sSFR versus centrals (Fig. 6), highlighting why satellites have higher gas fractions at high sSFR (Fig. 4); at low stellar masses, the gradient of mean gas fractions in satellites is steep (Fig. 3).

To close the loop on the connection between stellar mass, H i fraction, and specific star formation rate, we compare the quiescent fraction of the full-model Dark Sage galaxies against the observations in Fig. 5$a$. Matching the observed quiescent fraction of galaxies as a function of stellar mass has proven challenging for semi-analytic models, especially for satellites (see, e.g., Font et al. 2008; Guo et al. 2011, 2013, 2016; Croton et al. 2016; Luo et al. 2016; Henriques et al. 2017). For the centrals, Dark Sage exhibits the common feature of overproducing low-mass quenched galaxies and underproducing high-mass quenched galaxies. For satellite galaxies, at least at low masses, the results are more promising. Only for $10^{10}\!\lesssim\!m_{*}/{\rm M}_{\odot}\!\lesssim\!10^{11}$ is there an overabundance of quiescent satellites. Overall, the quiescent fraction of satellites is notably improved from recent results from both semi-analytic models and hydrodynamic simulations (cf. Guo et al. 2016).

4.2 Environmental processes

In panels $c$–$e$ of Figs. 3 & 4, we present the mean H i fraction of Dark Sage galaxies after altering the prescriptions of gas stripping in subhaloes in three ways. The results in panel $c$ of these figures exclude ram-pressure stripping of cold gas entirely (but maintain all other aspects of the model, including hot-gas stripping). Conversely, those in panel $d$ exclude any consideration of hot-gas stripping (but maintain cold-gas stripping). Panel $e$ includes both forms of stripping, but uses the alternative ram-pressure prescription for hot gas (Section 2.1.3). Similarly, panels $b$–$d$ of Fig. 5 present the quiescent fraction of galaxies for these model variants. In principle, despite these stripping processes only directly affecting satellites, the hierarchical nature of galaxy formation means that centrals could be causally affected as well. We do not find any note-worthy changes to the $z\!=\!0$ population of Dark Sage central galaxies here, however.

4.2.1 The overall effect of cold-gas stripping

In Fig. 3$c$, we see that, for Dark Sage galaxies, without cold-gas stripping, there is little to separate the H i fraction of satellites and centrals as a function of stellar mass, highlighting that this process is far more effective at reducing H i fractions of satellites in the model than removal of hot gas. That said, at least for $m_{*}\!\lesssim\!10^{10.2}\,{\rm M}_{\odot}$, the satellite H i fractions are closer to the observational data than the complete, default model (cf. panels $a$ and $c$ of Fig. 3).

In Fig. 4$c$, we see that without cold-gas stripping, the model galaxies now display similar qualitative behaviour to the observed galaxies, in terms of the H i fraction as a function of sSFR. That is, at low sSFR, the Dark Sage satellites and centrals show the same H i fraction on average, but for increasing sSFR, the satellites begin to show higher H i fractions than the centrals. Although, with this run of the model, this splitting occurs at $\sim$1 dex lower sSFR than the observations. Without any fast-acting processes to remove gas from satellites, there are two reasons for satellites showing higher gas fractions than centrals in Fig. 4$c$. First, the bias of satellites having lower stellar mass than centrals for fixed sSFR, as shown in Fig. 6, is still present in this run. Second, even if we were to consider only satellites and centrals of the same stellar mass for a given sSFR, the satellite would have a suppressed supply of fresh gas due to the stripping of hot gas from its subhalo. The satellite would then have to rely on its instantaneous cold-gas reservoir more heavily for the formation of stars, and hence would need to have a greater mass of H i at an instant than an equivalent central.

Our results suggest the direct stripping of cold gas in galaxies should have a large impact on their H i content (which we reaffirm in Section 5), but that the current implementation in Dark Sage is systematically too efficient at removing H i. Potential solutions could include: (i) non-instantaneous stripping of the cold gas within the annuli, (ii) a more physically strict criterion for when cold gas stripping is allowed (rather than the current case that $m_{\rm cold}\!+\!m_{*}\!>\!m_{\rm hot}$), (iii) adding inclination effects to equation (2), or (iv) the observed H i masses include gas that the model would not consider ‘part of the galaxy’ (which is hard to define, even when one has full three-dimensional information like in a hydrodynamic simulation – see Stevens et al. 2014). However, each of these is unlikely to be satisfactory by itself, as this will have consequences for the relative H i fractions of satellites in different environments (as the strength of stripping is dependent on local environment), which, as we show in Section 5, is an area the model performs well in. We discuss this further in Section 6.

Where cold-gas stripping has a negligible effect is on the quiescent fraction of galaxies (cf. panels $a$ and $c$ of Fig. 5). As indicated by equation (2), ram-pressure stripping of cold gas is more dominant in the low-surface-density areas of satellites’ discs, towards their outskirts. The hydrogen in annuli with low $\Sigma_{\rm gas}$ will also be predominantly atomic (regardless of the specific choice of $f_{\rm H_{2}}$ prescription – see Section 2.1.1 of this paper and section 3.4 of Stevens et al. 2016). Star-forming gas is predominantly molecular and hence deeper in the potential well. As a result, star-forming gas is less susceptible to ram-pressure stripping, and thus this process has little impact on the quiescent fraction of satellites.

We note that even without cold-gas stripping, our quiescent-fraction results are a significant improvement over sage (Croton et al. 2016, from which Dark Sage was developed) for both centrals and satellites, especially at low masses. To highlight this, we have included results from sage in Fig. 5$b$ (shown in panel $b$ as sage does not include cold-gas stripping).

4.2.2 The overall effect of hot-gas stripping

Fig. 3$d$ shows that when there is no given environmental mechanism to remove the hot-gas reservoir of subhaloes, the mean H i fraction of Dark Sage galaxies is remarkably close to observations for a given stellar mass. The maintenance of that reservoir allows gas to continue to cool onto satellite galaxies in the model for longer. Based on the severity of the difference seen in the model when cold-gas stripping is removed, it is safe to say the greatest reason removing hot-gas stripping from the model has an effect on the H i fraction of galaxies is because hot gas shields the cold gas from being stripped (see Section 2.1.2). This minimises the difference between the gas fractions of satellites and centrals of the same mass, which is favourable next to observational data.

Excluding hot-gas stripping does, however, lead to satellites having systematically lower H i fractions than centrals as a function of sSFR, for all sSFR (Fig. 4$d$). This is at odds with observations. Without hot-gas stripping, satellites can continue to accrete plenty of fresh gas, which, by design of the model, settles onto the disc in an approximately exponential form. Because cold-gas stripping is still present, the outer gas that cools is quickly removed, whereas the inner gas enhances star formation. Hot-gas stripping therefore affects the molecular-to-atomic gas ratio of satellites, where, without it, galaxies would have higher sSFR at fixed H i fraction, as seen in Fig. 4$d$.

The impact of hot-gas stripping on the star-forming gas content of satellites is highlighted further in Fig. 5$d$. Without hot-gas stripping, the fraction of quiescent galaxies drops dramatically for satellites with $m_{*}\!\gtrsim 10^{10}\,{\rm M}_{\odot}$. This result would be unsurprising for any semi-analytic model (see, e.g., the difference in gradual and instantaneous hot-gas stripping from Font et al. 2008), but the point here is that even when cold-gas stripping is included (which is not the case in many models), hot-gas stripping remains the dominant process for regulating the quiescent fraction of satellites. Moreover, we see that low-mass galaxies in Dark Sage are not plagued by being over-quenched from hot-gas stripping, which has been a persistent issue with semi-analytic models in the past, as raised in Section 4.1.

Finally, in Figs. 3$e$, 4$e$, and 5$e$, we show the results of Dark Sage with the alternative, ram-pressure-like prescription for hot gas (Section 2.1.3). With this implementation, high-sSFR satellites begin to exhibit slightly higher H i fractions than their central counter-parts, qualitatively more in line with the observations, whereas the opposite is true for high-mass galaxies. There is also a minimal increase in the fraction of quiescent satellites. By and large though, the two implementations of hot-gas stripping are relatively consistent when it comes to the mean differences in gas fractions of satellites and centrals.

4.3 Evolutionary processes

In this subsection, we assess processes that affect all galaxies (regardless of whether they are centrals or satellites). We focus on the effects of (i) disc instabilities, (ii) stellar feedback and evolution, and (iii) AGN feedback, which all appear to play some role in regulating the relations between H i fraction, stellar mass, and sSFR. We present results from Dark Sage with each of these features independently switched off in panels $f$–$h$ of Figs. 3 & 4.

It should be noted that unlike for the previous runs where we altered stripping processes, simply removing feedback or instabilities from the model has a significant impact on mass functions and scaling relations, which the model is usually constrained to meet. These runs have not been recalibrated; they are, therefore, an attempt to gauge what the H i fractions of galaxies would be if each physical phenomenon did not exist, rather than attempting to recover realistic galaxies in the absence of said phenomena. This provides a more direct measure of how each process independently helps to shape the H i fractions of galaxies.

We also note that, despite using the metallicity-based prescription for determining the fraction of gas in each annulus in the form of H i and H2 (equation 1) as the new default for Dark Sage, we have assumed the pressure-based prescription for the runs without feedback or instabilities. As noted in Section 2.1.1, passive star formation from H2 will only ignite after initial metal enrichment, which happens when a galaxy first experiences a starburst from either an instability or merger. Without instabilities in the model then, only galaxies that experience a merger will have the passive H2 star formation channel open to them. This channel of star formation is the dominant mechanism by which stars are formed in low-mass galaxies. Because the majority of galaxies in the Millennium simulation merger trees do not have a recorded merger, one would end up with a non-physically large number of star-less galaxies with large gas reservoirs. A similar problem arises when stellar evolution is removed from the model, as this removes metal enrichment entirely. Both cases are resolved by switching to the mid-plane pressure prescription for $f_{\rm H_{2}}$, which is unaffected by gas metallicity. For consistency, we also use this prescription for the run without AGN feedback (but this does not make a significant difference in this case).

4.3.1 The effect of disc instabilities

As shown in Fig. 1, most star formation in low-mass Dark Sage systems comes from the passive H2 channel, whereas the high-mass systems formed most of their stars through the instability channel. By eliminating any consideration of disc instabilities from the model, star formation, that would have otherwise more directly harnessed available H i reservoirs, becomes progressively less efficient towards the high-mass end. As a result, come $z\!=\!0$, not only are there fewer galaxies of high stellar mass, but the gradient of the mean H i fraction as a function of stellar mass has become shallower, which we show in Fig. 3$f$.

With less integrated star formation having occurred in the lead up to $z\!=\!0$ without instabilities in the model, galaxies are a lot more gas rich on average. Also, because instabilities are unresolved, most of the gas remains gravitationally unstable in dense clumps, and hence there is more H2 than H i. This means that the (instantaneous) star formation rates are actually higher at $z\!=\!0$ than they were for the full model. This is seen in Fig. 4$f$, where the average H i fraction for a given sSFR has decreased.

Stevens et al. (2016) showed that instabilities regulated the stellar specific angular momentum of Dark Sage galaxies as a function of stellar mass. It is no coincidence that they also regulate the gas fractions of galaxies; recent analytic work by Obreschkow et al. (2016) and results from hydrodynamic simulations (Lagos et al. 2017) have shown that the gas fraction and specific angular momentum of a galaxy are inherently related. We will address this connection in more specific detail with Dark Sage in future work.

4.3.2 The effect of stellar feedback and evolution

Without stellar feedback or stellar evolution, an equivalent amount of gas can be used to form a greater mass of stars. This is because the instantaneous recycling approximation no longer needs to be upheld, nor does a fraction of the cold gas in a star-forming annulus need to be reserved for being reheated out of the disc. Hence, the run of Dark Sage without stellar feedback (or evolution) shows a deficit in the mean H i fraction for fixed sSFR, as in Fig. 4$g$.

Provided galaxies continue to accrete fresh gas, the lack of gas reheating raises the frequency of disc instabilities. There also tends to be greater cold-gas reservoirs in galaxies during mergers. These both not only mean that starbursts become stronger in the model, but also that black holes are able to accrete gas faster. Larger black holes lead to stronger AGN feedback, which shuts down fresh gas accretion at lower redshift. This leads to a deficit in the stellar mass function in the range $10^{10.5}\!\lesssim\!m_{*}/{\rm M}_{\odot}\!\lesssim\!10^{11.5}$ and an excess of galaxies with lower masses. Resultantly, at $z\!=\!0$, the gas fractions of galaxies at lower masses are systematically reduced, as seen in Fig. 3$g$. In this situation, only galaxies that experience many mergers can accumulate enough cold baryons to grow to larger stellar masses. Because all mergers are more gas rich, the highest-mass galaxies have larger H i fractions than in the full model. This is also seen in Fig. 3$g$.

4.3.3 The effect of active galactic nuclei

Without an AGN engine to regulate the rate at which galaxies accrete gas, all galaxies, but especially the high-mass ones, become extra gas-rich. This is seen in Fig. 3$h$. Where an AGN makes a smaller difference in Dark Sage is with how gas fraction varies with sSFR. Fig. 4$h$ shows that the most star-forming galaxies have almost unchanged H i fractions when AGNs are removed from the picture. On the other hand, the more quiescent galaxies have lower H i fractions for fixed sSFR in this case. This comes back to the extra accretion galaxies receive; much of this extra gas can quickly be consumed in star formation, leading to higher sSFRs (which are averaged over a finite time-scale) for a given instantaneous gas fraction. This is equivalent to getting lower instantaneous gas fractions for a given sSFR.

5 Satellites in various environments

In this section, we investigate how the H i fractions of satellite galaxies vary as a function of parent halo mass, which serves as a measure of environment. We consider three bins of halo mass: (i) $M_{\rm vir}\!<\!10^{12}\,{\rm M}_{\odot}$, (ii) $10^{12}\!\leq\!M_{\rm vir}/{\rm M}_{\odot}\!<\!10^{13.5}$, and (iii) $10^{13.5}\!\leq\!M_{\rm vir}/{\rm M}_{\odot}\!<\!10^{15}$. For the Dark Sage galaxies, virial masses of haloes are directly given from the Millennium simulation. As discussed in Section 3, the halo mass measurements for the observed galaxies follow the group-finding algorithm of Yang et al. (2007).

There is an argument to be made that, for a fairer comparison, one should run the same group-finding algorithm on the Dark Sage output as was used for the observations, and use those returned halo masses. We have tested that by rank ordering the stellar and virial masses of haloes, and reassigning halo masses by the ranked virial mass, any changes to our results are entirely negligible. Of course, this only mimics part of the method used in estimating the halo masses of the observed galaxies. Any further differences would come as a result of cross-contamination of centrals and satellite populations (as raised in Section 4.1). We have opted for the simpler, more digestible comparison for our main results (presenting the true model output with the observational data as they are), but we discuss the effect cross-contamination could have in Section 6.

The main results of this section are presented in Fig. 7 & 8. In Fig. 7, instead of presenting the mean H i fractions of satellite galaxies for the various halo mass bins on the $y$-axis (e.g., as in Brown et al. 2017), we present $\Delta\log_{10}(m_{\mathrm{H}\,\textsc{i}}/m_{*})$, which is the separation between the mean H i fraction curve for that halo mass bin and the mean curve for all satellite galaxies (the curves and points in Figs. 3 & 4). This means we can directly compare how environment splits the H i content of satellite galaxies in the observations and model without being concerned about any normalisation issues between the two (which have been discussed in Section 4). Similarly, in Fig. 8, rather than presenting the raw quiescent fraction of satellites in different environments, we show the difference in those fractions from the overall quiescent fraction of satellites (given in Fig. 5).

5.1 Full model versus observations

In Fig. 7, we show how halo mass (environment) affects the H i fraction of satellite galaxies in the full Dark Sage model and the observations with stellar mass (panel $a$) and sSFR (panel $e$). For $m_{*}\!\gtrsim\!10^{9.5}\,{\rm M}_{\odot}$, the observations and model are in remarkable agreement: better than 0.1 dex for the mean H i fraction in all halo mass bins. Only at lower stellar masses do we begin to see a discrepancy, where the model shows less environmental splitting than the observations. This is true regardless of whether we only consider subhaloes from the $z\!=\!0$ snapshot of Millennium with a historical maximum number of particles of $N_{\rm p,max}\!\geq\!100$, or use all subhaloes with $\geq$20 particles. In Fig. 8$a$, we show that Dark Sage also recovers the relative role environment plays on the quiescent fraction of satellites to an encouraging degree. In this case, results at low-$m_{*}$ are more in line with observations when we apply the cut of $N_{\rm p,max}\!\geq\!100$.

An important consideration is whether the distribution of satellites in the various halo mass bins is consistent between the observations and model data. While environmental processes will not inherently care about the number of galaxies they act upon, the differences in H i fractions and quiescent fraction of satellites in a given environment from the overall mean will depend on the contribution of satellites in that environment to the overall mean. As such, we list the number of satellites in each halo mass bin for the observed and model samples in Table 1.

Another important consistency check between the observations and models is the probability distribution of halo masses (within each halo mass bin). Because the group finder used to associate halo masses with observed galaxies drew from a halo mass function produced by $N$-body simulations whose cosmology is very similar to Millennium’s (cf. Springel et al. 2005; Warren et al. 2006), we should expect these to be consistent. However, because the observed sample of galaxies is limited in its volume and stellar mass range, differences in the distributions of halo masses can arise. We have tried extracting, as precisely as possible, a set of haloes from the Millennium simulation that matches the halo catalogue from the observations + group finder. While the smaller sample size introduces noise, no significant conclusions from our figures come about by making this selection (not shown here). This is in spite of the fact that the specific cases of the Coma and Hercules clusters have a large influence on the largest halo mass bin for the observational data.

5.2 The relative impact of cold-gas stripping

In Fig. 3, we demonstrated that cold-gas stripping in Dark Sage was almost entirely responsible for the systematic separation in H i fractions between satellite and central galaxies of the same stellar mass. Similarly, removing cold-gas stripping from the model almost eliminates any environmentally driven splitting in the H i fractions of satellite galaxies. This is seen not only as a function of stellar mass in Fig. 7$b$, but also as a function of sSFR in Fig. 7$f$.

The results of Dark Sage presented here are in contrast to those of the Gonzalez-Perez et al. (2014) version of the galform semi-analytic model presented alongside the same observations in Brown et al. (2017). That is, despite the model being run on the same simulation (Millennium; Springel et al. 2005) and lacking a prescription for ram-pressure stripping of cold gas entirely, galform still displayed H i fraction splitting with halo mass for satellite galaxies. There are many differences between Dark Sage and galform in terms of the specifics of how galaxies are evolved (for details on galform, see Cole et al. 2000; Benson & Bower 2010; Lacey et al. 2016). One vital aspect, which we consider here, is how the H i and H2 surface density profiles of discs are ultimately derived. We note that a modified set of merger trees (Jiang et al. 2014) for Millennium is employed for galform, but this algorithm is only expected to have a significant impact on semi-analytic galaxy properties for higher-resolution $N$-body simulations (see Merson et al. 2013; Jiang et al. 2014).

The key feature of Dark Sage is that the full two-phase angular-momentum structure of discs (and hence their radial structure also) is numerically evolved self-consistently. Similar to many other semi-analytic models (cf. Hatton et al. 2003; Somerville et al. 2008; Guo et al. 2011), galform only evolves the total angular momentum of discs (gas and stars together), and then assumes that discs always carry an exponential surface density profile with radius. Although discs are seeded as approximately exponential in Dark Sage, they develop strong cusps come $z\!=\!0$, partly caused by instability-induced radial flows of gas (Stevens et al. 2016). These cuspy gas discs are qualitatively in agreement with both observations and hydrodynamic simulations (see Bigiel & Blitz 2012; Stevens et al. 2017, respectively). A cuspy disc will have higher $R_{\rm H_{2}}(r)$ at small $r$ and relatively lower $R_{\rm H_{2}}(r)$ at large $r$ (e.g., see Equation 1). Hence, a Dark Sage galaxy will have more H2 than a galform galaxy with the same stellar mass, cold gas mass, and halo properties. Some processes that affect the H i content of galform galaxies might, therefore, affect the H2 content of Dark Sage galaxies as well/instead. We quantify the impact of this further in Section 5.4.

Within Dark Sage, cold-gas stripping only affects the star formation activity of low-$m_{*}$ satellites in the most massive haloes. This is evidenced by Fig. 8$b$, where only in this regime is there a notable change to the relative quiescent fraction when cold-gas stripping is removed. The combination of these low-mass galaxies having shallow local potential wells and the high-mass haloes in which they reside inflicting strong ram pressure means that even the molecular gas in the satellites is susceptible to being directly stripped.

We also include results from the regular sage model (Croton et al. 2016) in Fig. 8$b$. The fact that sage does not include cold-gas stripping, but is still able produce a decent result for the relative role of environment on quiescent fraction, is testament to cold-gas stripping only being of secondary importance for satellites’ quiescence. Differences seen between sage and Dark Sage here are the result of detailing the evolution of discs’ structure.

5.3 The relative impact of hot-gas stripping

As discussed in Section 4.2.2, any changes to the mean H i fraction of satellite galaxies from removing hot-gas stripping from Dark Sage were mostly caused by a reduction in the opportunity for cold-gas stripping to take place. As a result, the splitting by environment of satellites’ H i fractions as a function of stellar mass is also reduced from the full model, as shown in Fig. 7$c$. Despite the fact that removing hot gas from the model brought the difference in the mean H i fractions of satellites and centrals more in line with observations (Fig. 3$c$), the splitting of satellites’ H i fractions by halo mass has become weaker than the observations for this run. This implies that while stripping processes in Dark Sage are systematically too strong (Fig. 3$a$), they are largely of the correct relative strength in haloes of different mass.

Hot-gas stripping plays a major role in regulating the quiescent fraction of satellites, especially at higher stellar masses (Section 4.2.2). As shown in Fig. 8$c$, removing hot-gas stripping also reduces the environmental splitting of the quiescent fraction of satellite galaxies with $m_{*}\!\lesssim\!10^{10}\,{\rm M}_{\odot}$. Satellites of equivalent stellar mass (or subhalo mass at infall) experience greater depletion of their hot-gas reservoir when falling into more massive haloes. This is true regardless of whether we prescribe hot gas to be stripped at the same rate as dark matter or impose ram pressure on the hot gas (cf. panels $a$ and $d$ of Fig. 8). However, Fig. 8 also highlights that a relative environmental dependence would still be present (and consistent with observations) for the quiescent fraction of satellites without any stripping processes. This is because the hot-gas reservoirs of satellites are not allowed to grow through cosmological accretion within Dark Sage (we come back to this in Section 6). As a result, the longer a galaxy exists as a satellite, the less gas it will have available to accrete, and thus the more quiescent it will become. Satellites in larger haloes have longer orbital and merging time-scales, and therefore will have greater opportunity to consume their gas and become quiescent.

5.4 The impact of disc structure modelling

As an experiment in post-processing, we reset the specific angular momentum of the Dark Sage discs so that $j_{\rm disc}\!=\!j_{\rm halo}$, and we redistribute both gas and stars in the disc such that they follow $\Sigma_{\rm gas}(r)\!\propto\!{\rm e}^{-r}$ with a common scale radius. This meets the disc model employed in the standard version of sage (Croton et al. 2016). We then recalculate the H i and H2 masses of each galaxy according to the $f_{\rm H_{2}}(P)$ prescription (which matches the method of galform – see Lagos et al. 2011). We further recalculate their SFR based on their adjusted H2 content, assuming $\Sigma_{\rm SFR}(r)\!=\!\epsilon_{\rm SF}\Sigma_{\rm H_{2}}(r)$, maintaing $\epsilon_{\rm SF}\!=\!1.3\!\times\!10^{-4}\,{\rm Myr}^{-1}$. This allows us to estimate what differences we would find regarding the relative influence of environment on the H i fractions of satellites if Dark Sage did not include a detailed treatment of the angular momentum of discs and did not allow them to be cuspy. We present results for this for a run of the model that uses the ram-pressure prescription for hot-gas stripping, and excludes cold-gas stripping entirely. These choices give us the most directly comparable results to those from galform version GP14+GRP (Gonzalez-Perez et al. 2014) published in figure 7 of Brown et al. (2017). Our results are given in panels $d$ and $h$ of Fig. 7.

Fig. 7$d$ shows that an environmental splitting of H i fractions as a function of stellar mass caused entirely by hot-gas stripping can be recovered for a more rudimentary model of galaxy discs. The strength of this splitting is still significantly less than what is observed in real galaxies. It is also less than half the strength seen in galform (cf. Brown et al. 2017). Meanwhile, the splitting with sSFR has become flat, now too weak at low sSFR and too strong at high sSFR versus observations (Fig. 7$h$). Of course, because the gas redistribution is done in post-processing, the evolution of the galaxies is not affected. It is an obvious statement that if aspects of galaxy evolution within one model were changed to match that of another model, the two models would end up giving very similar results. Such an exercise is hence not particularly informative. What the results of panels $d$ and $h$ of Fig. 7 do tell us is that the method for determining $f_{\rm H_{2}}$ only has a limited potential contribution to the dominance of cold-gas stripping in driving the environmental dependence of satellites’ H i fractions in Dark Sage. In other words, these results still favour ram-pressure stripping of cold gas to be the main mechanism for curbing the H i content of satellite galaxies.

6 Caveats, discussion and conclusions

We have shown that the relative role that environment is observed to have on the H i fractions of satellites (Brown et al. 2017) can be recovered within a semi-analytic model of galaxy formation with a self-consistent evolution of disc structure, namely Dark Sage (Stevens et al. 2016). However, this model also predicts the mean H i fractions of satellites as a whole to be too low. In the current design of the model, these two results are inherently tied together. One possible explanation is that the model is missing a physical process that would systematically (and exclusively) raise the gas content of satellites, irrespective of their environment.

One process hypothesised to help explain why low-mass galaxies in denser environments are observed to have brighter X-ray haloes (Mulchaey & Jeltema 2010) is ‘confinement pressure’. The idea here is that the temperature differential between the hot gas of a satellite and the hotter intracluster medium (ICM) creates an inward pressure, thus helping that satellite to retain its hot gas and dampening the effect of hot-gas stripping. However, it has been shown with both a simple analytic model and the results of hydrodynamic simulations that confinement pressure only has a strong dampening effect on hot-gas stripping once a satellite reaches pericentre, at which point most of its hot gas has been removed already (Bahé et al. 2012). The satellite galaxies in the Dark Sage results without hot-gas stripping presented in this paper were also prescribed to maintain any reheated gas from feedback within their own hot reservoir (Section 2.1.3). In effect, this run of the model shows an extreme for what confinement pressure could achieve.

More recently, using a cosmological hydrodynamic simulation focussed on a galaxy cluster, Quilis, Planelles & Ricciardelli (2017) have suggested that the hot gas of satellite galaxies can be replenished by accretion from the ICM. This helps to counter-act the stripping of hot gas, where cold-gas stripping is shown by Quilis et al. to have a greater net impact on these galaxies. Complementary results from the EAGLE simulations have explicitly shown that the gas accretion rates through a fixed 30-kpc aperture around satellite galaxies still show a strong dependence on environment (van de Voort et al. 2017). Note, though, that gas accretion rates through a fixed aperture can be wildly different to accretion rates onto a (sub)halo (Stevens et al. 2014).

Unlike central galaxies, which accrete hot gas cosmologically, satellites cannot acquire hot gas from an external source within the current Dark Sage model. This feature remains a shared standard for semi-analytic models in general (see the latest versions of other models, e.g. Gargiulo et al. 2015; Henriques et al. 2015; Lacey et al. 2016). If subhaloes were able to accrete hot gas from the ICM within a semi-analytic framework, this should systematically raise H i fractions and lower quiescent fractions. Both of these outcomes would be favourable for Dark Sage, potentially offering a solution (at least partially) for simultaneously recovering the absolute and relative impacts of environmental stripping processes.

Another possible explanation for the discrepancy in the mean H i fractions of satellites in the observations and model lies with the group-finding algorithm applied to the observational data. Campbell et al. (2015) studied how well several group-finding algorithms performed in correctly identifying satellites and centrals by applying them to a mock catalogue from simulated data. The authors found that the Yang et al. (2005, 2007) technique returns a population of satellites with a purity of $\sim\!0.6$–$0.7$ for all haloes with $M_{\rm vir}\!>\!10^{12}\,{\rm M}_{\odot}$ (see fig. 6 of Campbell et al. 2015). This means it is possible that 30–40 per cent of galaxies classed as satellites in our observational sample (Section 3) might actually be centrals. This could more than account for the difference in the fraction of satellites in the observed and model samples (38.2 and 29.7 per cent, respectively – Section 4.1). While repeating the efforts of Campbell et al. for Dark Sage is beyond the scope of this paper, we can easily test what the model results would be after contaminating the satellite population with misidentified centrals by taking a weighted mean of the true central and satellite galaxy properties. Using respective weights of 0.6 and 0.4 for satellites and centrals gives a mean H i fraction as a function of stellar mass that is consistent with a satellite population with a purity of 0.6. We have plotted this for Dark Sage in Fig. 9$a$, where we have also included central-galaxy results after applying a purity of 0.8. The contaminated satellite population meets the relation for observed satellites to a remarkable degree. While this could offer an easy out for the too-gas-poor nature of Dark Sage satellites, cross contamination of satellites and centrals leads to a reduction in the relative gas fractions of satellites in bins of halo masses, as seen in Fig. 9$c$. This would, therefore, simply trade one good result for a different one. Either way, the treatment of satellites in the semi-analytic model still requires further attention before it can reproduce all the observational results we have shown.

Our results with Dark Sage suggest the effects of hot- and cold-gas stripping on galaxies are almost entirely separable. Cold-gas stripping removes the highest-$j$ gas in the disc, which is at low surface density and is dominated in mass by H i. This generates a difference in the mean H i fractions between central and satellite galaxies at fixed stellar mass (Fig. 3), and nicely recovers observed trends in the relative H i fractions of satellites in different environments (probed by halo mass – Fig. 7). Molecular gas, which exists almost exclusively at the centres of these galaxies’ discs, is almost unaffected by this process. Recent CO observations of 3 satellite galaxies in the Virgo cluster by Lee et al. (2017) show that, while the morphology of molecular gas can be disturbed by ram pressure, there is no evidence that ram-pressure stripping affects the total molecular gas content. For Dark Sage, only low-$m_{*}$ galaxies in high-mass haloes have their molecular gas and star formation activity directly impacted by ram-pressure stripping (Fig. 8).

On the other hand, hot-gas stripping shuts down galaxy accretion, thereby affecting the evolution of gas across the entire disc. The replenishment of molecular gas in the galaxy is then suppressed, causing the satellites to become quiescent after consuming their available H2 in star formation. This is a major contributor behind the observed difference in quiescent fractions between satellites and centrals at fixed stellar mass (Fig. 5). Based on the relative colours, sizes, and masses of centrals and satellites, observational evidence supports a picture where hot-gas stripping is the responsible mechanism for satellite quenching (van den Bosch et al. 2008). While showing improvement relative to earlier results from semi-analytic models, Dark Sage still overproduces the fraction of quiescent satellites in the range $10^{10}\!\lesssim\!m_{*}/{\rm M}_{\odot}\!\lesssim\!10^{11}$. Cross-contamination of centrals and satellites in the observational data from the group finder could also be playing a role here, mind; after cross-contaminating the Dark Sage galaxies, we find the excess of quiescent satellites at mid masses goes away (Fig. 9$b$). Again though, this reduces the relative strength of environment on quiescence in the model (Fig. 9$d$).

We note that Stevens (2016) has shown that the stellar half-mass radii of disc-dominated Dark Sage galaxies are systematically lower than the half-light radii of galaxies in the GAMA††††††Galaxy And Mass Assembly survey by $\lesssim\!0.1$ dex (cf. Lange et al. 2016). Galaxies that never suffered an instability in their existence, however, meet the observed relation quite precisely. Instabilities drive mass inwards and angular momentum outwards in Dark Sage discs (Stevens et al. 2016). Because the model lacks a detailed consideration of radial dispersion support at the centres of discs, their surface density profiles grow an exaggerated cusp. This leads to smaller stellar sizes and higher molecular-to-atomic hydrogen mass ratios, especially at the discs’ centres. Because the integrated H i content of the model galaxies is calibrated to meet observations (Brown et al. 2015), the average radius of an H i element will be greater than in real galaxies (see fig. 4 of Stevens et al. 2016). This does not prevent the nominal H i radii of the galaxies (where $\Sigma_{\mathrm{H}\,\LARGE\textsc{i}}\!=\!1\,{\rm M}_{\odot}\,{\rm pc}^{-2}$) from agreeing with observations for fixed H i mass, however (see Lutz et al. in preparation). We compare the molecular-to-atomic ratio as a function of stellar mass for Dark Sage galaxies against data presented by Boselli et al. (2014b) from the Herschel Reference Survey (Boselli et al. 2010, 2014a) in Fig. 10. Indeed, the ratios of the model galaxies are systematically higher than the observational data, regardless of whether the conversion between CO luminosity and H2 mass is taken as variable or a constant. While the Boselli et al. data are not volume-limited and only account for 74 late-type galaxies with detections in both H i and CO, they are designed to be representative of average galaxies in the local Universe. With that in mind, we encourage that while the information in Fig. 10 is important to consider in interpreting our results, it is unlikely to be what drives them. Future work in improving how the centres of Dark Sage discs are modelled, as well as larger H2 surveys, will help to clarify our results.

We should be optimistic about soon converging on a detailed description of how environmental and secular processes affect the gas content and star formation activity of galaxies in a cohesive manner that explains our observations. Already, data from MUSE‡‡‡‡‡‡Multi Unit Spectroscopic Explorer (Bacon et al. 2010) have been used to show the optical effects of ram pressure in exquisite detail for a small sample of galaxies (e.g. Fumagalli et al. 2014; Poggianti et al. 2017), which provide clues for this picture. Upcoming surveys such as the ASKAP******Australian Square Kilometre Array Pathfinder H i All-Sky Survey (also known as Wallaby; Koribalski 2012) will complement these with large-number statistics of the H i content, structure, and kinematics of environmentally influenced satellite galaxies. On the theory side, where detailed cosmological hydrodynamic simulations have been limited in their volume ($\sim\!10^{6}\,{\rm Mpc}^{3}$ – e.g. Dubois et al. 2014; Vogelsberger et al. 2014; Khandai et al. 2015; Schaye et al. 2015), follow-up simulations with the same well-tested physical prescriptions are now targeting Local Group analogues and massive clusters, providing more opportunity to learn about the impact of stripping processes (Sawala et al. 2016; Bahé et al. 2017, respectively). Based on the level of agreement with observations we have presented in this paper, and without any concerns regarding volume or computational efficiency, semi-analytic models are bound to continue to play an important role in this quest as well.

Acknowledgements

We thank Claudia Lagos for comments on this paper, as well as Luca Cortese, Violeta Gonzalez-Perez, and Manodeep Sinha for helpful discussion regarding this work. ARHS also thanks Darren Croton for continued academic support during the writing of this manuscript.

The Dark Sage codebase is publicly available at https://github.com/arhstevens/DarkSage. This was developed from the sage codebase, which is also publicly available at https://github.com/darrencroton/sage. Catalogues from both models are available at https://tao.asvo.org.au/tao/.