# The new semi-analytic code GalICS 2.0 – reproducing the galaxy stellar mass function and the Tully–Fisher relation simultaneously

A. Cattaneo, J. Blaizot, J. Devriendt, A. Mamon, E. Tollet, A. Dekel, B.

Guiderdoni, M. Kucukbas, A. Thob

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

A. Cattaneo, J. Blaizot, J. Devriendt, A. Mamon, E. Tollet, et al.. The new semi-analytic code GalICS 2.0 – reproducing the galaxy stellar mass function and the Tully–Fisher relation simultaneously. Monthly Notices of the Royal Astronomical Society, 2017, 471 (2), pp.1401-1427. 10.1093/mn-ras/stx1597 . hal-02143138

HAL Id: hal-02143138

https://hal.science/hal-02143138v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# The new semi-analytic code GalICS 2.0 – reproducing the galaxy stellar mass function and the Tully–Fisher relation simultaneously

A. Cattaneo,$^{1,2*}$ J. Blaizot,$^3$ J. E. G. Devriendt,$^4$ G. A. Mamon,$^2$ E. Tollet,$^{1,5}$ A. Dekel,$^6$
B. Guiderdoni,$^3$ M. Kucukbas$^{1,5}$ and A. C. R. Thob$^7$

Accepted 2017 June 22. Received 2017 June 21; in original form 2016 October 18

## ABSTRACT

GalICS 2.0 is a new semi-analytic code to model the formation and evolution of galaxies in a cosmological context. N-body simulations based on a Planck cosmology are used to construct halo merger trees, track subhaloes, compute spins and measure concentrations. The accretion of gas on to galaxies and the morphological evolution of galaxies are modelled with prescriptions derived from hydrodynamic simulations. Star formation and stellar feedback are described with phenomenological models (as in other semi-analytic codes). GalICS 2.0 computes rotation speeds from the gravitational potential of the dark matter, the disc and the central bulge. As the rotation speed depends not only on the virial velocity but also on the ratio of baryons to dark matter within a galaxy, our calculation predicts a different Tully–Fisher relation from models in which $v_{\text{rot}} \propto v_{\text{vir}}$. This is why, GalICS 2.0 is able to reproduce the galaxy stellar mass function and the Tully–Fisher relation simultaneously. Our results are also in agreement with halo masses from weak lensing and satellite kinematics, gas fractions, the relation between star formation rate (SFR) and stellar mass, the evolution of the cosmic SFR density, bulge-to-disc ratios, disc sizes and the Faber–Jackson relation.

**Key words:** galaxies: evolution – galaxies: formation – galaxies: kinematics and dynamics – galaxies: luminosity function, mass function – galaxies: spiral – galaxies: structure.

## 1 INTRODUCTION

Semi-analytic models (SAMs) are a technique to model the formation and evolution of galaxies in a cosmological context. Pioneered by White & Frenk (1991) and Lacey & Cole (1993), this technique is based on the notion that galaxy formation is a two-stage process (White & Rees 1978). The gravitational instability of primordial density fluctuations in the dark matter (DM) forms haloes. The dissipative infall of gas within haloes forms luminous galaxies. SAMs follow these two stages separately. First, one constructs merger trees for the haloes in a representative cosmic volume. Then, the evolution of baryons within haloes is broken down into a number of elementary processes, which are modelled analytically.

This paper introduces the new SAM GalICS 2.0. An early version had already been presented in a comparison of all the main SAMs (Knebe et al. 2015). The models that participated in this comparison

are those by Bower et al. (2006), Font et al. (2008), Gonzalez-Perez et al. (2014), Croton et al. (2006), De Lucia & Blaizot (2007), Henriques et al. (2013), Benson (2012), Monaco, Fontanot & Taffoni (2007), Gargiulo et al. (2015), Somerville et al. (2008) and Lee & Yi (2013).

GalICS 2.0 builds on our previous experience with GalICS (Hatton et al. 2003; Cattaneo et al. 2006, 2008, 2013) but is more than a new version. The entire code has been rewritten from scratch. One of the reasons is to enable a more extensive use of the cosmological N-body simulation used to construct the merger trees. In GalICS 2.0, we use the information on DM substructures (merger rates are more accurate) and the density profiles of DM haloes (we can compute realistic rotation curves). Other advantages, besides an improved description of several physical processes, are a massive gain in computational speed and far greater modularity.

SAMs need complex baryon physics [radiative cooling, shock heating, active galactic nuclei, supernovae (SNe)] to explain why the galaxy stellar mass function (SMF) has a knee at $M_{\text{stars}} \sim 10^{11} M_{\odot}$ when the mass function of their host haloes is essentially a single

<footnote>$^1$Observatoire de Paris, GEPI, 77 av. Denfert-Rocherau, F-75014 Paris, France</footnote>

<footnote>$^2$Institut d'Astrophysique de Paris, 98bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote>$^3$Centre de Recherche Astrophysique de Lyon, 9 avenue Charles André, F-69561 Saint-Genis-Laval, France</footnote>

<footnote>$^4$Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>$^5$Université Paris VII, Sorbonne Paris Cité, 5 rue Thomas-Mann, F-75205 Paris Cedex 13, France</footnote>

<footnote>$^6$Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel</footnote>

<footnote>$^7$Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK</footnote>

<footnote>* E-mail: andrea.cattaneo@obspm.fr</footnote>
power law up to $M_{vir} \sim 10^{13} M_{\odot}$. Yet the assumption that the growth histories of DM haloes determine the properties of galaxies underpins the entire semi-analytic approach.

There are three main ways to measure halo masses and probe the galaxy–halo connection directly: from rotation curves, from satellite kinematics and from weak lensing data. The first method is the oldest. In fact, it is one of the ways DM was discovered. Hence, the importance of the relation between stellar mass and disc rotation speed (Tully–Fisher relation, TFR) as a key test for SAMs. However, reproducing the TFR and the SMF simultaneously has been a main challenge for SAMs since their inception. Either models are calibrated on the TFR and fail to fit the SMF/luminosity function (Kauffmann, White & Guiderdoni 1993) or they are calibrated on the SMF/luminosity function and fail to reproduce the TFR (Heyl, Hernquist & Spergel 1995). The discrepancy persists today, albeit to a lesser extent (e.g. Guo et al. 2011, also see Guo et al. 2010, although the latter is based on abundance matching rather than semi-analytic modelling). The need to compare the predictions of SAMs to direct probes of halo masses has played a major role in the development of GalICS 2.0. In this paper, we show that modelling the rotation curves of disc galaxies accurately is necessary for a meaningful comparison with Tully–Fisher data.

The structure of the paper is as follows. In Section 2, we describe GalICS 2.0 and we explain our strategy to set the values of the model parameters. In Section 3, we compare its predictions with the observations [SMFs, baryonic mass function, halo masses from weak lensing and satellite kinematics, relation of star formation rate (SFR) to gas and stellar mass, SFR function, evolution of the cosmic SFR density, gas fractions, bulge-to-disc ratios, disc sizes, stellar and baryonic TFR at $z \sim 0$, stellar TFR at $z \sim 1$, the relation between disc rotation speed and virial velocity from rotation-curve studies, and the Faber–Jackson relation]. Section 4 summarizes the conclusions of the paper.

## 2 THE MODEL

Luminous matter is composed of galaxies and the intergalactic medium. We could also say it is composed of gas and stars. The objects we use to describe the Universe depend on the scale we are looking at. GalICS 2.0 works the way we think. Different modules follow different objects, which capture the formation and evolution of galaxies on different scales. Each is written to be as self-contained as possible.

On the largest scale, the TREE module follows the hierarchical formation and merging of DM haloes. For TREE, a halo is just a point in a network of relationships (progenitor, descendant, host, subhalo). The flow of baryons in and out of haloes is followed by the HALO module. HALO follows the exchanges of matter between the cold gas, the hot gas and the galaxy (e.g. the rate at which gas accretes on to the galaxy) but not the galaxy's internal structure. The decision of what goes to the disc and what goes to the bulge is done in the GALAXY module, which computes all morphological, structural and kinematic properties. The relation between HALO and GALAXY can be compared to that between a mill and a baker. There is an exchange of matter both ways (inflows and outflows, flour and money) but the baker does not need to know if the flour has been ground with a water or a wind mill. Neither does the millman about the baker's recipes. This philosophy explains some practical choices, such as that of the time substeps in Section 3.4. A galaxy contains different components, such as a disc, a bulge or a bar (which we classify as a pseudo-bulge), but star formation and feedback

within a component are followed in the COMPONENT module. The lowest scale corresponds to the STAR (stellar evolution) and GAS (interstellar medium) modules.

GalICS 2.0 exists in both a FORTRAN 2003 and a C++ version. Their quantitative agreement to several significant digits on a galaxy-by-galaxy basis is one of the reasons why we are confident of the code's quality and reliability. Finally, GalICS 2.0 should not be confused with EGALICS (Cousin et al. 2015). EGALICS started from an early version of GalICS 2.0 but the two codes have been developed independently in the last few years and should be considered as different SAMs.

### 2.1 TREE

#### 2.1.1 Cosmology and analysis of the N-body simulation

GalICS 2.0 uses DM merger trees from cosmological N-body simulations. In this paper, we use a simulation with $\Omega_M = 0.308$, $\Omega_\Lambda = 0.692$, $\Omega_b = 0.0481$, $\sigma_8 = 0.807$ and $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Planck Collaboration XX 2014, Planck + WP + BAO). The simulation has a volume of $(100 \text{ Mpc})^3$ and contains $512^3$ particles (implying a particle mass of $2.9 \times 10^8 M_{\odot}$). As the Poisson equation is solved on a non-uniform mesh using a multi-grid method in RAMSES (Teyssier 2002 for details), the force resolution is not spatially uniform. The simulation reaches a maximum force resolution of 1.5 kpc (physical units) in the densest regions (the centres of DM haloes). Outputs have been saved at 265 time-steps equally spaced in the logarithm of the expansion factor between $z = 13.2$ and 0.

Haloes and subhaloes are identified with the halo finder HaloMaker, which is based on AdaptaHOP (Tweed et al. 2009). For each halo containing at least a hundred particles (corresponding to a minimum halo mass of $2.9 \times 10^{10} M_{\odot}$), we determine the inertia ellipsoid, centred on its centre of mass, after iteratively removing gravitationally unbound particles, and we keep rescaling it until the critical overdensity contrast inside the inertia ellipsoid equals the one that we compute with the fitting formulae of Bryan & Norman (1998) for a Planck cosmology.<fnref n="1" /> The virial mass $M_{vir}$ is the mass of the gravitationally bound N-body particles contained within the virial ellipsoid (i.e. the rescaled inertia ellipsoid). The virial radius $r_{vir}$ is that of a sphere whose volume equals that of the virial ellipsoid. At $z = 0$, Bryan & Norman (1998)'s formulae give $\Delta_c = 102$. Hence, a halo mass of $M_{vir} = 10^{12} M_{\odot}$ corresponds to $v_{vir} = \sqrt{GM_{vir}/r_{vir}} = 128 \text{ km s}^{-1}$.

Fig. 1 shows the halo mass function that we measure in our N-body simulation at $z = 0$. The resolution mass is clearly visible as the $M_{vir}$ below which the mass function drops (the peak is at $M_{vir} = 10^{10.625} M_{\odot}$ rather than at $M_{vir} \simeq 10^{10.5} M_{\odot}$ simply because this is the mid-point of the logarithmic mass bin $10^{10.5} M_{\odot} < M_{vir} < 10^{10.75} M_{\odot}$). The effects of mass variance/low-number statistics at high masses can be quantified by comparing our mass function (the grey shaded area) with the analytic fit of Sheth, Mo & Tormen (2001). This fit contains two free parameters, which the original paper adjusted on N-body simulations by Kauffmann et al. (1999) and which we have recalibrated on ours. The fit is consistent with the mass function of galaxy clusters from Wen et al. (2010, black squares in Fig. 1). The comparison of our halo mass

<footnote n="1">¹ $3M_{vir}/(4\pi\rho_c r_{vir}^3)$, where $\rho_c$ is the critical density of the Universe, cannot be exactly equal to $\Delta_c$ because $M_{vir}$ can only take values that are multiples of the N-body particle's mass. However, the only haloes for which the difference is significant are low-mass systems just above the detection threshold.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The halo mass function measured from our N-body simulation (shaded area) compared to a Sheth et al. (2001) fit for the same cosmology [the solid curve corresponds to $a = 0.6$ and $q = 0.354$, where the parameters $a$ and $q$ are defined as in Sheth et al. (2001)]. Here $\phi_h$ is the number density of haloes per dex of virial mass. If one defines the halo mass function $n$ so that $n(M_{vir}) dM_{vir}$ is the number density of haloes with mass between $M_{vir}$ and $M_{vir} + dM_{vir}$, then the relation between $\phi_h$ and $n$ is $\phi_h = n(M_{vir})M_{vir} \log(10)$. The thickness of the shaded curve is determined by the Poissonian error on the number of haloes in each mass bin. The black squares are data points from Wen, Han & Liu (2010) for the mass function of galaxy clusters. The vertical dashed line at $M_{vir} = 10^{10.625} M_{\odot}$ (~150 particles) provides a measure of the halo mass resolution.</figcaption></figure>

function (grey shaded area) with this fit (black solid curve) shows good agreement up to $M_{vir} \sim 10^{13.3} M_{\odot}$. However, at $M_{vir} \sim 10^{13.8} - 10^{14} M_{\odot}$, the centre of the grey shaded area is above the black solid curve by about ~40 per cent. Number densities for the central galaxies of haloes in this mass range may be overestimated by a similar factor.

For each halo, we measure the virial angular momentum $J_{vir}$, which we use to compute the spin parameter $\lambda$ (Section 2.3.2). Bett et al. (2007) argued that they needed at least 300 particles to measure halo spins accurately but, going from 300 to 100 particles, as we do here, the median spin parameter changes from $\lambda = 0.0425$ to 0.044 (fig. 7 of Bett et al. 2007). The difference (<4 per cent) is well within the uncertainties of the results presented in this work and is comparable to the uncertainties that derive from different unbinding procedures (e.g. Onions et al. 2013).

We also fit the mass distribution of each halo with an NFW profile (Navarro, Frenk & White 1997) to measure its concentration $c$. Neto et al. (2007) found that they needed at least $10^4$ particles to measure concentrations accurately. Indeed, our concentrations drop below the fitting formulae by Muñoz-Cuartas et al. (2011) and Dutton & Macciò (2014) for $M_{vir} < 10^{12} M_{\odot}$ (corresponding to ~3000 particles), suggesting that our measurements are affected by N-body resolution in a manner that may be significant for haloes with $M_{vir} < 10^{11.5} M_{\odot}$. Most of the figures shown in this paper are not sensitive to the value of $c$ but the TFR is. We could use the concentrations that we measure in our N-body simulation for haloes with more than 3000 particles and the fitting formulae of Dutton & Macciò (2014) for haloes with less than 3000 particles but this approach neglects the large scatter in measured concentration values and poses the problem that our $c-M_{vir}$ relation exhibits systematic differences with respect to that of Dutton & Macciò (2014) even for haloes with more than $10^5$ particles. There is also the problem of applying fitting formulae for haloes to subhaloes. In this paper, we do all calculations with concentrations from our N-body simulation for self-consistency. However, in Section 3.7,

we explore how the TFR varies when our concentration measurements are replaced with values from the fitting formulae of Dutton & Macciò (2014).

In the case of a subhalo, we use the same procedure applied to haloes to obtain a first estimate of $r_{vir}$. Then, we shrink the subhalo by peeling off its outer layers until the density at the recomputed virial radius is at least as large as the host density at the position where the subhalo is located. The concentration parameter $c$ is recomputed accordingly. The particles removed from the outer layers of a subhalo are reassigned to the host halo if they are gravitationally bound to it. The host halo masses used in GalICS 2.0 are exclusive, i.e. they do not include those of subhaloes. By construction, a host halo is always more massive than its most massive subhalo.

The TreeMaker algorithm (Tweed et al. 2009) is used to link haloes/subhaloes identified at different redshifts and to generate merger trees. A halo is identified as the descendant of another when it inherits more than half of its progenitor's particles. Because of this definition, a halo can have many progenitors but at most one descendant. The main progenitor is always the one with the largest virial mass. A halo/subhalo is found to have no descendant if it loses more than half its mass but no single halo accretes enough mass from it to qualify as its descendant. Haloes that disappear in this manner exist in our merger trees but they are so rare that they are not statistically significant (tidally stripped particles are normally accreted by the halo that causes the tide).

### 2.1.2 Halo representation

The TREE module reads the halo catalogues and the tree structure. The properties that are read and used for each halo are (1) virial mass, radius and angular momentum, $M_{vir}$, $r_{vir}$ and $J_{vir}$, (2) concentration $c$, (3) position and velocity of the centre of mass, (4) position in the hierarchy of substructures (host halo/subhaloes) and (5) position in the tree (progenitors/descendant).

The virial mass $M_{vir}$ measured in the N-body simulation is a total mass of DM and baryons, treated as if they were both collisionless. Assuming that the DM distribution is described by the NFW profile, the mass of DM enclosed within a sphere of radius $r$ is

$$M_{\text{dm}}(x) = \frac{\log(1 + cx) - cx/(1 + cx)}{\log(1 + c) - c/(1 + c)} \cdot \frac{\Omega_M - \Omega_b}{\Omega_M} M_{\text{vir}}, \quad (1)$$

where $x \equiv r/r_{vir}$. Here and throughout this paper, log is the natural logarithm. The decimal logarithm is Log.

### 2.1.3 Scheme to evolve baryons along merger trees

The code loops over all time-steps. At each time-step, it loops over all haloes and calls the routines that compute the evolution of the baryons within them. The transition from one time-step to the next is handled as follows. Let $t_s$ be the age of the Universe at time-step $s$. If a halo detected at time-step $s$ has one or more progenitors at time-step $s-1$, we compute a random merging time

$$t_{s-1} < t_m < t_s. \quad (2)$$

If the halo has no progenitors, we assume that $t_m = t_{s-1}$ is its formation time. The baryons are evolved using the DM properties at $t_{s-1}$ between $t_{s-1}$ and $t_m$. If there are more than one progenitor, they are merged at $t_m$. The merger remnant, or the halo if there is only one progenitor, is then evolved between $t_m$ and $t_s$ using the DM properties measured at $t_s$.

In the code's current version, there is a one-to-one correspondence between satellite galaxies and subhaloes. When the latter
merge, so do the former. However, we have developed a beta version that includes the possibility of delayed mergers. The beta version computes the delay using Jiang et al. (2008)'s formula for the dynamical friction time-scale (see Cattaneo et al. 2011 for a description of the method). Preliminary investigations with the beta version show little difference with respect to the conclusions of this paper.

## 2.2 HALO

The HALO module follows the accretion of gas on to a halo and the exchanges of matter between its baryonic components (the cold gas in the halo, the hot-gas halo and the central galaxy).

### 2.2.1 Accretion

The mass $M_b$ of the baryons within a halo is updated from its value at $t_m$ (with $t_s - 1 < t_m < t_s$) to its values at $t_s$ with the equation:

$$ \Delta M_b = \max[f_b M_{\text{vir}}(t_s) - M_b(t_m), 0], \quad (3) $$

where $M_b(t_m)$ is the sum of the baryon masses of all progenitors at the time of merging.

There is no accretion on to subhaloes or their descendants. This prescription models the physical phenomenon of strangulation through ram-pressure and tidal stripping (Gunn & Gott 1972; Abadi, Moore & Bower 1999; Balogh & Morris 2000; Balogh, Navarro & Morris 2000; Peng, Maiolino & Cochrane 2015). It is also imposed because a halo may become a subhalo, evade detection by the halo finder as it passes through the centre of its host, reappear as a subhalo one or two time-steps later, and finally be identified as a halo again when it comes out on the other side ('backsplash haloes'). When the subhalo disappears, the associated satellite galaxy merges with the central one unless the subhalo is not bound to the host, in which case the satellite galaxy and its baryons disappear from the model universe. However, this case is rare and has no impact on our statistical predictions. Backsplash haloes are much more frequent. Allowing them to accrete gas would duplicate baryons by reinstating galaxies that have just merged. Hence, in GalICS 2.0, backsplash haloes are haloes with no galaxies (we are working on a new version that improves the description of these systems).

In equation (3), $f_b = f_b(M_{\text{vir}}, z)$ is a function that models reionization feedback (gas will not accrete on to haloes with virial temperature lower than the temperature of the intergalactic medium). We assume that the intergalactic medium has a Maxwellian velocity distribution and that the escape speed is $v_{\text{esc}} = v_{\text{vir}}\sqrt{2}$, where $v_{\text{vir}}(M_{\text{vir}}, z)$ is the halo virial velocity. This assumption gives

$$ f_b = \frac{\Omega_b}{\Omega_M} \int_0^{v_{\text{vir}}\sqrt{2}} 4\pi \left( \frac{\mu}{2\pi k T_{\text{reio}}} \right)^{\frac{3}{2}} v^2 e^{-\frac{\mu v^2}{2k T_{\text{reio}}}} dv, \quad (4) $$

where $T_{\text{reio}}$ is the temperature at which the intergalactic medium is reionized, $\mu$ is the mean particle mass and $k$ is the Boltzmann constant. The cosmic baryon fraction $\Omega_b/\Omega_M$ appears in front of the right-hand side of equation (4) because the integral gives the fraction of the baryonic mass with $v < v_{\text{esc}}$, while $f_b$ is expressed in terms of the total mass $M_{\text{vir}}$. Let $\sigma_{\text{reio}}$ be the one-dimensional thermal velocity dispersion of the intergalactic medium, so that $\mu\sigma_{\text{reio}}^2 = kT_{\text{reio}}$. Then, equation (4) becomes

$$ f_b = \frac{\Omega_b}{\Omega_M} \left[ \text{erf} \left( \frac{v_{\text{vir}}}{\sigma_{\text{reio}}} \right) - \frac{2}{\sqrt{\pi}} \left( \frac{v_{\text{vir}}}{\sigma_{\text{reio}}} \right) e^{-\left( \frac{v_{\text{vir}}}{\sigma_{\text{reio}}} \right)^2} \right]. \quad (5) $$

Equation (5) implies that $f_b = 0.5\Omega_b/\Omega_M$ for $v_{\text{vir}}/\sigma_{\text{reio}} = 1.088$.

Our assumption for $v_{\text{esc}}$ is accurate only for a singular isothermal sphere. In the NFW model, $v_{\text{esc}} \sim 2.5v_{\text{vir}}$ (the exact value depends on concentration). However, the difference in $v_{\text{esc}}$ will simply be reabsorbed by $\sigma_{\text{reio}}$, which is a free parameter of the model. For this reason, $\sigma_{\text{reio}}$ may differ from the physical thermal velocity dispersion of the intergalactic medium by a factor of order unity. Finally, reionization feedback has been added with a view to running GalICS 2.0 on higher resolution N-body simulations because we expect that $\sigma_{\text{reio}} \sim 20-30 \text{ km s}^{-1}$ and the current simulation cannot resolve haloes with $v_{\text{vir}} < 40 \text{ km s}^{-1}$.

A fraction $f_{\text{hot}}$ of the accreted baryon mass $\Delta M_b$ is shock heated to the virial temperature while falling in and is added to the hot halo. The rest accretes on to the central galaxy through cold filamentary flows. The fraction $f_{\text{hot}}(M_{\text{vir}})$ is assumed to be $f_{\text{hot}} = 0$ for $M_{\text{vir}} \le M_{\text{shock}}$,

$$ f_{\text{hot}} = \frac{\text{Log } M_{\text{vir}} - \text{Log } M_{\text{shock}}}{\text{Log } M_{\text{shutdown}} - \text{Log } M_{\text{shock}}} \quad (6) $$

for $M_{\text{shock}} < M_{\text{vir}} < M_{\text{shutdown}}$, and $f_{\text{hot}} = 1$ for $M_{\text{vir}} \ge M_{\text{shutdown}}$. Here $M_{\text{shock}}$ is the halo mass at which the accreted gas begins to be shock heated, while $M_{\text{shutdown}}$ is the halo mass above which shock heating is complete and any residual cold gas in the filaments is evaporated by the hot phase. While $M_{\text{shock}}$ and $M_{\text{shutdown}}$ are in principle free parameters of the model, to be determined by fitting the galaxy SMF, the functional form in equation (6) is based on work by Ocvirk, Pichon & Teyssier (2008), who measured the flow rates of cold ($T < 10^5$ K) and hot ($T > 10^5$ K) gas through a spherical surface of radius $0.2r_{\text{vir}}$ in the Horizon–Mare Nostrum cosmological hydrodynamical simulation. The results of Ocvirk et al. (2008) for

$$ f_{\text{hot}} = \frac{\dot{M}_{\text{hot}}}{\dot{M}_{\text{cold}} + \dot{M}_{\text{hot}}} \quad (7) $$

at $0.2r_{\text{vir}}$ are very similar to a ramp between $\text{Log}(M_{\text{shock}}/M_{\odot}) = 10.7$ and $\text{Log}(M_{\text{shutdown}}/M_{\odot}) = 12.7$, at least in the redshift range $2 < z < 4$ (the Horizon–Mare Nostrum simulation stops at $z = 2$). In GalICS, Cattaneo et al. (2006) found a good fit to Sloan Digital Sky Survey (SDSS) data for a sharper cut-off of the form $M_{\text{shock}} = M_{\text{shutdown}} = 10^{12.3} M_{\odot}$.

The equations for the variations of the masses of the hot gas and the cold filaments are

$$ \Delta M_{\text{hot}} = f_{\text{hot}}\Delta M_{\text{b}}, \quad (8) $$

$$ \Delta M_{\text{fil}} = (1 - f_{\text{hot}})\Delta M_{\text{b}} - \frac{M_{\text{fil}}}{t_{\text{ff}}} \Delta t, \quad (9) $$

where $\Delta M_b$ is the accreted mass calculated in equation (3). The second term on the right-hand side of equation (9) is the accretion rate from the filaments on to the galaxy, which we assume to take place on a free-fall time-scale $t_{\text{ff}} = r_{\text{vir}}/v_{\text{vir}}$.

Following Dekel & Birnboim (2006) and Dekel et al. (2009), we assume that the accretion of cold gas is the main mode of galaxy formation and that hot gas never cools. Hence, there is no cooling term in equation (8). The assumption of a total shutdown of gas accretion in massive haloes is extreme (see e.g. Bildfell et al. 2008), but we know from previous work (Cattaneo et al. 2006) that its predictions are in good agreement with the galaxy colour–magnitude distribution, while letting the hot gas cool leads to results in clear disagreement with the observations. Introducing cooling makes sense only if one has a physical model of how AGN feedback mitigates it (see Cattaneo et al. 2009 for a review). Attempts in this direction have been made, starting with Croton et al. (2006), Bower et al. (2006) and Somerville et al. (2008). Following Benson &
Babul (2009), Fanidakis et al. (2011) have gone as far as to compute the mechanical luminosity of the jets from the accretion rate and the spin of the black holes that power them (an approach pioneered in SAMs by Cattaneo 2002). However, the physics of these models are uncertain. Hence, we have considered that it would be premature to include AGN feedback in our SAM, especially in a paper focused on spiral galaxies, which tend to live in haloes with $M_{vir} \lesssim 10^{12} M_{\odot}$.

In these lower mass systems, however, the assumption that hot gas never cools may be even more extreme because it prevents the re-accretion of ejected gas (see the discussion in Section 2.2.2) and thus the possibility of substantially delaying star formation with respect to gas accretion on to the halo. It is important to realize that, in this paper, we are not arguing for the absence of cooling on physical ground. This is the simplest possible assumption within the cold-flow paradigm, and we want to explore how far it can take us. In Section 3, we shall show that the results are reasonably good, although we have not compared them to all possible observations. [Bower et al. (2006) have suggested that the gradual re-accretion of ejected gas is necessary in order to predict a large enough fraction of galaxies on the blue sequence at low stellar masses].

The effects of introducing cooling and therefore re-accretion will be explored in a future publication. However, we note that, to an extent, equation (6) already includes some of the effects of cooling implicitly because the cold gas fraction $1 - f_{hot}$ that Ocvirk et al. (2008) measure in their simulation at $0.2r_{vir}$ includes any gas that may have been shock heated and cooled before reaching $0.2r_{vir}$ (the outer boundary of the H1 disc; Pichon, private communication), although it does not account for the possibility that a cooling flow may develop inside the galaxy itself and for the re-accretion of ejected gas.

### 2.2.2 Outflows

The gas that is blown out of the galaxy is either mixed to the hot gas in the halo or expelled from the halo altogether. In the second case, we store it in an outflow component, so that we always know how much gas and how many metals have been expelled from each halo.

The baryon mass $M_b$ that enters equation (3) is

$$M_b = M_{\text{hot}} + M_{\text{fil}} + M_{\text{gal}} + M_{\text{out}}, \quad (10)$$

where $M_{gal}$ is the total baryonic mass of the galaxy. This definition of $M_b$ includes the ejected mass $M_{out}$ to prevent its re-accretion. Hence, the actual halo baryon fraction $(M_{hot} + M_{fil} + M_{gal})/M_{vir}$ may be lower than $M_b/M_{vir}$. Had we defined $M_b$ as $M_b = M_{hot} + M_{fil} + M_{gal}$, $M_{out}$ would have been available for immediate re-accretion on to the halo.

In this paper, it makes no difference whether galactic winds mix with the hot gas or escape from the halo because hot gas is not allowed to cool. Physically, however, cooling can only be important if the mass of hot gas is at least comparable to the mass of cold gas in the filaments, in which case the material blown out of the galaxy will almost certainly mix with it.

## 2.3 GALAXY

The GALAXY module deals with the structural properties of galaxies (morphologies, scalelengths, kinematics) and is organized around two main routines. GALAXY EVOLVE follows the evolution of an individual galaxy over an interval of time. It computes disc radii, speeds and the formation of pseudo-bulges through disc instabilities. GALAXY MERGE models morphological transformations induced by mergers.

### 2.3.1 Galaxy structure

A galaxy is modelled as the sum of three components: a disc, built through gas accretion and minor mergers, a pseudo-bulge, built through disc instabilities, and a classical bulge, built through major mergers.

The disc has an exponential surface-density profile:

$$\Sigma(r) = \Sigma_0 \exp\left(-\frac{r}{r_d}\right). \quad (11)$$

The pseudo-bulge originates from the buckling and bar instability of the disc within radius $r_{pseudo}$. We therefore assume that, while instabilities transfer matter from the disc to the pseudo-bulge, the mass redistribution is mainly vertical and azimuthal, and that the radial exponential profile of the disc plus pseudo-bulge system is largely unaffected by this process.

We assume that the bulge is spherical and that its density distribution is described by a Hernquist (1990) profile. For this model, the bulge mass within radius $r$ is

$$M_{\text{bulge}}(r) = \left( \frac{r/r_{\text{bulge}}}{1 + r/r_{\text{bulge}}} \right)^2 M_{\text{bulge}}, \quad (12)$$

where $r_{bulge}$ is the scale radius of the bulge.

Having described the components with which we model a galaxy, we are ready to enter the details of the physical processes that determine their formation and characteristic quantities.

### 2.3.2 Disc radii and rotation speeds

Any gas that accretes on to a galaxy is automatically added to the disc component. Its scale radius $r_d$ is determined by solving the disc angular momentum equation

$$J_d = 2\pi\Sigma_0 \int_0^\infty e^{-r/r_d} r^2 v_c dr = M_d r_d \int_0^\infty e^{-x} x^2 v_c dx, \quad (13)$$

where $v_c(r)$ is the disc rotation curve (see equation 17 below) and $M_d$ is the total baryonic mass (stars and gas) of the disc and the pseudo-bulge combined (this is what $M_d$ stands for throughout this paper).

Following Fall & Efstathiou (1980) and Mo, Mao & White (1998), SAMs compute disc sizes by assuming that baryons and DM have the same initial angular momentum distribution and that specific angular momentum is conserved, i.e. that

$$\frac{J_d}{M_d} = \frac{J_{vir}}{M_{vir}}. \quad (14)$$

Cosmological hydrodynamic simulations have shown that the angular momentum of the gas is not conserved during infall (Kimm et al. 2011; Danovich et al. 2015). However, even though equation (14) can give incorrect results when applied to individual galaxies, it retains a statistical validity because the distribution of specific angular momentum of discs is similar to that of DM haloes (Danovich et al. 2015). This statistical validity is backed by observations both in the local Universe (Tonini et al. 2006) and at high redshift (Burkert et al. 2016).

Equation (13) can be rewritten as

$$r_d = \frac{\lambda r_{vir}}{\int_0^\infty e^{-x} x^2 \frac{v_c}{v_{vir}} dx}, \quad (15)$$

where

$$\lambda = \frac{J_h}{M_h r_{vir} v_{vir}} \quad (16)$$
is the halo spin parameter as defined by Bullock et al. (2001). This definition differs from the usual one by Fall & Efstathiou (1980) by a factor equal to $\sqrt{2}$ for a truncated singular isothermal sphere. If $v_c$ were equal to $v_{vir}$ at all radii, then the integral in equation (15) would make 2 and equation (15) would reduce to $r_d = \lambda r_{vir}/2$. However, $v_c(r)$ is not flat.

The disc rotation curve $v_c(r)$ is determined by the sum in quadrature of three terms:

$$v_c^2(r) = \frac{GM_{dm}(r)}{r} + \frac{GM_{bulge}(r)}{r} + v_d^2(r), \quad (17)$$

where $M_{dm}(r)$ is the DM mass with radius $r$ (equation 1), $M_{bulge}(r)$ is the bulge mass within radius $r$ (equation 12) and $v_d(r)$ is the contribution from the disc, which has a more complicated form because the disc has a cylindrical rather than spherical symmetry. We compute $v_d(r)$ exactly by using Bessel functions as in Freeman (1970). The dependence of $v_d(r)$ on $r_d$ is the reason why solving equation (13) is not trivial. The pseudo-bulge does not appear as a fourth contribution to $v_c(r)$ in equation (17) because it is simply the inner disc that buckles up (Section 2.3.3). There is no radial migration in our model.

### 2.3.3 Adiabatic contraction

Adiabatic contraction is the contraction of the DM halo in response to the infall of the baryons. Blumenthal et al. (1986) estimated it from the adiabatic invariance of the specific angular momentum $rv_c(r) \propto \sqrt{rM(r)}$, where $M(r)$ is the total mass enclosed in a sphere of radius $r$. If $r_i$ is the initial radius of a DM shell that contracts to radius $r$ and $M_{NFW}(r_i)$ is the initial DM profile (assumed to be described by the NFW model), then this assumption gives the equation

$$r[M_{NFW}(r_i) + M_d(r) + M_{bulge}(r)] = r_i \frac{M_{NFW}(r_i)}{1 - \frac{\Omega_b}{\Omega_M}}, \quad (18)$$

from which $r_i(r)$ and thus $M_{dm}(r) = M_{NFW}(r_i)$ can be computed.

Equation (18) is the standard description of adiabatic contraction in models of galaxy formation (e.g. Mo et al. 1998; Cole et al. 2000; Somerville et al. 2008). However, cosmological hydrodynamic simulations have shown that it overestimates its importance (Gnedin et al. 2004; Abadi et al. 2010). It is also likely that adiabatic contraction may be compensated by adiabatic expansion during massive outflows (Pontzen & Governato 2012; Teyssier et al. 2013; Tollet et al. 2016) because the haloes of dwarf galaxies have shallow cores rather than the central cusps predicted by cosmological simulations of dissipationless hierarchical clustering (Moore 1994 and Flores & Primack 1996, but also see Swaters et al. 2003). Hence, in the standard version of GallICS 2.0, there is no adiabatic contraction.

A version with adiabatic contraction (computed with equation 18) has however been explored. Its results will be briefly discussed in Section 3.7, when we talk about the possible effects of adiabatic contraction on the TFR.

### 2.3.4 Disc instabilities

Following Efstathiou, Lake & Negroponte (1982), Christodoulou, Shlosman & Tohline (1995) and van den Bosch (1998), we assume that discs are unstable when their self-gravity contributes more than a critical fraction of the circular velocity (see however, Athanassoula 2008 for a criticism of this model). In formulae, our instability condition is

$$v_d > \epsilon_{\text{inst}} v_c, \quad (19)$$

where $\epsilon_{inst} \le 1$ is a free parameter of the model that sets the instability threshold ($\epsilon_{inst} = 1$ corresponds to the unphysical assumption that all discs are always stable, that is, to turning off disc instabilities).

The pseudo-bulge radius $r_{pseudo}$ is the largest radius at which the instability condition (equation 19) is satisfied. Its value is used to increment the pseudo-bulge mass with the algorithm

$$\Delta M_{\text{pseudo}} = \max[M_d(r_{\text{pseudo}}) - M_{\text{pseudo}}, 0], \quad (20)$$

where $M_{pseudo}$ is the pseudo-bulge mass before incrementation. Equation (20) guarantees that the pseudo-bulge mass never decreases (except in major mergers, where all components form one large classical bulge). Gas and stars transferred from the disc to the pseudo-bulge contribute to $\Delta M_{pseudo}$ in a ratio equal to the disc gas-to-stellar mass ratio. The pseudo-bulge characteristic speed is $v_{pseudo} = v_c(r_{pseudo})$.

We note that, in our model, the pseudo-bulge is any structure formed by disc instabilities, be it a peanut-shaped pseudo-bulge, a bar or an oval.

### 2.3.5 Mergers

A merger of two galaxies is major if $M_1 < \epsilon_m M_2$ with $M_1 > M_2$, where $\epsilon_m$ is a parameter of the model. Here $M_1$ and $M_2$ are total masses of the two galaxies (baryons and DM) within their baryonic half-mass radii $r_1$ and $r_2$, which we compute numerically from the profiles of their discs and bulges.

In a major merger, the two galaxies are destroyed and all their baryons are put into one large bulge, the size of which is determined by energy conservation (Cole et al. 2000; Hatton et al. 2003; Shen et al. 2003):

$$\frac{1}{2} U_{\text{bulge}} = \frac{1}{2} U_1 + \frac{1}{2} U_2 + E_{12}, \quad (21)$$

where $U_{bulge}$, $U_1$ and $U_2$ are the gravitational potential energies of the resulting bulge, galaxy 1 and galaxy 2, respectively, while $E_{12}$ is the interaction energy of the two-body system. The $1/2$ coefficient in front of the potential-energy terms comes from the virial theorem, since, for each galaxy, the total energy equals half the gravitational potential energy. Equation (21) can be rewritten as

$$\frac{1}{2} C_{\text{remn}} \frac{G(M_1 + M_2)^2}{r_{1/2}} = \frac{1}{2} C_1 \frac{G M_1^2}{r_1} + \frac{1}{2} C_2 \frac{G M_2^2}{r_2} + \frac{1}{2} C_{12} \frac{G M_1 M_2}{r_1 + r_2}. \quad (22)$$

Here $r_{1/2}$ is the merger remnant's half-mass radius. $C_{remn}$, $C_1$ and $C_2$ are form factors that relate the gravitational potential energy of each system to its mass and half-mass radius. $C_{12} = 2E_{12}/U_{12}$, where $U_{12} = -G M_1 M_2 / (r_1 + r_2)$ is the two-body system gravitational binding energy. The factor of 2 in the definition of $C_{12}$ comes, once again, from the virial theorem. $C_{12} = 1$ corresponds to merging from circular orbits. $C_{12} = 0$ corresponds to a motion that starts with zero speed at infinity.

The form factor that relates the gravitational potential energy of a system to its mass and half-mass radius is about 0.4 for a Hernquist (1990) profile and 0.49 for the purely academic case of a self-gravitating exponential disc. Real form factors are complicated by the presence of several components but can be simplified by assuming that they are all equal. Covington et al. (2008) have shown that equation (22) with $C_{remn} = C_1 = C_2 = 0.5 C_{12}$ is in good agreement with the sizes of the remnants of dissipationless mergers
in hydrodynamic simulations, though the best fit for $C_{\text{remn}}/C_{12}$ depends on the simulations' assumptions for the initial orbits. Here, we use $C_{12} = 0$ because Shankar et al. (2013, 2014) showed that this assumption is in better agreement with the size evolution of spheroids. Covington et al. (2008) have also shown that energy dissipation through radiation in gas-rich mergers can result in final values of $r_{1/2}$ smaller than those obtained from energy conservation (equation 22) by up to a factor of 2.

The scale radius of the bulge $r_{\text{bulge}}$ is computed from the relation (Hernquist 1990)

$$r_{1/2} = (1 + \sqrt{2})r_{\text{bulge}}. \qquad (23)$$

This equation assumes that the half-mass radius for the stars is equal to the half-mass radius for the stars and the DM within the galaxy (see the definition of $M_1$ and $M_2$ at the beginning of this subsection). This assumption is reasonable, inasmuch as we know that DM makes a negligible contribution to the stellar dynamics of elliptical galaxies.

The radius $r_{1/2}$ computed with equation (23) is the half-mass radius in three dimensions. The radius that contains half of the mass in a two-dimensional projection is $R_e = r_{1/2}/1.33 = 1.8153 r_{\text{bulge}}$ assuming a Hernquist profile (Hernquist 1990).

The one-dimensional stellar velocity dispersion averaged over a cylindrical aperture on the sky of projected radius $R_e$ is given by

$$\sigma_e^2 = G \frac{M_1 + M_2}{C R_e}, \qquad (24)$$

where C is a structure coefficient. For a Hernquist profile with isotropic velocity dispersion, $C = 3.31$ (Courteau et al. 2014, chapter 5, appendix B). However, the Hernquist profile is just a phenomenological model. It has the advantage that its total mass, gravitational potential and velocity dispersion can be computed analytically, but there is no physical reason why galaxies should follow it. In fact, systematic departures from C = 3.31 are observed (Cappellari et al. 2013; Courteau et al. 2014). We therefore treat C as a parameter of the model to be constrained by observations.

A fraction $\epsilon_*$ of the total gas mass of the merging galaxies falls to the centre and feeds the growth of a supermassive black hole. The post-merger black hole mass is computed with the formula

$$M_* = M_{*1} + M_{*2} + \epsilon_* (M_{\text{gas1}} + M_{\text{gas2}}), \qquad (25)$$

where $M_{*1}$, $M_{*2}$, $M_{\text{gas1}}$, $M_{\text{gas2}}$ are the black hole and gas masses of the merging galaxies prior to the merging event.

In a minor merger, the gas and the stars in the disc of the smaller galaxy are added to the disc of the larger galaxy, while the gas and the stars in the bulge of the smaller galaxy are added to the bulge of the larger galaxy.

## 2.4 COMPONENT

A galactic component (a disc, a pseudo-bulge or bulge) is composed of a stellar population and its interstellar medium. The COMPONENT module follows the exchanges of matter between gas and stars, as well as the ejection of gas from a component, i.e. the processes of star formation and feedback.

The equation that governs the evolution of the gas mass within a component is

$$\dot{M}_{\text{gas}} = \dot{M}_{\text{accr}} + \dot{M}_{\text{sml}} - \text{SFR} - \dot{M}_{\text{out}}. \qquad (26)$$

Here $\dot{M}_{\text{accr}}$ is the accretion rate on to the galaxy, which is entirely due to cold flows since there is no cooling in our model (Section 2.2.1). $\dot{M}_{\text{sml}}$ is the gas deposited into interstellar medium by the later stages

of stellar evolution (stellar mass-loss), SFR is the star formation rate and $\dot{M}_{\text{out}}$ is the rate at which gas is blown out by stellar feedback (outflow rate).

The time separation between two output time-steps $t_s - 1$ and $t_s$ in the N-body simulation used to construct the merger trees ranges from 6 Myr at $z \simeq 13$ to 144 Myr at $z \simeq 0$, but we follow star formation and feedback on smaller substeps of $\Delta t = 1$ Myr. We have chosen this value because it is short compared to the time-scale on which a stellar population evolves, i.e. the time-scale on which $\dot{M}_{\text{sml}}$ varies (even massive OB stars spend >6 Myr on the main sequence). $\dot{M}_{\text{accr}}$ varies on a time-scale (set by $t_{\text{dyn}}$) that is even longer. Hence, without loss of generality, we can write SFR = $M_{\text{gas}}/t_{\text{sf}}$ and $\dot{M}_{\text{out}} = \eta$SFR, where our only assumption about the star formation time-scale $t_{\text{sf}}$ and the mass-loading factor $\eta$ is that they are constant during a $\Delta t = 1$ Myr substep, even though they vary on longer time-scales in response to stellar evolution or changes in the structural properties of the component (computed in GALAXY).

In conclusion, the equations for $M_{\text{gas}}$ and $M_{\text{stars}}$ take the form

$$\dot{M}_{\text{gas}} = \dot{M}_{\text{accr}} + \dot{M}_{\text{sml}} - (1 + \eta) \frac{M_{\text{gas}}}{t_{\text{sf}}}, \qquad (27)$$

$$\dot{M}_{\text{stars}} = \frac{\dot{M}_{\text{gas}}}{t_{\text{sf}}} - \dot{M}_{\text{sml}}. \qquad (28)$$

These equations have an analytic solution (Cole et al. 2000; Lilly et al. 2013; Dekel & Mandelker 2014; Peng & Maiolino 2014),<fnref n="2" /> which we use to evolve $M_{\text{gas}}$ and $M_{\text{stars}}$ on substeps of $\Delta t = 1$ Myr (a direct numerical integration would require $\Delta t \ll 1$ Myr to provide the required accuracy).

In this paper, we assume instantaneous recycling (a fraction R of the gas that forms stars is immediately returned to the interstellar medium). This assumption is coded in the STAR module (Section 2.5). Since $\dot{M}_{\text{sml}} = R M_{\text{gas}}/t_{\text{sf}}$, the presence of an analytic solution means that there is no need for any substepping whatsoever. However, we retain the substepping with $\Delta t = 1$ Myr because, in the future, we may want to replace the instantaneous recycling approximation with a stellar evolution model in which $\dot{M}_{\text{sml}}$ does depend on the age of the stellar population (COMPONENT assumes no knowledge of what is assumed in STARS).

This description hides the complexity of star formation and feedback in the values of $t_{\text{sf}}$ and $\eta$, which may depend on several galaxy properties. Their calculation is done by the star formation model (Section 2.4.1) and the feedback model (Section 2.4.2), respectively.

### 2.4.1 The star formation model

The star formation model computes the star formation time-scale $t_{\text{sf}}$. For discs, we assume that $t_{\text{sf}}$ is proportional to $t_{\text{dyn}}$, where $t_{\text{dyn}}$ is the orbital time at the disc scale radius $r_a$ ($t_{\text{sf}} = t_{\text{dyn}}/\epsilon_{\text{sf}}$; see e.g. Kauffmann et al. 1993 and Kennicutt 1998). Hence

$$\text{SFR} = \epsilon_{\text{sf}} \frac{M_{\text{gas}}}{t_{\text{dyn}}}. \qquad (29)$$

The star formation efficiency $\epsilon_{\text{sf}} = 1/25$ adopted in this paper is calibrated on a local galaxy sample by Boselli et al. (2014). With this efficiency, our model reproduces the statistical relation between $M_{\text{gas}}$ and SFR in local galaxies (Fig. 2, left) and the SFRs of individual galaxies within a factor of 2 (Fig. 2, right). The gas

<footnote n="2">² The analytic solution has the interesting feature that $M_{\text{gas}} \rightarrow \dot{M}_{\text{accr}} t_{\text{sf}}/(1 + \eta)$ for $t \rightarrow \infty$ when $\dot{M}_{\text{sml}} = 0$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relation between total gas mass (H₁ plus H₂; left) and goodness of $M_{gas}/t_{dyn}$ as an SFR estimator (right) in a local sample from Boselli et al. (2014, black squares). Here, $t_{dyn} = 2\pi v_d/v_{rot}$, where $r_d$ and $v_{rot}$ are observational measurement from Boselli et al. (2014)'s sample, from which we have retained only galaxies classified as S or Irr for which there is a measurement of the exponential scalelength $r_d$. In the left-hand panel, the observed SFR-$M_{gas}$ relation is compared to what one would expect for SFR = $M_{gas}/25t_{dyn}$ (red circles). The red dashed line corresponds to a constant star formation time-scale of 3.5 Gyr. The right-hand panel shows that Boselli et al. (2014)'s data (black squares) follow the SFR = $M_{gas}/25t_{dyn}$ relation (red solid line) within a scatter of a factor of 2 (red dashed lines).</figcaption></figure>

masses plotted in Fig. 2 are total masses of cold neutral gas. Had we used H₂ masses instead, the observed $M_{gas}$–SFR relation would have had a higher normalization and would have been tighter because the molecular gas forms stars. However, deciding what fraction of the cold gas is molecular would add another layer of complexity and uncertainty to our model, possibly larger than the difference in scatter between the SFR-$M_{gas}$ relation and the SFR-$M_{H_2}$ relation.

Explaining why $t_{sf}$ is so long compared to the local free-fall time of the gas is a hot topic in star formation theory (e.g. Renaud, Kraljic & Bournaud 2012; Hopkins et al. 2014; Kraljic et al. 2014; Forbes et al. 2016; Gatto et al. 2017). By calibrating equation (29) on observational data, we are effectively short-circuiting these complex physics, whose outcome is summarized in the value of $\epsilon_{sf}$.

It is also important to remark that our star formation model is based on observations of disc galaxies. We assume that we can generalize it to pseudo-bulges and bulges by simply redefining the dynamical time. In GalICS 2.0, $t_{dyn}$ is the orbital time at $r_{pseudo}$ for pseudo-bulges and the half-crossing time of the starbursting region $t_{dyn} = r_{starburst}/\sigma$ for bulges.

Mergers are the only mechanism through which bulges can accrete gas. When they are gas-rich, they induce intense starbursts. The SFR has a first peak at the first pericentric passage and another later when the two galaxies coalesce (e.g. Di Matteo et al. 2008). At the first pericentric passage, the system's morphology is completely irregular and star formation is concentrated in a series of knots along the galaxies' spiral arms. By the time the galaxies coalesce, most of the gas has sunk to the centre of the merger remnant. The scalelength of its distribution is of the order of $r_{starburst} \sim 0.1r_{bulge}$ (Cattaneo et al. 2005a). In SAMs, mergers are instantaneous and galaxies jump from their initial morphologies directly to this final state. However crude, this assumption is in line with observational evidence that the star formation time-scale $t_{sf} = M_{gas}/\dot{M}_{sf}$ for star-burst galaxies is about 10 times shorter than it is for normal galaxies (Bigiel et al. 2008).

Observationally, galaxies begin to depart from the mean Schmidt–Kennicutt law (Kennicutt 1998) between SFR surface density $\Sigma_{SFR}$ and gas surface density $\Sigma_{gas}$ for $\Sigma_{gas} < \Sigma_{th} \sim 9 M_{\odot} \text{pc}^{-2}$ (Bigiel et al. 2008), where $\Sigma_{gas}$ is the mean gas surface density

(H₁ plus H₂) within the optical radius $r_{opt} = 3.2r_d$. However, there are galaxies on the relation (including some of the black squares in Fig. 2) with values of $\Sigma_{gas}$ as low as $\sim 2 M_{\odot} \text{pc}^{-2}$. Therefore, the threshold is not sharp. In GalICS 2.0, we set SFR = 0 for $\Sigma_{gas} < \Sigma_{th}$, where $\Sigma_{th}$ is a parameter of the model. We set it to the relatively low value $\Sigma_{th} = 2 M_{\odot} \text{pc}^{-2}$ because higher values suppress star formation too much in low-mass haloes, leading to galaxies that are all gas and no stars, though this may be a resolution artefact. The surface area on which we spread the gas to compute $\Sigma_{gas}$ is $2\pi(r_{opt}^2 - r_{pseudo}^2)$ for discs, $2\pi r_{pseudo}^2$ for pseudo-bulges and $2\pi r_{starburst}^2$ for bulges.

### 2.4.2 The feedback model

Feedback is a generic term for the effects that star formation and black hole accretion exert on the surrounding gas. These effects influence the processes that cause them and can regulate their rates. This section is on stellar feedback but even that is multifaceted because it results from the synergy of different processes (SNe, radiation pressure, photoionization and photoelectric heating) that act on different scales.

Mathews & Baker (1971) and Larson (1974) were the first to suggest that gas may be strongly heated by SN blastwaves and driven out of galaxies in hot winds. While SNe have certainly the energy to this, and have become for this reason a standard ingredient of galaxy formation theory, their efficiency and the mass scale at which they become important are affected by the fraction of SN energy that is radiated (Dekel & Silk 1986). If SNe explode inside dense molecular clouds, most of their energy will be quickly lost to X-rays. Radiation pressure and stellar winds from massive OB stars must disperse giant molecular clouds rapidly, after they have turned just a few per cent of their mass into stars, for this not to occur (Hopkins et al. 2013 and references therein). Photoelectrons extracted from dust grains by ultraviolet radiation are the primary source of heating for the neutral interstellar medium and suppress star formation by preventing its overheating and overcondensation into dense molecular clouds (Forbes et al. 2016).
These complex physics are beyond the scope of our feedback model, whose purpose is to computes the mass-loading factor $\eta = \dot{M}_{out}/SFR$, i.e. the rate at which cold gas is removed from galaxies. Any feedback mechanism that regulates star formation without removing gas from galaxies is already incorporated phenomenologically in our star formation efficiency $\epsilon_{sf}$ (Section 2.4.1). Similarly, the fraction $\epsilon_{SN}$ of the power output from SN explosions that is converted into wind kinetic energy and/or thermalized in the hot atmosphere is chosen to reproduce the observation and therefore includes the effects of all the other processes (e.g. radiation pressure, stellar winds, photoionization, photoelectric heating) that may affect the outflow rate.

If $\Psi_{SN} \simeq 1/(140 M_{\odot})$ is the number of SNe per unit stellar mass formed [assuming a Chabrier (2003) initial mass function in the stellar mass range 0.1–100 $M_{\odot}$ and a minimum mass for core-collapse SNe of $8 M_{\odot}$] and $E_{SN} \sim 10^{51}$ erg is the energy released by one SN, then the power output from SNe will be $E_{SN}\Psi_{SN}SFR$. If a fraction $\epsilon_{SN}$ of this power is used to drive a wind with speed $v_w$, then the outflow rate $\dot{M}_{out}$ from the component will satisfy

$$ \frac{1}{2} \dot{M}_{out} v_w^2 \sim \epsilon_{SN} E_{SN} \Psi_{SN} SFR \quad (30) $$

(Silk 2003). Expulsion from the gravitational potential well of the DM requires $v_w \ge v_{esc} \sim 2.5 v_{vir}$ (the numerical coefficient in front of $v_{vir}$ depends on halo concentration), but here we make no assumption as to whether the wind escapes from the halo or settles into a hot circumgalactic medium. We therefore reabsorb the uncertainty on $v_w/v_{vir}$ into the free parameter $\epsilon_{SN}$ and define mass-loading factor $\eta$ so that

$$ \eta \equiv \frac{\dot{M}_{out}}{SFR} = \frac{2\epsilon_{SN}E_{SN}\Psi_{SN}}{v_{vir}^2}. \quad (31) $$

The only inconvenient of this definition is that $\epsilon_{SN}$ underestimates the real SN efficiency required to produce the mass-loading factors assumed by our model. The difference is small (a factor of $\sim 1.5$) if the gas blown out of the galaxy is heated to the virial temperature and mixed with the hot atmosphere (for a singular isothermal sphere, $\frac{3}{2}kT_{vir} = \frac{3}{2} \cdot \frac{1}{2}\mu v_{vir}^2$, White & Frenk 1991). Much larger energetic efficiencies ($\gtrsim 2.5^2 \epsilon_{SN}$) are required if the gas expelled from galaxies is also blown out of the halo.

The problem of this simple scaling with $v_{vir}^{-2}$ is that it cannot reproduce the shallow slope of the low-mass end of the galaxy SMF (unless we use merger trees from a low-resolution N-body simulation that misses low-mass haloes, as in Cattaneo et al. 2006, but there we focused on massive galaxies). A phenomenological solution is to introduce an SN efficiency $\epsilon_{SN}$ that depends on both $v_{vir}$ and redshift $z$, and to impose a plausible maximum $\epsilon_{max}$ to the values that $\epsilon_{SN}$ can take

$$ \epsilon_{SN} = \min \left[ \left( \frac{v_{vir}}{v_{SN}} \right)^{\alpha_v} (1+z)^{\alpha_z}, \epsilon_{max} \right], \quad (32) $$

where $v_{SN}$, $\alpha_v$ and $\alpha_z$ are free parameters of the model to be determined by fitting the galaxy SMF. The speed $v_{SN}$ corresponds to the virial velocity for which $\epsilon_{SN} = 1$ at $z=0$ if no maximum efficiency is imposed.

As the laws of physics do not vary with time, one could argue that a physical model should not contain any explicit dependence on $z$. A simple answer is that this objection does not apply to a phenomenological model (see Peirani et al. 2012 for evidence from cosmological hydrodynamic simulations supporting more efficient feedback at high $z$). We also remark that the values with which we fit

the data ($\alpha_v = -4$, $\alpha_z = 3$, $v_{SN} = 24 \text{ km s}^{-1}$; Table 1 and Section 3) give a simple relation between mass-loading factor and halo mass:

$$ \eta \simeq 3.8 \left( \frac{M_{vir}}{10^{11} M_{\odot}} \right)^{-2}, \quad (33) $$

since $v_{vir} \propto M_{vir}^{1/3}(1+z)^{1/2}$. $M_{vir}$ is a physical quantity, though it is not clear why the outflow rate should scale with $M_{vir}$ rather than $v_{vir}$. The mass resolution of the N-body simulation used to construct the merger trees is $M_{vir} \sim 3 \times 10^{10} M_{\odot}$. In this paper, we formulate our model in terms of $v_{vir}$ and $z$ rather than $M_{vir}$ to ease comparison with previous work, for the sake of greater generality and because, with approach, it is easier to check that our feedback model is energetically plausible.

Equation (33) corresponds to a very strong dependence of the mass-loading factor on the virial velocity ($\eta \propto v_{vir}^{-6}$). For comparison, the exponents used by other SAMs are $-5.5$ (Cole et al. 1994), $-3.5$ (Guo et al. 2011), $-2.5$ (Somerville et al. 2012), $-0.92$ (Henriques et al. 2013) and $-3.2$ with an allowable range between 0 and $-5.5$ (Lacey et al. 2016), though the details of how stellar feedback is implemented vary from one model to another (see Hirschmann, De Lucia & Fontanot 2016 for a discussion of the mass loading in different SAMs and simulations).

Our normalization of $\eta$ at $M_{vir} = 10^{11} M_{\odot}$, $\eta = 3.8$, is comparable to those of Guo et al. ($\eta \sim 1$) and Henriques et al. ($\eta \sim 2.5$), but much lower than that of Lacey et al. (2016). As the mass $M_{vir} = 10^{11} M_{\odot}$ is only a factor of 3 larger than our resolution limit, our normalization combined to our much steeper dependence on $v_{vir}$ implies that our mass-loading factors are lower than those assumed by Guo et al. (2011), Henriques et al. (2013) and Lacey et al. (2016) at all but the smallest halo masses probed in this paper. It is possible that we fit the observations with lower mass-loading factors for a given halo mass because our current model neglects the re-accretion of ejected gas.

Physically, $\eta$ is limited by the maximum energetic efficiency of SNe $\epsilon_{max}$. Without such maximum, equation (32) implies $\epsilon_{SN} \to \infty$ for $v_{vir} \to 0$, which is absurd (the wind cannot contain more energy than it is available). In the most generous case, $\epsilon_{max} = 1$. The real efficiency will probably be much lower. In practice, $\epsilon_{SN}$ is limited by the mass resolution of the N-body simulation, $M_{vir} \sim 3 \times 10^{10} M_{\odot}$. Inserted into equation (33), this mass gives a maximum mass-loading factor of $\eta = 30-40$. As $M_{vir} \sim 3 \times 10^{10} M_{\odot}$ corresponds to $v_{vir} \sim 40 \text{ km s}^{-1}$ at $z=0$, our default parameter values (Table 1) imply $\epsilon_{SN} \lesssim 0.1$ at $z \sim 0$ for all haloes that we can resolve. At high $z$, however, $\epsilon_{SN}$ can take much larger values if no maximum efficiency is prescribed.

## 2.5 STAR

The STAR module follows the evolution of a component’s stellar population. In the code’s current version, this evolution is computed based on the instantaneous recycling approximation. Stellar evolution is, therefore, entirely described by two parameters: the returned fraction $R$ and the metal yield $y$. The explicit equations for the stellar mass-loss rate $\dot{M}_{sml}$ and the mass-loss rate in metals $\dot{M}_{sml, z}$ are

$$ \dot{M}_{sml} = R \cdot SFR \quad (34) $$

and

$$ \dot{M}_{sml, z} = y(1-R) \cdot SFR. \quad (35) $$

In equation (35), $1-R$ is the fraction of the star-forming gas that remains in stars and contributes to the final stellar masses
<table><caption>**Table 1.** Model parameters: symbols, units (for dimensional quantities), default values and alternative models (where they differ).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Units</th><th>Default</th><th>Model 1</th><th>Model 2</th><th>Model 3</th></tr></thead><tbody><tr><td colspan="7"><strong>Cosmology</strong></td></tr><tr><td>Matter density</td><td>Ω<sub>M</sub></td><td></td><td>0.308</td><td></td><td></td><td></td></tr><tr><td>Baryon density</td><td>Ω<sub>b</sub></td><td></td><td>0.0481</td><td></td><td></td><td></td></tr><tr><td>Cosmological constant</td><td>Ω<sub>Λ</sub></td><td></td><td>0.692</td><td></td><td></td><td></td></tr><tr><td>Hubble constant</td><td>H<sub>0</sub></td><td>100 km s<sup>-1</sup> Mpc<sup>-1</sup></td><td>0.678</td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>N-body simulation</strong></td></tr><tr><td>Box size</td><td>L<sub>box</sub></td><td>Mpc</td><td>100</td><td></td><td></td><td></td></tr><tr><td>Resolution</td><td>N<sub>part</sub></td><td></td><td>512<sup>3</sup></td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>Dimensional parameters</strong></td></tr><tr><td>Thermal velocity dispersion of IGM</td><td>σ<sub>reio</sub></td><td>km s<sup>-1</sup></td><td>25</td><td></td><td></td><td></td></tr><tr><td>Star formation threshold</td><td>Σ<sub>th</sub></td><td>M<sub>☉</sub> pc<sup>-2</sup></td><td>2</td><td></td><td></td><td></td></tr><tr><td>Minimum shock heating mass</td><td>M<sub>shock</sub></td><td>M<sub>☉</sub></td><td>10<sup>10.7</sup></td><td>10<sup>11.3</sup></td><td></td><td>10<sup>11.0</sup></td></tr><tr><td>Shutdown mass</td><td>M<sub>shutdown</sub></td><td>M<sub>☉</sub></td><td>10<sup>12.7</sup></td><td>10<sup>12.4</sup></td><td></td><td>10<sup>12.3</sup></td></tr><tr><td>SN feedback saturation scale</td><td>ν<sub>SN</sub></td><td>km s<sup>-1</sup></td><td>24</td><td></td><td></td><td>43</td></tr><tr><td>SN energy</td><td>E<sub>SN</sub></td><td>erg</td><td>10<sup>51</sup></td><td></td><td></td><td></td></tr><tr><td>SN rate</td><td>Ψ<sub>SN</sub></td><td>M<sub>☉</sub><sup>-1</sup></td><td>1/140</td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>Adimensional efficiency factors</strong></td></tr><tr><td>Star formation</td><td>ε<sub>sf</sub></td><td></td><td>0.04</td><td></td><td></td><td></td></tr><tr><td>Disc instabilities</td><td>ε<sub>inst</sub></td><td></td><td>0.9</td><td></td><td></td><td></td></tr><tr><td>Mass ratio for major mergers</td><td>ε<sub>m</sub></td><td></td><td>4</td><td></td><td></td><td>3</td></tr><tr><td>Structure coefficient of bulges</td><td>C</td><td></td><td>2.5</td><td></td><td></td><td></td></tr><tr><td>Black hole accretion</td><td>ε<sub>•</sub></td><td></td><td>0.01</td><td></td><td></td><td></td></tr><tr><td>Maximum SN feedback efficiency</td><td>ε<sub>max</sub></td><td></td><td>1</td><td></td><td>0.12</td><td>0.12</td></tr><tr><td>Returned fraction</td><td>R</td><td></td><td>0.45</td><td></td><td></td><td></td></tr><tr><td>Metal yield</td><td>y</td><td></td><td>0.06</td><td></td><td></td><td></td></tr><tr><td colspan="7"><strong>SN feedback scaling exponents</strong></td></tr><tr><td>ν<sub>vir</sub> scaling</td><td>α<sub>v</sub></td><td></td><td>-4</td><td></td><td></td><td>-6.2</td></tr><tr><td>z scaling</td><td>α<sub>z</sub></td><td></td><td>3</td><td></td><td></td><td></td></tr></tbody></table>

of galaxies, while y is the metal mass ejected into the interstellar medium per unit mass locked into stars.

Metal enrichment has been included in GalICS 2.0 to pave the way for future developments but has no effect whatsoever on any of the results presented in this paper because we are not computing cooling or any properties that depend on the spectral energy distribution of galaxies, such as magnitudes and colours.

## 2.6 GAS

The GAS module defines what composes a gas. Currently, an object of type gas has only two attributes: total mass and metal mass. The metal yield y in STAR determines the metallicity of gas returned to the interstellar medium. This is the only place where metals enter GalICS 2.0 explicitly outside the GAS module. The reason is that, whenever an object of type gas is transferred from one gas component to another, its metals are transferred with it automatically in a manner completely transparent to the other modules.

## 2.7 Summary of parameters and models explored

Table 1 summarizes the parameters of the models considered in this paper. The first two sets of parameters (cosmology and N-body simulation) are set by the N-body simulation used to build the halo catalogues and merger trees. They are not free parameters of the SAM.

There are 16 free parameters in GalICS 2.0. We have separated them into three groups: dimensional parameters, efficiency factors and scaling exponents. The dimensional parameters set the

characteristic surface density, velocity and mass scales for star formation, stellar feedback and shock heating. Efficiency factors and scaling exponents are dimensionless. The former are multiplicative factors that set the efficiency of a process (star formation, bar formation, bulge formation, black hole accretion). The latter determine the exponent of the power law with which a quantity depends on another (in our case, how the energetic efficiency of SN feedback scales with virial velocity and redshifts).

Some parameters, such as those related to SN feedback, are highly uncertain. Others are reasonably well constrained by observations, previous models and simulations.

$E_{SN}$, $\Psi_{SN}$, R and y are determined by stellar evolution. We have used the values for a Chabrier (2003) initial mass function with Romano et al. (2010)'s stellar yields (see Vincenzo et al. 2016), and we have not allowed them to vary.

Okamoto, Gao & Theuns (2008) find that cosmic reionization suppresses gas accretion on to haloes up to $M_{reio} = 6.5 \times 10^9 h_0^{-1} M_{\odot}$ at $z=0$. This mass corresponds to $v_{vir} \sim 25$ km s$^{-1}$. We therefore assume $\sigma_{reio} = 25$ km s$^{-1}$ but notice that this assumption will have little consequence on our results as this scale is below the resolution of our N-body simulation.

Cosmological hydrodynamic simulations by Ocvirk et al. (2008) find $M_{shock} \sim 10^{10.7} M_{\odot}$ and $M_{shutdown} \sim 10^{12.7} M_{\odot}$.

The two parameters $\epsilon_m$ and $\epsilon_{inst}$ affect morphology only. For the formation of bulges, simulations of galaxy mergers have consistently found that the critical mass ratio that separates major and minor mergers is of $\epsilon_m = 3-5$. We choose $\epsilon_m = 4$ as our default value. For disc instabilities, Efstathiou et al. (1982) used N-body simulations and showed that $\epsilon_{inst} \simeq 0.91$. Christodoulou et al. (1995)
argued for a lower value ($\epsilon_{\text{inst}} \simeq 0.83$) in stellar discs but remarked that the presence of gas could raise $\epsilon_{\text{inst}}$. The values reported above are the inverse of those contained in the original papers because of our different definition of $\epsilon_{\text{inst}}$. Furthermore, the original papers used a global instability criterion at the radius $r = 2.15r_d$ at which the rotation curve of a self-gravitating exponential disc peaks. Hence, there is no reason why their values should apply to our SAM, in which equation (19) is applied at each radius to find the one, if any, at which an instability develops. However, in Section 3.5, we shall show that assuming $\epsilon_{\text{inst}} = 0.9$ leads to morphologies in reasonable agreement with observations.

The structure coefficient of bulges (defined so that $\sigma_e$ is the average one-dimensional velocity dispersion within an aperture $R_e$; Section 2.3.5) is a parameter of the model but not a free one. Observations of early-type galaxies find $C = 2.5$ on average (Cappellari et al. 2006, 2013). Cappellari et al. (2006) remarked that, in a self-consistent model, this value corresponds to a profile with Sérsic (1963) index $n = 5.5$, but warned that this conclusion is based on assuming a spherical isotropic system with uniform mass-to-light ratio.

We have set the gas fraction that accretes on to the central black hole in major mergers to $\epsilon_* = 0.01$ because we know from experience with different codes (Cattaneo 2001; Cattaneo et al. 2005b) that this value is in good agreement with the black hole–bulge mass relation (Magorrian et al. 1998; Marconi & Hunt 2003; Häring & Rix 2004).

The greatest uncertainty is in the parameters that control the efficiency of SN feedback ($v_{\text{SN}}$, $\alpha_v$, $\alpha_z$, $\epsilon_{\text{max}}$). They are the only true free parameters of our model and they have been calibrated on the galaxy SMF in the local Universe. Our default assumption is $\epsilon_{\text{max}} = 1$. When all the other parameters are kept fix to their default values, the best fit to the local SMF is found for $v_{\text{SN}} = 24 \text{ km s}^{-1}$, $\alpha_v = -4$ and $\alpha_z = 3$ (Section 3.1).

In this paper, we have explored three models in addition to our default parameter combination. The corresponding parameters are listed in Table 1 when they differ from the default values. Model 1 corresponds to a more abrupt shutdown of cold accretion ($M_{\text{shock}}$ has been increased to $M_{\text{shock}} \sim 10^{11.3} M_{\odot}$ and $M_{\text{shutdown}}$ has been lowered to $M_{\text{shutdown}} \sim 10^{12.4} M_{\odot}$). Model 2 limits the efficiency of SN feedback to $\epsilon_{\text{max}} = 0.12$ (models with $\epsilon_{\text{max}} > 0.2$ are indistinguishable from the default model because there are few galaxies with an energetic efficiency of SNe $> 20$ per cent). Model 3 produces an SMF that is flat around $M_{\text{stars}} \sim 10^{10} M_{\odot}$ and rises steeply at lower masses in agreement with observations by Baldry, Glazebrook & Driver (2008) and Baldry et al. (2012). The default model predicts a shape of the SMF in better agreement with that of Bernardi et al. (2013, Section 3.1).

In addition to the four main models in Table 1, we have run four other models to test our sensitivity to specific assumptions. Their results are shown only in connection with the relevant figures. They are a model without disc instabilities (Section 3.5, morphologies), a model in which all haloes have the same spin parameter $\lambda = 0.05$ (Section 3.6, disc sizes), a model in which the halo concentration parameter $c$ is computed with the fitting formulae of Dutton & Macciò (2014) rather than by using the values measured in our N-body simulation and the same model when we also include adiabatic contraction (Section 3.7, TFR).

## 3 COMPARISON WITH OBSERVATIONS

In this section, we compare the models in Table 1 with observations. However, let us start with a foreword on our plotting conventions.

In most figures (SMFs, SFR function, cosmic SFR density, early-type fraction, disc sizes, stellar and baryonic TFR, Faber–Jackson relation), model predictions are shown by curves and observations by points with error bars. In these figures, black solid curves correspond to the default model, red curves to model 1, black dashed curves to model 2 and green curves to model 3. Some correlations ($M_{\text{stars}}$ versus $M_{\text{vir}}$, gas-to-stellar mass ratio versus $M_{\text{stars}}$, SFR versus $M_{\text{stars}}$) are shown as scatter plots. These figures are shown for the default model and model 3 only, the former with a black point cloud, the latter with a green one.

In this paper, we compare our predictions to derived data (stellar masses, SFRs) rather than to primary data (magnitudes, colours) because GalICS 2.0 has not been interfaced with stellar population synthesis models yet. All data have been corrected for a Hubble constant of $H_0 = 67.8 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and a Chabrier (2003) initial mass function.

### 3.1 Mass functions

Fig. 3 compares the galaxy SMFs predicted by the models in Table 1 (curves) to observations at different redshifts in the range $0 < z < 2.5$ (data points with error bars). We begin our analysis from the local Universe ($z \simeq 0.1$), where there are also data for the baryonic mass function (Fig. 4; the baryonic mass is the total mass of stars and cold neutral gas). In Fig. 4, we also show once and for all which line corresponds to each model.

Without adjusting any parameter besides $v_{\text{SN}}$, $\alpha_v$ and $\alpha_z$, which are otherwise completely undetermined, the default model is in good agreement with the local SMFs by Baldry et al. (2012), Bernardi et al. (2013), Moustakas et al. (2013) and Yang, Mo & van den Bosch (2009) in the mass range $10^9 M_{\odot} < M_{\text{stars}} < 10^{11.8} M_{\odot}$.

The SMFs observed by different authors are overall fairly similar. The most noteworthy difference is that the SMF of Bernardi et al. (2013) contains a larger number of massive galaxies. The reason is that Bernardi et al. (2013) did not use the photometry from the SDSS pipeline. They fitted the surface brightness profiles of galaxies with the combination of a Sérsic (1963) and an exponential profile, from which they computed magnitudes by extrapolating it to infinity.

At $M_{\text{stars}} > 10^{11.8} M_{\odot}$, we overpredict galaxy number densities by a factor of 2 even with respect to Bernardi et al. (2013). This could be an effect of cosmic variance (the halo mass function shows an excess of objects by $\sim 40$ per cent at $M_{\text{vir}} \sim 3 \times 10^{13} M_{\odot}$; Fig. 1), exacerbated by overmerging because we have assumed that subhalo mergers result in immediate galaxy mergers, though tests based on a beta version with delayed merging show that this effect cannot be large.

Model 1 ($M_{\text{shock}} = 10^{11.3} M_{\odot}$ and $M_{\text{shutdown}} = 10^{12.4} M_{\odot}$; red curves) corresponds to a more abrupt shutdown of cold accretion, reflected in a sharp change of slope of the SMF just below $M_{\text{stars}} = 10^{11} M_{\odot}$ (red curve at $z \simeq 0.1$). This model is in better agreement with the shape of the SMF measured by Yang et al. (2009) and Bernardi et al. (2013) at $10^9 M_{\odot} < M_{\text{stars}} < 10^{11.6} M_{\odot}$.

Model 2 is identical to the default model except that we have limited the efficiency of SN feedback to $\epsilon_{\text{max}} = 0.12$. Capping the efficiency of SN feedback raises the slope of the galaxy SMF at low masses. This improves the agreement with $z \simeq 0.1$ data (compare the black dashes with the data at $M_{\text{stars}} \lesssim 10^9 M_{\odot}$) but makes things worse at higher $z$.

The SMF of Baldry et al. (2008, 2012) has a different shape from that of Bernardi et al. (2013), even though they are within each other’s error bars everywhere except at the highest masses. While the SMF of Bernardi et al. is consistent with a double power law,
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Model SMFs at $0 < z < 2.5$ (curves) compared with observations (data points with error bars) from several groups (Yang et al. 2009; Baldry et al. 2012; Bernardi et al. 2013; Ilbert et al. 2013; Moustakas et al. 2013; Muzzin et al. 2013; Tomczak et al. 2014). Here, $\phi$ is the number density of galaxies per dex of stellar mass. Black solid curves, red curves, black dashed curves and green curves correspond to the default model, model 1, model 2 and model 3, respectively. The parameter values for each model are in Table 1. The magenta curve at $z \simeq 0.1$ is the SMF in the EAGLE simulation (Schaye et al. 2015, Ref-L100N1504). All models have been convolved with a Gaussian random error of $0.04(1+z)$ dex on stellar masses.</figcaption></figure>

that of Baldry et al. is steeper at $M_{\text{stars}} < 10^{9.5} M_{\odot}$, almost flat at $10^{9.5} M_{\odot} < M_{\text{stars}} < 10^{10.5} M_{\odot}$ and then drops more rapidly at higher masses. Model 3 corresponds to a combination of parameters that was chosen to reproduce this behaviour. We make the

dependence of feedback on $v_{\text{vir}}$ stronger (we pass from $\alpha_v = -4$ to $-6.2$) so that the SMF becomes almost flat at intermediate masses but then we cap the efficiency of SN feedback at $\epsilon_{\text{max}} = 0.12$ to increase the low-mass slope, like in model 2. We also raise $M_{\text{shock}}$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The baryonic mass function in GalICS 2.0 (curves) and in the observations by Papastergis et al. (2012, grey shaded area). Here, $\phi$ is the number density of galaxies per dex of baryonic mass. The black solid curve corresponds to our default model. The red curve, the black dashed curve and the green curve correspond to models 1, 2 and 3 of Table 1, respectively. In the observations, $M_{\text{gal}} = M_{\text{stars}} + 1.4M_{\text{H1}}$. The lower boundary of shaded area is based on measured H1 masses (ALFALFA). The upper boundary is computed based on the maximum H1 masses that the galaxies could contain given the survey's detection limit.</figcaption></figure>

and lower $M_{\text{shutdown}}$ (a bit like in model 1) to make the change of slope around $M_{\text{stars}} \sim 10^{11} M_{\odot}$ more pronounced.

The turnovers seen at low masses in model 3 are caused by the limited resolution of the *N*-body simulation and provide a measure of the real stellar mass up to which resolution effects can propagate. This mass (of nearly $10^9 M_{\odot}$) is an order of magnitude larger than the formal resolution limit (defined, in Section 3.2, as the $M_{\text{stars}}$ corresponding to the minimum halo mass resolved by the *N*-body simulation).

At $0.5 < z < 0.8$, the agreement is still good. The default model and model 1 fit better the SMFs by Ilbert et al. (2013) and Muzzin et al. (2013). Models 2 and 3 reproduce better the steep low-mass slope found by Tomczak et al. (2014). The data of Ilbert et al. are closer to those of Muzzin et al. than to those of Tomczak et al. However, Ilbert et al. and Tomczak et al. find the same behaviour at $0.5 < z < 0.8$ that Baldry et al. find in the local Universe: their SMFs flatten around $M_{\text{stars}} \sim 10^{10} M_{\odot}$ and steepen again at lower masses. The only difference between the SMFs of Ilbert et al. and Tomczak et al. is that, at $M_{\text{stars}} < 10^{11} M_{\odot}$, the one of Tomczak et al. is shifted to higher masses by 0.3 dex on average. In contrast, the SMF of Muzzin et al. displays a single slope at $M_{\text{stars}} < 10^{11} M_{\odot}$, like that of Bernardi et al. in the local Universe.

Once a major challenge for SAMs, reproducing the number density of massive galaxies at high $z$ is no longer a problem when we convolve our theoretical predictions with the observational errors to account for the Eddington bias. Ilbert et al. (2013) quote an error of 0.04(1 + $z$) dex on stellar masses, which they model with a Lorentzian distribution. If we apply this assumption to our results, we find a small tail of galaxies, the masses of which are overestimated by orders of magnitude. We therefore make the conservative assumption that the errors are Gaussian. While an error of 0.04(1 + $z$) dex may not apply to the data of other authors, who have not always stated how their errors vary with redshift, we assume that the errors in the other data sets are of comparable magnitude.

The main discrepancy with the observations is below the knee of the galaxy SMF. At $z > 1$, models begin to overestimate the number density of galaxies with respect to all data sets. The discrepancy is more severe when we limit the efficiency of SN feedback to $\epsilon_{\text{max}} = 0.12$ (models 2 and 3) and is a general problem of all SAMs (Fontanot et al. 2009; Guo et al. 2011; Henriques et al. 2012; Asquith et al., in preparation). **LGALAXIES** (Henriques et al. 2013) is the only model that is marginally consistent with the observations because it combines high ejection rates with a re-accretion time-scale that is inversely proportional to halo mass. In GalICS 2.0, there is no re-accretion because there is no cooling, but simply reintroducing cooling, without a gradual return of gas to the halo, would not solve the problem because the re-accretion time-scale required to make this picture work (6 Gyr time for a halo with $M_{\text{vir}} \sim 3 \times 10^{10} M_{\odot}$) is much longer than the radiative cooling time-scale. Alternative explanations are overefficient star formation in dwarf galaxies in SAMs (but see Section 3.4) or that the observations are missing faint galaxies at high $z$.

Comparing GalICS 2.0 to the baryonic mass function is useful because we can see the extent to which our mass functions are affected by our star formation law. Unfortunately, these data are only available for the local Universe. Papastergis et al. (2012) determined the baryonic mass function from a sample for which both optical (SDSS) and H1 (ALFALFA) data were available. The baryonic mass was assumed to be $M_{\text{gal}} = M_{\text{stars}} + 1.4M_{\text{H1}}$, where $M_{\text{stars}}$ is the stellar mass derived from optical data, $M_{\text{H1}}$ is the H1 mass from radio data and the factor of 1.4 accounts for the presence of helium. This is a lower limit for $M_{\text{gal}}$ because gas could be present and not be detected. One can find an upper limit by giving to galaxies not detected in H1 the maximum H1 mass consistent with their non-detection. The grey shaded area in Fig. 4 shows the region between these limits.

Massive galaxies have low gas fractions (Section 3.3). Therefore, the baryonic mass function is essentially identical to the SMF at high masses. Papastergis et al. (2012)'s SMF is intermediate between Baldry et al. (2008)'s and Yang et al. (2009)'s but closer to the former than to the latter. Model 3 fits the baryonic mass function of Papastergis et al. at $10^{11} M_{\odot} < M_{\text{gal}} < 3 \times 10^{11} M_{\odot}$ better than the other three models because it is calibrated on Baldry et al.'s data. None of the models fits the baryonic mass function at $M_{\text{gal}} > 3 \times 10^{11} M_{\odot}$ but neither do they fit the SMF of Baldry et al. at $M_{\text{stars}} > 3 \times 10^{11} M_{\odot}$.

The agreement of the default model, model 1 and model 2 with the baryonic mass function of Papastergis et al. (2012) is good down to $M_{\text{gal}} \sim 10^{10} M_{\odot}$. At $10^9 M_{\odot} < M_{\text{gal}} < 10^{10} M_{\odot}$, baryonic masses are underestimated by 0.2 dex on average. $M_{\text{gal}} \sim 10^9 M_{\odot}$ is our resolution in baryonic mass, which produces the turnover at low masses seen in all models. Hence, it would make no sense to extend the comparison with the data to lower masses. In model 3, the baryonic mass function is reproduced correctly around $M_{\text{gal}} \sim 10^9 M_{\odot}$ but is underpredicted around $M_{\text{gal}} \sim 10^{10} M_{\odot}$. Overall, the comparison with the baryonic mass function of Papastergis et al. (2012) seems to suggest that, in GalICS 2.0, low-mass galaxies do not contain enough gas for their stellar masses. See however, Section 3.3 for a direct comparison with gas fractions.

## 3.2 Halo masses

Fig. 5 compares the $M_{\text{stars}}-M_{\text{vir}}$ relation predicted by the default model (left, black point cloud) and model 3 (right, green point cloud) with data from weak lensing (Reyes et al. 2012) and satellite kinematics (Wojtak & Mamon 2013). The excellent agreement
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The $M_{vir}-M_{stars}$ relation predicted by the default model (black point cloud, left) and model 3 (green point cloud, right) compared with observational determinations of halo masses from weak lensing (Reyes et al. 2012, cyan squares) and the kinematics of satellite galaxies (Wojtak & Mamon 2013). The results from the latter are shown separately for late-type galaxies (blue triangles) and early-type galaxies (red circles). The yellow point cloud shows the relation for the baryonic rather than the stellar mass. The black and yellow arrows show our resolution in $M_{stars}$ and $M_{gal}$, respectively. The dot-dashed line corresponds to $M_{stars} = (\Omega_b/\Omega_M)M_{vir}$. The vertical dashed line shows our halo mass resolution.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

(particularly with weak lensing data) provides additional and independent evidence that local haloes harbour galaxies with sensible stellar masses, at least for systems with $M_{stars} \gtrsim 10^{10} M_{\odot}$. On the other hand, Fig. 5 shows that these data are not very powerful to discriminate between models.

The dot-dashed diagonal line in Fig. 5 corresponds to $M_{stars} = (\Omega_b/\Omega_M)M_{vir}$. It has been plotted to emphasize that the $M_{stars}-M_{vir}$ relation is steeper than linear at $M_{vir} < 10^{12} M_{\odot}$ and shallower at higher masses. The logarithmic slope of the relation varies from ~2 at low mass to ~1/2 at high mass. This change is the reason why we observe a knee in the galaxy SMF (Fig. 3). In our model, it results from the combined effects of SN feedback (at low masses) and shock heating (at high masses), both of which limit the baryon mass that can be converted into stars.

This point is illustrated in Fig. 6: $1/(1+\eta)$ is the retained (not ejected) gas fraction (blue curve); $1-f_{hot}$ is the gas fraction that is able to accrete on to galaxies (red curve). Both are shown as a function of $M_{vir}$. The gas that makes stars is the one that is able to accrete and to avoid ejection. Its mass fraction, $(1-f_{hot})/(1+\eta)$, is shown by the black curve. The black curve peaks for $M_{vir} \simeq 10^{11.6} M_{\odot}$ in the default model and $M_{vir} \simeq 10^{11.5} M_{\odot}$ in model 3. Fig. 6 shows that the formation of galaxies is efficient only in a narrow range of halo masses, between $10^{11}$ and a few times $10^{12} M_{\odot}$, in agreement with previous studies by Bouché et al. (2010), Guo et al. (2010), Cattaneo et al. (2011), Behroozi, Wechsler & Conroy (2013) and Birrer et al. (2014). Even in this mass range, it is very difficult for haloes to convert more than a third of the baryons into stars (Fig. 5).

The yellow points in Fig. 5 show the $M_{gal}-M_{vir}$ relation for the baryonic rather than the stellar mass. The vertical dashed lines correspond to the resolution of the N-body simulation, while the black and the yellow arrows mark the formal resolution masses for the stellar mass and the baryonic mass, respectively. We define them as the masses $M_{stars}$ and $M_{gal}$ that correspond to our halo mass resolution ($M_{vir} = 10^{10.5} M_{\odot}$) on a galaxy mass-halo mass diagram (Fig. 5). Both the default model and model 3 have formal resolution $M_{stars} \lesssim 10^8 M_{\odot}$. Hence, the SMFs in Fig. 3 should be well resolved over the entire plotting range, while the baryonic mass function in Fig. 4 is expected to be only at $M_{gal} > 10^9 M_{\odot}$. However, resolution effects can trickle above the formal resolution mass (e.g. Cattaneo et al.

2011). The turnovers in SMFs between $10^8$ and $10^9 M_{\odot}$ are a clear sign of that. We therefore only trust results above $10^9 M_{\odot}$ for both stellar and baryonic masses.

### 3.3 Gas fractions

The default model's predictions for the gas-to-stellar mass ratio as a function of stellar mass are in good agreement with the measurements of $(M_{H1} + M_{H2})/M_{star}$ in a volume-limited, K-band-selected sample of nearby late-type galaxies (Boselli et al. 2014; Fig. 7, left). Unsurprisingly, these measurements find slightly higher gas-to-stellar mass ratios than studies based on H1 data only (Garnett 2002; Swaters & Balcells 2002; Noordermeer et al. 2005; Zhang et al. 2009; also shown in Fig. 7), but the differences are small.

The H1 measurements of Swaters & Balcells (2002), Garnett (2002) and Noordermeer et al. (2005) are all for spiral or irregular galaxies, but their results are not systematically different from those of Zhang et al. (2009), who did not operate any morphological selection [Zhang et al. did not measure the H1 mass of each individual galaxy in their SDSS sample but inferred it from its colour and luminosity using an empirical relation calibrated on 800 galaxies with optical photometry from the SDSS and H1 masses from the HyperLeda catalogue of Paturel et al. (2003)].

While the point clouds that show the results of GalICS 2.0 in Fig. 7 are composed of galaxies on which no selection has been performed, galaxies with so little gas that would not be detected in H1 have not been included in the calculation of the mean gas fraction (shown by the black curve for the default model and the green curve for model 3; averages are logarithmic). This selection, based on assuming a minimum detectable gas mass of $M_{gas} = 10^{7.5} M_{\odot}$ for a galaxy with $M_{stars} = 10^8 M_{\odot}$ and $M_{gas} = 10^{7.5} M_{\odot}$ for a galaxy with $M_{stars} = 10^{11} M_{\odot}$ (Boselli, private communication), is effectively equivalent to a morphological selection.

The direct measurements in Fig. 7 should be compared to indirect constraints on gas-to-stellar mass ratios from the baryonic mass function (Papastergis et al. 2012) and the TFR relation in dwarf galaxies (Papastergis, Adams & van der Hulst 2016). Fig. 4 suggests that, at $M_{gal} \sim 3 \times 10^9 M_{\odot}$, the default model underestimates
<figure><img src="image_10.png" /><figcaption>**Figure 6.** The baryon fraction available for star formation (black curve) is the product of the fraction that can accrete on to the galaxy ($1 - f_{\text{hot}}$; red curve) times the one that is not ejected ($1/(1+\eta)$; blue curve).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** The local $M_{\text{gas}}/M_{\text{star}}$–stellar mass relation in the default model (black point cloud, left), in model 3 (green point cloud, right) and in the observations (data points with error bars, identical in both panels). The cyan triangles and squares are $H_1+H_2$ data from Boselli et al. (2014). The triangles are for a constant CO-to-$H_2$ conversion factor. The squares are for a conversion factor that depends on luminosity. The other symbols are observational determinations of $M_{\text{gas}}/M_{\text{star}}$ based on $H_1$ data only (blue: Swaters & Balcells 2002; green: Garnett 2002; red: Noordermeer et al. 2005; magenta: Zhang et al. 2009). The curves (black to the left, green to the right) show the mean value of Log (Mgas/Mstars) in bins of Log Mstars for galaxies that contain enough gas to be detected in $H_1$. The detection limit is about $M_{\text{gas}} = 10^{7.5} M_{\odot}$ at $M_{\text{stars}} = 10^8 M_{\odot}$ and $M_{\text{gas}} = 10^{7.5} M_{\odot}$ at $M_{\text{stars}} = 10^{11} M_{\odot}$ (Boselli, private communication). The $H_1$ detection limit for stellar masses intermediate between those above has been determined by logarithmic interpolation. Only one model galaxy out of five has been shown not to overcrowd the plot.</figcaption></figure>

$M_{\text{gal}}$ by ~0.2 dex on average, possibly because the gas content is underpredicted, since the SMF by Papastergis et al. (2012) is intermediate between those of Baldry et al. (2008, 2012) and Yang et al. (2009). TFR informs us on gas fractions because there are data sets for which we have both $M_{\text{stars}}$ and $M_{\text{gal}}$ as a function of rotation speed but the interpretation of these data is not straightforward and will be discussed at length in Section 3.7.

## 3.4 Star formation rates

Both the default model and model 3 are broadly consistent with the slope and the normalization of the $M_{\text{star}}$–SFR relation in the local Universe (Elbaz et al. 2007; Salim et al. 2007; Wuyts et al. 2011; Fig. 8). On a closer inspection, however, both have their shortcomings. The slope of the main sequence of star-forming galaxies is correctly reproduced by the standard model but appears tilted in model 3. In contrast, model 3 reproduces correctly the characteristic mass $M_{\text{stars}} \sim 10^{10.7}-10^{10.8} M_{\odot}$ at which galaxies migrate from the star-forming population to the passive one according to Wuyts et al.

(2011). In the default model, this mass is overestimated because galaxies in the mass range $M_{\text{stars}} \sim 10^{10.7}-10^{10.8} M_{\odot}$ have still got plenty of gas (Fig. 7).

In both cases, one remarks a small but non-negligible tail of star-forming galaxies at $M_{\text{stars}} > 3 \times 10^{11} M_{\odot}$. This tail is evidence that shutting down cold accretion may not be enough and that an additional quenching mechanism (e.g. quasar feedback) is probably needed [but see Bildfell et al. (2008), who argued for reactivated star formation in cD galaxies].

Fig. 9 compares the local ($z \simeq 0.1$) SFR functions in our four models with the data by Gruppioni et al. (2015) at $0 < z < 0.3$. The most interesting differences are between the default model and model 3. Overall, model 3 predicts lower SFRs than the default model, which fits the SFR function better at all but the highest SFRs.

Having examined star formation in the local Universe, we move our attention to the evolution of the cosmic SFR density across the Hubble time (Fig. 10). Measured SFR densities were taken from the compilation by Behroozi et al. (2013), which includes $H\alpha$ (Shim
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The local SFR-stellar mass relation in the default model (black point cloud, left) and model 3 (green point cloud, right) compared to observations. The observations are the same in both panels. The thick dashed and solid lines delimit the main sequence of star-forming (SF) galaxies as determined observationally by Elbaz et al. (2007) and Salim et al. (2007), respectively. The contours show the distribution of SDSS galaxies in the SFR-$M_{\text{stars}}$ plane according to Wuyts et al. (2011). Only one model galaxy out of five has been shown not to overcrowd the plot. The red dots are galaxies with zero SFR. We have assigned them an arbitrary SFR to be able to show them on the diagram.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** The SFR function in the local Universe. The curves show the prediction of GalICS 2.0 at $z \simeq 0.1$ (black solid: default model; red: model 1; black dashed: model 2; green: model 3). The data points with error bars are the observations by Gruppioni et al. (2015) at $0 < z < 0.3$.</figcaption></figure>

et al. 2009; Tadaki et al. 2011; Sobral et al. 2012), radio (1.4 GHz; Dunne et al. 2009; Smolčić et al. 2009; Karim et al. 2011), combined ultraviolet/infrared (Salim et al. 2007; Zheng et al. 2007; Kajisawa et al. 2010; Cucciati et al. 2012; Gruppioni et al. 2015), ultraviolet data only (Yoshida et al. 2006; van der Burg et al. 2010; Robotham & Driver 2011; Bouwens et al. 2012) and far-infrared data only (Le Borgne et al. 2009; Rujopakarn et al. 2010). The default model is in good agreement with observation at $z > 2$ (especially with UV data) because it has strong feedback that suppresses star formation but forms too many stars at low $z$ (Figs 8 and 9). In model 3, feedback is capped to $\epsilon_{\text{max}} = 0.12$. Therefore, even though the efficiency scales as $\epsilon_{\text{SN}} \propto (1+z)^3$, its increase at high $z$ cannot be as strong as in the default model and the cosmic SFR density is overpredicted. However, the lower value of $M_{\text{shutdown}}$ curbs star formation in massive galaxies much more effectively and this improves the fit at low $z$. In fact, the agreement of model 3 with the observations at $0 < z < 2$ is unprecedented.

## 3.5 Morphologies

The data points with error bars in Fig. 11 show the fraction of elliptical galaxies. They are based on galaxies visually classified by Conselice (2006), who separated them into E, S and Irr types. In order to compare the results of GalICS 2.0 to these observational data, it is necessary to determine to which bulge-to-total mass ratios ($B/T$) this classification corresponds. Weinzirl et al. (2009) have shown that S-type galaxies have $B/T$ in the range 0–0.5 (two thirds have $B/T < 0.2$). The range extends to 0–0.7 if we broaden our definition of S-type galaxies to include S0s (Laurikainen et al. 2010). We therefore follow Wilman et al. (2013) and Fontanot et al. (2015) in using $B/T = 0.7$ as the watershed $B/T$ that separates S- and E-type galaxies.

GalICS 2.0 computes $B/T$ assuming that the bulge mass is the total stellar mass of the classical bulge and the pseudo-bulge. The total mass is the sum of the stellar masses of the disc, the classical bulge and the pseudo-bulge. The fraction of galaxies with $B/T > 0.7$ increases with $M_{\text{stars}}$, first more gently at $M_{\text{stars}} < 10^{11} M_{\odot}$, then more rapidly at $M_{\text{stars}} > 10^{11} M_{\odot}$, where mergers become the main mechanism of galaxy growth (e.g. Cattaneo et al. 2011) and bulge formation. If mergers are the only mechanism to form bulges, a significant population of galaxies with $B/T > 0.7$ appears only for $M_{\text{stars}} > 1-3 \times 10^{11} M_{\odot}$ (the dot-dashed curve in Fig. 11 corresponds to the default model without disc instabilities; also see Lacey et al. 2016).

When disc instabilities are activated alongside morphological transformation in major mergers, the default model (black solid curve) and model 2 (black dashed curve) reproduce a good agreement with the data points in Fig. 11 for $\epsilon_m = 4$ and $\epsilon_{\text{inst}} = 0.9$. In contrast, model 1 (red curve) and model 3 (green curve) overestimate the fraction of massive ellipticals even when the critical mass ratio for major mergers is lowered to $\epsilon_m = 3$. These differences come from the gas mass that accretes on to galaxies because the merger rate is set by the DM and is the same in all models (in models with a larger value of $M_{\text{shutdown}}$, even galaxies as large as the Milky Way have a chance to regrow a disc after a merger).

Fig. 11 proves that our morphologies are reasonable. However, as the association of visually classified elliptical galaxies with a
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Evolution of the cosmic SFR density across the Hubble time in GalICS 2.0 (curves) and the observations (points with error bars from the compilation by Behroozi et al. 2013). Yellow symbols show Hα data (squares: Sobral et al. 2012; triangles: Tadaki et al. 2011; circles: Shim et al. 2009), black symbols 1.4 GHz data (triangles: Karim et al. 2011; squares: Dunne et al. 2009; down-pointed triangles: Smolčić et al. 2009), red symbols combined UV/IR data (stars: Gruppioni et al. 2015; squares: Cucciati et al. 2012; circles: Kajisawa et al. 2010; down-pointed triangles: Salim et al. 2007; triangles: Zheng et al. 2007), blue symbols UV data only (squares: Bouwens et al. 2012; down-pointed triangles: Robotham & Driver 2011; triangles: van der Burg, Hildebrandt & Erben 2010; circles: Yoshida et al. 2006) and green symbols FIR/IR data only (circles: Rujopakarn et al. 2010; squares: Le Borgne et al. 2009). The line styles for the models are the same as in Fig. 4.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** The fraction of galaxies with bulge-to-total mass ratio $B/T > 0.7$ as a function of stellar mass in GalICS 2.0 (curves) compared to the fraction of elliptical galaxies in the observations of Conselice (2006). As in Figs 3, 4 and 9, the black solid curve refers to the default model, while the red curve, the black dashed curve and the green curve refer to models 1, 2 and 3, respectively. Here we also show a fourth model, identical to the default one, except that we turn off disc instabilities (black dot-dashed curve).</figcaption></figure>

critical bulge-to-total mass ratio of $B/T = 0.7$ is somewhat arbitrary, a quantitative comparison with an observational sample with measured $B/T$ is highly desirable. Nevertheless, the only purpose of morphologies in this paper is to select spiral galaxies when

comparing to observations for disc sizes (Section 3.6) and the TFR (Section 3.7). As we have verified that neither disc sizes nor the TFR were sensitive to the critical $B/T$ used to select spiral galaxies (choosing $B/T < 0.3$ or $B/T < 0.7$ does not change the TFR significantly), we have decided to defer the comparison with quantitative morphologies to a future publication.

## 3.6 Disc sizes

In the last section, we checked that our morphologies are reasonable. Now, we focus on spiral galaxies and, in particular, on the disc mass-size relation. The solid curves in Fig. 13 show the predictions of the default model (±1 standard deviation from the mean) for galaxies with $B/T < 0.3$ at $0 < z < 2.5$ (see Section 3.5 for a discussion of $B/T$ in GalICS 2.0).

The criterion $B/T < 0.3$ has been chosen to match van der Wel et al. (2014)'s morphological selection by Sérsic (1963) index (Huertas-Company, private communication). Concerning the other data sets used for this comparison, Shibuya et al. (2015) selected late-type galaxies based on star formation, while Bruce et al. (2014) at high $z$ and Bernardi et al. (2014) in the local Universe performed a bulge/disc decomposition. Their results are therefore more directly comparable to ours. Also notice that both van der Wel et al. and Bruce et al. based their investigations on CANDELS data. The local data from Boselli et al. (2014) are individual spiral galaxies from the Herschel Reference Survey.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** The distribution for the spin parameter $\lambda$ in our N-body simulation (solid curve) compared to a lognormal distribution with $\bar{\lambda} = 0.049$ and $\sigma_{\ln\lambda} = 0.57$ (dashed curve).</figcaption></figure>

The data points sit comfortably in the range predicted by the default model, with the only possible exception of the most massive discs in the local Universe (at high $M_d$, theoretical predictions are affected by poor statistics because of the decline of the galaxy SMF combined to the increase of $B/T$ with $M_{\text{stars}}$). However, as de Jong & Lacey (2000) had already found in an earlier SAM that computed disc radii from the halo spin distribution measured in N-body simulations, the scatter is much larger in GalICS 2.0 than in the observations. In fact, it is so large that it covers any difference between models. Hence, in Fig. 13, only the default model has been shown.

Most of the scatter comes from the halo spin parameter. To prove it, we have rerun the default model using $\lambda = 0.05$ for all haloes rather than the values measured in the N-body simulation. The results are shown by the dot-dashed curves in Fig. 13. Some scatter is still present because galaxies differ in halo concentration, $M_{\text{gal}}/M_{\text{vir}}$ (Fig. 5) and $B/T$ (if any of these quantities increases, the rotation speed will increase, too; so, $r_d$ has to shrink if specific angular momentum is to be conserved). However, the disc size-mass relation is much tighter when the scatter in $\lambda$ is removed.

This finding is puzzling because (a) there are many processes and sources of errors that could contribute to the observational scatter and that our model does not include, and (b) the spin distribution in our N-body simulation is in agreement with previous studies. We fit our distribution for $\lambda$ with a lognormal distribution with $\bar{\lambda} = 0.049$ and $\sigma_{\ln\lambda} = 0.57$ (Fig. 12), in agreement with Muñoz-Cuartas et al. (2011), who find $\bar{\lambda} = 0.044$ and $\sigma_{\ln\lambda} = 0.57$, and Burkert et al. (2016), who find $\bar{\lambda} = 0.052$ and $\sigma_{\ln\lambda} = 0.46$.

The most likely explanation is that specific angular momentum is not conserved during infall. Sharma & Steinmetz (2005) found that, in adiabatic cosmological simulations, the specific angular momentum distribution is narrower for gas than DM, owing to the presence of counter-rotating material in the latter but not in the former. Kimm et al. (2011) confirmed this finding for the specific angular momentum distribution of the gas and the stars in a Milky Way-type galaxy, which they simulated both with and without feedback. At the end of their simulations ($z = 3$), the galaxy and the halo had the same specific angular momentum on a global scale despite their different internal distributions, but that was not true at all times and nobody has performed a systematic study on a representative volume to determine whether the scatter in the specific angular momentum distributions of discs and haloes is the same. Fig. 13 suggests that

for some yet unexplained reason there is less scatter for discs than for haloes.

We remark that the model with $\lambda = 0.05$ for all haloes (corresponding to the dot-dashed lines in Fig. 13) contains much fewer early-type galaxies than the default model (less than half). Many of the galaxies that become unstable and develop a pseudo-bulge in the default model live in haloes with $\lambda \ll 0.05$. They are systems in which the disc is very concentrated with respect to the DM. Turning off disc instabilities causes the lower envelope of the disc mass-size relation (the lower solid line in the six panels of Fig. 13) to drop considerably because the spin distribution in Fig. 12 implies a tail of galaxies with very small radii. When disc instabilities are activated, these galaxies cease to contribute to the disc mass-size relation because they develop massive pseudo-bulges and therefore no longer satisfy the $B/T < 0.3$ selection criterion.

## 3.7 The Tully–Fisher relation

In this section, we compare the TFR with observations at $z \sim 0$ (Fig. 14) and $z \sim 1$ (Fig. 16), and we explain the physical reason why its shape is reproduced correctly by GalICS 2.0, while it was not by previous models in which the rotation speed was proportional to the virial velocity (e.g. Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999), although the Durham model had already relaxed this assumption in Cole et al. (2000).

The TFR links the stellar mass $M_{\text{stars}}$ or the baryon mass $M_{\text{gal}}$ (i.e. the total mass of stars and cold neutral gas) of a spiral galaxy to its rotation speed $v_{\text{rot}}$. In this paper, $v_{\text{rot}}$ is the circular velocity $v_c$ at the optical radius $r_{\text{opt}} = 3.2r_d$, which contains 83 per cent of the mass of an exponential disc. We choose this definition for consistency with previous work (Cattaneo et al. 2014) and because measuring $v_{\text{rot}}$ at the outer edge of the disc makes our results less sensitive to the real form of the DM density profile, which is likely to differ at the centre from the NFW model assumed in GalICS 2.0 (Moore 1994; Flores & Primack 1996; Persic, Salucci & Stel 1996, but also see Swaters et al. 2003).

Observations of both the stellar (Fig. 14, left) and the baryonic (Fig. 14, right) TFR show that each individual data set is consistent with a single power law within its intrinsic scatter, although different data sets differ in both slope and normalization.

In GalICS 2.0, $v_{\text{vir}} \propto M_{\text{vir}}^{1/3}$. Hence, the $M_{\text{stars}}-v_{\text{vir}}$ relation is entirely determined by the $M_{\text{stars}}-M_{\text{vir}}$ relation, which changes slope at $M_{\text{stars}} \sim 2 \times \sim 10^{10} M_{\odot}$ (Fig. 5). There is no such feature in the TFR. Hence, understanding the TFR comes down to understanding how $v_{\text{rot}}$ depends on $v_{\text{vir}}$.

The models presented in Fig. 14 (curves) include the presence of bulges but let us focus on pure discs to make the interpretation simpler. In this case:

$$v_{\text{rot}}^2 = v_c^2(r_{\text{opt}}) = \frac{GM_{\text{dm}}(r_{\text{opt}})}{r_{\text{opt}}} + 1.11 \frac{GM_d}{r_{\text{opt}}}, \quad (36)$$

where the first addend is the halo contribution and the second addend is the square of the rotation speed at $r_{\text{opt}} = 3.2r_d$ for a self-gravitating exponential disc (Freeman 1970). By using equation (1) for the mass $M_{\text{dm}}(r_{\text{opt}})$ of the DM within $r_{\text{opt}}$ and by defining $x \equiv r_{\text{opt}}/r_{\text{vir}}$, equation (36) can be rewritten as

$$\frac{v_{\text{rot}}^2}{v_{\text{vir}}^2} = \frac{\frac{\log(1+\text{cx})}{x} - \frac{c}{1+\text{cx}}}{\log(1+c) - \frac{c}{1+c}} + \frac{1.1}{x} \frac{M_d}{M_{\text{vir}}}. \quad (37)$$

For $\lambda = 0.05$ and $c = 8$ (the mean values of the spin and concentration parameters), equation (15) gives values of $x$ that range from $x = 0.083$ for a dwarf galaxy ($M_d/M_{\text{vir}} = 0.004$) to $x = 0.063$ for a
<figure><img src="image_17.png" /><figcaption>**Figure 13.** The stellar mass–size (exponential scalelength) relation for discs (galaxies with $B/T < 0.3$) at different $z$ (curves: ±1σ from the mean) in the default model (solid curves) and a simplified version in which all discs have the same spin parameter $\lambda = 0.05$ (dot–dashed curves). The data points with error bars show the mean observational value of $r_d$ in a bin of stellar mass, while the circles without error bars (the local data by Boselli et al. 2014) are measurements for individual galaxies and are colour-coded according to morphology (blue: Sd or Sc; cyan: Sb; yellow: Sa). The data of van der Wel et al. (2014) and Shibuya, Ouchi & Harikane (2015) are for galaxies classified as late type, while Bruce et al. (2014) and Bernardi et al. (2014) performed a bulge/disc decomposition. The solid curves stop at $10^{10.75} M_{\odot}$ because, in the default model, there are not enough discs with $M_d > 10^{10.75} M_{\odot}$ to compute meaningful averages.</figcaption></figure>

Milky Way ($M_d/M_{vir} = 0.04$). For comparison, the standard model $r_d = \lambda r_{vir}/2$ gives $x = 0.08$. The range of typical values extends to $0.03 \lesssim x \lesssim 0.16$ when the scatter in $\lambda$ is considered. In this range, the first term on the right-hand side of equation (37) is a function

that grows slowly from 0.74 at $x = 0.03$ to 1.1 at $x = 0.16$ (Fig. 17, red curve). The second term ranges from being a small correction for dwarf galaxies to being comparable to the first term for galaxies with masses comparable to the Milky Way (blue dashed and solid
<figure><img src="image_18.png" /><figcaption>**Figure 14.** The stellar (left) and baryonic (right) TFR for disc galaxies in the local Universe. The lines show the predictions of GalICS 2.0 (upper and lower quartiles). They are colour-coded as in Fig. 4: black solid, red, black dashed and green curves correspond to the default model and models 1, 2, 3, respectively. We have selected disc galaxies using the criterion $B/T < 0.7$, where $B/T$ is computed considering both the classical and pseudo-bulge mass, as in Fig. 11. The data points correspond to observed galaxies. Those for the stellar TFR (left) are from Bradford, Geha & van den Bosch (2016), Papastergis et al. (2016), Reyes et al. (2011), Pizagno et al. (2005) and Rhee et al. (2004). In the case of the baryonic TFR (right), we have also shown data from Hall et al. (2012), McGaugh (2012), Gurovich et al. (2010) and Yegorova & Salucci (2007). The data points by Gurovich et al. have converted for their different estimate of $v_{rot}$ (they inferred $v_{rot}$ from the H I linewidth at 20 per cent rather than 50 per cent of the peak height).</figcaption></figure>

curves in Fig. 17, respectively). The dependence on $M_d/M_{vir}$ is the reason why GalICS 2.0 reproduces the $v_{rot}-v_{vir}$ relation from studies of rotation curves (Fig. 18, the default model and model 3 are indistinguishable in this respect) and therefore the TFR (Fig. 14).

The dependence of $v_{rot}/v_{vir}$ on $M_{gal}/M_{vir}$ is also the reason why, in GalICS 2.0, the TFR is very tight despite the large scatter in disc radii. The second term on the right-hand side of equation (37) decreases with $x$ and, summed to first one, conspires to produce a rotation curve that is nearly flat at the optical radius (dashed and solid black curves in Fig. 14). Hence, $v_{rot}/v_{vir}$ is not very sensitive to $x$.

Having established that the dependence of $v_{rot}/v_{vir}$ on $M_{gal}/M_{vir}$ is the physical reason why GalICS 2.0 is able to reproduce the shape of the TFR, we can now look in closer details at the predictions of our different models. Model 1 is very similar to the default model and model 2 is very similar to model 3. Hence, we shall focus our discussion of Fig. 14 on the default model (black solid curves) and model 3 (green solid curves).

We start by noting that the TFR in the two models is very similar down to $M_{stars} \sim 10^9 M_{\odot}$. As our results below this mass may be affected by halo mass resolution (discussion in Sections 3.1 and 3.2), the differences between the default model and model 3 should not be overinterpreted. Nevertheless, these differences make sense when one considers the SMFs predicted by these models. In the default model, the SMF has a constant shallow slope in the mass range $10^8 M_{\odot} < M_{stars} < 10^{10} M_{\odot}$ and the stellar TFR has a small concavity below $M_{stars} \sim 10^9 M_{\odot}$. In model 3, the SMF has a steeper faint-end slope in better agreement with the observations ($M_{stars}/M_{vir}$ decreases less rapidly at low masses) and the concavity is absent.

The baryonic TFR becomes slightly concave below $M_{gal} \sim 10^{10} M_{\odot}$ but the concavity turns to a convexity below $M_{gal} \sim 10^9 M_{\odot}$ because, below this mass, the gas surface density is lower than the threshold $\Sigma_{th}$ for star formation (Section 2.4.1 and Table 1). Hence, there is neither star formation nor ejection of gas. It is possible that these concavities and convexities would not be noticeable if measurement errors on masses and speeds were included in our analysis.

The greatest uncertainty in our predictions for the TFR derives from our incapacity to measure concentrations accurately for haloes with $M_{vir} < 10^{12} M_{\odot}$ (Section 2.1.1). To understand the impact that errors on concentration measurements could have on our results, we have rerun the default model replacing our measurements with concentrations computed with the fitting formulae of Dutton & Macciò (2014, Fig. 15, magenta curves). The difference with respect to the TFR obtained using our concentrations (black curves) is negligible at $M_{stars} > 10^{9.5} M_{\odot}$. Below this mass, the higher concentrations from the fitting formulae of Dutton & Macciò (2014) cause galaxies with $M_{stars} > 10^{8.5} M_{\odot}$ to turn with rotation speeds ~15 per cent higher. Hence, the concavity of the stellar TFR in the default model becomes more pronounced.

All the models considered until now assume no adiabatic contraction (our standard assumption throughout this paper). However, the role of adiabatic contraction in the zero-point of the TFR has been a matter of discussion in previous SAM studies (Cole et al. 2000; Guo et al. 2011). We have therefore considered another variant of the default model, in which not only do we use concentrations computed with the fitting formulae of Dutton & Macciò (2014), but also we include adiabatic contraction modelled as in Blumenthal et al. (1986, Section 2.3.3). The TFR for this model is shown by the region within the gold curves in Fig. 15. Adiabatic contraction produces measurable effects only in galaxies with $M_{gal} > 10^{10} M_{\odot}$. At lower masses, $M_{gal}/M_{vir}$ is too low for the baryons to produce any noticeable effect on the DM, even when the higher concentration of the former with respect to the latter is accounted for. Above $M_{gal} > 10^{10} M_{\odot}$, adiabatic contraction increases the DM mass within the galaxy and thus the disc rotation speed. The increase is stronger at higher masses. Hence, adiabatic contraction tilts the slope of the TFR (it causes it to be less steep at high masses).

Assessing the agreement with observations is not straightforward because the answer depends on the data used for the comparison. The measurements by Papastergis et al. (2016) at low masses provide the main observational hint for a concave stellar TFR when combined with those by Bradford et al. (2016) and Reyes et al. (2011) at higher masses, but the same data are perfectly
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Sensitivity of the TFR to errors on concentration measurements and to adiabatic contraction. The default model (the region between the black curves) and the data points are the same as in Fig. 14. The magenta curves show how the stellar and the baryonic TFR vary when the halo concentrations that we measure in our *N*-body simulation are replaced by values obtained from our halo masses by using the fitting formulae of Dutton & Macciò (2014). The gold curves show how adding adiabatic contraction (modelled with equation 18) modifies the magenta curves.</figcaption></figure>

consistent with a single power law if taken in conjunction with those by Rhee et al. (2004), instead. This stresses the danger of combining different data sets without fully understanding their systematics.

The TFRs by Bradford et al. (2016), Papastergis et al. (2016), McGaugh (2012) and Gurovich et al. (2010) are based on HⅠ linewidths (McGaugh uses the rotation speed in the outermost region, where the rotation curve is approximately flat) and are all more or less consistent with one another. In contrast, Hall et al. (2012, HⅠ linewidths), Yegorova & Salucci (2007, resolved rotation curves from Hα and HⅠ data), Pizagno et al. (2005, Hα linewidths) and Rhee et al. (2004, Hα linewidths) find a normalization that is higher by ~0.3 dex in mass [the systematic shift in $M_{\text{stars}}$ between the SMFs of Baldry et al. (2008) and Bernardi et al. (2014)]. Yegorova & Salucci (2007) are the only ones among the aforementioned authors to measure $M_{\text{gal}}$ dynamically by performing a halo/disc decomposition (the other authors measured $M_{\text{stars}}$ by assuming a mass-to-light ratio or by using a stellar population synthesis model). By using resolved rotation curves, they could measure $v_{\text{rot}}$ at exactly $r_{\text{opt}} = 3.2r_{\text{d}}$. The Hα measurements by Reyes et al. (2011) are intermediate but closer to the results of other Hα studies. Therefore, for the same $M_{\text{stars}}$ or $M_{\text{gal}}$, the rotation speeds from HⅠ studies tend to be systematically higher than those from Hα studies.

This finding has a simple explanation. The Hα emission probes star formation, which is mainly concentrated in the inner regions of galaxies ($r < r_{\text{opt}}$), while the HⅠ emission extends beyond the limit $r_{\text{opt}}$ of the stellar disc. The rotation curves of discs with $M_{\text{stars}} \lesssim 10^{10.5} M_{\odot}$ are usually rising (Persic et al. 1996). Hence, the gas seen in HⅠ has higher rotation speeds than the gas seen in Hα. This also explains why the stellar TFR in GalICS 2.0 is intermediate between the HⅠ determination by Bradford et al. (2016) and the Hα determination by Pizagno et al. (2005), though all our models are consistent with the Hα measurements by Reyes et al. (2011).

An important question in relation to the dwarf-galaxy data by Papastergis et al. (2016) is why GalICS 2.0 appears to be in good agreement with them for the baryonic TFR but not for the stellar one. A straightforward interpretation of this finding is that the baryonic masses are computed correctly but star formation has not been efficient enough. This interpretation would suggest that gas-to-stellar mass ratios are overestimated, not underestimated, as the baryonic mass function seemed to suggest (direct measurements of

$M_{\text{gas}}/M_{\text{stars}}$ suggested that gas fractions are reproduced correctly, particularly in the default model). Another explanation is that the discrepancy is a selection effect. Papastergis et al. (2016) selected their galaxies to be heavily gas dominated. Thus, it is normal that, for a fixed baryon mass, their galaxies contain less stars than average. The difficulty with this explanation is that we see a similar behaviour in the data by Bradford et al. (2016), which should be less biased in this sense. A third explanation is that dwarf galaxies have rising rotation curves. Therefore, HⅠ data overestimate the rotation speed at $r_{\text{opt}}$ (see above). If one followed this explanation, one should apply the same consideration to the baryonic TFR and conclude that the gas fractions of dwarf galaxies are indeed underestimated [a possibility that we had previously argued based on the baryonic mass function by Papastergis et al. (2012)].

A final point to consider when comparing models to observations is the implicit assumption that discs are infinitely thin and cold, so that the rotation speed $v_{\text{rot}}$ is equal to the circular velocity $v_c$ required to support a particle on a circular orbit against gravity. However, $v_{\text{rot}}$ can be lower than $v_c$ if random motions contribute to the support against gravity. Therefore, a model like the one with adiabatic contraction, which appears to get the wrong slope for the TFR (gold curves in Fig. 15), might not be incompatible with the data if there were a lesser degree of rotational support in the discs of massive spiral/S0 galaxies.

Fig. 16 shows the evolution of the stellar TFR in the default model from $z \sim 0$ (region between the solid curves) to $z \sim 1$ (grey shaded area) and $z \sim 2-2.5$ (region between the dashed curves). The evolution in the other models is similar. Only the detailed shape of the relation at low masses differs. This evolution is driven not only by a decrease in $M_{\text{stars}}/M_{\text{vir}}$ when moving to higher redshifts (Behroozi et al. 2013; Birrer et al. 2014; Coupon et al. 2015; McCracken et al. 2015; Skibba et al. 2015), which, in GalICS 2.0, is not strong enough to explain the evolution of the GSMF with $z$ (Fig. 3), but also by the increase of $v_{\text{vir}}$ with $z$ for a fixed $M_{\text{vir}}$.

The stellar TFR at $z \sim 1$ can be compared to Hα data from the KMOS Redshift One Spectroscopic Survey (KROSS; Tiley et al. 2016). To quantify the evolution of the TFR from $z \sim 0$ to $z \sim 1$, Tiley et al. fitted the local data by Pizagno et al. (2005), Rhee et al. (2004) and Reyes et al. (2011) with a slope of 3.68, which they then used to fit their $z \sim 1$ data and measure the change
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Predicted evolution of the stellar TFR from $z \sim 0$ to $z \sim 2$ (regions between lines/grey shaded area for $z \sim 1$) and comparison to $z \sim 1$ data (Tiley et al. 2016, points with error bars).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Dependence of $v_{\text{rot}}/v_{\text{vir}}$ on $r_{\text{opt}}/r_{\text{vir}}$ and $M_{\text{d}}/M_{\text{vir}}$. The total value (black curves) is split in the contributions of the halo (red dot-dashed curve, computed for a concentration parameter of $c = 8$) and the disc (blue curves). The disc contribution depends on $M_{\text{d}}/M_{\text{vir}}$ and has been shown for two values: $M_{\text{d}}/M_{\text{vir}} = 0.04$ (blue solid curve) and $M_{\text{d}}/M_{\text{vir}} = 0.004$ (dashed blue curve). The black solid curve and the black dashed curve are obtained by summing the red dot-dashed curve in quadrature with the blue solid curve and the blue dashed curve, respectively. Notice that this figure is similar but not identical to the rotation curve ($v_{\text{rot}}/v_{\text{vir}}$ as a function of $r/r_{\text{vir}}$) because the factor of 1.11 in equation (37) is specific to the optical radius $r_{\text{opt}} = 3.2r_{\text{d}}$.</figcaption></figure>

in normalization. Any conclusion from this comparison should be tempered by the scatter in the observations. However, it is encouraging that when we fit our default model at $z \sim 1$ with a slope of 3.68 over the $v_{\text{rot}}$ range covered by the KROSS data, we find exactly the same normalization as Tiley et al. (2016), possibly because their procedure to measure $v_{\text{rot}}$ is very similar to ours. They use the radius that contains 80 per cent of the light. We use the radius that contains 83 per cent of the mass.

## 3.8 The Faber–Jackson relation

The Faber–Jackson relation (hereafter FJR; Faber & Jackson 1976) links the luminosity (in our case, the stellar mass) of an elliptical galaxy to its stellar velocity dispersion. The FJR is to elliptical galaxies what the TFR is to spirals, although the scatter of the FJR relation is larger than that of the TFR. Our analysis and discussion

<figure><img src="image_22.png" /><figcaption>**Figure 18.** The FJR for galaxies with $B/T > 0.7$ in GallICS 2.0 and in the data points of Cappellari et al. (2013, red circles). $M_{\text{stars}}$ is the total stellar mass and $\sigma_e$ is the stellar velocity dispersion within an aperture $R_e$. The predictions of the default model, model 1, model 2 and model 3 correspond to the region of the diagram between the pair of black solid curves, red curves, black dashed curves and green curves, respectively. For each model, the curves correspond to upper and lower quartiles for $M_{\text{stars}}$ in bins of $\sigma_e$. The arrows correspond to a displacement of 0.15 dex. They are based on the default model and they show the maximum increase in $\sigma_e$ that can be reasonably attributed to dissipation in gas-rich mergers.</figcaption></figure>

of elliptical galaxies will not be as in depth as for spiral galaxies because elliptical galaxies are not the focus of this paper and because potential errors in the modelling of elliptical galaxies have no repercussions on the properties of spirals (the converse is not true). We nevertheless include this paragraph on the FJR for completeness and as a sanity check to demonstrate that our modelling of elliptical galaxies is reasonable.

Fig. 18 compares the FJR predicted by the four models in Table 1 with the data points of Cappellari et al. (2013). The structure coefficient $C$ (equation 24) used to pass from the mass and the effective $R_e$ of an elliptical galaxy to the velocity dispersion $\sigma_e$ within an aperture $R_e$ has been calibrated on the data of Cappellari et al. (2006, 2013). They find $C = 2.5$. The structure coefficient predicted by the Hernquist model assuming isotropic velocity dispersion returns values of $\sigma_e$ that are systematically lower than those found using $C = 2.5$ by a factor of 0.87 (Courteau et al. 2014, chapter 5, appendix B). Therefore, what our comparison with the FJR really probes is the dependence of $R_e$ on $M_{\text{stars}}$.

GallICS 2.0 computes bulge radii from energy conservation in major mergers (equation 22). This assumption gives results in reasonable agreement with the FJR (Fig. 18). The agreement is better for the default model, which allows the accretion of gas on to the central galaxies of haloes with masses up to $M_{\text{shutdown}} = 10^{12.7} M_{\odot}$, than it is for models 1 and 3, where gas accretion shuts down at a lower halo mass.

However, Fig. 18 shows that our SAM tends to overestimate the lower envelope of the FJR for $M_{\text{stars}} < 10^{11} M_{\odot}$. This tendency has a simple physical explanation. Gas can radiate its internal energy, invalidating the assumption of energy conservation. Hydrodynamic simulations by Covington et al. (2008) show that, in gas-rich mergers, dissipation can cause the radii of bulges to contract by up to a factor of 2. At constant mass, shrinking radii translate into higher velocity dispersions.

Dissipation will be more significant for ellipticals with $M_{\text{stars}} \lesssim 10^{11} M_{\odot}$, the cuspy profiles of which are consistent with a
dissipational origin in gas-rich ('wet') mergers. In contrast, ellipticals with $M_{\text{stars}} \gg 10^{11} M_{\odot}$ have cores that are consistent with their formation through dissipationless ('dry') mergers (Kormendy et al. 2009; Bernardi et al. 2011). Cattaneo et al. (2011) confirmed this picture by showing that the transition from a dissipative regime to a dissipationless one at $M_{\text{stars}} \sim 10^{11} M_{\odot}$ is a direct consequence of the shutdown of gas accretion above a critical halo mass. The arrows in Fig. 18 show the maximum increase in $\sigma_e$ that dissipation could plausibly cause in bulges with $M_{\text{stars}} \lesssim 10^{11} M_{\odot}$. They correspond to a contraction in radius by a factor of 2 (the maximum value allowed by Covington et al. 2008). Fig. 18 shows that dissipation could easily explain any systematic difference between the default model and the data points.

Incidentally, we note that the data points of Cappellari et al. (2013) overlap with those for the stellar TFR relation if we plot $M_{\text{stars}}$ on the y-axis and $\sqrt{3}\sigma_e$ on the x-axis. The overlap of the FJR with the TFR after correcting $\sigma_e$ by a factor of $\sqrt{2}$ or $\sqrt{3}$ was first remarked by Kassin et al. (2007).

## 4 CONCLUSION

GalICS 2.0 is a new semi-analytic code that we run on merger trees from an N-body simulation in a *Planck* cosmology. The simulation is used to follow the evolution of DM haloes in mass, position, angular momentum and concentration. In the version used for this paper, GalICS 2.0 assumes a one-to-one correspondence between galaxies and haloes/subhaloes. A beta version that includes the possibility of a delay between halo and galaxy mergers shows no difference with respect to our conclusions.

The masses of luminous galaxies within haloes are determined by (i) the rate at which gas flows to the centre and accretes on to galaxies, (ii) the rate at which the accreted gas is converted into stars and (iii) the rate at which gas is blown out of galaxies.

In GalICS 2.0, galaxies grow through accretion of cold flows (Kereš et al. 2005; Dekel & Birnboim 2006; Dekel et al. 2009). Shock-heated gas never cools and ejected gas is never re-accreted. Hence, the shock-heated fraction determines the accretion rate (reionization heating is included but negligible on the scales resolved by our N-body simulation). In the Horizon/Mare Nostrum cosmological hydrodynamic simulation (Ocvirk et al. 2008), shock heating begins above $M_{\text{vir}} = 10^{10.7} M_{\odot}$ and is complete when $M_{\text{vir}} = 10^{12.7} M_{\odot}$. This is the default assumption to compute the accretion rate on to galaxies in GalICS 2.0. A sharper transition from cold accretion to shock heating (model 1) improves the fit to the SMFs by Yang et al. (2009) and Bernardi et al. (2013) but not to other data sets (Fig. 3, $z \simeq 0.1$). Assuming no cooling and no re-accretion is admittedly extreme but we deliberately intended to explore how far it can take us. The effects of reintroducing these processes will be explored in a future publication, where we shall also discuss which observations can constrain their importance.

The assumption SFR = $\epsilon_{\text{sf}}M_{\text{gas}}/t_{\text{dyn}}$ is standard in SAMs and has been calibrated on observational data (Boselli et al. 2014, Fig. 2). In contrast, outflow rates are considerably uncertain. We have computed them assuming that a fraction $\epsilon_{\text{SN}}$ of the power output from SNe is used to heat gas to $T_{\text{vir}}$ and/or to blow it out of haloes. This fraction depends on the depth of gravitational potential well (measured by $v_{\text{vir}}$) with an exponent tuned to fit the local SMF.

In agreement with previous studies (Cole et al. 1994; Somerville et al. 2008; Benson & Bower 2010; Guo et al. 2011; Benson 2012, but see Henriques et al. 2013), we find that $\eta_{\text{SN}}$ must depend strongly on $v_{\text{vir}}$ to reproduce the shallow slope of the local SMF at $10^9 M_{\odot} < M_{\text{stars}} < 10^{11} M_{\odot}$.

Below $M_{\text{stars}} \sim 10^9 M_{\odot}$, the slope of the observed SMF steepens again (Baldry et al. 2008, 2012). This is unlikely for a reionization signature because 33 out 38 dwarf galaxies in the Local Group formed the bulk of their stars after cosmic reionization (Weisz et al. 2014). In contrast, GalICS 2.0 shows that this feature arises naturally when we require that $\epsilon_{\text{SN}}$ cannot exceed a maximum value. The efficiency of SN feedback increases as $\epsilon_{\text{SN}} \propto v_{\text{vir}}^{-4}$ to keep the slope of the SMF shallow until it reaches a saturation value $\epsilon_{\text{max}}$, after which the slope starts rising again. For our default choice $\epsilon_{\text{max}} = 1$, the saturation scale is lower than the N-body resolution. However, lowering the maximum feedback efficiency<fnref n="3" /> to $\epsilon_{\text{max}} = 0.12$ (model 2) changes the slope at $M_{\text{stars}} \lesssim 10^{9.5} M_{\odot}$, which becomes steeper (Fig. 3). This is beneficial for the agreement with the local SMF but worsens the problem at high $z$, where even the default model overestimates the number density of low-mass galaxies. This is a general problem of SAMs (Guo et al. 2011; Henriques et al. 2012; Asquith et al., in preparation, but see Henriques et al. 2013). In contrast, GalICS 2.0 has no difficulty with the number density of massive galaxies that are already in place at $z \sim 2-2.5$, which used to be a major challenge for SAMs until a few years ago.

If we look at the SFR density integrated over all galaxy masses, then, despite overpredicting the number densities of low-mass galaxies, the default model is consistent with the observations within the errors at $2 < z < 8$ (Fig. 10). In contrast, model 2 is clearly above the data points because the low-mass galaxy excess is too large this time. Overall, our analysis favours strong feedback at high $z$ (the default model is the best in relation to both the SMF and the cosmic SFR density) and a lower efficiency at low $z$ to explain why the slope of the SMF rises again below $M_{\text{stars}} \sim 10^9 M_{\odot}$.

The main problem of the default model with respect to the cosmic SFR density is that it does not seem to be capable of reproducing its decline by a factor of $\sim 20$ from $z \sim 2$ to $z \sim 0$. This is another general problem of SAMs (see Knebe et al. 2015 for an extensive comparison of all major codes). In GalICS 2.0, the local SFR density excess (a factor of $\sim 2$ with respect to the observations) is due to an excess in the number density of galaxies with SFR > $30 M_{\odot} \text{yr}^{-1}$ (Fig. 9). This excess is also visible on the SFR-$M_{\text{stars}}$ diagram (Fig. 8). We can curb it by lowering the shutdown mass, above which shock heating is complete, as in model 3, a parameter combination with $M_{\text{shutdown}} = 10^{12.3} M_{\odot}$ and $\epsilon_{\text{max}} = 0.12$ designed to fit the local SMF by Baldry et al. (2008, 2012). This choice is in better agreement with the number density of galaxies with SFR > $30 M_{\odot} \text{yr}^{-1}$, though it underpredicts the number density of galaxies with SFR = $10-20 M_{\odot} \text{yr}^{-1}$, which the default model reproduces correctly (Fig. 9). Model 3 also reproduces correct bimodality scale at which the transition from the main sequence of star-forming galaxies to the quiescent population occurs (Fig. 8) and the evolution of the cosmic SFR density in the redshift range $0 < z < 2$. In fact, we are not aware of another SAM that reproduces so well the decline of the cosmic SFR density between $z \sim 2$ and $z \sim 0$.

The assumption that the masses and growth histories of DM haloes determine the properties of galaxies is at the heart and foundation of SAMs. Halo quantities determine the accretion rate on to galaxies, the dynamical time, which sets the star formation

<footnote n="3">³ Our efficiency parameter is likely to underestimate the real feedback efficiency required by our model (Section 2.4.2) but the difference is only a factor of $\sim 1.5$ if the ejected gas is heated to the virial temperature and mixed with the hot atmosphere.</footnote>
time-scale, and the feedback efficiency. Hence, direct probes of the galaxy-halo connection are particularly significant. GalICS 2.0 is in good agreement with the halo mass estimates from weak lensing for late-type galaxies (Reyes et al. 2012) and satellite kinematics for early-type galaxies (Wojtak & Mamon 2013, Fig. 5). We also find an unprecedented fit to the baryonic and stellar TFR at $z \sim 0$ (Fig. 14), as well as the correct evolution when passing from $z \sim 0$ to $z \sim 1$ (Fig. 16; recent data by Tiley et al. 2016).

Reproducing the SMF and the TFR simultaneously has been a major challenge for SAMs since their inception (Kauffmann et al. 1993; Heyl et al. 1995). This difficulty has a simple explanation. The characteristic knee of the galaxy SMF around $M_{\text{stars}} \sim 10^{10.5} M_{\odot}$ arises because the slope of the $M_{\text{stars}}-M_{\text{vir}}$ relation changes from approximately $M_{\text{stars}} \propto M_{\text{vir}}^2$ to approximately $M_{\text{stars}} \propto M_{\text{vir}}^{1/2}$ around $M_{\text{vir}} \sim 7 \times 10^{11} M_{\odot}$ (Fig. 5). Therefore, assuming $v_{\text{rot}} \propto v_{\text{vir}} \propto M_{\text{vir}}^{1/3}$ will necessarily predict a bend in the TFR if a model is in agreement with the SMF.

GalICS 2.0 is not affected by this problem because we do not assume $v_{\text{rot}} \propto v_{\text{vir}}$. We compute the disc rotation speed from the disc self-gravity, the gravity of the central bulge and the gravitational potential of the DM halo (we fit NFW profiles to the haloes identified in the N-body simulation). Interestingly, Cole et al. (2000) followed an approach similar to ours (they computed the rotation speeds of spiral galaxies at the disc half-mass radius taking the gravity of the disc, the bulge and the halo properly into account), and they did find a much better agreement with the shape of the TFR than previous studies. In those days, however, the offset of the zero-point received more attention. For the same I-band magnitude, the disc rotation speeds predicted by Cole et al. (2000) are systematically higher than those of the Sd-Sc galaxies of de Jong & Lacey (2000) by ~30 per cent on average.

Cole et al. (2000)'s model for adiabatic contraction (based on Blumenthal et al. 1986) is partially responsible for this offset because it overestimates the phenomenon (Gnedin et al. 2004; Abadi et al. 2010) and thus the DM mass with the disc. Besides the limitations of the model used to describe adiabatic contraction, it is possible that the phenomenon itself may be largely compensated by adiabatic expansion during massive outflows, the extent of which is not easily quantified.

The models presented in this paper neglect adiabatic contraction but we have made the experiment to run our default model with adiabatic contraction modelled as in Blumenthal et al. (1986). The only noticeable differences were in the TFR (Fig. 14), which shifted to higher rotation speeds but not by a large amount. The differences that adiabatic contraction makes are larger at high masses, where $M_{\text{gal}}/M_{\text{vir}}$ is larger and thus we expect a stronger effect of the baryons on the central density of the DM in the model with adiabatic contraction.

In recent versions (e.g. Guo et al. 2011), the Munich model, too, has relaxed the assumption $v_{\text{rot}} \simeq v_{\text{vir}}$ by assuming that $v_{\text{rot}}$ is the maximum circular velocity of the DM halo. However, even with this improvement, the agreement with the shape of the TFR is not as good as the one in Fig. 14 because the model does not consider the dependence of $v_{\text{rot}}/v_{\text{vir}}$ on $M_{\text{gal}}/M_{\text{vir}}$, and this is the essential point to reproduce the galaxy SMF and the TFR simultaneously.

Fig. 17 shows that $v_{\text{rot}} \simeq v_{\text{vir}}$ at $r_{\text{opt}} = 3.2\lambda v_{\text{vir}}/2 \sim 0.08 v_{\text{vir}}$ for dwarf galaxies with $M_{\text{gal}}/M_{\text{vir}} \sim 10^{-3}-10^{-2}$ (black dashed curve), while $v_{\text{rot}} \simeq 1.5 v_{\text{vir}}$ at the same radius for galaxies with masses comparable to the Milky Way, in which $M_{\text{gal}}/M_{\text{vir}} \sim 0.04$ (Fig. 5). The relation between $v_{\text{rot}}$ and $v_{\text{vir}}$ that the dependence of $v_{\text{rot}}/v_{\text{vir}}$ on $M_{\text{gal}}/M_{\text{vir}}$ produces (Fig. 19) is in agreement with previous findings from rotation-curve studies (Cattaneo et al. 2014) and compensates

<figure><img src="image_23.png" /><figcaption>**Figure 19.** The local disc rotation speed-virial velocity relation for model galaxies (point cloud). Disc rotation speeds are given at the optical radius ($3.2r_d$). The black diagonal line is $v_{\text{disc}} = v_{\text{vir}}$. The red dashed curve is the relation that Cattaneo, Salucci & Papastergis (2014) infer from rotation-curve studies. The plot excludes galaxies with $B/T > 0.3$ (the bulge-to-disc ratio $B/T$ is computed considering both the bulge and the pseudo-bulge mass).</figcaption></figure>

the relation between $M_{\text{stars}}$ and $v_{\text{vir}}$, producing a nearly straight TFR.

Fig. 14 shows that the compensation is not perfect because the stellar TFR deviates from a single power law at $M_{\text{stars}} < 10^9 M_{\odot}$. In the default model and model 1, which correspond to a shallow low-mass end of the galaxy SMF (as in Bernardi et al. 2013), the predicted TFR is slightly concave below $v_{\text{rot}} \sim 80 \text{ km s}^{-1}$. In models 2 and 3, which correspond to the steep low-mass slope of Baldry et al. (2008, 2012), the predicted TFR is slightly convex below $v_{\text{rot}} \sim 80 \text{ km s}^{-1}$. Three caveats must be taken into account when discussing these differences. First, observational errors on $M_{\text{stars}}$ and $v_{\text{rot}}$ may blur these differences. When one compares theoretical predictions to observations at $v_{\text{rot}} < 80 \text{ km s}^{-1}$, one must keep in mind that the data for dwarf galaxies are all based on H I measurements. The H I disc may extend well beyond the optical radius at which we measure $v_{\text{rot}}$ in GalICS 2.0. This concern is more significant for dwarf galaxies because of their rising rotation curves. Last but not least, the turnovers at low masses in model SMFs (Fig. 3) show that resolution effects can propagate up to $M_{\text{stars}} \sim 10^9 M_{\odot}$ (discussion in Sections 3.1 and 3.2). Hence, any result at $M_{\text{stars}} \lesssim 19^9 M_{\odot}$ (below the horizontal dash-dotted line in Fig. 14) should be taken with caution to avoid the risk of overinterpretation.

Fig. 19 differs from the relation between $v_{\text{rot}}$ and $v_{\text{vir}}$ that Sales et al. (2017) measure in the **EAGLE** and **APOSTLE** cosmological hydrodynamic simulations. They find $v_{\text{rot}} \simeq 1.15 v_{\text{vir}}$ over a broad range of virial velocities ($v_{\text{vir}} = 30-200 \text{ km s}^{-1}$) possibly because their $M_{\text{gal}}/M_{\text{vir}}$ ratios do not vary with $M_{\text{vir}}$ as much as ours (a natural consequence of less extreme feedback). The proportionality of $v_{\text{rot}}$ with $v_{\text{vir}}$ (fig. 5 of Sales et al. 2017) is the reason why the **EAGLE** SMF (Schaye et al. 2015, simulation Ref-L100N1504 at $z=0.1$) is steeper than the observations at low $M_{\text{stars}}$ and shallower than the observations at high $M_{\text{stars}}$ (Fig. 3). Although the **EAGLE** simulation was calibrated on the SMF and its agreement with the data of Baldry et al. (2012) is quite reasonable considering that a hydrodynamic simulation does not have the same freedom as a SAM, this is the typical behaviour of models calibrated on the TFR.

The optical radius $r_{\text{opt}} = 3.2r_d$ at which we measure $v_{\text{rot}}$ in GalICS 2.0 depends on scalelengths computed from angular momentum conservation. This assumption may be inaccurate on a
galaxy-by-galaxy basis because hydrodynamic simulations have shown that angular momentum is not conserved for individual fluid elements (Sharma & Steinmetz 2005; Kimm et al. 2011). However, the assumption is consistent with the observed spin distribution of disc galaxies both at $z \sim 0$ (Tonini et al. 2006) and $z \sim 1-3$ (Burkert et al. 2016), and the predictions of GalICS 2.0 are in good agreement with disc sizes over the entire redshift range $0 < z < 2.5$, though we tend to see more scatter in GalICS 2.0 than in the observations (Fig. 13). In any case, rotation curves around the optical radius are usually nearly flat, so the dependence of $v_{\text{rot}}$ on $r_{\text{opt}}$ is weak (Fig. 17) and any error on disc sizes has limited impact on our predictions for the TFR.

Fig. 5 implies a formal resolution of $M_{\text{stars}} \lesssim 10^8 M_{\odot}$ but we know from experience (e.g. Cattaneo et al. 2011) that resolution effects could trickle up to scales an order of magnitude larger. Therefore, it is possible that the SMF at $M_{\text{stars}} < 10^9 M_{\odot}$ could have contained more galaxies if we had used merger trees from an N-body simulation with $1024^3$ rather than $512^3$ particles. As the difference would come from the number of low-mass haloes, we expect it to have limited impact on the TFR. We intend to make a resolution study and to present it in a future publication but we are confident that our results are robust at $M_{\text{stars}} > 10^9 M_{\odot}$.

The FJR is the equivalent of the TFR for elliptical galaxies. The focus of this paper is on spiral galaxies but we added a section on the FJR for completeness. In the same way that the radii of discs are computed by assuming that the infall of gas conserves angular momentum, the radii of bulges are computed by assuming that mergers conserve the total mechanical energy of the merging galaxies. Accounting for systematic departures of observed galaxies from the Hernquist model (Cappellari et al. 2006; Courteau et al. 2014), this assumption is in reasonably good agreement with observations of the FJR (Cappellari et al. 2013), although it tends to underestimate the velocity dispersions of lower mass ellipticals. This tendency can be explained as an effect of dissipation, which plays an important role in major mergers at $M_{\text{stars}} < 10^{11} M_{\odot}$ but not in major mergers at larger masses (Kormendy et al. 2009; Bernardi et al. 2011; Cattaneo et al. 2011).

## ACKNOWLEDGEMENTS

We thank P. S. Behroozi, A. Boselli, J. D. Bradford, J. Freundlich, M. Huertas-Company, O. Ilbert, E. Papastergis, P. Salucci and A. L. Tiley for providing electronic data that have assisted us in the comparison with observations. We also thank C. Pichon for useful conversation.

## REFERENCES

* Abadi M. G., Moore B., Bower R. G., 1999, MNRAS, 308, 947

* Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435

* Athanassoula E., 2008, MNRAS, 390, L69

* Baldry I. K., Glazebrook K., Driver S. P., 2008, MNRAS, 388, 945

* Baldry I. K. et al., 2012, MNRAS, 421, 621

* Balogh M. L., Morris S. L., 2000, MNRAS, 318, 703

* Balogh M. L., Navarro J. F., Morris S. L., 2000, ApJ, 540, 113

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Benson A. J., 2012, New Astron., 17, 175

* Benson A. J., Babul A., 2009, MNRAS, 397, 1302

* Benson A. J., Bower R., 2010, MNRAS, 405, 1573

* Bernardi M., Roche N., Shankar F., Sheth R. K., 2011, MNRAS, 412, L6

* Bernardi M., Meert A., Sheth R. K., Vikram V., Huertas-Company M., Mei S., Shankar F., 2013, MNRAS, 436, 697

* Bernardi M., Meert A., Vikram V., Huertas-Company M., Mei S., Shankar F., Sheth R. K., 2014, MNRAS, 443, 874

* Bett P., Eke V., Frenk C. S., Jenkins A., Helly J., Navarro J., 2007, MNRAS, 376, 215

* Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

* Bildfell C., Hoekstra H., Babul A., Mahdavi A., 2008, MNRAS, 389, 1637

* Birrer S., Lilly S., Amara A., Paranjape A., Refregier A., 2014, ApJ, 793, 12

* Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27

* Boselli A., Cortese L., Boquien M., Boissier S., Catinella B., Lagos C., Saintonge A., 2014, A&A, 564, A66

* Bouché N. et al., 2010, ApJ, 718, 1001

* Bouwens R. J. et al., 2012, ApJ, 754, 83

* Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

* Bradford J. D., Geha M. C., van den Bosch F. C., 2016, ApJ, 832, 11

* Bruce V. A. et al., 2014, MNRAS, 444, 1660

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240

* Burkert A. et al., 2016, ApJ, 826, 214

* Cappellari M. et al., 2006, MNRAS, 366, 1126

* Cappellari M. et al., 2013, MNRAS, 432, 1709

* Cattaneo A., 2001, MNRAS, 324, 128

* Cattaneo A., 2002, MNRAS, 333, 353

* Cattaneo A., Combes F., Colombi S., Bertin E., Melchior A.-L., 2005a, MNRAS, 359, 1237

* Cattaneo A., Blaizot J., Devriendt J., Guiderdoni B., 2005b, MNRAS, 364, 407

* Cattaneo A., Dekel A., Devriendt J., Guiderdoni B., Blaizot J., 2006, MNRAS, 370, 1651

* Cattaneo A., Dekel A., Faber S. M., Guiderdoni B., 2008, MNRAS, 389, 567

* Cattaneo A. et al., 2009, Nature, 460, 213

* Cattaneo A., Mamon G. A., Warnick K., Knebe A., 2011, A&A, 533, A5

* Cattaneo A., Woo J., Dekel A., Faber S. M., 2013, MNRAS, 430, 686

* Cattaneo A., Salucci P., Papastergis E., 2014, ApJ, 783, 66

* Chabrier G., 2003, ApJ, 586, L133

* Christodoulou D. M., Shlosman I., Tohline J. E., 1995, ApJ, 443, 551

* Cole S., Aragon-Salamanca A., Frenk C. S., Navarro J. F., Zepf S. E., 1994, MNRAS, 271, 781

* Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168

* Conzelice C. J., 2006, MNRAS, 373, 1389

* Coupon J. et al., 2015, MNRAS, 449, 1352

* Courteau S. et al., 2014, Rev. Mod. Phys., 86, 47

* Cousin M., Lagache G., Bethermin M., Guiderdoni B., 2015, A&A, 575, A33

* Covington M., Dekel A., Cox T. J., Jonsson P., Primack J. R., 2008, MNRAS, 384, 94

* Croton D. J. et al., 2006, MNRAS, 365, 11

* Cucciati O. et al., 2012, A&A, 539, A31

* Danovich M., Dekel A., Hahn O., Ceverino D., Primack J., 2015, MNRAS, 449, 2087

* de Jong R. S., Lacey C., 2000, ApJ, 545, 781

* De Lucia G., Blaizot J., 2007, MNRAS, 375, 2

* Dekel A., Birnboim Y., 2006, MNRAS, 368, 2

* Dekel A., Mandelker N., 2014, MNRAS, 444, 2071

* Dekel A., Silk J., 1986, ApJ, 303, 39

* Dekel A. et al., 2009, Nature, 457, 451

* Di Matteo P., Bournaud F., Martig M., Combes F., Melchior A.-L., Semelin B., 2008, A&A, 492, 31

* Dunne L. et al., 2009, MNRAS, 394, 3

* Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

* Efstathiou G., Lake G., Negroponte J., 1982, MNRAS, 199, 1069

* Elbaz D. et al., 2007, A&A, 468, 33

* Faber S. M., Jackson R. E., 1976, ApJ, 204, 668

* Fall S. M., Efstathiou G., 1980, MNRAS, 193, 189
Fanidakis N., Baugh C. M., Benson A. J., Bower R. G., Cole S., Done C., Frenk C. S., 2011, MNRAS, 410, 53

Flores R. A., Primack J. R., 1996, ApJ, 457, L5

Font A. S. et al., 2008, MNRAS, 389, 1619

Fontanot F., De Lucia G., Monaco P., Somerville R. S., Santini P., 2009, MNRAS, 397, 1776

Fontanot F., Macciò A. V., Hirschmann M., De Lucia G., Kannan R., Somerville R. S., Wilman D., 2015, MNRAS, 451, 2968

Forbes J. C., Krumholz M. R., Goldbaum N. J., Dekel A., 2016, Nature, 535, 523

Freeman K. C., 1970, ApJ, 160, 811

Gargiulo I. D. et al., 2015, MNRAS, 446, 3820

Garnett D. R., 2002, ApJ, 581, 1019

Gatto A. et al., 2017, MNRAS, 466, 1903

Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16

Gonzalez-Perez V., Lacey C. G., Baugh C. M., Lagos C. D. P., Helly J., Campbell D. J. R., Mitchell P. D., 2014, MNRAS, 439, 264

Gruppioni C. et al., 2015, MNRAS, 451, 3419

Gunn J. E., Gott J. R., III, 1972, ApJ, 176, 1

Guo Q., White S., Li C., Boylan-Kolchin M., 2010, MNRAS, 404, 1111

Guo Q. et al., 2011, MNRAS, 413, 101

Gurovich S., Freeman K., Jerjen H., Staveley-Smith L., Puerari I., 2010, AJ, 140, 663

Hall M., Courteau S., Dutton A. A., McDonald M., Zhu Y., 2012, MNRAS, 425, 2741

Hatton S., Devriendt J. E. G., Ninin S., Bouchet F. R., Guiderdoni B., Vibert D., 2003, MNRAS, 343, 75

Henriques B. M. B., White S. D. M., Lemson G., Thomas P. A., Guo Q., Marleau G.-D., Overzier R. A., 2012, MNRAS, 421, 2904

Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R. E., Guo Q., Lemson G., Springel V., 2013, MNRAS, 431, 3373

Hernquist L., 1990, ApJ, 356, 359

Heyl J. S., Hernquist L., Spergel D. N., 1995, ApJ, 448, 64

Hirschmann M., De Lucia G., Fontanot F., 2016, MNRAS, 461, 1760

Hopkins P. F., Cox T. J., Hernquist L., Narayanan D., Hayward C. C., Murray N., 2013, MNRAS, 430, 1901

Hopkins P. F., Kereš D., Ōńorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Häring N., Rix H., 2004, ApJ, 604, L89

Ilbert O. et al., 2013, A&A, 556, A55

Jiang C. Y., Jing Y. P., Faltenbacher A., Lin W. P., Li C., 2008, ApJ, 675, 1095

Kajisawa M., Ichikawa T., Yamada T., Uchimoto Y. K., Yoshikawa T., Akiyama M., Onodera M., 2010, ApJ, 723, 129

Karim A. et al., 2011, ApJ, 730, 61

Kassin S. A. et al., 2007, ApJ, 660, L35

Kauffmann G., White S. D. M., Guiderdoni B., 1993, MNRAS, 264, 201

Kauffmann G., Colberg J. M., Diaferio A., White S. D. M., 1999, MNRAS, 303, 188

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2

Kimm T., Devriendt J., Slyz A., Pichon C., Kassin S. A., Dubois Y., 2011, preprint (arXiv:1106.0538)

Knebe A. et al., 2015, MNRAS, 451, 4029

Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216

Kraljic K., Renaud F., Bournaud F., Combes F., Elmegreen B., Emsellem E., Teyssier R., 2014, ApJ, 784, 112

Lacey C., Cole S., 1993, MNRAS, 262, 627

Lacey C. G. et al., 2016, MNRAS, 462, 3854

Larson R. B., 1974, MNRAS, 169, 229

Laurikainen E., Salo H., Buta R., Knapen J. H., Comerón S., 2010, MNRAS, 405, 1089

Le Borgne D., Elbaz D., Ocvirk P., Pichon C., 2009, A&A, 504, 727

Lee J., Yi S. K., 2013, ApJ, 766, 38

Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119

McCracken H. J. et al., 2015, MNRAS, 449, 901

McGaugh S. S., 2012, AJ, 143, 40

Magorrian J. et al., 1998, AJ, 115, 2285

Marconi A., Hunt L. K., 2003, ApJ, 589, L21

Mathews W. G., Baker J. C., 1971, ApJ, 170, 241

Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319

Monaco P., Fontanot F., Taffoni G., 2007, MNRAS, 375, 1189

Moore B., 1994, Nature, 370, 629

Moustakas J. et al., 2013, ApJ, 767, 50

Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584

Muzzin A. et al., 2013, ApJ, 777, 18

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Neto A. F. et al., 2007, MNRAS, 381, 1450

Noordermeer E., van der Hulst J. M., Sancisi R., Swaters R. A., van Albada T. S., 2005, A&A, 442, 137

Ocvirk P., Pichon C., Teyssier R., 2008, MNRAS, 390, 1326

Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920

Onions J. et al., 2013, MNRAS, 429, 2739

Papastergis E., Cattaneo A., Huang S., Giovanelli R., Haynes M. P., 2012, ApJ, 759, 138

Papastergis E., Adams E. A. K., van der Hulst J. M., 2016, A&A, 593, A39

Paturel G., Theureau G., Bottinelli L., Gouguenheim L., Coudreau-Durand N., Hallet N., Petit C., 2003, A&A, 412, 57

Peirani S., Jung I., Silk J., Pichon C., 2012, MNRAS, 427, 2625

Peng Y., Maiolino R., Cochrane R., 2015, Nature, 521, 192

Peng Y.-j., Maiolino R., 2014, MNRAS, 443, 3643

Persic M., Salucci P., Stel F., 1996, MNRAS, 281, 27

Pizagno J. et al., 2005, ApJ, 633, 844

Planck Collaboration XX, 2014, A&A, 571, A20

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Renaud F., Kraljic K., Bournaud F., 2012, ApJ, 760, L16

Reyes R., Mandelbaum R., Gunn J. E., Pizagno J., Lackner C. N., 2011, MNRAS, 417, 2347

Reyes R., Mandelbaum R., Gunn J. E., Nakajima R., Seljak U., Hirata C. M., 2012, MNRAS, 425, 2610

Rhee G., Valenzuela O., Klypin A., Holtzman J., Moorthy B., 2004, ApJ, 617, 1059

Robotham A. S. G., Driver S. P., 2011, MNRAS, 413, 2570

Romano D., Karakas A. I., Tosi M., Matteucci F., 2010, A&A, 522, A32

Rujopakarn W. et al., 2010, ApJ, 718, 1171

Sales L. V. et al., 2017, MNRAS, 464, 2419

Salim S. et al., 2007, ApJS, 173, 267

Schaye J. et al., 2015, MNRAS, 446, 521

Sérsic J. L., 1963, Bol. Asociacion Argentina Astron. La Plata Argentina, 6, 41

Shankar F., Marulli F., Bernardi M., Mei S., Meert A., Vikram V., 2013, MNRAS, 428, 109

Shankar F. et al., 2014, MNRAS, 439, 3189

Sharma S., Steinmetz M., 2005, ApJ, 628, 21

Shen S., Mo H. J., White S. D. M., Blanton M. R., Kauffmann G., Voges W., Brinkmann J., Csabai I., 2003, MNRAS, 343, 978

Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1

Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15

Shim H., Colbert J., Teplitz H., Henry A., Malkan M., McCarthy P., Yan L., 2009, ApJ, 696, 785

Silk J., 2003, MNRAS, 343, 249

Skibba R. A. et al., 2015, ApJ, 807, 152

Smolčić V. et al., 2009, ApJ, 690, 610

Sobral D., Best P. N., Matsuda Y., Smail I., Geach J. E., Cirasuolo M., 2012, MNRAS, 420, 1926

Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087

Somerville R. S., Hopkins P. F., Cox T. J., Robertson B. E., Hernquist L., 2008, MNRAS, 391, 481

Somerville R. S., Gilmore R. C., Primack J. R., Domínguez A., 2012, MNRAS, 423, 1992

Swaters R. A., Balcells M., 2002, A&A, 390, 863

Swaters R. A., Verheijen M. A. W., Bershady M. A., Andersen D. R., 2003, ApJ, 587, L19
Tadaki K.-I., Kodama T., Koyama Y., Hayashi M., Tanaka I., Tokoku C., 2011, Publ. Astron. Soc. Japan, 63, 437

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Tiley A. L. et al., 2016, MNRAS, 460, 103

Tollet E. et al., 2016, MNRAS, 456, 3542

Tomczak A. R. et al., 2014, ApJ, 783, 85

Tonini C., Lapi A., Shankar F., Salucci P., 2006, ApJ, 638, L13

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

van den Bosch F. C., 1998, ApJ, 507, 601

van der Burg R. F. J., Hildebrandt H., Erben T., 2010, A&A, 523, A74

van der Wel A. et al., 2014, ApJ, 788, 28

Vincenzo F., Matteucci F., Belfiore F., Maiolino R., 2016, MNRAS, 455, 4183

Weinzirl T., Jogee S., Khochfar S., Burkert A., Kormendy J., 2009, ApJ, 696, 411

Weisz D. R., Dolphin A. E., Skillman E. D., Holtzman J., Gilbert K. M., Dalcanton J. J., Williams B. F., 2014, ApJ, 789, 148

Wen Z. L., Han J. L., Liu F. S., 2010, MNRAS, 407, 533

White S. D. M., Frenk C. S., 1991, ApJ, 379, 52

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Wilman D. J., Fontanot F., De Lucia G., Erwin P., Monaco P., 2013, MNRAS, 433, 2986

Wojtak R., Mamon G. A., 2013, MNRAS, 428, 2407

Wuyts S. et al., 2011, ApJ, 738, 106

Yang X., Mo H. J., van den Bosch F. C., 2009, ApJ, 693, 830

Yegorova I. A., Salucci P., 2007, MNRAS, 377, 507

Yoshida M. et al., 2006, ApJ, 653, 988

Zhang W., Li C., Kauffmann G., Zou H., Catinella B., Shen S., Guo Q., Chang R., 2009, MNRAS, 397, 1243

Zheng X. Z., Bell E. F., Papovich C., Wolf C., Meisenheimer K., Rix H.-W., Rieke G. H., Somerville R., 2007, ApJ, 661, L41

This paper has been typeset from a TeX/LATEX file prepared by the author.