Semi-analytic galaxies – I. Synthesis of environmental and star-forming regulation mechanisms

ABSTRACT

We present results from the semi-analytic model of galaxy formation sagapplied on the MultiDark simulation MDPL2. sag features an updated supernova (SN) feedback scheme and a robust modelling of the environmental effects on satellite galaxies. This incorporates a gradual starvation of the hot gas halo driven by the action of ram pressure stripping (RPS), that can affect the cold gas disc, and tidal stripping (TS), which can act on all baryonic components. Galaxy orbits of orphan satellites are integrated providing adequate positions and velocities for the estimation of RPS and TS. The star formation history and stellar mass assembly of galaxies are sensitive to the redshift dependence implemented in the SN feedback model. We discuss a variant of our model that allows to reconcile the predicted star formation rate density at z ≳ 3 with the observed one, at the expense of an excess in the faint end of the stellar mass function at z= 2. The fractions of passive galaxies as a function of stellar mass, halo mass, and the halo-centric distances are consistent with observational measurements. The model also reproduces the evolution of the main sequence of star-forming central and satellite galaxies. The similarity between them is a result of the gradual starvation of the hot gas halo suffered by satellites, in which RPS plays a dominant role. RPS of the cold gas does not affect the fraction of quenched satellites but it contributes to reach the right atomic hydrogen gas content for more massive satellites (M_⋆ ≳ 10¹⁰ M_⊙).

methods: numerical, catalogues, galaxies: clusters: general, galaxies: evolution, galaxies: formation

1 INTRODUCTION

The properties of galaxies we observe at the present epoch are the result of a variety of physical mechanisms acting throughout their assembly history. Galaxy evolution leads to a bimodal colour distribution that features a red and a blue branch in the colour–stellar mass diagram (e.g. Baldry et al. 2004). Interestingly, each colour branch is also associated with a different galaxy morphology: early-type galaxies are located in the well-defined red sequence, whereas late-type galaxies populate the blue cloud (Conselice 2006). The fraction of galaxies in the red sequence increases as a function of local environmental density (Baldry et al. 2006), especially at high stellar masses. However, the environmental dependence of the morphology fractions is milder than that of colour fractions for low masses (Bamford et al. 2009). Outliers of these two broad colour branches include galaxies in a so-called intermediate green valley. These galaxies are interpreted as experiencing a decrease of their star formation rate (SFR) (Schawinski et al. 2014; Bait, Barway & Wadadekar 2017; Powell et al. 2017; Coenda, Martínez & Muriel 2018). Star formation (SF) quenching is higher for more massive galaxies, and increases at fixed stellar mass with environmental density (given by either local density or halo mass) as becomes evident from the strong correlations of the fractions of quiescent galaxies with stellar mass, environment, and halo-centric radius, found both locally (Wetzel, Tinker & Conroy 2012) and at high redshift (Peng et al. 2010; Muzzin et al. 2012; Lin et al. 2014; Jian et al. 2017; Kawinwanichakij et al. 2017). Besides, the way in which red and quiescent fractions vary with environment and stellar mass is different for central and satellite galaxies (Weinmann et al. 2006; Wetzel et al. 2012; Kovač et al. 2014; Darvish et al. 2017; Smethurst et al. 2017; Wang et al. 2018).

Central and satellite galaxies of a given stellar mass are subject to very different processes. Central galaxies are believed to reside at the centres of their host dark matter (DM) haloes. They can receive new gas via cooling flows, and they cannibalize satellite galaxies whose orbits in the group potential decay due to dynamical friction. Dry mergers (both minor and major) can explain the shallower slope of the massive end of the red sequence (Jiménez et al. 2011). On the other hand, major mergers contribute to produce massive, passive galaxies (van der Wel et al. 2009), doubling the mass of the massive brightest cluster galaxies since z ∼ 1 (e.g. Shankar et al. 2015). While orbiting in the main group, satellite galaxies are also affected by galaxy encounters that may result in morphological transformations (Kannan et al. 2015). However, the key processes that are believed to quench the SF in satellite galaxies are triggered by gas disruption due to tidal stripping (TS, Merritt 1983) and ram pressure stripping (RPS, Gunn & Gott 1972; Abadi, Moore & Bower 1999).

During TS, the material is pulled from a galaxy by the global tidal field of the host halo or a larger neighbour galaxy, the former being more relevant than the latter in low-velocity close encounters taking place in galaxy groups (Villalobos, De Lucia & Murante 2014). Tidal shocks at the pericentres of the orbit of satellite galaxies induce dynamical instabilities and impulsive tidal heating of the stellar distribution that are particularly relevant for dwarf galaxies in environments similar to our Local Group (Kazantzidis et al. 2011). Besides, the impact of TS (from negligible effect to complete galaxy disruption) also depends on the satellite morphology (Chang, Macciò & Kang 2013).

RPS is a result of the ram-pressure (RP) exerted by the hot diffuse gas of a group or cluster of galaxies on satellites moving through it at velocities that could be close to supersonic. This physical process may account for the gradual removal of both the hot gas halo and the cold gas disc (Boselli & Gavazzi 2006; Jaffé et al. 2013; Poggianti et al. 2017), giving place to ‘jellyfish’ galaxies in cases of extreme RPS (Bellhouse et al. 2017) where the trailing stripped tails show signs of shock heating and gas compression. Hydrodynamical simulations show that, on one hand, RP contributes to SF quenching through gas stripping and, on the other, can also produce a temporary and moderate enhancement of SF (more likely to occur in Milky Way-type disc galaxies) as a result of gas compression that takes place at pericentre passage (Bekki 2014; Yozin & Bekki 2015; Steinhauser, Schindler & Springel 2016; Ruggiero & Lima Neto 2017).

The instantaneous complete removal of the hot diffuse gas halo of a galaxy after its infall into a larger halo, that takes place regardless of which physical mechanism may actually be responsible for it, has been standard ingredient in semi-analytic models (SAMs). This effect is generally referred to as strangulation (Larson, Tinsley & Caldwell 1980; Balogh, Navarro & Morris 2000). Although this simple prescription for halo stripping can account for the presence of passive galaxies in clusters at high redshift, it results in fractions of red satellite galaxies which are higher than those observed in groups of galaxies (e.g. Baldry et al. 2006; Weinmann et al. 2006; Kimm et al. 2009); this excess of red satellites constitutes the satellite overquenching problem. Based on these results, there is now general consensus that the instantaneous stripping of the hot gas of satellites due to shock-heating is a crude approximation of the process. This is additionally supported by the observational evidence that large fractions of near-IR-bright, early-type galaxies in groups (Jeltema, Binder & Mulchaey 2008) and also in clusters (Sun et al. 2007; Wagner, McDonald & Courteau 2017) have extended X-ray emission, indicating that they retain significant hot gas haloes even in these dense environments. Studies based on smoothed particle hydrodynamics simulations have found that the hot gas haloes of satellites are not stripped instantly (McCarthy et al. 2008), although RPS of halo gas is much more effective than that of disc gas (Bekki 2009). Moreover, the effect of RPS on the hot corona seems to be even milder than expected because of turbulence in the interface with the intracluster medium (ICM) and fuelling of gas from this component to the satellite hot gas halo (Quilis, Planelles & Ricciardelli 2017).

This mounting evidence motivated several authors to include prescriptions for a gradual removal of hot halo gas in SAMs, with the inherent limitations of this numerical technique. Some of them only consider the effect of RPS (e.g. Font et al. 2008; Kang & van den Bosch 2008; Weinmann et al. 2010; Gonzalez-Perez et al. 2014) while others also take into account the influence of TS (Guo et al. 2011; Kimm, Yi & Khochfar 2011; Henriques et al. 2017; Stevens & Brown 2017). If hot gas removal proceeds gradually by the action of RPS and/or TS, the time-scale of SF quenching becomes longer and the discrepancies in the fraction of red/quiescent galaxies between models and observations are alleviated although some disagreements remain. The prescriptions adopted are generally based on the analytic formulation of McCarthy et al. (2008). The effect of RPS on the cold gas of galaxy discs has been implemented only in few SAMs (Okamoto & Nagashima 2003; Lanzoni et al. 2005; Brüggen & De Lucia 2008; Tecce et al. 2010; Luo et al. 2016; Stevens & Brown 2017). The combination of suppressed cooling flows due to strangulation and RPS of cold disc gas as in Tecce et al. (2010) worsens the problem of excess of red satellites. The presence of hot gas halo during gradual removal acts as a shield that helps to prevent a too efficient action of RPS on the cold gas phase. Although this is taken into account by Stevens & Brown (2017), the condition used seems to still produce larger amounts of removed cold gas than expected.

The star formation history of galaxies and, consequently, the fractions of star-forming and passive galaxies are strongly affected by supernova (SN) feedback. Hirschmann, De Lucia & Fontanot (2016) tested the impact of different schemes for stellar feedback and gas recycling using the semi-analytic model gaea. They find that a model characterized by a strong ejection coupled to a mass-loading that is explicitly dependent on redshift is reasonably successful in reproducing the observed ‘antihierarchical’ trends in galaxy mass assembly, with a delayed metal enrichment and a realistic present-day stellar and gaseous metallicity. They adopt a redshift dependence of the mass-loading factor (ratio between the gas outflow rate from galaxies and the SFR) similar to the one provided by Muratov et al. (2015), who analyse the galaxy-scale gaseous outflows in the cosmological hydrodynamical zoom FIRE (Feedback in Realistic Environments) simulations.

In this work, we present an updated version of the semi-analytic model sag (acronym for semi-analytic galaxies, Cora 2006; Lagos, Cora & Padilla 2008; Orsi et al. 2014; Gargiulo et al. 2015; Muñoz Arancibia et al. 2015), which includes an improved treatment of environmental effects through the implementation of the gradual removal of hot gas in satellites by RPS and TS, allowing also the action of these processes on the cold gas disc, under certain conditions; the modelling of supernova (SN) feedback is also modified. This version of sag was used to generate one of the galaxy catalogues of the MultiDark Galaxies project (Knebe et al. 2018a), which is based on the Planck cosmology 1 h⁻¹ Gpc MultiDark simulation MDPL2 (Klypin et al. 2016). This catalogue is publicly available¹ in the cosmosim data base² together with those generated by the semi-analytic models sage (Croton et al. 2016) and galacticus (Benson 2012).

This paper presents a detailed description of the modifications introduced, emphasizing the advantages with respect to previous works. Details of the MDPL2 simulation and a brief summary of the basic aspects of sag are given in Section 2. Improvements introduced regarding environmental effects and SN feedback are explained in Sections 3 and 4, respectively. The method applied to tune the free parameters of the model and the comparison of model results with observational constraints imposed to calibrate the model are presented in Section 5. Section 6 includes an analysis of general galaxy properties predicted by sag (complementing the discussion already presented in Knebe et al. 2018a), and of those related with SF quenching: dependence of the fraction of quenched galaxies as a function of stellar mass, halo mass, and halo-centric distance. We show that the relations involving the fraction of passive galaxies are better reproduced by modifying a parameter introduced in the modelling of SN feedback which has impact on the SF history of galaxies. For this variant of the model, we analyse the role of environmental processes in Section 7, focusing on the stripped mass and atomic gas content. In Section 8, we compare details of our model and its main results with other SAMs. In Section 9, we present the summary and conclusions of this work. Appendix  A describes the estimation of the hot gas stripping radius dueto RP.

2 GALAXY FORMATION MODEL

Our semi-analytic model of galaxy formation and evolution sag uses DM haloes extracted from a cosmological DM simulation and their corresponding merger trees as the basic inputs to construct the galaxy population. We use the MultiDark simulation MDPL2, which is part of the CosmoSim data base.

2.1 MDPL2 simulation

MDPL2 simulation follows the evolution of 3840³ particles within a box of side-length 1 h⁻¹ Gpc, with a mass resolution m_p = 1.5 × 10⁹ h⁻¹ M_⊙ per DM particle. This simulation is an analogue of the MDPL simulation described in Klypin et al. (2016), but ran with a different realization of the initial conditions. It is consistent with a flat Lambda cold dark matter model characterized by Planck cosmological parameters: Ω_m = 0.307, |$\Omega _\Lambda$| = 0.693, Ω_B = 0.048, n_s= 0.96, and H₀ = 100 h⁻¹ km s⁻¹ Mpc⁻¹, where h= 0.678 (Planck Collaboration XVI 2014).

DM haloes have been identified with the rockstar halo finder (Behroozi, Wechsler & Wu 2013a), and merger trees were constructed with consistenttrees (Behroozi et al. 2013b). Overdensities with at least N_min = 20 DM particles were considered in the detection by the halo finder. Each halo is characterized by measuring the physical properties defined by the particle distribution, assuming spherical overdensity approximation. The virial mass is defined as the mass enclosed by a sphere of radius r_vir, so that the mean density reaches a constant factor Δ = 200 times the critical density of the Universe ρ_c, i.e.

\begin{equation*} M_{\rm vir}(\lt r_{\rm\!vir}) = \Delta \rho _{\rm c} \frac{4 \pi }{3} r_{\rm vir}^{3}. \end{equation*}

(1)

Moreover, the virial velocity of each halo is defined in terms of these properties as |$V_{\rm vir}= \sqrt{G M_{\rm vir}/r_{\rm vir}}$|⁠, where G is the gravitational constant. The detected DM haloes can be over the background density or lie within another DM haloes. To differentiate them, henceforth the former will be referred to as main host haloes, whereas the latter as subhaloes. The calculation of the physical properties of subhaloes considers only the bound particles of the substructure identified by the halo finder.

2.2 Semi-analytic model of galaxy formation SAG

The version of our semi-analytic model of galaxy formation and evolution sag presented here is a further development of the model described by Cora (2006), which is based on Springel et al. (2001), and later improved in Lagos et al. (2008), Tecce et al. (2010), Orsi et al. (2014), Muñoz Arancibia et al. (2015), and Gargiulo et al. (2015).

sag includes the effects of radiative cooling of hot gas, star formation, feedback from SN explosions, chemical enrichment with a scheme that tracks several chemical elements contributed by different sources (stellar winds and supernovae Type Ia and II) taking into account the lifetime of progenitors (Cora 2006), growth of supermassive black holes (BHs) in galaxy centres and the consequent active galactic nucleus (AGN) feedback, and starbursts triggered by disc instabilities or galaxy mergers (Lagos et al. 2008). These starbursts contribute to the formation of a bulge component, whose sizes are estimated as described in Muñoz Arancibia et al. (2015); the cold gas that has been transferred to the bulge is gradually consumed, thus starbursts are characterized by a time-scale (Gargiulo et al. 2015) instead of being instantaneous. Tecce et al. (2010) added the effect of RPS on the cold gas disc of satellite galaxies implementing the criterion from Gunn & Gott (1972).

Two important new features have been added to sag as described in this work. One is related to the inclusion of an explicit redshift dependence in the estimation of the reheated and ejected mass as a result of SN feedback, guided by one of the ejective feedback schemes implemented in the gaea semi-analytic model (Hirschmann et al. 2016). This redshift dependence becomes crucial to reproduce observed galaxy properties at high redshift, as we will show in Sections 5.2 and 6. The second major improvement replaces the strangulation scheme that removes the hot gas from satellites instantaneously (standard practice in SAMs until recently) by a gradual removal of the halo gas by the action of TS and/or RPS (gradual starvation³).

The features that the current version of sag inherits from its predecessors are described in Gargiulo et al. (2015), except for those related to the RPS of cold gas (Tecce et al. 2010), which will be described in Section 3 along with the other environmental processes introduced in this work. The new SN feedback scheme is described in Section 4. Furthermore, the model of radio mode AGN feedback, in which BHs generate jets and bubbles as they accrete gas from hot atmospheres suppressing gas cooling, has been replaced by that described in Henriques et al. (2015, see their equation S24 in the supplementary material). This means that the dependence of the BH accretion rate on the fraction of hot gas and the virial velocity of the halo (introduced by Lagos et al. 2008) are substituted by a dependence on the hot gas mass. The aim of this change is to make AGN feedback more efficient at late times.

The semi-analytic model assigns one galaxy to each new detected halo in the simulation, and it follows the halo merger trees to compute the evolution of the galaxy properties. Each considered system of haloes is constituted by only one group/cluster central galaxy, the one associated with the main host halo, so that the other galaxies act as satellites. When two haloes merge, the smaller one loses mass to tides as it orbits within the larger structure, until the satellite is no longer identified by the halo finder. We assume that the galaxy contained within this disappeared subhalo survives until it eventually merges with the central galaxy of its host halo. During this temporary stage, these galaxies are called orphan satellite galaxies.

3 ENVIRONMENTAL PROCESSES

As satellite galaxies orbit within their main host haloes, DM and the baryonic components in their own subhaloes are affected by numerous environmental processes. We model mass stripping produced by tidal forces and RP. TS will act on every component of the galaxy (stars, gas, and DM) whereas RPS is a hydrodynamical process which will only affect the gas. Thus, a two-stage model for gas removal is now included in sag. TS and RPS both act first on the hot gas halo: whichever effect is stronger determines the amount of gas stripped. Once a significant fraction of the hot halo is removed, RPS can start affecting the cold gas disc, following Tecce et al. (2010). The combination of the RPS of the cold gas phase with the gradual starvation scheme applied to the hot gas constitutes an improvement of our model with respect to previous works that ignore the RPS of the cold gas (e.g. Guo et al. 2011; Kimm et al. 2011; Gonzalez-Perez et al. 2014). Recently, this effect has been taken into account by Luo et al. (2016) and Stevens & Brown (2017); while the former authors do not consider the hot gas as a shield that regulates the action of RP on the cold gas, the latter take into account this possibility using a different criterion than the one adopted in sag.

3.1 Estimation of RP

There have been numerous studies using hydrodynamic simulations focused on the effects of RPS on cold gas discs (e.g. Abadi et al. 1999; Roediger & Brüggen 2007; Tonnesen, Bryan & van Gorkom 2007) and hot gas haloes (McCarthy et al. 2008; Quilis et al. 2017). RP is defined as the product of the ICM density at the location of the satellite galaxy and the square of its velocity relative to the ICM. The importance of the local variations of the ambient gas makes analytical calculations of the RP effect slightly inaccurate, particularly for the resulting size of the gas disc. Tonnesen & Bryan (2008) show that these variations can produce differences in the pressure of more than one order of magnitude at fixed cluster-centric distance. As demonstrated by Tecce et al. (2010) from the analysis of adiabatic hydrodynamic simulations of galaxy clusters, the values of RP present a radial profile mainly determined by the dependence on cluster-centric distance of the ICM⁴ density. Values of RP of the order of 10⁻¹², 10⁻¹¹, and 10⁻¹⁰ h² dyn cm⁻² are considered weak, medium, and strong, respectively (Roediger & Brüggen 2006). At z = 0, mean RP values range from weak to medium for clusters masses of ≈10¹⁴ and ≈10¹⁵ M_⊙, respectively. At z ≈ 2, most of satellite galaxies within progenitors of currently massive clusters are already experiencing medium-level RP, and at z ≲ 0.5 a significant fraction of galaxies (≈20 per cent) suffer strong levels of RP (see fig. 5 of Tecce et al. 2010). This evolution depends mainly on the build-up of the ICM density over time.

We take into account all these effects through the implementation in sag of fitting formulae that estimate the RP experienced by galaxies in haloes of different mass as a function of halo-centric distance and redshift, which are derived from a procedure similar to that described by Tecce, Cora & Tissera (2011). The new fit uses alternative formulae that capture more adequately the RP behaviour, as described in Vega-Martínez et al. (in preparation) using the self-consistent information provided by Tecce et al. (2010) from the analysis of hydrodynamic simulations of galaxy clusters. This procedure, called the ‘gas-particles’ method in Tecce et al. (2010), automatically takes into account any local variation of the density and/or the velocity field which may result from the dynamics of the gas in the hydrodynamic simulation. The fitting formulae generated from these data allow to include RP effects into SAMs applied on dark matter-only simulations which lack gas physics, as the one used in this work. This way of estimating RP values constitutes the main advantage of our treatment of RPS.

3.2 Orbits of orphan galaxies

For those satellites with identified host subhaloes, the effects of dynamical friction and TS are provided self-consistently by tracking the subhalo evolution in the base simulation. This approach breaks down at the subhalo resolution limit, when the galaxy host by a subhalo that is no longer identified becomes an orphan satellite. In previous versions of sag, positions and velocities of orphans were assigned using the current position and velocity of the most-bound DM particle of the galaxy’s host subhalo at the time it was last identified (Lagos et al. 2008; Tecce et al. 2010), or assuming a circular orbit with decaying radial distance estimated from the dynamical friction and initial halo-centric distance and velocity determined by the virial radius and velocity of the host (sub)halo, respectively (Gargiulo et al. 2015). In any case, the merging criterion applied to orphans compares the time elapsed since the time of infall with a dynamical friction time-scale; when the former is larger than the latter, an orphan galaxy is considered merged with a central galaxy regardless of the relative distance between them.

In the current version of sag, the phase space of orphans is obtained from the integration of the orbits of subhaloes that will give place to an orphan galaxy once sag is applied to the merger trees extracted from the underlying simulation, as described in detail by Vega-Martínez et al. (in preparation). When a subhalo is not longer tracked, the last known position, velocity, and virial mass are taken as initial conditions to integrate its orbit numerically. The orbital evolution is consistent with the potential well of the host halo, and takes into account mass-loss by TS and dynamical friction effects, following some aspects of the works by Gan et al. (2010) and Kimm et al. (2011). The integration is carried out until a merger event occurs, which is defined when the halo-centric distance becomes smaller than 10 per cent of the virial radius of the host halo, taken as an estimation of the radius of a central galaxy, since in this stage previous to the application of sag we do not have information on the properties of galaxies that will populate these haloes. This merging criterion allows to reach general good agreement with the redshift evolution of the merger rate of galaxies obtained by considering the merger time-scale given by Jiang et al. (2008), which is inferred from a high-resolution cosmological hydro/N-body simulation with star formation.

3.3 Assumptions for the action of RPS and TS

We assume that the time-scale for RPS and TS of any of the baryonic components is given by t_dyn,orb= 2  |$\pi$|ω⁻¹, where ω is the angular velocity of the satellite, according to Zentner et al. (2005). In each interval ΔT between simulation snapshots, the RPS and/or TS effect removes a fraction ΔT/t_dyn,orb of gas and stars outside the corresponding stripping radius, r_s, if this fraction is less than unity; otherwise, the whole estimated stripped mass, M_strip, is removed. In the following, we describe the way in which the stripping radius r_s is derived for hot gas, cold gas disc, disc stars, and bulge stars, referring to this quantity as r_s,hot, R_s,cold, R_s,disc, and r_s,bulge, respectively.

3.4 Gradual stripping of the hot gas haloes of satellites

The new version of sag presented in this paper allows galaxies to keep their hot gas haloes when they become satellites. These haloes are then gradually stripped by the action of RPS and/or TS, and as long as they survive they can replenish the satellite’s cold gas via gas cooling. This is the scheme we call gradual starvation.

In this new scenario, the mass of hot gas M_hot available for cooling to the central galaxy of each main host halo is calculated at the beginning of each simulation snapshot as

\begin{eqnarray*} M_\text{hot} &=& f_{\rm b} M_\text{vir} - M_{\star ,\text{cen}} - M_\text{cold,cen} - M_\text{BH,cen} \nonumber\\ &-& \sum _{i=1}^{N_\text{sat}} \left( M_{\star ,i} + M_{\text{cold,}i} + M_{\text{BH,}i} + M_{\text{hot,}i} \right), \end{eqnarray*}

(2)

where f_b= 0.1569 is the universal baryon fraction (Planck Collaboration XVI 2014) and M_⋆, M_cold, and M_BH are, respectively, the total masses of stars, of cold gas, and of the central BH of each of the galaxies contained within a given main host halo. Central and satellite galaxies are considered separately in order to discount the total mass in the remaining hot gas haloes of the N_sat satellite galaxies within the same main host halo.

Our prescription for RPS of the hot gas halo is based on the model from Font et al. (2008), which uses the criterion for RPS for a spherical distribution of gas determined by McCarthy et al. (2008) from the results of hydrodynamic simulations. The gas beyond a satellite-centric radius r_sat will be removed if the value P_ram of RP, which is given by the fitting formulae described in Section 3.1, meets the condition

\begin{equation*} P_{\rm ram} \gt \alpha _\text{RP} \frac{G M_\text{sat}(r_\text{sat}) \rho _\text{hot} (r_\text{sat})}{r_\text{sat}}, \end{equation*}

(3)

where ρ_hot is the density of the satellite’s hot gas halo, and α_RP is a geometrical constant of order unity chosen to match the results of hydrodynamic simulations. McCarthy et al. (2008) find that their simulations support α_RP = 2. They also test the effect of values as large as α_RP = 10, which produce milder effects. Since they do not consider gas cooling that can also reduce the hot gas content, we adopt an intermediate value of α_RP = 5 in our model. The total mass M_sat of a satellite is

\begin{eqnarray*} M_\text{sat}(r_\text{sat}) &=& M_\star + M_\text{cold}\\ &+& 4\pi \int _0^{r_\text{sat}} \left[ \rho _\text{hot}(r) + \rho _\text{DM}(r) \right] r^2 \,\mathrm{ d}r, \end{eqnarray*}

(4)

assuming that r_sat is large enough to contain all the stars and cold gas. All density profiles are represented using isothermal spheres, i.e. ρ_hot = M_hot/  (4|$\pi$|r_hotr²), where r_hot is the radius that contains all of M_hot. This radius initially adopts the value of the subhalo virial radius, r_vir; in the case of orphan satellites, r_vir preserves the value corresponding to the last time the subhalo was identified. Combining equations (3) and (4) one can numerically solve for the hot gas stripping radius due to RP, |$r_\text{s,hot}^\text{RPS}$| (see Appendix  A).

To determine the halo gas loss of satellite galaxies via TS we assume that the hot gas distributes parallel to the DM. The bounding radius for the DM, r_DM, is given by r_vir. This radius gives the tidal radius determined by TS, that is |$r_\text{s,hot}^\text{TS}=r_\text{DM}$|⁠, and is compared with the hot gas stripping radius due to RP, |$r_\text{s,hot}^\text{RPS}$|⁠. The smaller of the two will be the stripping radius r_s,hot which will contain the remaining hot gas mass; all hot gas beyond that radius can be stripped. Thus, the value of r_hot is updated such that |$r_\text{hot}^\text{new}=r_\text{s,hot}$|⁠.

Once the hot gas outside the new stripping radius is removed, we assume that the remaining gas quickly redistributes its mass and restores an isothermal profile, but truncated at r_hot, as in Font et al. (2008) and Kimm et al. (2011). The pressure of the ambient ICM could act to confine the satellite’s hot gas, preventing any outflow from going beyond the stripping radius (Mulchaey & Jeltema 2010). If in a subsequent simulation time-step the value obtained for r_hot is larger than the previous one, r_hot is not updated and no gas is lost in that time-step.

In previous versions of sag, based on the strangulation scheme, all feedback processes occurring in satellite galaxies were assumed to transfer cold gas and its associated metals from the satellite into the hot phase of the central galaxy of their main host halo (i.e. the ICM). In our new implementation, we assume that the hot gas halo of a satellite is depleted when the hot gas mass drops below a fraction f_{hot, sat} of the baryonic mass of the satellite, specified in Section 5.1. Before this situation takes place, all feedback processes will transport gas and metals to the hot phase of the same galaxy, proportionally to its content of hot gas; the remaining reheated mass and associated metals are transferred to the corresponding central galaxy. If the latter is a central galaxy of a subhalo (i.e. it is also a satellite) and its hot gas halo also drops below f_hot,sat times the baryon fraction, then the reheated gas and metals of the orphan satellite are transferred to the main host halo.

In any given simulation time-step, all those galaxies, either central or satellite, that still have a hot gas reservoir will proceed to cool gas. Gas cooling rates are calculated using the simple model presented in Springel et al. (2001), but considering the total radiated power per chemical element given by Foster et al. (2012).

3.5 Ram pressure and tidal stripping of cold gas disc

Ram pressure exerted by the intragroup/cluster medium can also affect the cold gas disc of satellite galaxies. Its action could be regulated by the presence of hot gas halo in the gradual starvation scheme. Hot gas haloes in satellite galaxies will be replenished by SNe feedback as long as there is some star formation and the ratio between the hot gas halo and the baryonic mass of the galaxy is larger than the fraction f_hot,sat (see Section 3.4). At a certain point, this ratio will become smaller than f_hot,sat, and the hot gas halo will be gradually reduced by gas cooling and/or stripping processes (RPS, TS) reaching very low values. We assume that as long as this ratio is larger than 0.1, the hot gas halo shields the cold gas disc from the action of RPS. Once this condition is not fulfilled anymore, the ambient RP starts affecting the cold gas. This threshold has been chosen small enough to allow the role of the hot gas as a shield for a sufficiently long time; the action of RPS on the cold disc gas becomes too effective without any restriction of this kind.

We consider the model for RPS of cold gas disc introduced in Tecce et al. (2010), which is based on the simple criterion proposed by Gunn & Gott (1972). The cold gas of the galactic disc located at a galactocentric radius R will be stripped away if the RP exerted by the ambient medium on the galaxy exceeds the restoring force per unit area due to the gravity of the disc,

\begin{equation*} P_{\rm ram} \gt 2\pi G \Sigma _\text{disc}(R) \Sigma _\text{cold}(R). \end{equation*}

(5)

Here Σ_disc and Σ_cold are the surface densities of the galactic disc (stars plus cold gas) and of the cold gas disc, respectively.

The discs of stars and gas are modelled by an exponential surface density profile given by Σ(R) = Σ₀exp (−R/R_d), where Σ₀ is the central surface density and R_d is the scale length of the disc. This scale length is estimated as |$R_\text{d} = (\lambda /\sqrt{2})R_\text{vir}$| (Mo, Mao & White 1998), where λ is the spin parameter of the DM halo in which the galaxy resides. Starting from condition (5) it can be shown that RPS will remove from a galaxy all cold gas beyond a stripping radius R_s,cold given by

\begin{equation*} R_\text{s,cold}^\text{RPS} = -0.5 R_\text{d} \ln \left( \frac{P_{\rm ram}}{2\pi G \Sigma _\text{0,disc} \Sigma _\text{0,cold}} \right), \end{equation*}

(6)

where Σ_0,disc and Σ_0,cold are the central surface densities of the stellar disc and of the cold gas disc, respectively. We assume for simplicity that both disc components have the same scale length.

To account for the effect of tides, at each simulation snapshot the value for |$R_\text{s,cold}^\text{RPS}$| obtained from equation (6) is compared to r_DM, which determines the stripping radius due to TS, |$R_\text{s,cold}^\text{TS}$|⁠; if the former is smaller than the latter we set |$R_\text{s,cold}=R_\text{s,cold}^\text{RPS}$|⁠. If the resulting R_s,cold is smaller than the current value of the cold gas disc radius, R_cold, all the gas beyond R_s,cold is stripped away. The stripped gas is added to the hot gas component of the central galaxy (either of a main halo or of a subhalo). If the latter is central galaxy of a subhalo already stripped of halo gas, the stripped gas goes to the ICM.

After a stripping event produced by RPS and/or TS, the remaining disc gas is assumed to form an exponential disc truncated at a radius |$R_\text{cold}^\text{new}=R_\text{s,cold}$| and with a new scale length defined as |$R_\text{d,cold}^\text{new}=R_\text{cold}^\text{new}/7$| (assuming that 99 per cent of the cold gas disc is contained within 7 × R_d,cold).

3.6 Tidal stripping of stars

Unlike RPS, TS may also affect the stellar components of a satellite galaxy. To process the stellar stripping we consider the disc and bulge separately. In the case of the disc, we compare the current value of r_DM, with the size of the stellar disc, R_disc. If R_disc > r_DM, we assume that the galaxy loses all stars beyond r_DM and that the remaining disc stars still form an exponential disc, truncated at R_s,disc = r_DM, but with the stellar mass redistributed with a new scale length R_d,disc = R_s,disc/7, as it is considered for the cold gas disc. The scale length of the galaxy disc composed by cold gas and stars is then defined from the mass-weighted scale length of their respective components as R_d = (M_cold R_s,cold + M_disc R_s,disc)/(M_cold + M_disc).

If a bulge component is present, to evaluate the TS of its stars we assume that they are distributed according to a Hernquist (1990) profile,

\begin{equation*} \rho (r) = \frac{M_\text{bulge}}{2\pi }\, \frac{a_\text{b}}{r}\, \frac{1}{(r + a_\text{b})^3}, \end{equation*}

(7)

where a_b is a scale length related to the bulge half-mass radius by |$r_\text{b,h}= (1 + \sqrt{2}) a_\text{b}$|⁠. To calculate r_b,h we follow the procedure outlined in Cole et al. (2000) and adapted to our model in Muñoz Arancibia et al. (2015). If the radius that contains the 99 per cent of the bulge mass, given by r_bulge= 198.5 × a_b, satisfies the condition r_bulge > r_s,bulge with r_s,bulge = r_DM, we assume that the stars in the bulge beyond r_s,bulge are stripped and that the remaining bulge stars are redistributed still following a Hernquist profile truncated at r_s,bulge, with a new scale length defined as |$a_\text{b}^\text{new}=r_\text{s,bulge}/198.5$|⁠.

The stars removed by TS are assumed to orbit freely within the main host halo. For practical purposes they are assigned to a new ‘stellar halo’ component of the corresponding central galaxy. As these halo stars evolve and die, they inject gas and metals directly into the ICM.

4 SUPERNOVAE FEEDBACK

In previous versions of sag, the amount of reheated mass produced by the SNe arising in each star-forming event is assumed to be

\begin{equation*} \Delta M_{\rm reheated} = \frac{4}{3} \epsilon {\frac{\eta E_\text{SN}}{V_{\rm vir}^2}} \Delta M_{\star }, \end{equation*}

(8)

where η is the number of SNe generated from the stellar population of mass ΔM_⋆ formed, E_SN = 10⁵¹ erg is the energy released by a SN, V_vir is the virial velocity of the host (sub)halo, which is a measure of its potential well, and ε is the SNe feedback efficiency, i.e. a free parameter that controls the amount of cold gas reheated by the energy generated by SNe. The number of SNe depends on the initial mass function (IMF) adopted and for core collapse supernovae (SNe CC) is estimated as

\begin{equation*} \eta = \frac{\int ^{\infty }_8 \phi (m)\ {\rm d}m}{\int _0^{\infty } \phi (m)\, m\ {\rm d}m}, \end{equation*}

(9)

where m is the stellar mass and ϕ(m) is the IMF. This quantity is constant for an universal IMF, like the Chabrier IMF (Chabrier 2003) adopted here. The energy and metals generated by SNe and stellar winds produced by progenitors in different mass ranges are released with different time-scales depending on the lifetime of the progenitors, as given by Padovani & Matteucci (1993).

In order to satisfy observational constraints at high redshifts, we found it necessary to modify this model of SNe feedback. Hirschmann et al. (2016) find that the observed trends in galaxy assembly are reproduced by replacing the ‘fiducial’ feedback scheme implemented in gaea model with a parametrization of the mass-loading factor similar to the one inferred by Muratov et al. (2015) from the analysis of the FIRE simulations. We use this information to modify the estimation of the reheated mass (equation 8) by simply adding new factors that take into account the dependence on redshift and an additional modulation with virial velocity. Thus, the feedback scheme in the current version of sag produces a reheated gas mass given by

\begin{equation*} \Delta M_{\rm reheated} = \frac{4}{3} \epsilon {\frac{\eta E_\text{SN}}{V_{\rm vir}^2}} \,(1+z)^{\beta }\,\left(\frac{V_{\rm vir}}{60\,{\rm km}\,{\rm s}^{-1}}\right)^{\alpha _\text{F}}\Delta M_{\star }, \end{equation*}

(10)

where the exponent α_F takes the values −3.2 and −1.0 for virial velocities smaller and larger than 60 km s⁻¹, respectively. Although the fit provided by Muratov et al. (2015) indicates that the value of β is 1.3, we consider this exponent to be a free parameter of the model sag, allowing it to absorb aspects of physical processes that are not properly captured by the model and helping to reproduce the observational constraints imposed, detailed in Section 5.1. In Section 5.2, we show the results of this procedure and also discuss the impact of assigning to β the value found by Muratov et al. (2015).

The reheated gas is transferred from the cold to the hot phase, subsequently returning to the cold phase through gas cooling taking place in both central and satellite galaxies. However, to avoid an excess of stellar mass at high redshifts, some of the hot gas must be ejected out of the halo reducing the hot gas reservoir available for gas cooling (Guo et al. 2011; Henriques et al. 2013; Hirschmann et al. 2016). Hence, we also consider the energy conservation argument presented by Guo et al. (2011) to calculate the ejected hot gas mass

\begin{equation*} \Delta M_{\rm ejected}= \frac{\Delta E_\text{SN} - 0.5\,\Delta M_{\rm reheated}\,V_{\rm vir}^2}{0.5\,V_{\rm vir}^2}, \end{equation*}

(11)

where ΔE_SN is the energy injected by massive stars which we model in a way similar to the modified reheated mass, as also done by Hirschmann et al. (2016),

\begin{equation*} \Delta E_{\rm SN} = \frac{4}{3} \epsilon _\text{ejec} {\frac{\eta E_\text{SN}}{V_{\rm vir}^2}} \,(1+z)^{\beta }\,\left(\frac{V_{\rm vir}}{60\,{\rm km}\,{\rm s}^{-1}}\right)^{\alpha _{\rm F}}\Delta M_{\star }\,0.5\,V_\text{SN}^2. \end{equation*}

(12)

Here ε_ejec is the corresponding efficiency, considered as another free parameter of the model, and |$0.5\,V_\text{SN}^2$| is the mean kinetic energy of SN ejecta per unit mass of stars formed. Following the analysis of Muratov et al. (2015), we adopt the fit for 95th percentile wind velocity as a function of the virial velocity of the halo (see their fig. 8 and eq. 10), such that |$V_\text{SN}=1.9\,V_\text{vir}^{1.1}$|⁠.

In the gradual starvation scheme implemented in sag, galaxies keep their hot gas halo when they become satellites. These haloes are reduced by gas cooling and environmental effects (RPS, TS), but they also can be reconstructed by the injection of reheated gas or the reincorporation of ejected gas. The ejected gas mass is assumed to be re-incorporated back on to the (sub)halo from which it was expelled within a time-scale that depends on the inverse of (sub)halo mass, M_vir, as assumed by Henriques et al. (2013) in order to reproduce the observed evolution of the galaxy stellar mass function (SMF). Thus, the reincorporated mass is given by

\begin{equation*} \Delta M_{\rm reinc}= \gamma \,\Delta M_\text{ejected}\,\frac{M_{\rm vir}}{10^{10}\,{\rm M}_{{\odot }}}, \end{equation*}

(13)

where the parameter γ regulates the efficiency of the process and is also a free parameter of sag model.

In all the above parametrizations, we have added the energy injected by SNe Ia. We adopt the single degenerate model in which an SN Ia occurs by carbon deflagration in C–O white dwarfs in binary systems whose components have masses between 0.8 and 8 M_⊙ (Greggio & Renzini 1983). We implement the formalism presented by Lia, Portinari & Carraro (2002), choosing a fixed fraction of binary systems, A_bin = 0.05. As for SNe CC, we also consider the lifetime of the progenitors.

SNe CC and SNe Ia, together with low- and intermediate-mass stars, contribute with metals that pollute the cold and hot gas, affecting the cooling rates and regulating also in this way the subsequent events of star formation. Details on the chemical model implemented in sag are given in Cora (2006), with the latest updates on chemical yields presented in Gargiulo et al. (2015).

5 CALIBRATION OF THE MODEL SAG

The physical processes included in sag model, as occurs for all SAMs, involve several free parameters which regulate the modelling of these processes. Proper values of these free parameters are found by imposing constraints given by observed galaxy properties that the SAM should satisfy.

5.1 Calibration process

The calibration process of sag model is performed by implementing the particle swarm optimization (PSO) technique presented in Ruiz et al. (2015). This numerical tool allows us to tune the free parameters via the exploration of the sag parameter space by random walks of a set of ‘particles’ (the swarm, equivalent to the chains in Monte Carlo Markov techniques) that share information between them, thus determining the evolution of the exploration from both their individual and collective experience. By comparing the model results against a given set of observables, the PSO method yields a set of best-fitting values for the free parameters.

For this current version of sag, we consider nine free parameters for the calibration process. The corresponding equations that involve these free parameters and a brief description of the related physical processes can be found in Ruiz et al. (2015), except those related to the new SN feedback and ejection/reincorporation scheme implemented in the current version of sag, presented in Section 4. The free parameters are the star formation efficiency (α), the efficiency of SN feedback from stars formed in both the disc and the bulge (ε, equation 10), the efficiency of ejection of gas from the hot phase (ε_ejec, equation 12) and of its reincorporation (γ, equation 13), the exponent that regulates the redshift evolution of the mass-loading factor of the reheated and ejected mass (β, equations 10 and 12), the growth of super massive BHs and efficiency of AGN feedback (f_BH and κ_AGN, respectively), the factor involved in the distance scale of perturbation to trigger disc instability events (f_pert), and the fraction that determine the destination of the reheated cold gas (f_hot,sat) introduced in Section 3.4.

The set of observables used for the calibration process is the SMFs at z = 0 and z = 2, the star formation rate distribution function (SFRF), the fraction of mass in cold gas as a function of stellar mass (CGMF) and the relation between bulge mass and the mass of the central supermassive BH (BHB). For the SMFs, we adopt the compilation data used by Henriques et al. (2015), which for z = 0 is a combination of the SMF of the SDSS from Baldry, Glazebrook & Driver (2008) and Li & White (2009), and of the Galaxy And Mass Assembly from Baldry et al. (2012), while for z = 2 is a combination of the data of the Cosmic Evolution Survey (COSMOS) from Domínguez Sánchez et al. (2011), the Ultra Deep Survey from Muzzin et al. (2013) and Ilbert et al. (2013), and the FourStar Galaxy Evolution Survey from Tomczak et al. (2014). For the SFRF, which is the number density of galaxies in a certain SFR interval, we use data from a flux-limited sample of galaxies observed with the Herschel satellite which gives the total (IR + UV) instantaneous SFR for the redshift interval z ∈ [0.0, 0.3], as presented in Gruppioni et al. (2015). This observable is compared with the sag output at z = 0.14, which is the closest output of the model to z = 0.15, the mean value of the redshift range of observations. For the CGMF, we adopt observational data from Boselli et al. (2014), which is based on a volume limited sample, within the range |$\rm {log}(\mathrm{ M}_{\star }[{\rm \mathrm{ M}}_{{\odot }}]) \in [9.15,10.52]$| in stellar mass. Finally, for the BHB relation, we combine the data sets from McConnell & Ma (2013) and Kormendy & Ho (2013). This particular set of constraints and the observational data adopted were defined as a common set for calibration during the Cosmic CARNage workshop aimed to compare different calibrated galaxy formation models based on the same cosmological DM-only simulation (Knebe et al. 2018b). This comparison project was initiated in the nIFTy Cosmology workshop (Knebe et al. 2015).

Due to the large size of the simulation, the computing time needed to run the model several times becomes quickly prohibitively long. Hence, the calibration is carried out by running sag with the merger trees extracted from a smaller box of the whole MDPL2 simulation that constitute a representative sample of all the merger trees contained in the 1 h⁻¹ Gpc side-length box. We divide the full volume of the MDPL2 simulation in 9³ parallelepiped sub-boxes, each of them with roughly the same number density of haloes of the whole volume. Despite each of these samples constitutes only ∼0.137 per cent of the whole simulation, they are large enough to be considered themselves simulations from which a statistically significant set of merger trees can be extracted for the further application of the sag model. Therefore, we compared the halo mass function for all these sub-boxes with the one obtained for the whole simulation, choosing the most similar one simply by visual inspection, selecting for this case a sub-box with an effective length side of 111.19h⁻¹ Mpc.

The values of the parameters obtained from the calibration of the model to match the aforementioned set of observables are shown in Table 1.

Table 1.

Best-fitting values of the free parameters of sag model obtained with the PSO technique. This set of values is obtained from the application of sag to the merger trees of the sub-box selected from the MDPL2 simulation.

Parameter	Best-fitting value
α	0.04
ε	0.33
εejec	0.022
γ	0.055
β	1.99
fBH	0.06
κAGN	3.02 × 10−5
fpert	14.56
fhot,sat	0.277

5.2 Functions and relations used as constraints

The new features included in sag allow the model to reach a good match with the observational constraints imposed to calibrate the model, i.e. the SMF at z = 0 and z = 2, the SFRF at z = 0.14, and the CGMF and BHB and relations at z = 0, as we detail below. In addition to the results of the calibrated model, we also evaluate the impact on these functions and relations of fixing the value of the parameter β involved in the redshift dependence of the reheated and ejected mass. For this alternative model, we adopt the value given by the fit of Muratov et al. (2015), i.e. β = 1.3, leaving the values of the remaining free parameters of sag unchanged. We refer to this slightly modified version of the model as sag_β1.3.

It is worth noting that the small value of the parameter γ that regulates the amount of ejected gas that is reincorporated makes the results of the model insensitive to the fate of the reincorporated gas, i.e. subhalo from which it was ejected or main host halo where the satellite resides. We tested the latter possibility and found no change in the galaxy properties shown in this work.

5.2.1 Stellar mass functions

In Fig. 1 both the local SMF (top panel) and the SMF at z= 2 (bottom panel) generated by the model are shown, where the break and low-mass end given by observations are reproduced. This is due to the new implementation of SN feedback. The classical estimation of the reheated mass considered in previous versions of sag (equation 8) allows a good match of the SMF only at z = 0. The addition of ejection acting on the hot phase and the redshift dependence of both reheated and ejected mass (equations 10 and 12) become key ingredients to avoid an excess of star formation overproducing stellar mass at high redshifts. This good agreement is obtained by a stronger redshift dependence than the one found by Muratov et al. (2015), whose parametrization is characterized by a parameter β = 1.3. From Table 1, which shows the best-fitting values of the free parameters of sag, we can see that the exponent that regulates this redshift dependence takes a value β = 1.99. On the other hand, model sag_β1.3 produces an excess in the low-mass end of the SMF at z = 2, as shown by the dashed line in Fig. 1. This is due to the higher levels of SFR this model features at higher redshifts, as shown and discussed in Section 6.1. However, the SMF at z = 0 is not significantly affected by the change in the value of the parameter β.

Figure 1.

SMFs used as constraints to calibrate the sag model. Results obtained from the calibrated model sag with the best-fitting values given in Table 1 are represented by red solid lines. Top panel: SMF at z = 0. Bottom panel: SMF at z = 2. Observational data compiled by Henriques et al. (2015) are represented by filled circles with error bars. The agreement between model and observations is rather good, except for the excess at the high-mass end at z = 0. This excess becomes less evident when comparing with estimations compiled and analysed by Bernardi et al. (2017, grey shaded area), which include data from Chen et al. (2012) and Mendel et al. (2014). In both panels, results obtained from model sag_β1.3 are shown by a blue dashed line. The SMF at z = 0 is not significantly affected by the change in the value of the parameter β, whereas this gives rise to an excess in the low-mass end of the SMF at z = 2.

None of the improvements included in sag leads to a better fit of the massive end of the SMF at z = 0, which shows an excess of galaxies with stellar mass M_⋆ ≳ 2 × 10¹¹ M_⊙. This excess is not avoided by the implementation of the more locally efficient radio mode AGN feedback (Henriques et al. 2015). In any case, AGN feedback allows to recover the break of the SMF which characteristic stellar mass highly depends on the feedback efficiency regulated by the two free parameters involved in the AGN feedback (κ_AGN and f_BH). However, the slope of the high-mass end of the SMF is insensitive to the AGN feedback efficiency, as we have checked by manually exploring values of the parameters related with this process around the best-fitting values.

It is important to take into account the uncertainties in the observational data used as constraints, considered in the comparison plot of the SMF at z = 0. As specified in Section 5.1, the compilation made by Henriques et al. (2015) includes data from Li & White (2009) for z = 0. As discussed by Bernardi et al. (2017), this particular set of data gives the lowest values of the comoving number density of galaxies for stellar masses ≳10¹¹ M_⊙ in comparison with other results presented in the literature because of the use of inappropriate algorithms for estimating the observed flux of the most luminous galaxies and inadequate assumptions about the stellar population modelling. While systematic effects on the SMF related with photometry account for differences of only 0.1 dex, systematics arising from different treatments of the stellar population can be as large as 0.5 dex. The estimations compiled and analysed by Bernardi et al. (2017), which include data from Chen et al. (2012) and Mendel et al. (2014), are represented by a grey shaded area. When comparing with these sets of data, the excess in the comoving number density of massive galaxies predicted by sag becomes less pronounced.

5.2.2 Star formation rate function

In addition to the SMF at different redshifts, the SFRF at z = 0.14 helps to recover the observed evolution of SFR and mass growth in galaxies. Fig. 2 presents the SFRF at z = 0.14 for the galaxy population generated by the calibrated model sag (solid line). The agreement with the observations by Gruppioni et al. (2015) is very good for intermediate values of the SFR, in the range 1 M_⊙ yr⁻¹ ≲ SFR ≲ 20 M_⊙ yr⁻¹, being slightly underpredicted (still within the dispersion of observed data) for lower values. sag gives higher number density of galaxies for SFR ≈ 60 M_⊙ yr⁻¹, reaching again a good agreement for the highest SFR bin (≈200 M_⊙ yr⁻¹). This particular behaviour of our model is quite similar to the predictions at low redshifts provided by other SAMs considered by Gruppioni et al. (2015) in a direct comparison with their data. The agreement shown here for SFR above 100 M_⊙ yr⁻¹ is better than the one reached by the ‘FIRE feedback model’ analysed by Hirschmann et al. (2016) from which the current feedback scheme implemented in sag is inspired. This better agreement could have been achieved simply because the SFRF was used to calibrate the model. However, it is not always possible to recover the complete behaviour of an observational constraint, as it is the case of the massive end of the SMF at z = 0. This tension that emerges during the calibration procedure denotes the need of improvement in the physical processes modelled and/or inconsistencies between different sets of observed galaxy properties. On the other hand, model sag_β1.3 produces a decrease in the number density of high SFR values, as shown by the dashed line.

Figure 2.

SFRF at z = 0.14 used as constraint to calibrate the sag model. Results obtained from the calibrated model sag with the best-fitting values given in Table 1, and from the sag_β1.3 model are shown by red solid and blue dashed lines, respectively. The latter produces a decrease in the number density of high SFR values. Model results are compared with observations by Gruppioni et al. (2015) represented by filled circles with error bars.

5.2.3 Cold gas mass fraction

The fact that galaxy growth and star formation quenching are well captured by sag becomes evident in the cold gas content of model galaxies (Fig. 3). The CGMF obtained from sag also behaves as expected within the stellar mass range considered for the calibration of the model (log(M_⋆[M_⊙]) ∈ [9.15, 10.52]). Their mean values (thick solid line) are in good agreement with observational data from Boselli et al. (2014, filled circles with error bars); thin solid lines represent the standard deviation around the mean and the shaded contours show the distribution of model values. For high stellar masses (log(M_⋆[M_⊙]) > 10.52), sag predicts lower fractions of cold gas with respect to those expected from a simple extrapolation of the trend found at lower masses. However, these values are still contained within the range of cold gas fractions inferred from the atomic and molecular gas fractions presented by Saintonge et al. (2016, see their fig. 3), obtained from data of the Arecibo Legacy Fast ALFA (ALFALFA), the GALEX Arecibo SDSS Survey (GASS), and the CO Legacy Database for GASS (COLD GASS) surveys. The same trend is obtained with model sag_β1.3.

Figure 3.

Fraction of mass in cold gas as a function of stellar mass at z = 0 used as constraints to calibrate the sag model. Results obtained from the calibrated model sag with the best-fitting values given in Table 1 are represented by red solid lines and reddish filled contours. Red thin solid lines denote the standard deviation around the mean (thick red solid line). In both cases, the contour levels are: [0.01, 0.19, 0.26, 0.38, 0.68, 0.95, 0.997] in terms of the maximum number density of points of each sample. CGMF is almost unaffected by adopting β = 1.3, as shown by results from the sag_β1.3 model represented by a blue dashed line and empty dotted contours. This relation is compared with observations by Boselli et al. (2014, filled circles with error bars).

5.2.4 Black hole–bulge mass relation

The BHB relation displayed by the calibrated sag model (Fig. 4, filled contours) is in good general agreement with the observational trend denoted by data from McConnell & Ma (2013) and Kormendy & Ho (2013). However, black holes seem to have higher masses than expected for large bulge masses, but still within the allowed ranges defined by the dispersions of the observational data. This trend is preserved in model sag_β1.3 (empty dotted contours). This is a result of the restriction imposed by both the high-mass end of the SMF at z = 0 and the BHB relation. During the calibration process, the parameter that controls the growth of super massive BHs, f_BH, takes a high value in order to make the AGN feedback as effective as possible. Due to the existing degeneration between the f_BH and κ_AGN parameters, a similar behaviour of the high-mass end of the SMF at z = 0 can be achieved by taking lower values of f_BH at the expense of increasing the feedback efficiency, κ_AGN. Nevertheless, large values of f_BH are preferred as a direct consequence of the statistical test used in the calibration process. The chi-square calculated for each constraint (Ruiz et al. 2015, see their equation 27) gives heavier weights to the mass bulge ranges where the bulk of the galaxies are located, which decreases the statistical significance of the high-mass galaxies in the selection of the best-fitting parameters.

Figure 4.

Relation between black hole and bulge mass at z = 0 used as constraints to calibrate the sag model. Results obtained from the calibrated model sag with the best-fitting values given in Table 1 are represented by reddish filled contours, while results obtained from the sag_β1.3 model are shown by empty dotted contours. In both cases, the contour levels are the same as those of Fig. 3. The BHB relation is almost no affected when adopting β = 1.3. Observational data are taken from Kormendy & Ho (2013, black circles) and McConnell & Ma (2013, black triangles).

5.2.5 Origin of the excess in the high-mass end of the SMF at z = 0

The AGN feedback together with additional modifications introduced in the model, that is, the redshift dependence of the reheated and ejected material and the scalings adopted for the reincorporation of the latter (equations 10, 12, and 13), allow to recover downsizing in the SFRs, as shown later in Fig. 6. Hence, there is not an excess of SF in high-mass galaxies that could be responsible for the excess of the high-mass end of the SMF at z = 0.

A test made with the semi-analytic model sagediscussed in Knebe et al. (2018a) shows that stars added to an intracluster component as a result of tidal disruption of satellite galaxies would enhance the massive end of the SMF if they would have merged with the central galaxy. In sag, satellite galaxies are not tidally disrupted when becoming orphans as in sage but suffer TS. The stripped mass is added to the stellar halo of the corresponding central galaxy (see Section 3.6) and represents the intracluster stars. The stellar mass density of this component is ∼2 and ∼3 orders of magnitude lower than the one characteristic of the stellar mass of the whole galaxy population, for redshifts z = 0 and 2, respectively. This means that the excess of the high-mass end of the SMF at z = 0 could be attributed to inefficient stripping or disruption of satellites by effects of tides in the model.

In order to evaluate the efficiency of TS in our model, we estimate the fraction of mass in stars contained in the stellar halo of central spiral galaxies and compare them with results from Merritt et al. (2016). They estimate the stellar mass in excess of a disc + bulge beyond five half-mass radius of eight spiral galaxies in the Dragonfly Nearby Galaxies Survey and obtain an average halo fraction of 0.009 ± 0.005. In our model, spiral galaxies are those characterized by a ratio between stellar bulge and total stellar mass B/T < 0.85; this cut allows to reproduce the observed morphological distribution (Conselice 2006). Median values of the halo fractions (estimated with respect to the total stellar mass of central galaxies, including the stellar halo) increase with stellar mass and are comprised within the range ≈3 × 10⁻⁵ to 7 × 10⁻⁴ for stellar masses M_⋆ ≈ 10¹⁰−10¹¹ M_⊙. These median values are at least one order of magnitude smaller than the average halo fraction obtained by Merritt et al. (2016). Although the difference between the methods used to obtain the mass of the stellar halo in the observations and in the model might contribute to such a discrepancy, it seems quite plausible that TS is not efficient enough in our model, as hinted by the excess in the high-mass end of the SMF at z = 0. In other words, the mass of merging satellites might be higher than expected because they do not suffer enough TS.

Another possible cause of the excess in the high-mass end of the SMF at z = 0 might be related to dry mergers at z ≲ 2 of low-mass galaxies with massive ones. In this context, it is worth mentioning a caveat regarding the merger condition adopted for orphan galaxies. As we briefly explain in Section 3.2, it is based on a criterion that demand an assumption about the radius of the central galaxy (assumed to be 10 per cent of the virial radius of the host halo). The drawback of this approximate estimation is that it can be very large for cluster-size haloes. Therefore, the integrated orbit can easily take the orphan satellite within the sphere determined by those radii, thus merging with its central. In order to test this possibility, we estimate the conditional SMFs for central and satellite galaxies within different ranges of halo mass in sag, which are shown in Fig. 5. Following Kang (2014), they are compared with data from SDSS DR4 (Yang, Mo & van den Bosch 2009). There are less model satellites than observed in all the halo mass ranges analysed. This deficiency is more pronounced in more massive haloes (log (M_halo[M_⊙]) ∈ [13.5, 13.8], [14.1, 14.4]), and it is accompanied by a shift to higher stellar masses of the population of central galaxies.

Figure 5.

Conditional SMF for central (solid line) and satellite (dashed line) galaxies in model sag. Galaxies are selected according to the mass of the main host halo they reside in as indicated in the legend of the different panels. Model results are compared with data from the group catalogues constructed by Yang et al. (2009) from SDSS DR4; central galaxies are identified by filled circles and satellites are represented by open squares.

Hence, both an inefficient TS and overmerging of satellite galaxies are responsible for the excess of high-mass galaxies found in the SMF at z= 0. These drawbacks in the model are mainly related to aspects of the integration of orphan galaxies, like the treatment of TS and the merging criterion, which will be revisited (Vega-Martínez et al. in preparation).

6 MODEL PREDICTIONS

We present galaxy properties predicted by the model, that are not used as constraints to calibrate it. We first focus on the evolution of the cosmic SFR density (SFRD) and of the specific SFR (sSFR). The former is given by the volume-averaged sum of SF of all galaxies at any given time, and the latter is defined as the ratio between the SFR and the stellar mass of the galaxy. We connect the trends of the sSFR of satellites and centrals with the stellar mass content of galaxies populating DM haloes of different mass through the stellar-to-halo mass relation (SHM). According to the value of the sSFR, we classify galaxies as active and passive. We analyse the dependences of the fraction of passive galaxies with stellar mass, halo mass and halo-centric distance.

6.1 Star formation rate density

The evolution of the cosmic SFRD predicted by sag is presented in Fig. 6 (thick solid line) and is compared with data compiled by Behroozi, Wechsler & Conroy (2013c). Although model results show consistency with observations for z< 4, they underpredict the SFRD at higher redshifts. Besides, the decline in the SFRD towards low redshifts is less pronounced in the model than in observations, thus resulting in an excess of SFRD at z = 0. The peak and normalization of the SFRD vary significantly among different galaxy formation models (e.g. Guo et al. 2016).

Figure 6.

Evolution of the cosmic SFR density for all galaxies (red thick solid line) compared with observational data compiled by Behroozi et al. (2013c, fiiled circles with error bars). Contribution from galaxies lying in different mass ranges at z = 0 are represented by thin black lines of different style, as indicated in the legend. Downsizing in stellar mass assembly becomes evident. Intermediate-mass galaxies log (M_⋆[M_⊙]) ∈ [10.1, 10.5] are the main responsible of the excess in the SFRD at z = 0. The global SFRD obtained from the calibrated model sag but adopting β = 1.3 (model sag_β1.3) is added for comparison, being represented by a blue thick dashed line. Model sag_β1.3 is characterized by higher (lower) levels of SFR at high (low) redshifts with respect to the calibrated model sag, in better agreement with observations.

The low SFRD predicted at high redshifts is a consequence of the redshift dependence of the reheated and ejected gas fractions in the new feedback scheme implemented in sag (equations 10 and 12), an effect that becomes even stronger in our model because of the preferred high value of β in the calibration process, as discussed in Section 5.2. The lower value suggested by the fit found by Muratov et al. (2015), e.g. β = 1.3, allows to reconcile model predictions with observed SFRD at high redshifts, as shown in Fig. 6 by model sag_β1.3 (thick dashed line). This better match takes place at the expense of increasing the number density of low-mass galaxies at z = 2 (dashed line in bottom panel of Fig. 1). Therefore, it is not possible to satisfy both observational constraints simultaneously. This tension between the observed redshift evolution of the SFRD and z = 2 galactic SMF has been also noted by Hirschmann et al. (2016) from their analysis of the evolution of the correlation of sSFR with the stellar mass (main sequence of star-forming galaxies, see their fig. 9); they find that ejective models that successfully reproduce the measured evolution of SMF tend to underestimate the sSFR at high redshifts. In the same line, the recent work of Rodrigues, Vernon & Bower (2017), based on the exploration of the parameter space of a version of the semi-analytic model, shows that the evolution of the SMF is recovered by a particular set of parameters that favours larger SN feedback efficiency at higher redshifts, although the rise of the SFRD with redshift in the range 0.5 <z < 1 is not as steep as observed.

The difficulty in satisfying simultaneously the observed evolution of the SFRD and the SMF at high redshift also emerges in the analysis of observational results. It is well known that the star formation history inferred from the observed evolution of the stellar mass differs from instantaneous indicators of star formation, being the former ≈0.6 dex smaller than the latter at z = 3 (Wilkins, Trentham & Hopkins 2008). Thus, the integrated star formation history implies a local stellar mass density in excess of that measured. According to recent studies (Madau & Dickinson 2014), this discrepancy in the stellar mass density can be ≈0.2 dex for z ≲ 3 depending on the data considered, being smaller than previously estimated. In any case, these discrepancies can arise for several reasons such as inaccurate dust extinction corrections, underestimation of stellar masses due to the outshining of old stellar populations in star-forming galaxies, or the evolution of the integrated stellar IMF in galaxies. Both the star formation history implied by instantaneous indicators, which are typically dominated by very massive stars, and the star formation history inferred from the evolution of the average stellar mass density are affected by the IMF assumed. The discrepancies between them might be mitigated when considering an evolving IMF that is top heavy at high redshifts (Wilkins et al. 2008).

The stronger suppression of SF at high redshifts is evident in galaxies of all masses being more pronounced in low-mass galaxies, thus allowing to reach the antihierarchical assembly of stellar mass. This becomes evident from the evolution of the galactic SMF presented by Hirschmann et al. (2016, see their fig. 1); for ejective feedback models, high-mass galaxies (M_⋆ ≳ 10¹⁰ M_⊙) are already in place at z ≈ 2−3, while the population of low-mass galaxies keeps increasing towards lower redshifts, in agreement with the trend denoted by observational measurements. Downsizing in galaxy assembly is also produced by sag, as can be appreciated from the contribution to the global SFRD of galaxies in different mass ranges selected at z = 0, as indicated by thin black lines of different style in Fig. 6. The peak of the SFRD of galaxies with stellar mass M_⋆ > 10¹¹ M_⊙ is located at z ≈ 2, while intermediate-mass galaxies (log (M_⋆[M_⊙]) ∈ [10.5, 10.9]) reach the peak at z ≈ 1. Lower mass galaxies are still in the regime of increasing SFRD while approaching to z = 0. Downsizing is also present in model sag_β1.3 but with a change in the SF history of galaxies which start to form stars earlier. As expected, this general trend is consistent with the predictions of ejective feedback models in Hirschmann et al. (2016, see bottom row of their fig. 4).

Galaxies with stellar masses M_⋆ ≲ 3 × 10¹⁰ M_⊙ are responsible for the excess of the local SFRD with respect to observations, especially those in the mass range log (M_⋆[M_⊙]) ∈ [10.1, 10.5] that dominate the contribution to the SFRD at z ≲ 0.3. Hirschmann et al. (2016) note this problem for the most massive galaxies (M_⋆ ∼ 10¹² M_⊙) and attribute this failure in reproducing the SFRD at z ≲ 0.5 to a radio-mode AGN feedback not efficient enough in suppressing SF in these galaxies. Although this could be a possible explanation for galaxies within the stellar mass range log (M_⋆[M_⊙]) ∈ [10.1, 10.5] in our model, the restrictions imposed to calibrate sag with the PSO method could also favour this excess. Both the high-mass end of the SMF at z = 0 and the local BHB relation help finding appropriate values of the free parameters involved in the AGN feedback scheme. However, the SFRF at z = 0.14 also plays a role. We have shown that the tuned parameters allow to achieve a highly satisfactory agreement with the data presented by Gruppioni et al. (2015) (see Fig. 2). However, the IR + UV SFRD estimated from these data through the integration of the best-fitting modified Schechter function to the IR+UV SFRFs down to  log(SFR)=−1.5 is higher than the optical SFRD presented by Behroozi et al. (2013c) at low redshift (z < 0.5), which is reflected in the excess shown by the model in Fig. 6. We evaluate the impact of the constraint used in the calibration process by imposing the evolution of the SFRD up to z = 2 as an observational restriction in place of the SFRF at z = 0.14. We consider the values for the SFR at several redshifts within that range taken from the fit to observational data presented in Behroozi et al. (2013c). We find that the excess of the SFRD at z = 0 is not avoided even with this restriction. This aspect deserves more investigation, in order to disentangle what particular combination of physical processes is causing this effect, or what process is still not well captured by the modelling. The analysis done in Section 5.2.5 points to TS, tidal disruption, and overmerging as processes that deserve special attention.

6.2 Specific star formation rate: main sequence

Fig. 7 shows the relation between the sSFR and stellar mass of galaxies generated by the calibrated model sag at different redshifts (z= 0, 0.7, 2.2, and 3.5) by the density map. At z= 0 (left panel), this distribution resembles very much the one presented by Salim et al. (2007), both in shape and normalization. They estimate the SFR of 50 000 optically selected galaxies in the local Universe from gas-rich dwarfs to massive ellipticals using ultraviolet and optical data from the Galaxy Evolution Explorer (GALEX) and SDSS Data Release 4, respectively. Comparing with their fig. 15, it is evident that the dependence of sSFR with stellar mass for model galaxies is shallower than observed. This can be better quantified by considering the mean values of sSFR for different stellar mass bins for star-forming galaxies. The tight correlation followed by them is known as the main sequence of galaxies. Following Brown et al. (2017), who investigate environment driven gas depletion in satellite galaxies using a sample from SDSS, we classify galaxies as star-forming when they have sSFR > 10^−10.7 yr⁻¹. This limit allows a better separation between active and passive galaxies in our model than the cut sSFR = 10⁻¹¹ yr⁻¹, commonly used in the literature (e.g. Wetzel et al. 2012, W12 hereafter). Although this is inferred from the bimodality that emerges for massive galaxies within groups and clusters (see Section 6.3), we apply this cut to galaxies in all mass ranges, following W12. This criterion also comprises the star-forming sample of galaxies selected by Salim et al. (2007, see their fig. 17) through optical emission lines and the BPT diagram.

Figure 7.

sSFR as a function of stellar mass for all galaxies in the calibrated model sag (filled contours, considering the same levels as the Fig. 3) at different redshifts (z= 0, 2.2, and 3.5). Mean values of star-forming galaxies (sSFR > 10^−10.7 yr⁻¹) are estimated for central and satellite galaxies (red solid thick and thin lines, respectively); error bars represent 1σ-standard deviation around the mean. They are compared with redshift dependent fits to observational data compiled by Behroozi et al. (2013c, filled circles with error bars). The gradual removal of hot gas in satellites through RPS and TS allows satellite galaxies to have a behaviour similar as central ones. Results of model sag_β1.3 for both central and satellite galaxies are represented by blue dashed thick and thin lines, respectively.

Mean values of sSFR of star-forming galaxies are shown for central and satellite galaxies separately. These main sequences are compared with redshift dependent fits to observational data collected by Behroozi et al. (2013c), all corrected to a Chabrier IMF (see their table 8 in appendix F). The data set, specified in their table 5, includes the sSFR measured by Salim et al. (2007) described above. The general agreement is rather good at all redshifts, especially regarding the normalization. Reproducing the right normalization of this correlation has been challenging; galaxy formation models have failed in reproducing this feature underpredicting it, especially at high redshifts (Daddi et al. 2007; Weinmann et al. 2012; Xie et al. 2017).

The similarity between the main sequence of central and satellite galaxies for any of the redshifts considered highlights theimportance of modelling environmental effects through RPS and TS, with the consequent gradual removal of hot halo gas. The strangulation scheme considered in previous versions of sag quenches star formation in satellite galaxies too early, leading to values of SFR one order of magnitude lower than those characterizing central galaxies, thus generating a complete separate main sequence for satellites. Satellite galaxies in the current model lie on a main sequence that is systematically lower than the one traced by centrals at all redshifts, but differences are within the 1σ-standard deviation around the mean. Although the normalization is quite good for both central and satellite galaxies, the model predicts flatter sequences than those inferred from observations. Mean values of sSFR of central galaxies decrease for stellar masses M_⋆ ≳ 10¹⁰ M_⊙, being this trend more pronounced for centrals than for satellite galaxies; the main sequence for satellites remains almost flat at z = 3.5. This leads to more similar values of mean sSFR for these two galaxy populations at high stellar masses. For low-mass galaxies (M_⋆ ≲ 10¹⁰ M_⊙), the agreement between sag predictions and observational fits is better at z = 0 than at higher redshifts, while the opposite situation occurs for high-mass galaxies. The higher values of sSFR at z = 0 for galaxies with stellar masses within the range log (M_⋆[M_⊙]) ∈ [10.1, 10.5] are consistent with the fact that these galaxies are the main responsible of the excess of the cosmic SFRD at z = 0 (see Fig. 6).

Both the normalization and shape of the main sequence of central and satellite galaxies are affected when considering model sag_β1.3, as it is shown in Fig. 7 by dashed thick and thin lines, respectively. Central and satellite galaxies from model sag_β1.3 achieve smaller sSFR values than those characterizing the same galaxy types generated by model sag in all the redshifts considered, being these differences larger for smaller stellar masses. They originate in the large amount of stellar mass acquired by low mass galaxies at high redshift in model sag_β1.3 (see bottom panel of Fig. 1) because the SFRD in this model is higher than in sag for z ≳ 1 (see Fig. 6). Besides, the SFRD in sag_β1.3 is slightly lower at z = 0, a pattern that is common to galaxies within different z = 0 stellar mass ranges, contributing to reduce the sSFR. Moreover, the decrease of sSFR with decreasing stellar mass for M_⋆ ≲ 10¹⁰ M_⊙, which is a trend opposite to the one typical of observational data, becomes more evident in model sag_β1.3. This particular behaviour is similar to the one obtained by Xie et al. (2017, see their fig. 11) who use an updated version of the semi-analytic model gaea (Hirschmann et al. 2016) considering the same ejective model of SN feedback in which the new feedback scheme implemented in our model is inspired; the similarity with results from model sag_β1.3 is not surprising since they adopt a parameter β = 1.25.

Xie et al. (2017) argue that central galaxies with underestimated SFRs are responsible for the decreasing trend of the sSFR for low stellar masses. In fact, from Fig. 7, we can clearly see that both central and satellite galaxies contribute to the decrease of sSFR with decreasing stellar mass. The responsible of this trend is the additional modulation with the virial velocity introduced in the estimation of reheated and ejected mass (equations 10 and 12), following the broken power law suggested by Muratov et al. (2015). The exponent α_F = −3.2 adopted for low-mass galaxies (V_vir < 60 km s⁻¹) makes the SN feedback too strong for them, which then have too low SFR. We have verified this by running a variant of model sag_β1.3 in which the slope α_F is fixed to −1 (the value corresponding to V_vir > 60 km s⁻¹), for any virial velocity. In this case, the main sequence of central galaxies become flatter at low stellar masses, while satellites present a mild trend in which sSFR increases monotonically with decreasing stellar mass.

The differences in normalization and shape of the main sequence of central and satellite galaxies are not explained by the particular treatment of environmental effects introduced in sag. Avoiding the action of RPS on the hot and cold gas phases only produces a higher number of active satellite galaxies, without changing the mean values of the sSFR of star-forming ones. In order to understand the origin of these differences, we examine the relationship between stellar mass and DM halo mass. We consider the mass of the main host haloes for central galaxies and the subhaloes masses for satellite galaxies, excluding orphans. Fig. 8 shows such relations, representing galaxies for the calibrated model sag with a coloured contour map on top of which the mean values of the stellar mass for each bin of halo mass (depicted by squares) are superimposed. Mean values of model sag_β1.3 are also added (represented by triangles), showing that the SHM relation is not appreciably affected by a change in the parameter β. Mean values are compared with the fitting function for the stellar-to-halo mass relation at z = 0 derived by Moster et al. (2010); such parametrization was estimated by requiring that the observed SMF obtained from SDSS DR3 (Panter, Heavens & Jimenez 2004) is reproduced when populating with galaxies the haloes and subhaloes in an N-body simulation. Note that this comparison is not fair for satellite galaxies, but it is done with the purpose of highlighting the differences in the relation for centrals and satellites. The general trend is very well reproduced by model central galaxies. For halo mass M_vir ≲ 10¹² M_⊙, differences between centrals and satellites become more evident. For a given stellar mass, central galaxies inhabit more massive haloes, while satellites prefer less massive ones, being the dispersion of the relation (vertical error bars) larger for smaller stellar masses. Systematic differences in the halo virial mass between main host and satellite haloes translate into differences in the hot halo mass and associated cooling rates, which impact on the SF activity leaving different imprints in the main sequence of central and satellite galaxies.

Figure 8.

Stellar mass as a function of DM halo mass for galaxies generated by both the calibrated model sag and model sag_β1.3 at z = 0 compared with the parametrization given by Moster et al. (2010). Top panel: Relation for central galaxies of main host haloes represented by a coloured contour map, considering the same contour levels as those of Fig. 4. Mean values of the stellar mass for each bin of halo mass for models sag and sag_β1.3 are depicted by red squares and blue triangles, respectively; vertical error bars are 1σ-standard deviation around the mean. Bottom panel: Same as top panel but for satellite galaxies within DM subhaloes. Differences between this relation for central and satellite galaxies translate into the different normalization and shape of their respective main sequences (see Fig. 7).

6.3 Specific star formation rate: distribution

We analyse the distribution of sSFR for central and satellite galaxies within different stellar mass ranges, discriminating the latter according to the mass of their main host DM haloes.⁵ We select galaxies within the same stellar and halo mass ranges as W12, who construct galaxy group catalogues using a group-finding algorithm based on the one presented by Yang et al. (2005). W12 derive stellar masses and sSFR for galaxies within groups and clusters from the spectroscopic information provided by the SDSS Data Release 7; the median redshift of the sample is z= 0.045. Left-hand, middle, and right-hand panels of Fig. 9 show, respectively, the sSFR distributions⁶ of galaxies with stellar masses within the ranges log(M_⋆[M_⊙]) ∈ [9.7, 10.1], log (M_⋆[M_⊙]) ∈ [10.1, 10.5] and log (M_⋆[M_⊙]) ∈ [10.5, 11.0]. Satellite galaxies are discriminated according to the mass of their main host haloes. In all cases, both central and satellite galaxies have a well distinguished peak at high values of sSFR (sSFR > 10⁻¹⁰ yr⁻¹). A bimodal distribution emerges only for the most massive galaxies considered. The break of the bimodality takes place at sSFR = 10^−10.7 yr⁻¹, a value in good agreement with the one adopted by Brown et al. (2017) and slightly higher than the one identified from the observational catalogue constructed by W12. While the sSFR distribution of galaxies in their observational sample exhibits a clear bimodality at all stellar masses, regardless of their classification in centrals or satellites (see their fig. 1), the authors clarify that the strong sharpness of the peak near sSFR = 10⁻¹² yr⁻¹ is a result of an artefact in spectral reductions. They also emphasize that the distributions exhibit a tail to much lower sSFR, like the sSFR distributions of our model galaxies do, which are in better agreement with the behaviour shown by Salim et al. (2007), as it was described in Section 6.2. Moreover, Kauffmann (2014) concludes that such a distribution is expected for galaxies characterized by a bursty mode of star formation with times of formation between 1 and 10 Gyr. This is consistent with the results we obtain from a detailed analysis of quenching time-scales of satellites (Cora et al. 2018). Consistent with the discussion also presented in Section 6.2 when analysing the main sequence of central and satellite galaxies at z = 0, the high sSFR peak of the distributions of satellite galaxies is slightly shifted towards lower sSFR values than for centrals, while such a trend is not evident in the observational sample.

Figure 9.

sSFR distributions for galaxies with stellar masses in different mass bins at z = 0: log (M_⋆[M_⊙]) ∈ [9.7, 10.1] (left-hand panel), log (M_⋆[M_⊙]) ∈ [10.1, 10.5] (middle panel), and log (M_⋆[M_⊙]) ∈ [10.5, 11.0] (right-hand panel). Coloured lines of different style represent satellite galaxies residing in main host DM haloes within different mass ranges, as indicated in the legend. Central galaxies are identified with black thin solid lines; all central galaxies are considered, without any restriction regarding their main host DM haloes. A bimodal distribution emerges only for the most massive galaxies considered, with a break at sSFR = 10^−10.7 yr⁻¹, a value in good agreement with the one adopted by Brown et al. (2017) and slightly higher than the one identified from the observational catalogue constructed by W12.

Despite these differences between our results and those of W12, there are several similarities in relevant aspects of the distributions. As in the observational sample, more massive central galaxies have larger probability of having low sSFR at the expense of a reduction in the peak height at high sSFR. However, satellites only show such a trend for those cases where the bimodality becomes evident in the model. Both in the simulated and observed sample, there are more satellites with low sSFR than centrals while the opposite occurs for high sSFR, regardless of the stellar mass bin considered. Besides, the number of satellite galaxies with low(high) sSFR increases(decreases) systematically as we consider more massive DM haloes. These general results remain valid for galaxies generated by model sag_β1.3. A quantitative comparison with observational results is attained through the estimation of the fraction of quenched galaxies for different stellar mass and halo massranges.

6.4 Quenched fractions of central and satellite galaxies from the model sag

From a given population of model galaxies selected according to their stellar mass and/or main host halo mass, we estimate the quenched fraction by considering those galaxies that are passive at a given time of interest, that is, those that satisfy the criterion sSFR < 10^−10.7 yr⁻¹, as emerges from the distribution of sSFR analysed in the previous section. We compare model results with those obtained by W12. Although these authors consider a galaxy as passive when its sSFR < 10⁻¹¹ yr⁻¹, we prefer to estimate the quenched fraction taking into account the separation between active and passive galaxies that emerge from the model at the expense of accepting a shift in the sSFR cut with respect to the one proposed by W12. Stellar mass and main host halo mass dependence of the quenched fraction at z = 0 (fq_z0) are shown in the left-hand and right-hand panels of Fig. 10, respectively. In both cases, the same mass bins as those chosen by W12 are used. Model results are identified with lines. Error bars show the 68 per cent Bayesian confidence interval estimated following Cameron (2011); their method is applied to estimate errors of fractions in all subsequent plots. The overwhelmingly numerous galaxy population obtained with the MDPL2 simulation makes a very good statistics in these results, so most of these error bars are too small to be distinguished. For the stellar mass dependence, we present the quenched fraction of both central (thin dashed line) and satellite (different line styles) galaxies; the latter are grouped according to the halo mass they reside in. For the halo mass dependence, only satellites within different mass ranges are considered. Symbols depict the results of W12 (see their figs 3a and b). From this comparison, we can see that the model underpredicts the fraction of quenched galaxies. Only the most massive galaxies, those in the stellar mass range log (M_⋆[M_⊙]) ∈ [10.9, 11.3], residing in haloes with masses M_halo ≳ 10¹⁴ M_⊙ are characterized by a quenched fraction similar to the observedone.

Figure 10.

Fraction of quenched galaxies (sSFR < 10^−10.7 yr⁻¹) in sag model at z= 0 as a function of stellar mass (top panel) and main host halo mass (bottom panel). Satellite galaxies are binned according to their main host halo or stellar mass, and are represented by coloured lines of different style, as indicated in the legend. Central galaxies, represented by a thin black dashed line, are only included in the top panel without making any distinction to the mass of the halo they reside in. Error bars show the 68 per cent Bayesian confidence interval estimated following Cameron (2011); they are hardly visible for the satellite population. These fractions are compared to those obtained by W12 (different symbols associated with different line styles). The calibrated sag model underpredicts the fraction of quenched galaxies. Only the most massive galaxies, those in the stellar mass range log (M_⋆[M_⊙]) ∈ [10.9, 11.3] within DM host haloes with masses M_halo ≳ 10¹⁴ M_⊙ achieve a quenched fraction similar to the observed one. This lack of general agreement is related to the underpredicted SFRD at high redshifts, as discussed in Section 6.1. The value β = 1.99 that regulates the redshift dependence of the reheated and ejected mass in sag model is too high and produces a shift of the SF activity to lower redshifts not leaving enough time for quenching.

This drawback of the model is the result of the high value of the parameter β achieved during the calibration process when trying to satisfy the constraint imposed by the SMF at z= 2, giving rise to an underprediction of the SFRD at high redshifts, as discussed in Section 6.1. Thus, the star formation activity is shifted to later epochs having less time to be quenched. We have also shown that fixing β in the lower value suggested by Muratov et al. (2015) allows to obtain a SFRD in better agreement with observations at high redshift, at the expense of having an excess in the number density of low-mass galaxies at z= 2. Results on the quenched fractions are in favour of choosing a model that prioritizes the right evolution of the SFRD. In the following, we show that quenched fractions are considerably improved for galaxies generated from model sag_β1.3.

6.5 Quenched fractions of central and satellite galaxies from the model sag_β1.3

Fig. 11 shows the stellar and main host halo mass dependence of the quenched fraction of galaxies generated by model sag_β1.3. The fractions of quenched galaxies at z = 0 are larger than those obtained with the calibrated model sag characterized by the parameter β = 1.99. The general agreement with values inferred by W12 is rather good, reproducing the observed trends for both central and satellite galaxies. This fact highlights the importance of an adequate efficiency of SN feedback at high redshifts in determining the passive fraction of galaxies. Despite this general good agreement, some differences in a certain range of stellar mass for satellites still remain. Namely, the quenched fractions of satellite galaxies with stellar masses log (M_⋆[M_⊙]) ∈ [10.1, 10.5] are clearly underpredicted for any halo mass considered. This mass range is responsible of the excess of the SFRD at z = 0 (see Fig. 6). However, satellites do not contribute significantly to the SFRD excess. This underprediction is caused by the lack of satellite galaxies in this stellar mass range as a result of overmerging, as it is evident from the depletion in the conditional SMF of satellites at log (M_⋆[h⁻² M_⊙]) ≈ 10 (see Fig. 5).

Figure 11.

Same as Fig. 10 for galaxies generated with model sag_β1.3. The general agreement with observations is rather good except for the underprediction of the quenched fraction for galaxies within the mass range log (M_⋆[M_⊙]) ∈ [10.1, 10.5].

As noted by W12, both central and satellite galaxies are more likely to be quenched as their stellar masses increase, being this dependence stronger for the former and mainly produced by self-regulating processes such as AGN and SN feedback, i.e. mass quenching (e.g. Peng et al. 2010; Henriques et al. 2017). Instead, satellites have a milder dependence of the quenched fraction on stellar mass, and its increment is more gradual for satellites hosted by more massive haloes, as a result of the larger effect of environmental processes on satellites with lower stellar masses. Hence, for a given stellar mass, the fraction of passive satellites increases with increasing halo mass, being the dependence on halo mass much milder for higher mass galaxies. This is more clearly shown in the right-hand panel of Fig. 11, where we can also appreciate that, for a given halo mass, more massive satellites are more likely to be quenched.

Fig. 12 shows how the fraction of z = 0 quenched satellite galaxies in model sag_β1.3 varies with the projected halo-centric distance, normalized by the virial radius of the main host DM halo in which the galaxies reside. Satellites are grouped according to their main host halo mass (top panel) or stellar mass (bottom panel), and their quenched fractions are depicted by different line styles. In the first case, only galaxies within the stellar mass range log (M_⋆[M_⊙]) ∈ [9.7, 10.5] are considered. In the second case, all galaxies within haloes with masses M_halo > 3.16 × 10¹² M_⊙ are taken into account. These selections allow to make a direct comparison with the quenched fractions of the sample analysed by W12, represented by different symbols according to the stellar or halo mass considered. In general terms, the radial gradient of the observed quenched fraction is well recovered by the model for both the stellar mass and halo mass selected sets.

Figure 12.

Fraction of quenched satellite galaxies (sSFR < 10^−10.7 yr⁻¹) in model sag_β1.3 at z = 0 as a function of the projected halo-centric distance normalized with the virial radius of the main host halo. Satellites are binned according to their main host halo mass (top panel) and stellar mass (bottom panel). Model results are shown by different lines; error bars show the 68 per cent Bayesian confidence interval estimated following Cameron (2011). Only satellites within the stellar mass range log (M_⋆[M_⊙]) ∈ [9.7, 10.5] are considered in the top panel and with main host halo mass log (M_halo[M_⊙]) > 12.5 in the bottom one. Such selection is done following W12 for comparison purposes. Their data are represented by different symbols according to the different stellar or halo mass ranges.

When the halo mass is varied (top panel), the radial profiles are a bit steeper than observed. The model gives a good match of fq_z0 for R_proj/r_vir ≲ 0.2 and R_proj/r_vir ≳ 0.5 for galaxies within haloes of mass log (M_halo[M_⊙]) ∈ [12.3, 13.2] and log (M_halo[M_⊙]) ∈ [13.2, 14.1], respectively. Quenched fractions for the highest halo mass bin (log (M_halo[M_⊙]) ∈ [14.1, 15.0]) are below W12 results for all halo-centric distances, and similar to those corresponding to the immediately smaller halo mass range. This result is consistent with the rather flat behaviour of the predicted quenched fractions as a function of halo mass for masses M_halo ≳ 4 × 10¹³ M_⊙, shown in the bottom panel of Fig. 11.

The discrepancies found when the halo mass is varied are reflected in the mismatch between model predictions and values inferred from observations for the two lowest stellar mass bins considered in the bottom panel of Fig. 12. As we can see, predicted values of fq_z0 are comprised within the range delimited by the corresponding observed quenched fractions, with galaxies within the stellar mass range log (M_⋆[M_⊙]) ∈ [10.1, 10.5] only achieving a good agreement near the cluster centre; this highlights once again the inability of the model in making adequate predictions for this particular mass range because of overmerger, as already noted. This might have affected the fraction of quenched low-mass satellites, since most of surviving orphans are those on less eccentric orbits that keep far from the central galaxy avoiding the merging region and preventing orphans from experiencing medium to strong-level RP.

On the other hand, the predicted quenched fractions of satellites within the two more massive stellar mass bins (log (M_⋆[M_⊙]) > 10.5) follow the observed radial profiles very well. Since SF in more massive galaxies is mainly suppressed by mass quenching processes at both low and high redshifts (Lin et al. 2014; Kawinwanichakij et al. 2017), both the radial and halo mass dependence of the fractions of quenched massive galaxies predicted by model sag_β1.3 support the modelling of those mass quenching processes that are particularly relevant for galaxies with stellar masses log (M_⋆[M_⊙]) > 10.5, such as AGN feedback and disc instabilities. While the former reduces the amount of cooled gas, the latter produces starbursts leading to a rapid exhaustion of the cold gas reservoir. These self-regulating physical processes explain the SF quenching in massive central galaxies but are still active as quenching mechanisms even when galaxies become satellites, as demonstrated by Peng et al. (2012), playing a dominant role with respect to environmental effects Cora et al. (2018). This does not mean that environmental processes are irrelevant for high-mass galaxies. Indeed, these galaxies suffer the effects of RPS but they are smaller than those experienced by low-mass satellites, as shown in the next Section.

7 ROLE OF ENVIRONMENT: STRIPPED MASS AND ATOMIC GAS CONTENT

In order to understand the role of environment on satellite galaxies of different mass, we show in Fig. 13 mean values of the fraction of hot gas mass stripped by RP, f_strip,hot, at each snapshot of the simulation as a function of redshift, for model sag_β1.3. This quantity is estimated considering the ratio between the stripped mass and the mass of hot gas just before being stripped; the latter is given by the sum of the stripped mass and the mass of hot gas at the redshift considered. Satellites are grouped according to their current stellar mass (different lines) and the virial mass of their current main host DM haloes (different panels). As expected, RPS is more efficient in removing hot gas from galaxies with smaller stellar masses as demonstrated from both the analysis of observational data (e.g. Fillingham et al. 2016) and hydrodynamical simulations (Bahé & McCarthy 2015). For any galaxy mass and halo mass, the fraction of hot gas mass stripped by RP increases with decreasing redshift. This is explained by the mass growth of the host DM halo and the consequent higher densities in the ICM, and high-redshift infalling galaxies that have more time to approach the cluster core, thus suffering stronger RP (Jaffé et al. 2015, 2016; Vega-Martínez et al., in preparation). Local values of the removed fraction at each stripping episode are quite similar for any halo mass, being of the order of ≈0.03−0.05, with lower values corresponding to lower mass haloes. The dependence of this fraction with halo mass becomes more evident at higher redshifts; the effect of RP is more pronounced in clusters of M_vir ≈ 10¹⁵ M_⊙ since z = 1.5 than in less massive ones, consistent with the distribution of RP values presented by Tecce et al. (2010, see their fig. 5).

Figure 13.

Mean values of the fraction of hot gas mass stripped by RP at each snapshot of the simulation as a function of redshift for model sag_β1.3. All satellites (passive and active) at z= 0 are included. They are grouped according to their local stellar mass (different line styles) and main host halo mass (different panels), as indicated in the legends. For any galaxy mass and halo mass, the fraction of hot gas mass stripped by RP increases with decreasing redshift.

Fig. 14 shows the mean values of the cumulative stripped hot gas fraction as a function of stellar mass of satellite galaxies binned according to the mass of the main host halo they inhabit. This fraction is estimated as the ratio between the accumulated stripped hot gas mass since the first RP stripping event suffered by the satellite and the sum of its current hot gas mass and the accumulated stripped hot gas. Again, we see that the mass stripped by RP is larger for less massive satellites residing in more massive haloes. For satellites within low-mass haloes (log (M_halo[M_⊙]) ∈ [12.3, 13.2]), the cumulative stripped hot gas fractions decrease from ≈50 per cent for M_⋆ ≈ 3 × 10⁹ M_⊙ to ≈20 per cent for M_⋆ ≈ 10¹¹ M_⊙. For high-mass haloes (log (M_halo[M_⊙]) ∈ [14.1, 15.0]), these fractions increase to ≈70 and ≈40 per cent, respectively.

Figure 14.

Mean values of the fraction of accumulated stripped hot gas by RP as a function of stellar mass for model sag_β1.3. All satellites (passive and active) at z= 0 are included. They are grouped according to their local main host halo mass (different line styles), as indicated in the legend. l

In our model, both RP and TS act on the hot gas halo, but in general the former process gives rise to a smaller stripping radius than the latter, indicating that RPS dominates TS. Thus, TS is considered as a secondary effect (Font et al. 2008; McCarthy et al. 2008; Bahé & McCarthy 2015). This is also consistent with inferences from observational results, like the efficient SF quenching detected in Virgo cluster (Boselli et al. 2016).

We demonstrate that low-mass satellites are the galaxies mainly affected by the way in which environmental processes regulate the content of the hot gas reservoir by making a test in which satellites are allowed to keep their hot gas halo at infall but neither RPS nor TS are activated, so that the reduction of hot gas is a result of ejection and/or gas cooling. Only low-mass satellites (log (M_⋆[M_⊙]) ∈ [9.5, 10.5]) experience a reduction in the quenched fraction which is as high as ≈0.15 for galaxies residing within the most massive haloes (log (M_halo[M_⊙]) ∈ [14.1, 15.0]). However, the quenched fractions of more massive satellites remain unchanged. The lack of significant change in the fraction of quenched massive satellites when their hot gas reservoir suffers all type of possible effects (strangulation, gradual removal, suppression of environmental processes) is a strong evidence that environmental processes do not play a significant role on the SF activity in massive satellites although these galaxies do suffer them as shown in Figs. 13 and 14. This was already deduced from the good match between model results and observations for the radial distribution of the fraction of quenched massive galaxies (see Fig. 12).

The cold gas reservoir of satellite galaxies can also be reduced by the action of RPS, as inferred from images of highly asymmetric H i distribution well within the stellar disc in cluster spirals (e.g. Abramson et al. 2011; Kenney, Abramson & Bravo-Alfaro 2015; Bellhouse et al. 2017). Fig. 15 shows the mean values of the cumulative stripped cold gas fraction as a function of stellar mass of satellite galaxies in model sag_β1.3 which interstellar medium is not longer shielded by the hot gas halo. They are grouped according to the mass of their main host haloes. The dependences of these fractions with stellar mass and halo mass are quite similar to those characterizing the stripped hot gas (Fig. 13), although the values of the fractions of cumulative stripped cold gas are one order of magnitude smaller than for the hot gas.

Figure 15.

Mean values of the fraction of accumulated stripped cold gas by RP as a function of stellar mass for satellites in model sag_β1.3 which interstellar medium is not longer shielded by the hot gas halo. They are grouped according to the mass of their local main host haloes (different line styles), as indicated in the legend.

The small fraction of cold gas removed by RPS has no impact on the fractions of quenched galaxies, which remain unchanged if RPS of cold gas is suppressed. However, this slight reduction of the cold gas content allows to obtain a better agreement with the observed relation of the atomic hydrogen (H i) gas content as a function of stellar mass than in the case where RPS of the cold disc is deactivated. This good agreement is shown in Fig.16 for satellites within main host haloes of different mass in model sag_β1.3. H i mass is estimated assuming that total cold gas is the sum of the contribution of molecular (H2) and atomic hydrogen content and 30 per cent helium (Boselli et al. 2014, equation 2). The molecular hydrogen content is estimated adopting the scaling relation of molecular gas-to-stellar mass ratio (Boselli et al. 2014, table 3), in which the coefficients correspond to molecular hydrogen masses estimated from CO intensities using the H-band luminosity-dependent CO-to-H₂ conversion factor. We compare against measurements by Brown et al. (2017) from a spectral stacking technique applied to a multiwavelength sample of satellite galaxies selected from the SDSS, with H i data from ALFALFA survey, with stellar masses larger than 10⁹ M_⊙ (different symbols correspond to averagevalues for different bins of halo mass). The match between model results and observations is good for satellites with stellar mass M_⋆ ≳ 10¹⁰ M_⊙ residing in haloes more massive than ≈10¹² M_⊙, consistent with the good behaviour of the radial dependence of the quenched fraction for this galaxy population, giving additional support to the modelling of mass quenching processes, as discussed in Section6.5. For smaller stellar masses, the mean H i-to-stellar mass fraction is underestimated by sag_β1.3 regardless of the host halo mass. This is not related to the modelling of RPS since this general trend persists even when the action of RPS on the cold gas is deactivated, in which case there is a systematic increase of this fraction by ≈0.25 dex for all stellar masses. Hence, the low-mass content of atomic gas might be explained by levels of star formation higher than expected in low-mass galaxies as becomes evident from the excess of the cosmic SFRD at z = 0 (see Fig. 6).

Figure 16.

Atomic hydrogen gas fraction as a function of stellar mass for satellite galaxies in model sag_β1.3 within different ranges of main host halo mass, as indicated in the legend. Different lines and associated shaded areas represent the corresponding median values and 10 and 90 percentiles, respectively. Median values are compared with the average atomic gas fraction presented by Brown et al. (2017), who consider H i data from ALFALFA survey (different symbols). The match between model results and observations is rather good for satellites with stellar mass larger than 10¹⁰ M_⊙ residing in haloes more massive than ≈10¹² M_⊙, consistent with the good behaviour of the radial dependence of the fraction of quenched massive satellites.

Results from sag_β1.3 regarding the levels of current SF and atomic gas content call for a more refined treatment of quiescent SF, like distinguishing the neutral and molecular phases of the cold gas as (e.g. Lagos et al. 2011). Molecular gas is not significantly affected by RPS as arises from both observations (e.g. Abramson et al. 2011) and hydrodynamical simulations (Tonnesen & Bryan 2009); this cold gas phase could be highly perturbed by high-speed tidal interactions (Scott et al. 2015).

8 COMPARISON WITH OTHER SAMS

In a recent work, Henriques et al. (2017) study environmental and mass quenching using the model described in Henriques et al. (2015), which is an updated version of the Munich semi-analytic model l-galaxies (Guo et al. 2011). Since we adopt the prescription used by Henriques et al. (2017) to model AGN feedback, the discrepancies between the predicted fractions of quenched galaxies from l-galaxies and those inferred from observations seem to originate in the modelling of other baryonic (e.g. SN feedback, SF, recycling) and environmental processes, which affect the complex baryon cycle.

The new features implemented in sag, including the new SN feedback scheme and treatment of environmental processes, differ considerably with respect to the corresponding implementation in l-galaxies. Regarding gas disruption processes, the most relevant difference is that sag allows gradual removal of hot gas reservoir by TS and RPS in both types of satellite galaxies (those that keep their DM subhaloes and orphans), while in l-galaxies orphan galaxies lose immediately their hot gas reservoir by TS, as this process is considered responsible of the disruption of their DM subhaloes (in general, hot gas is assumed to be tidally stripped at the same rate as dark matter). Their assumption might be too strong, not taking into account the fact that the mass resolution of the simulation prevents subhalo detection in some cases. Therefore, the assumption adopted in l-galaxies prevents orphans from the gradual removal of hot gas through RPS. The action of RPS is only allowed in satellites that keep their substructures. As described in Guo et al. (2011), the stripping radius is estimated from the balance of self-gravity and ram pressure. After the stripping event, the remaining hot gas expands out to the subhalo r_vir recorded at infall. Guo et al. (2011) argue that feedback could provide the required energy to redistribute the gas in this manner, which would require outflow velocities comparable to the circular velocity of the host group. This assumption on the distribution of the remaining hot gas results in an increased efficiency of RPS compared to our procedure, since there is always hot gas beyond the current stripping radius. Besides, the halo becomes increasingly more diffuse which leads to reduced cooling rates. Therefore, their model results in an increased fraction of passive low-mass satellites (Guo et al. 2011; Henriques et al. 2013; Hirschmann et al. 2014).

With the aim of reducing the excess of passive satellites predicted by l-galaxies, Henriques et al. (2015) define a minimum threshold halo mass (10¹⁴ M_⊙) for RP to be effective as an stripping mechanism of the hot gas halo of satellites that keep their substructure, below which the action of RP is suppressed. However, despite of this additional condition, discrepancies between predicted and observed fractions of quenched galaxies still remain. Their condition for the action of RP not only differs from the one we adopt (equation 3), which action is self-regulated without any restriction, but it also considers that the orbital velocity of the satellite is approximated by the virial circular velocity of the main halo. This aspect is improved in our implementation by taking the information of the actual orbits of satellites, including those of orphans. Furthermore, model l-galaxies does not include RPS of the cold gas disc.

Gonzalez-Perez et al. (2014) present a variant of galform semi-analytic model with a molecular-based SF law that also assumes gradual stripping of the hot gas by the action of RP adopting the same parametrization for the gravitational restoring as in sag (see equation 3) but with the value α_RP = 2 suggested by McCarthy et al. (2008). The main difference with our implementation is that RP is set to its maximum value for the whole orbit of the satellite galaxy since it is calculated at its pericentre. Therefore, the stripped hot gas is overestimated leading to lower fractions of atomic gas than those inferred from observational data (see fig. 7b of Brown et al. 2017). This disagreement would become even larger if RPS of cold gas were included in their model.

The impact of both cold gas and hot gas stripping on the H i fractions of local galaxies is investigated by Stevens & Brown (2017) using an updated version of the dark sage semi-analytic model, which evolves the one-dimensional structure of galactic discs in annuli of fixed specific angular momentum. Quiescent SF depends on the molecular gas fraction of each annulus; they consider both metallicity- and pressure-based prescriptions for determining the ratio of H₂. The action of RP on both hot gas and cold gas phases is regulated by conditions similar to those applied in our model sag, that is, those given by McCarthy et al. (2008) and Gunn & Gott (1972), respectively. An additional important similarity to note is that dark sage also considers that the hot gas halo protects the cold gas disc from the action of RP. However, the condition adopted to allow RPS differs from ours, since RP is able to remove cold gas from the disc when the total baryonic mass of the galaxy (cold gas and stars) exceeds the mass of hot gas halo. Stevens & Brown (2017) show that the satellite H i fractions are underpredicted by their full model even when it has been calibrated using the observed H i fraction of galaxies as a constraint with the greatest weight. These fractions become closer to the observational data when cold gas stripping is suppressed. Among the possible explanations of these results offered by the authors, there is one directly related to the limitations of SAMs, namely, the fact that the hot gas of satellite galaxies is not replenished by accretion from the ICM.

The effect of RPS of cold gas is also taken into account by Luo et al. (2016) in a branch of l-galaxies (Fu et al. 2013) which also assumes galaxy discs divided into multiple rings, and an SFR dependent on the local surface density of the molecular gas. This model is unable to reproduce the observed trends of quenched fractions of satellites and centrals. The number of low-mass passive galaxies are overpredicted despite gradual stripping of the hot gas is considered. The action of RPS of the cold gas increases the discrepancies with observational data. However, it is difficult to quantify the effect of this environmental process because the passive fraction of central galaxies is underpredicted for any stellar mass hinting to problems not related with environment.

Even with the aforementioned drawback affecting the hot gas halo of satellites and a simple modelling of galaxy discs and quiescent star formation, our model sag_β1.3 is able to achieve better agreement with observed atomic gas content and quiescent fractions than other current SAMs.

9 SUMMARY AND CONCLUSIONS

The updated version of our semi-analytic model of galaxy formation sag has been used to generate one of the galaxy catalogues of the MultiDark Galaxies project (Knebe et al. 2018a), which is based on the Planck cosmology 1 h⁻¹ Gpc MultiDark simulation MDPL2. The model has been calibrated using the PSO technique (Ruiz et al. 2015) combined with a particular set of constraints and the observational data defined in Knebe et al. (2018b). This catalogue is publicly available in the CosmoSim data base. We also consider a second galaxy population generated by a slightly different version of the model, characterized by a change in one of the parameters involved in the modelling of SN feedback (model sag_β1.3). We analyse the properties of the galaxies in these two catalogues, both of star-forming and quiescent galaxies, splitting the population in centrals and satellites, and also sampling them according to the mass of the main host haloes they inhabit. The detailed analysis performed allowed us to evaluate the impact of the improved treatment of environmental effects implemented in the model, and the modification introduced in the SN feedback scheme, given by the addition of an explicit redshift dependence in the estimation of the reheated and ejected mass.

We implement the gradual removal of hot gas in satellites from RPS and TS, allowing also the action of these processes on the cold gas disc under certain conditions. This is an improvement of our model with respect to previous studies that ignore the RPS of the cold gas (Guo et al. 2011; Kimm et al. 2011; Gonzalez-Perez et al. 2014; Henriques et al. 2017). The advantage of our implementation of RPS with respect to previous works (e.g. Henriques et al. 2017; Stevens & Brown 2017) resides in the use of a fitting formulae for RP experienced by galaxies in haloes of different mass as a function of halo-centric distance and redshift instead of analytic estimations. The effect of TS is considered as an additional mechanism contributing to the removal of the hot and cold gas. It can also affect the stellar components. Another important improvement with respect to other SAMs is the integration of orphan satellites (Vega-Martínez et al., in preparation).

In the following, we summarize the main results of this work obtained from the analysis of the two galaxy catalogues generated with the calibrated model sag and model sag_β1.3:

• SN feedback plays an important role in the star formation history of galaxies. In order to avoid the excess of the faint end of the SMF at z = 2 and to give rise to downsizing in the stellar mass assembly, it is necessary to reduce the availability of cold gas for star formation at high redshifts. Thus, both the reheating of the cold gas phase and the ejection of part of the hot gas reservoir, with the possibility of being reincorporated later, must be stronger at earlier epochs. This has been demonstrated by the inclusion of an explicit redshift dependence of both the reheated and ejected mass in the new feedback scheme implemented in sag. Results are sensitive to the value of the power-law slope of such dependence, denoted by the parameter β in our model. When β = 1.99, as emerged from the calibration process in which the free parameters are restricted by imposing the SMF at z = 2 as constraint, the predicted SFRD becomes too low with respect to observational data at high redshifts. When fixing that parameter in the value suggested by the fit of Muratov et al. (2015) obtained from the FIRE hydrodynamical simulations (β = 1.3), the trend of SFRD is reconciled with the observed one (Fig. 6) at the expense of an excess in the faint end of the SMF at z = 2 (Fig. 1).

• The gradual removal of hot halo gas through a robust treatment of environmental effects, as the one implemented in our model, becomes crucial to produce main sequences that evolve consistently with observational data for both central and satellite galaxies (Fig. 7). Cooling flows keep replenishing the cold gas reservoir of satellites after first infall during gradual starvation, delaying the beginning of the quenching of SF. The differences in normalization and shape of the main sequence of central and satellite galaxies that still remain arise because, for a given stellar mass, central galaxies inhabit more massive haloes than satellites (Fig. 8), which impact on the mass of the hot halo and associated cooling rates. The fact that galaxies start to form stars earlier in model sag_β1.3 is reflected in the lower normalization of the main sequence of star-forming galaxies achieved by this model with respect to model sag. In both models, the main sequences of central and satellite galaxies are characterized by a decreasing trend of sSFR with decreasing stellar mass for M_⋆ ≲ 10¹⁰ M_⊙. The strong dependence with the virial velocity for low-mass galaxies introduced in the estimation of the reheated and ejected mass is responsible for this effect.

• The higher efficiency of SN feedback at high redshifts, required to recover the downsizing in stellar mass assembly, must be characterized by a mild redshift dependence of the reheated and ejected mass, as in model sag_β1.3, in order to achieve the expected behaviour of the fractions of quenched galaxies (those with sSFR < 10^10.7 yr⁻¹) as a function of stellar mass, halo mass, and the halo-centric distances (Figs. 11 and 12). The right fractions of quenched massive galaxies (log (M_⋆[M_⊙]) > 10.5) as a function of both halo mass and halo-centric distance support the modelling of those self-regulating physical processes that are particularly relevant for massive galaxies (AGN feedback; disc instabilities), which are still active as quenching mechanisms even when galaxies become satellites, dominating over environmental quenching; the effect exerted by the latter is milder in high-mass galaxies than in low-mass ones.

• RPS plays a dominant role among the environmental processes considered in our model and contributes to regulate adequately the mass of the hot gas halo and the cold gas disc (when it is no longer shielded by the hot halo), being more efficient in removing gas from galaxies with smaller stellar masses, as expected (Figs. 13, 14, and 15). For any galaxy mass, the fraction of hot and cold gas mass stripped by RP increases with decreasing redshift, as a result of higher densities achieved by the ICM, associated with the mass growth of the main host DM haloes, and the fact that high redshift infalling galaxies have more time to approach the cluster core where they suffer strong RP. The total stripped mass is higher in more massive clusters (M_vir ≈ 10¹⁵ M_⊙) than in less massive ones, consistent with the distribution of RP values inherent to groups and clusters (Tecce et al. 2010). The cumulative stripped hot gas fraction increases from ≈20 to ≈50 per cent for satellites with M_⋆ ≈ 3 × 10⁹ M_⊙ and M_⋆ ≈ 10¹¹ M_⊙, respectively, residing in low-mass haloes (log (M_halo[M_⊙]) ∈ [12.3, 13.2]). For high-mass haloes (log (M_halo[M_⊙]) ∈ [14, 115.]), these fractions increase to ≈40 and ≈70 per cent, respectively.

• The cold gas disc is less affected by RP than the hot gas, with cumulative stripped cold gas fractions being one order of magnitude smaller for the former than for the latter. Thus, RPS on the cold gas does not affect the fraction of quenched galaxies but it definitely contributes to reach the right atomic hydrogen (H i) gas content for satellites with different stellar mass, especially for more massive ones (M_⋆ ≳ 10¹⁰ M_⊙) residing in haloes with masses M_vir ≳ 10¹² M_⊙ (Fig. 16). Distinguishing the neutral and molecular phases of the cold gas could help to reproduce the observed levels of H i gas content of low-mass satellites.

Our results highlight the impact of specific aspects of mass and environmental quenching on galaxy evolution. The stronger effect of SN feedback at higher redshifts and the physics of environmental processes captured by the RP fitting formulae combined with the orbital evolution of orphan galaxies make our model sag_β1.3 able to generate a galaxy population with fractions of gas and levels of SF that allow to achieve better agreement with observed atomic gas content and quiescent fractions than other current SAMs. These attainments of our model are reached even though it does not include a detailed treatment of galaxy discs as in dark sage(Stevens & Brown 2017), or a molecular-based SF law as galform (Gonzalez-Perez et al. 2014). Differences still found for particular ranges of stellar and halo mass give hints for further improvement of themodel.

Footnotes

1

http://dx.doi.org/10.17876/cosmosim/mdpl2/007

2

https://www.cosmosim.org

3

The process of removal of the hot gas haloes of satellite galaxies is usually indistinctly called ‘strangulation’ or ‘starvation’, regardless of which physical mechanism may actually be responsible for it. In this study, we will use the term ‘strangulation’ to refer to the instantaneous removal of hot gas from satellites, and we will use ‘gradual starvation’ to refer to the gradual removal of halo gas by any process.

4

We refer to the hot gas associated with a main host halo using the term ‘intracluster medium’, even for host halo masses smaller than that of a galaxy cluster.

5

From hereafter, all the sampling of satellite galaxies is done according to the masses of their main host haloes, even if it is simply referred as halo.

6

Probability density function at the bin, normalized such that the integral over the range is equal to unity.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, http://www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, http://www.lrz.de). The MDPL2 simulation has been performed under grant pr87yi. Our collaboration has been supported by the DFG grant GO 563/24-1. The CosmoSim database used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (IAP). We acknowledge in particular the support from the eSciene team at the AIP, Kristin Riebe, Harry Enke and Anastasia Galkin. This work was done in part using the Geryon computer at the Center for Astro-Engineering UC, part of the BASAL PFB-06, which received additional funding from QUIMAL 130008 and Fondequip AIC-57 for upgrades. We thank the referee for useful comments and suggestions that have contributed to improve this work. We acknowledge Tomás E. Tecce for his valuable contribution in initial stages of the development of the current version of sag and Mario G. Abadi for useful discussion. SAC acknowledges funding from Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, PIP-0387), Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT, PICT-2013-0317), and Universidad Nacional de La Plata (G11-124), Argentina. CAVM, TH, FC, and IDG acknowledge CONICET, Argentina, for their supporting fellowships. ANR acknowledges funding from ANPCyT (PICT-2014-2862) and from Secretaría de Ciencia y Tecnología de la Universidad Nacional de Córdoba (PID 30720450100484). AO acknowledges support from project AYA2015-66211-C2-2 of the Spanish Ministerio de Economia, Industria y Competitividad. AMMA acknowledges support from CONICYT-PCHA/Doctorado Nacional 2011-21110870, BASAL PFB-06 and FONDECYT grant 3160776. GY acknowledges financial support from the Ministerio de Economía y Competitividad and the Fondo Europeo de Desarrollo Regional (MINECO/FEDER, UE) in Spain through grant AYA2015-63810-P.

APPENDIX A: RAM PRESSURE STRIPPING RADIUS OF THE HOT GAS

Assuming that the DM and hot gas in a satellite distribute following an isothermal sphere density profile, the total mass of the satellite from equation (4) can be written as

\begin{equation*} M_\text{sat}(r_\text{sat}) = M_\text{gx} + \left( \frac{M_\text{hot}}{r_\text{hot}} + \frac{M_\text{DM}}{r_\text{DM}} \right) r_\text{sat}, \end{equation*}

(A1)

where we have defined M_gx ≡ M_* + M_cold and r_sat is the current radius of the satellite. We assume that r_DM has already been determined by the effect of TS. Using equation (A1) in the RPS condition for a spherically symmetrical distribution of gas from McCarthy et al. (2008, equation 3 in this paper) and rearranging terms a bit, we can obtain

\begin{equation*} P_\text{ram} = \frac{\alpha _\text{RP} M_\text{hot}}{4\pi r_\text{hot} r_\text{sat}^3} \left[ G M_\text{gx} + \left( \frac{G M_\text{hot}}{r_\text{hot}} + \frac{G M_\text{DM}}{r_\text{DM}} \right) r_\text{sat} \right]. \end{equation*}

(A2)

The equality actually holds when r_sat equals the gas stripping radius r_s, which represents the radius |$r_{\rm s,hot}^{\rm RPS}$| defined in Section 3.4. Notingthat terms of the form GM/r have units of velocity squared, we can define

\begin{equation*} V_h^2 \equiv \frac{G M_\text{hot}}{r_\text{hot}} + \frac{G M_\text{DM}}{r_\text{DM}} \equiv V_\text{hot}^2 + V_\text{DM}^2 \end{equation*}

(A3)

and also

\begin{equation*} V_\text{gx}^2 \equiv G M_\text{gx} / r_\text{DM}. \end{equation*}

(A4)

We then have

\begin{equation*} P_\text{ram} r_{\rm s}^3 = \frac{\alpha _\text{RP} M_\text{hot}}{4\pi r_\text{hot}} \left[ V_\text{gx}^2 r_\text{DM} + V_h^2 r_{\rm s} \right]. \end{equation*}

(A5)

Rearranging terms once again, we obtain that the stripping radius r_s must satisfy the following condition:

\begin{equation*} r_{\rm s}^3 - \frac{\alpha _\text{RP} M_\text{hot}}{4\pi r_\text{hot} P_\text{ram}} \left[ V_\text{gx}^2 r_\text{DM} + V_h^2 r_{\rm s} \right] = 0. \end{equation*}

(A6)

This is a third-order polynomial |$a r_{\rm s}^3 + b r_{\rm s}^2 + c r_{\rm s} + d$|⁠, with a = 1, b = 0 and both c and d are always negative. A polynomial of this form has only one real and positive root, which corresponds to the RP stripping radius; equation (A6) must be solved numerically to obtain r_s.