The new semianalytic code GalICS 2.0 - Reproducing the galaxy stellar mass function and the Tully-Fisher relation simultaneously

TL;DR

GalICS 2.0 is a semianalytic galaxy formation model that successfully reproduces key galaxy properties, including the stellar mass function and Tully-Fisher relation, by incorporating detailed halo and baryonic physics.

Contribution

The paper introduces GalICS 2.0, a novel semianalytic code that models galaxy evolution with improved rotation speed calculations, enabling simultaneous reproduction of multiple galaxy observables.

Findings

Reproduces galaxy stellar mass function and Tully-Fisher relation.

Aligns with halo mass measurements from weak lensing and satellite kinematics.

Matches observed gas fractions, star formation rates, and galaxy structural properties.

Abstract

GalICS 2.0 is a new semianalytic 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 onto 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 semianalytic 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 the rotation speed is proportional to the virial velocity. 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.

Tables

Table 1. Table 1: Model parameters: symbols, units (for dimensional quantities), default values and alternative models (where they differ).

Parameter	Symbol	Units	Default	Model1	Model 2	Model 3
Cosmology
Matter density	${\Omega }_M$		$0.308$
Baryon density	${\Omega }_b$		$0.0481$
Cosmological constant	${\Omega }_{\Lambda }$		$0.692$
Hubble constant	$H_0$	$100\mathrm{km}{s}^{-1}{\mathrm{Mpc}}^{-1}$	$0.678$
N-body simulation
Box size	$L_{\mathrm{box}}$	Mpc	$100$
Resolution	$N_{\mathrm{part}}$		${512}^3$
Dimensional parameters
Thermal velocity dispersion of IGM	${\sigma }_{\mathrm{reio}}$	km s-1	$25$
Star formation threshold	${\Sigma }_{\mathrm{th}}$	$M_{\odot }{\mathrm{pc}}^{-2}$	$2$
Minimum shock heating mass	$M_{\mathrm{shock}}$	$M_{\odot }$	${10}^{10.7}$	${10}^{11.3}$		${10}^{11.0}$
Shutdown mass	$M_{\mathrm{shutdown}}$	$M_{\odot }$	${10}^{12.7}$	${10}^{12.4}$		${10}^{12.3}$
SN feedback saturation scale	$v_{\mathrm{SN}}$	$\mathrm{km}{s}^{-1}$	$24$			$43$
SN energy	$E_{\mathrm{SN}}$	erg	${10}^{51}$
SN rate	${\Psi }_{\mathrm{SN}}$	$M_{\odot }^{-1}$	$1/140$
Adimensional efficiency factors
Star formation	${ϵ}_{\mathrm{sf}}$		$0.04$
Disc instabilities	${ϵ}_{\mathrm{inst}}$		$0.9$
Mass ratio for major mergers	${ϵ}_m$		$4$			$3$
Structure coefficient of bulges	$C$		$2.5$
Black hole accretion	${ϵ}_{\bullet }$		$0.01$
Maximum SN feedback efficiency	${ϵ}_{\max }$		$1$		$0.12$	$0.12$
Returned fraction	$R$		$0.45$
Metal yield	$y$		$0.06$
SN feedback scaling exponents
$v_{\mathrm{vir}}$ scaling	${\alpha }_v$		$-4$			$-6.2$
$z$ scaling	${\alpha }_z$		$3$

Equations

$$M_{\mathrm{dm}}(x)=\frac{\log (1+cx)-cx/(1+cx)}{\log (1+c)-c/(1+c)}\cdot \frac{{\Omega }_M-{\Omega }_b}{{\Omega }_M}{M}_{\mathrm{vir}},$$

$$t_{s-1}<t_{\mathrm{m}}<t_s.$$

$$\Delta M_b=\max [f_b{M}_{\mathrm{vir}}(t_s)-M_b(t_m),0],$$

$$f_b=\frac{{\Omega }_b}{{\Omega }_M}{\int }_0^{v_{\mathrm{vir}}\sqrt{2}}4\pi {\left(\frac{\mu }{2\pi k{T}_{\mathrm{reio}}}\right)}^{\frac{3}{2}}{v}^2{e}^{-\frac{\mu v^2}{2k{T}_{\mathrm{reio}}}}\mathrm{d}v,$$

$$f_b=\frac{{\Omega }_b}{{\Omega }_M}\left[\mathrm{erf}\left(\frac{v_{\mathrm{vir}}}{{\sigma }_{\mathrm{reio}}}\right)-\frac{2}{\sqrt{\pi }}\left(\frac{v_{\mathrm{vir}}}{{\sigma }_{\mathrm{reio}}}\right)e^{-{\left(\frac{v_{\mathrm{vir}}}{{\sigma }_{\mathrm{reio}}}\right)}^2}\right].$$

$$f_{\mathrm{hot}}=\frac{\mathrm{Log}{M}_{\mathrm{vir}}-\mathrm{Log}{M}_{\mathrm{shock}}}{\mathrm{Log}{M}_{\mathrm{shutdown}}-\mathrm{Log}{M}_{\mathrm{shock}}}$$

$$f_{\mathrm{hot}}=\frac{{\dot{M}}_{\mathrm{hot}}}{{\dot{M}}_{\mathrm{cold}}+{\dot{M}}_{\mathrm{hot}}}$$

$$\Delta M_{\mathrm{hot}}=f_{\mathrm{hot}}\Delta M_b,$$

$$\Delta M_{\mathrm{fil}}=(1-f_{\mathrm{hot}})\Delta M_b-\frac{M_{\mathrm{fil}}}{t_{\mathrm{ff}}}\Delta t,$$

$$M_b=M_{\mathrm{hot}}+M_{\mathrm{fil}}+M_{\mathrm{gal}}+M_{\mathrm{out}},$$

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

$$M_{\mathrm{bulge}}(r)={\left(\frac{r/r_{\mathrm{bulge}}}{1+r/r_{\mathrm{bulge}}}\right)}^2{M}_{\mathrm{bulge}},$$

$$J_d=2\pi {\Sigma }_0{\int }_0^{\infty }{e}^{-r/r_d}{r}^2{v}_c\mathrm{d}r=M_d{r}_d{\int }_0^{\infty }{e}^{-x}{x}^2{v}_c\mathrm{d}x,$$

$$\frac{J_d}{M_d}=\frac{J_{\mathrm{vir}}}{M_{\mathrm{vir}}}.$$

$$r_d=\frac{\lambda r_{\mathrm{vir}}}{{\int }_0^{\infty }{e}^{-x}{x}^2\frac{v_c}{v_{\mathrm{vir}}}\mathrm{d}x},$$

$$\lambda \equiv \frac{J_h}{M_h{r}_{\mathrm{vir}}{v}_{\mathrm{vir}}}$$

$$v_c^2(r)=\frac{G{M}_{\mathrm{dm}}(r)}{r}+\frac{G{M}_{\mathrm{bulge}}(r)}{r}+v_d^2(r),$$

$$r[M_{\mathrm{NFW}}(r_i)+M_d(r)+M_{\mathrm{bulge}}(r)]=r_i\frac{M_{\mathrm{NFW}}(r_i)}{1-{\Omega }_b/{\Omega }_M},$$

$$v_d>{ϵ}_{\mathrm{inst}}{v}_c,$$

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

$$\frac{1}{2}{U}_{\mathrm{bulge}}=\frac{1}{2}{U}_1+\frac{1}{2}{U}_2+E_{12},$$

$$\frac{1}{2}{C}_{\mathrm{remn}}\frac{\mathrm{G}({\mathcal{M}}_1+{\mathcal{M}}_2{)}^2}{r_{1/2}}=$$

$$=\frac{1}{2}{C}_1\frac{\mathrm{G}{\mathcal{M}}_1^2}{r_1}+\frac{1}{2}{C}_2\frac{\mathrm{G}{\mathcal{M}}_2^2}{r_2}+\frac{1}{2}{C}_{12}\frac{\mathrm{G}{\mathcal{M}}_1{\mathcal{M}}_2}{r_1+r_2}.$$

$$r_{1/2}=(1+\sqrt{2})r_{\mathrm{bulge}}.$$

$${\sigma }_e^2=\mathrm{G}\frac{{\mathcal{M}}_1+{\mathcal{M}}_2}{C{R}_e},$$

$$M_{\bullet }=M_{\bullet 1}+M_{\bullet 2}+{ϵ}_{\bullet }(M_{gas1}+M_{gas2}),$$

$${\dot{M}}_{\mathrm{gas}}={\dot{M}}_{\mathrm{accr}}+{\dot{M}}_{\mathrm{sml}}-\mathrm{SFR}-{\dot{M}}_{\mathrm{out}}.$$

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

$${\dot{M}}_{\mathrm{stars}}=\frac{{\dot{M}}_{\mathrm{gas}}}{t_{\mathrm{sf}}}-{\dot{M}}_{\mathrm{sml}}.$$

$$\mathrm{SFR}={ϵ}_{\mathrm{sf}}\frac{M_{\mathrm{gas}}}{t_{\mathrm{dyn}}}.$$

Full text

††pagerange: The new semianalytic code GalICS 2.0

• Reproducing the galaxy stellar mass function and the Tully-Fisher relation simultaneously–The new semianalytic code GalICS 2.0
• Reproducing the galaxy stellar mass function and the Tully-Fisher relation simultaneously††pubyear: 2016

1 Introduction

Semianalytic 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 article 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 to 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 et al. (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 re-written from scratch. One of the reasons is to enable a more extensive use of the cosmological N-body simulation used to contruct 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) to explain why the galaxy stellar mass function (SMF) has a knee at $M_{\rm stars}\sim 10^{11}\,M_{\odot}$ when the mass function of their host haloes is essentially a single power law up to $M_{\rm vir}\sim 10^{13}\,M_{\odot}$. Yet the assumption that the growth histories of DM haloes determine the properties of galaxies underpins the entire semianalytic 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 et al., 1993) or they are calibrated on the SMF/luminosity function and fail to reproduce the TFR (Heyl et al., 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 semianalytic 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 article, 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 article 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 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 summarises the conclusions of the article.

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 onto 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 recipies. 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 pseudobulge), 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 article, 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{\rm\,km\,s}^{-1}{\rm Mpc}^{-1}$ (Planck Collaboration et al., 2014, Planck + WP + BAO). The simulation has a volume of $(100{\rm\,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 dark matter haloes). Outputs have been saved at $265$ timesteps equally spaced in the logarithm of the expansion factor between $z=13.2$ and $z=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 cosmology111$3M_{\rm vir}/(4\pi\rho_{c}r_{\rm vir}^{3})$, where $\rho_{c}$ is the critical density of the Universe, cannot be exactly equal to $\Delta_{c}$ because $M_{\rm 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.. The virial mass $M_{\rm 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_{\rm 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_{\rm vir}=10^{12}\,M_{\odot}$ corresponds to $v_{\rm vir}=\sqrt{{\rm G}M_{\rm vir}/r_{\rm vir}}=128{\rm\,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_{\rm vir}$ below which the mass function drops (the peak is at $M_{\rm vir}=10^{10.625}\,M_{\odot}$ rather than at $M_{\rm vir}\simeq 10^{10.5}\,M_{\odot}$ simply because this is the midpoint of the logarithmic mass bin $10^{10.5}\,M_{\odot}<M_{\rm 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 gray shaded area) with the analytic fit of Sheth et al. (2001). This fit contains two free parameters, which the original article 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 function (gray shaded area) with this fit (black solid curve) shows good agreement up to $M_{\rm vir}\sim 10^{13.3}\,M_{\odot}$. However, at $M_{\rm vir}\sim 10^{13.8}-10^{14}\,M_{\odot}$, the centre of the gray shaded area is above the black solid curve by about $\sim 40\%$. Number densities for the central galaxies of haloes in this mass range may overestimated by a similar factor.

For each halo, we measure the virial angular momentum ${\bf J}_{\rm vir}$, which we use to compute the spin parameter $\lambda$ (Section 2.3.2). Bett et al. (2007) argued that they needed at least three hundred particles to measure halo spins accurately but, going from three hundred to one hundred particles, as we do here, the median spin parameter changes from $\lambda=0.0425$ to $\lambda=0.044$ (Fig. 7 of Bett et al., 2007). The difference ($<4\%$) 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 et al., 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_{\rm vir}<10^{12}M_{\odot}$ (corresponding to $\sim 3000$ particles), suggesting that our measurements are affected by N-body resolution in a manner that may be significant for haloes with $M_{\rm vir}<10^{11.5}\,M_{\odot}$. Most of the figures shown in this article 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_{\rm 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 article, 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 concentrations 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_{\rm 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 peeled off 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 descendent 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 descendent. The main progenitor is always the one with the largest virial mass. A halo/subhalo is found to have no descendent if it looses more than half its mass but no single halo accretes enough mass from it to qualify as its descendent. 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_{\rm vir}$, $r_{\rm vir}$ and ${\bf J}_{\rm 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_{\rm 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

[TABLE]

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

2.1.3 Scheme to evolve baryons along merger trees

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

[TABLE]

If the halo has no progenitors, we assume that $t_{\rm 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_{\rm 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 timescale (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 article.

2.2 halo

The halo module follows the accretion of gas onto 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_{\rm m}$ (with $t_{s-1}<t_{\rm m}<t_{s}$) to its values at $t_{s}$ with the equation:

[TABLE]

where $M_{b}(t_{m})$ is the sum of the baryon masses of all progenitors at the time of merging.

There is no accretion onto subhaloes or their descendents. This prescription models the physical phenomenon of strangulation through ram-pressure and tidal stripping (Gunn & Gott, 1972; Abadi et al., 1999; Balogh et al., 2000; Balogh & Morris, 2000; Peng et al., 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 timesteps 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 Eq. (3), $f_{b}=f_{b}(M_{\rm vir},z)$ is a function that models reionisation feedback (gas will not accrete onto 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_{\rm esc}=v_{\rm vir}\sqrt{2}$, where $v_{\rm vir}(M_{\rm vir},z)$ is the halo virial velocity. This assumption gives:

[TABLE]

where $T_{\rm reio}$ is the temperature at which the intergalactic medium is reionised, $\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 Eq. (4) because the integral gives the fraction of the baryonic mass with $v<v_{\rm esc}$, while $f_{b}$ is expressed in terms of the total mass $M_{\rm vir}$. Let $\sigma_{\rm reio}$ be the one-dimensional thermal velocity dispersion of the intergalactic medium, so that $\mu\sigma_{\rm reio}^{2}=kT_{\rm reio}$. Then, Eq. (4) becomes:

[TABLE]

Eq. (5) implies that $f_{b}=0.5\Omega_{b}/\Omega_{M}$ for $v_{\rm vir}/\sigma_{\rm reio}=1.088$.

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

A fraction $f_{\rm 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 onto the central galaxy through cold filamentary flows. The fraction $f_{\rm hot}(M_{\rm vir})$ is assumed to be $f_{\rm hot}=0$ for $M_{\rm vir}\leq M_{\rm shock}$,

[TABLE]

for $M_{\rm shock}<M_{\rm vir}<M_{\rm shutdown}$, and $f_{\rm hot}=1$ for $M_{\rm vir}\geq M_{\rm shutdown}$. Here $M_{\rm shock}$ is the halo mass at which the accreted gas begins to be shock-heated, while $M_{\rm 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_{\rm shock}$ and $M_{\rm shutdown}$ are in principle free parameters of the model, to be determined by fitting the galaxy SMF, the functional form in Eq. (6) is based on work by Ocvirk et al. (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_{\rm vir}$ in the Horizon-Mare Nostrum cosmological hydrodynamical simulation. The results of Ocvirk et al. (2008) for

[TABLE]

at $0.2r_{\rm vir}$ are very similar to a ramp between ${\rm Log}(M_{\rm shock}/M_{\odot})=10.7$ and ${\rm Log}(M_{\rm 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 SDSS data for a sharper cut-off of the form $M_{\rm shock}=M_{\rm shutdown}=10^{12.3}\,M_{\odot}$.

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

[TABLE]

[TABLE]

where $\Delta M_{b}$ is the accreted mass calculated in Eq. (3). The second term on the right hand side of Eq. (9) is the accretion rate from the filaments onto the galaxy, which we assume to take place on a freefall timescale $t_{\rm ff}=r_{\rm vir}/v_{\rm 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 Eq. (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 an article focused on spiral galaxies, which tend to live in haloes with M_{\rm vir}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}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 reaccretion 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 onto the halo. It is important to realize that, in this article, 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 reaccretion 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 reaccretion will be explored in a future publication. However, we note that, to an extent, Eq. (6) already includes some of the effects of cooling implicitly because the cold gas fraction $1-f_{\rm hot}$ that Ocvirk et al. (2008) measure in their simulation at $0.2r_{\rm vir}$ includes any gas that may have been shock-heated and cooled before reaching $0.2r_{\rm vir}$ (the outer boundary of the Hi disc; C. Pichon, private communication), although it does not account for the possibility that a cooling-flow may develop inside the galaxy itself and for the reaccretion 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 Eq. (3) is

[TABLE]

where $M_{\rm gal}$ is the total baryonic mass of the galaxy. This definition of $M_{b}$ includes the ejected mass $M_{\rm out}$ to prevent its re-accretion. Hence, the actual halo baryon fraction $(M_{\rm hot}+M_{\rm fil}+M_{\rm gal})/M_{\rm vir}$ may be lower than $M_{b}/M_{\rm vir}$. Had we defined $M_{b}$ as $M_{b}=M_{\rm hot}+M_{\rm fil}+M_{\rm gal}$, $M_{\rm out}$ would have been available for immediate reaccretion onto 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, scale lengths, 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 pseudobulges 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 pseudobulge, built through disc instabilities, and a classical bulge, built through major mergers.

The disc has an exponential surface-density profile:

[TABLE]

The pseudobulge originates from the buckling and bar instability of the disc within radius $r_{\rm pseudo}$. We therefore assume that, while instabilities transfer matter from the disc to the pseudobulge, the mass redistribution is mainly vertical and azimuthal, and that the radial exponential profile of the disc plus pseudobulge 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

[TABLE]

where $r_{\rm 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 onto a galaxy is automatically added to the disc component. Its scale radius $r_{d}$ is determined by solving the disc angular momentum equation

[TABLE]

where $v_{c}(r)$ is the the disc rotation curve (see Eq. 17 below) and $M_{d}$ is the total baryonic mass (stars and gas) of the disc and the pseudobulge combined (this is what $M_{d}$ stands for throughout this article).

Following Fall & Efstathiou (1980) and Mo et al. (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

[TABLE]

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 Eq. (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).

Eq. (13) can be re-written as

[TABLE]

where

[TABLE]

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_{\rm vir}$ at all radii, then the integral in Eq. (15) would make $2$ and Eq. (15) would reduce to $r_{d}=\lambda r_{\rm 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:

[TABLE]

where $M_{\rm dm}(r)$ is the DM mass with radius $r$ (Eq. 1), $M_{\rm bulge}(r)$ is the bulge mass within radius $r$ (Eq. 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 Eq. (13) is not trivial. The pseudobulge does not appear as a fourth contribution to $v_{c}(r)$ in Eq. (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_{\rm NFW}(r_{i})$ is the initial DM profile (assumed to be described by the NFW model), then this assumption gives the equation:

[TABLE]

from which $r_{i}(r)$ and thus $M_{\rm dm}(r)=M_{\rm NFW}(r_{i})$ can be computed.

Eq. (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 GalICS 2.0, there is no adiabatic contraction.

A version with adiabatic contraction (computed with Eq. 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 et al. (1982), Christodoulou et al. (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:

[TABLE]

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

The pseudobulge radius $r_{\rm pseudo}$ is the largest radius at which the instability condition (Eq. 19) is satisfied. Its value is used to increment the pseudobulge mass with the algorithm

[TABLE]

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

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

2.3.5 Mergers

A merger of two galaxies is major if ${\cal M}_{1}<\epsilon_{m}{\cal M}_{2}$ with ${\cal M}_{1}>{\cal M}_{2}$, where $\epsilon_{m}$ is a parameter of the model. Here ${\cal M}_{1}$ and ${\cal 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):

[TABLE]

where $U_{\rm bulge}$, $U_{1}$ and $U_{2}$ are the gravitational potential energies of the resulting bulge, galaxy $1$ and galaxy $2$, respectively, while $E_{\rm 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. Eq. (21) can be re-written as:

[TABLE]

[TABLE]

Here $r_{1/2}$ is the merger remnant’s half-mass radius. $C_{\rm 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}\equiv 2E_{12}/U_{12}$ where $U_{12}=-{\rm G}{\cal M}_{1}{\cal M}_{2}/(r_{1}+r_{2})$ is the two-body system gravitational binding energy. The factor of two in the definition of $C_{12}$ comes, once agains, 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 they are all equal. Covington et al. (2008) have shown that Eq. (22) with $C_{\rm remn}=C_{1}=C_{2}=0.5C_{12}$ is in good agreement with the sizes of the remnants of dissipationless mergers in hydrodynamic simulations, though the best fit for $C_{\rm 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 (Eq. 22) by up to a factor of two.

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

[TABLE]

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 ${\cal M}_{1}$ and ${\cal 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 Eq. (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.8153r_{\rm 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:

[TABLE]

where $C$ is a structure coefficient. For a Hernquist profile with isotropic velocity dispersion, $C=3.31$ (Courteau et al., 2014, chapter 5, section 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_{\bullet}$ 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

[TABLE]

where $M_{\bullet 1}$, $M_{\bullet 2}$, $M_{\rm gas1}$, $M_{\rm 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 pseudobulge 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:

[TABLE]

Here $\dot{M}_{\rm accr}$ is the accretion rate onto the galaxy, which is entirely due cold flows since threre is no cooling in our model (Section 2.2.1). $\dot{M}_{\rm 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}_{\rm out}$ is the rate at which gas is blown out by stellar feedback (outflow rate).

The time separation between two output timesteps $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 timescale on which a stellar population evolves, i.e., the timescale on which $\dot{M}_{\rm sml}$ varies (even massive OB stars spend $>6\,$Myr on the main sequence). $\dot{M}_{\rm accr}$ varies on a timescale (set by $t_{\rm dyn}$) that is even longer. Hence, without loss of generality, we can write ${\rm SFR}=M_{\rm gas}/t_{\rm sf}$ and $\dot{M}_{\rm out}=\eta{\rm SFR}$, where our only assumption about the star formation timescale $t_{\rm 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 timescales in response to stellar evolution or changes in the structural properties of the component (computed in galaxy).

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

[TABLE]

[TABLE]

These equations have an analytic solution (Cole et al., 2000; Lilly et al., 2013; Dekel & Mandelker, 2014; Peng & Maiolino, 2014)222The analytic solution has the interesting feature that $M_{\rm gas}\rightarrow\dot{M}_{\rm accr}t_{\rm sf}/(1+\eta)$ for $t\rightarrow\infty$ when $\dot{M}_{\rm sml}=0$., which we use to evolve $M_{\rm gas}$ and $M_{\rm 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 article, 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}_{\rm sml}=RM_{\rm gas}/t_{\rm 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}_{\rm 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 star formation and feedback in the values of $t_{\rm 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 timescale $t_{\rm sf}$. For discs, we assume that $t_{\rm sf}$ is proportional to $t_{\rm dyn}$, where $t_{\rm dyn}$ is the orbital time at the disc scale radius $r_{d}$ ($t_{\rm sf}=t_{\rm dyn}/\epsilon_{\rm sf}$; see, e.g., Kauffmann et al., 1993 and Kennicutt, 1998). Hence:

[TABLE]

The star formation efficiency $\epsilon_{\rm sf}=1/25$ adopted in this article is calibrated on a local galaxy sample by Boselli et al. (2014). With this efficiency, our model reproduces the statistical relation between $M_{\rm gas}$ and SFR in local galaxies (Fig. 2, left) and the SFRs of individual galaxies within a factor of two (Fig. 2, right). The gas masses plotted in Fig. 2 are total masses of cold neutral gas. Had we used H2 masses instead, the observed $M_{\rm gas}$ - SFR relation would have had a higher normalization and woud 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_{\rm gas}$ relation and the SFR - $M_{{\rm H}_{2}}$ relation.

Explaining why $t_{\rm sf}$ is so long compared to the local freefall time of the gas is a hot topic in star formation theory (e.g., Renaud et al., 2012; Hopkins et al., 2014; Kraljic et al., 2014; Forbes et al., 2016; Gatto et al., 2017). By calibrating Eq. (29) on observational data, we are effectively short-circuiting these complex physics, whose outcome is summarised in the value of $\epsilon_{\rm 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 pseudobulges and bulges by simply redefining the dynamical time. In GalICS 2.0, $t_{\rm dyn}$ is the orbital time at $r_{\rm pseudo}$ for pseudobulges and the half-crossing time of the starbursting region $t_{\rm dyn}=r_{\rm 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 pericentic 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 scale-length of its distribution is of order $r_{\rm starburst}\sim 0.1r_{\rm bulge}$ (Cattaneo et al., 2005). 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 timescale $t_{\rm sf}=M_{\rm gas}/\dot{M}_{\rm sf}$ for starburst galaxies is about ten 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_{\rm SFR}$ and gas surface density $\Sigma_{\rm gas}$ for $\Sigma_{\rm gas}<\Sigma_{\rm th}\sim 9\,M_{\odot}{\rm\,pc}^{-2}$ (Bigiel et al., 2008), where $\Sigma_{\rm gas}$ is the mean gas surface density (Hi plus H2) within the optical radius $r_{\rm opt}=3.2r_{d}$. However, there are galaxies on the relation (including some of the black squares in Fig. 2) with values of $\Sigma_{\rm gas}$ as low as $\sim 2\,M_{\odot}{\rm\,pc}^{-2}$. Therefore, the threshold is not sharp. In GalICS 2.0, we set SFR $=0$ for $\Sigma_{\rm gas}<\Sigma_{\rm th}$, where $\Sigma_{\rm th}$ is a parameter of the model. We set it to the relatively low value $\Sigma_{\rm th}=2\,M_{\odot}{\rm\,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 artifact. The surface area on which we spread the gas to compute $\Sigma_{\rm gas}$ is $2\pi(r_{\rm opt}^{2}-r_{\rm pseudo}^{2})$ for discs, $2\pi r_{\rm pseudo}^{2}$ for pseudobulges and $2\pi r_{\rm 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 supernova (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 percent 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 overcooling 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}_{\rm out}/{\rm 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_{\rm sf}$ (Section 2.4.1). Similarly, the fraction $\epsilon_{\rm 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_{\rm 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_{\rm SN}\sim 10^{51}\,$erg is the energy released by one SN, then the power output from SNe will be $E_{\rm SN}\Psi_{\rm SN}{\rm SFR}.$ If a fraction $\epsilon_{\rm SN}$ of this power is used to drive a wind with speed $v_{\rm w}$, then the outflow rate $\dot{M}_{\rm out}$ from the component will satisfy:

[TABLE]

(Silk, 2003). Expulsion from the gravitational potential well of the DM requires $v_{\rm w}\geq v_{\rm esc}\sim 2.5v_{\rm vir}$ (the numerical coefficient in front of $v_{\rm 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 reabsorbe the uncertainty on $v_{\rm w}/v_{\rm vir}$ into the free parameter $\epsilon_{\rm SN}$ and define mass-loading factor $\eta$ so that:

[TABLE]

The only inconvenient of this definition is that $\epsilon_{\rm 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$) is the gas blown out of the galaxy if heated to the virial temperature and mixed with the hot atmosphere (for a singular isothermal sphere, ${3\over 2}kT_{\rm vir}={3\over 2}\cdot{1\over 2}\mu v_{\rm vir}^{2}$; White & Frenk, 1991). Much larger energetic efficiencies (\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}2.5^{2}\epsilon_{\rm SN}) are required if the gas expelled from galaxies is also blown out of the halo.

The problem of this simple scaling with $v_{\rm 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 focussed on massive galaxies). A phenomenological solution is to introduce a SN efficiency $\epsilon_{\rm SN}$ that depends on both $v_{\rm vir}$ and redshift $z$, and to impose a plausible maximum $\epsilon_{\rm max}$ to the values that $\epsilon_{\rm SN}$ can take:

[TABLE]

where $v_{\rm SN}$, $\alpha_{v}$ and $\alpha_{z}$ are free parameters of the model to be determined by fitting the galaxy SMF. The speed $v_{\rm SN}$ corresponds to the virial velocity for which $\epsilon_{\rm 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_{\rm SN}=24{\rm\,km\,s}^{-1}$; Table 1 and Section 3) give a simple relation between mass-loading factor and halo mass:

[TABLE]

since $v_{\rm vir}\propto M_{\rm vir}^{1/3}(1+z)^{1/2}$. $M_{\rm vir}$ is a physical quantity, though it is not clear why the outflow rate should scale with $M_{\rm vir}$ rather than $v_{\rm vir}$. The mass resolution of the N-body simulation used to construct the merger trees is $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$. In this article, we formulate our model in terms of $v_{\rm vir}$ and $z$ rather than $M_{\rm 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.

Eq. (33) corresponds to a very strong dependence of the mass-loading factor on the virial velocity ($\eta\propto v_{\rm 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 [math] and $-5.5$ (Lacey et al., 2016), though the details of how stellar feedback is implemented vary from one model to another (see Hirschmann et al., 2016 for a discussion of the mass-loading in different SAMs and simulations).

Our normalization of $\eta$ at $M_{\rm 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_{\rm vir}=10^{11}\,M_{\odot}$ is only a factor of three larger that our resolution limit, our normalization combined to our much steeper dependence on $v_{\rm 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 article. It is possible that we fit the observations with lower mass-loading factors for a given halo mass because our current model neglects the reaccretion of ejected gas.

Physically, $\eta$ is limited by the maximum energetic efficiency of supernovae $\epsilon_{\rm max}$. Without such maximum, Eq. (32) implies $\epsilon_{\rm SN}\rightarrow\infty$ for $v_{\rm vir}\rightarrow 0$, which is absurd (the wind cannot contain more energy than it is available). In the most generous case, $\epsilon_{\rm max}=1$. The real efficiency will probably be much lower. In practice, $\epsilon_{\rm SN}$ is limited by the mass resolution of the N-body simulation, $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$. Inserted into Eq. (33), this mass gives a maximum mass-loading factor of $\eta=30-40$. As $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$ corresponds to $v_{\rm vir}\sim 40{\rm\,km\,s}^{-1}$ at $z=0$, our default parameter values (Table 1) imply \epsilon_{\rm SN}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.1 at $z\sim 0$ for all haloes that we can resolve. At high $z$, however, $\epsilon_{\rm 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}_{\rm sml}$ and the mass loss rate in metals $\dot{M}_{\rm sml,\,Z}$ are

[TABLE]

and

[TABLE]

In Eq. (35), $1-R$ is the fraction of the star-forming gas that remains in stars and contributes to the final stellar masses 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 future developments but has no effect whatsoever on any of the results presented in this article 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 article. 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 sixteen 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 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_{\rm SN}$, $\Psi_{\rm 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 et al. (2008) find that cosmic reionization suppresses gas accretion onto haloes up to $M_{\rm reio}=6.5\times 10^{9}h_{0}^{-1}M_{\odot}$ at $z=0$. This mass corresponds to $v_{\rm vir}\sim 25{\rm\,km\,s}^{-1}$. We therefore assume $\sigma_{\rm reio}=25{\rm\,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_{\rm shock}\sim 10^{10.7}\,M_{\odot}$ and $M_{\rm shutdown}\sim 10^{12.7}\,M_{\odot}$.

The two parameters $\epsilon_{\rm m}$ and $\epsilon_{\rm 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_{\rm m}=3-5$. We choose $\epsilon_{\rm m}=4$ as our default value. For disc instabilities, Efstathiou et al. (1982) used N-body simulations and showed that $\epsilon_{\rm inst}\simeq 0.91$. Christodoulou et al. (1995) argued for a lower value ($\epsilon_{\rm inst}\simeq 0.83$) in stellar discs but remarked that the presence of gas could raise $\epsilon_{\rm inst}$. The values reported above are the inverse of those contained in the original articles because of our different definition of $\epsilon_{\rm inst}$. Furthermore, the original articles 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 Eq. (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_{\rm 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; Cappellari et al., 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 onto the central black hole in major mergers to $\epsilon_{\bullet}=0.01$ because we know from experience with different codes (Cattaneo, 2001; Cattaneo et al., 2005) 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_{\rm SN}$, $\alpha_{v}$, $\alpha_{z}$, $\epsilon_{\rm 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_{\rm 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_{\rm SN}=24{\rm\,km\,s}^{-1}$, $\alpha_{v}=-4$ and $\alpha_{z}=3$ (Section 3.1).

In this article, 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_{\rm shock}$ has been increased to $M_{\rm shock}\sim 10^{11.3}\,M_{\odot}$ and $M_{\rm shutdown}$ has been lowered to $M_{\rm shutdown}\sim 10^{12.4}\,M_{\odot}$). Model 2 limits the efficiency of SN feedback to $\epsilon_{\rm max}=0.12$ (models with $\epsilon_{\rm max}>0.2$ are indistinguishable from the default model because there are few galaxies with an energetic efficiency of SNe $>20\%$). Model 3 produces a SMF that is flat around $M_{\rm stars}\sim 10^{10}\,M_{\odot}$ and rises steeply at lower masses in agreement with observations by Baldry et al. (2008); 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_{\rm stars}$ vs. $M_{\rm vir}$, gas-to-stellar mass ratio vs. $M_{\rm stars}$, SFR vs. $M_{\rm 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 article, 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{\rm\,km\,s}^{-1}{\rm Mpc}^{-1}$ and a Chabrier (2003) inital mass function.

3.1 Mass functions

Fig. 3 compare 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_{\rm 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 et al. (2009) in the mass range $10^{9}\,M_{\odot}<M_{\rm 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_{\rm stars}>10^{11.8}\,M_{\odot}$, we overpredict galaxy number densities by a factor of two 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\%$ at $M_{\rm 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_{\rm shock}=10^{11.3}\,M_{\odot}$ and $M_{\rm 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_{\rm 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_{\rm 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_{\rm 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_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{9}\,M_{\odot}) but makes things worse at higher $z$.

The SMF of Baldry et al. (2008); Baldry et al. (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, that of Baldry et al. is steeper at $M_{\rm stars}<10^{9.5}\,M_{\odot}$, almost flat at $10^{9.5}\,M_{\odot}<M_{\rm stars}<10^{10.5}$ 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_{\rm vir}$ stronger (we pass from $\alpha_{v}=-4$ to $\alpha_{v}=-6.2$) so that the SMF becomes almost flat at intermediate masses but then we cap the efficiency of SN feedback at $\epsilon_{\rm max}=0.12$ to increase the low-mass slope, like in model 2. We also raise $M_{\rm shock}$ and lower $M_{\rm shutdown}$ (a bit like in model 1) to make the change of slope around $M_{\rm 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_{\rm 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_{\rm 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_{\rm 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_{\rm 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 datasets are of comparable magnitude.

The main discreapancy 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_{\rm 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 reaccretion timescale that is inversely proportional to halo mass. In GalICS 2.0, there is no reaccretion 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 reaccretion timescale required to make this picture work ($6\,$Gyr time for a halo with $M_{\rm vir}\sim 3\times 10^{10}\,M_{\odot}$) is much longer than the radiative cooling timescale. 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 Hi (ALFALFA) data were available. The baryonic mass was assumed to be $M_{\rm gal}=M_{\rm stars}+1.4M_{\rm HI}$, where $M_{\rm stars}$ is the stellar mass derived from optical data, $M_{\rm HI}$ is the Hi mass from radio data and the factor of $1.4$ accounts for the presence of helium. This is a lower limit for $M_{\rm gal}$ because gas could be present and not be detected. One can find an upper limit by giving to galaxies not detected in Hi the maximum Hi mass consistent with their non-detection. The gray 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_{\rm 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_{\rm gal}>3\times 10^{11}\,M_{\odot}$ but neither do they fit the SMF of Baldry et al. at $M_{\rm 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_{\rm gal}\sim 10^{10}\,M_{\odot}$. At $10^{9}\,M_{\odot}<M_{\rm gal}<10^{10}\,M_{\odot}$, baryonic masses are underestimatee by $0.2\,$dex on average. $M_{\rm 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_{\rm gal}\sim 10^{9}\,M_{\odot}$ but is underpredicted around $M_{\rm 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_{\rm stars}$ - $M_{\rm 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 (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_{\rm stars}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{10}\,M_{\odot}. On the other hand, Fig. 5 shows that these data are not very powerful to discriminate between models.

The dotted-dashed diagonal line in Fig. 5 corresponds to $M_{\rm stars}=(\Omega_{b}/\Omega_{M})M_{\rm vir}$. It has been plotted to emphasize that the $M_{\rm stars}$ - $M_{\rm vir}$ relation is steeper than linear at $M_{\rm vir}<10^{12}\,M_{\odot}$ and shallower at higher masses. The logarithmic slope of the relation varies from $\sim 2$ at low mass to $\sim 1/2$ at high mass. This change is the reason why we observe a knee in the galaxy stellar mass function (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_{\rm hot}$ is the gas fraction that is able to accrete onto galaxies (red curve). Both are shown as function of $M_{\rm vir}$. The gas that makes stars is the one that is able to accrete and to avoid ejection. Its mass fraction, $(1-f_{\rm hot})/(1+\eta)$, is shown by the black curve. The black curve peaks for $M_{\rm vir}\simeq 10^{11.6}\,M_{\odot}$ in the default model and $M_{\rm 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 et al. (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_{\rm gal}$ - $M_{\rm 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_{\rm stars}$ and $M_{\rm gal}$ that corresponds to our halo-mass resolution ($M_{\rm 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_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}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_{\rm 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}\,M_{\odot}$ 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_{\rm HI}+M_{{\rm H}_{2}})/M_{\rm 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 Hi data only (Swaters & Balcells, 2002; Garnett, 2002; Noordermeer et al., 2005; Zhang et al., 2009; also shown in Fig. 7), but the differences are small.

The Hi 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 Hi 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 Hi 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 Hi 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_{\rm gas}=10^{7.5}\,M_{\odot}$ for a galaxy with $M_{\rm stars}=10^{8}\,M_{\odot}$ and $M_{\rm gas}=10^{7.5}\,M_{\odot}$ for a galaxy $M_{\rm stars}=10^{11},M_{\odot}$ (A. 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 et al. (2016). Fig. 4 suggests that, at $M_{\rm gal}\sim 3\times 10^{9}\,M_{\odot}$, the default model underestimates $M_{\rm gal}$ by $\sim 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); Baldry et al. (2012) and Yang et al. (2009). TFR informs us on gas fractions because there are data sets for which we have both $M_{\rm stars}$ and $M_{\rm 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_{\rm 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 is correctly reproduced by the standard model but appears titlted in model 3. In contrast, model 3 reproduces correctly the characteristic mass $M_{\rm 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_{\rm 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_{\rm 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$ (Sobral et al., 2012; Tadaki et al., 2011; Shim et al., 2009) radio ($1.4\,$GHz; Karim et al., 2011; Dunne et al., 2009; Smolčić et al., 2009), combined ultraviolet/infrared (Gruppioni et al., 2015; Cucciati et al., 2012; Kajisawa et al., 2010; Salim et al., 2007; Zheng et al., 2007), ultraviolet data only (Bouwens et al., 2012; Robotham & Driver, 2011; van der Burg et al., 2010; Yoshida et al., 2006), and far infrared data only (Rujopakarn et al., 2010; Le Borgne et al., 2009). 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_{\rm max}=0.12$. Therefore, even though the efficiency scales as $\epsilon_{\rm 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_{\rm 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 bulge-to-total mass ratio that separates S- and E-type galaxies.

GalICS 2.0 computes $B/T$ ratios assuming that the bulge mass is the total stellar mass of the classical bulge and the pseudobulge. The total mass is the sum of the stellar masses of the disc, the classical bulge and the pseudobulge. The fraction of galaxies with $B/T>0.7$ increases with $M_{\rm stars}$, first more gently at $M_{\rm stars}<10^{11}\,M_{\odot}$, then more rapidly at $M_{\rm 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_{\rm stars}>1-3\times 10^{11}\,M_{\odot}$ (the dotted-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 to the data points in Fig. 11 for $\epsilon_{\rm m}=4$ and $\epsilon_{\rm 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_{\rm m}=3$. These differences come from the gas mass that accretes onto 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_{\rm 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 to a 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$ ratios is highly desirable. Nevertheless, the only purpose of morphologies in this article 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 ($\pm 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 bulge-to-total mass ratios 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 (M. 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.

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_{\rm 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 dotted-dashed curves in Fig. 13. Some scatter is still present because galaxies differ in halo concentration, $M_{\rm gal}/M_{\rm 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 log-normal distribution with $\bar{\lambda}=0.049$ and $\sigma_{{\rm ln\,}\lambda}=0.57$ (Fig. 12), in agreement with Muñoz-Cuartas et al., 2011, who find $\bar{\lambda}=0.044$ and $\sigma_{{\rm ln\,}\lambda}=0.57$, and Burkert et al. (2016), who find $\bar{\lambda}=0.052$ and $\sigma_{{\rm 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 dotted-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 pseudobulge 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 (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 pseudobulges 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_{\rm stars}$ or the baryon mass $M_{\rm gal}$ (i.e., the total mass of stars and cold neutral gas) of a spiral galaxy to its rotation speed $v_{\rm rot}$. In this article, $v_{\rm rot}$ is the circular velocity $v_{c}$ at the optical radius $r_{\rm opt}=3.2r_{\rm d}$, which contains $83\%$ of the mass of an exponential disc. We choose this definition for consistency with previous work (Cattaneo et al., 2014) and because measuring $v_{\rm 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 et al., 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_{\rm vir}\propto M_{\rm vir}^{1/3}$. Hence, the $M_{\rm stars}$ - $v_{\rm vir}$ relation is entirely determined by the $M_{\rm stars}$ - $M_{\rm vir}$ relation, which changes slope at $M_{\rm 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_{\rm rot}$ depends on $v_{\rm 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:

[TABLE]

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

[TABLE]

For $\lambda=0.05$ and $c=8$ (the mean values of the spin and concentration parameters), Eq. (15) gives values of $x$ that range from $x=0.083$ for a dwarf galaxy ($M_{d}/M_{\rm vir}=0.004$) to $x=0.063$ for a Milky Way ($M_{d}/M_{\rm vir}=0.04$). For comparison, the standard model $r_{d}=\lambda r_{\rm vir}/2$ gives $x=0.08$. The range of typical values extends to 0.03\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}x\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}0.16 when the scatter in $\lambda$ is considered. In this range, the first term on the right hand side of Eq. (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 curves in Fig. 17, respectively). The dependence on $M_{d}/M_{\rm vir}$ is the reason why GalICS 2.0 reproduces the $v_{\rm rot}$ - $v_{\rm 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_{\rm rot}/v_{\rm vir}$ on $M_{\rm gal}/M_{\rm 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 Eq. (37) decreases with $x$ and, summed to first one, conspires to a produce a rotation curve that is nearly flat at the optical radius (dashed and solid black curves in Fig. 14). Hence, $v_{\rm rot}/v_{\rm vir}$ is not very sensitive to $x$.

Having established that the dependence of $v_{\rm rot}/v_{\rm vir}$ on $M_{\rm gal}/M_{\rm 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 $M_{\rm 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_{\rm stars}<10^{10}\,M_{\odot}$ and the stellar TFR has a small concavity below $M_{\rm stars}\sim 10^{9}\,M_{\odot}$. In model 3, the SMF has a steeper faint-end slope in better agreement with the observations ($M_{\rm stars}/M_{\rm vir}$ decreases less rapidly at low masses) and the concavity is absent.

The baryonic TFR becomes slightly concave below $M_{\rm gal}\sim 10^{10}\,M_{\odot}$ but the concavity turns to a convexity below $M_{\rm gal}\sim 10^{9}\,M_{\odot}$ because, below this mass, the gas surface density is lower than the threshold $\Sigma_{\rm 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_{\rm 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_{\rm stars}>10^{9.5}\,M_{\odot}$. Below this mass, the higher concentrations from the fitting formulae of Dutton & Macciò (2014) cause galaxies with $M_{\rm stars}>10^{8.5}\,M_{\odot}$ to turn with rototation speeds $\sim 15\%$ higher. Hence, the concavity of the stellar TFR in default model becomes more pronounced.

All the models considered until now assume no adiabatic contraction (our standard assumption throughout this article). However, the role of adiabatic contraction on 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 modeled 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_{\rm gal}>10^{10}\,M_{\odot}$. At lower masses, $M_{\rm gal}/M_{\rm 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_{\rm 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 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 Hi 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; Hi linewidths), Yegorova & Salucci (2007; resolved rotation curves from H$\alpha$ and Hi data), Pizagno et al. (2005; H$\alpha$ linewidths) and Rhee et al. (2005; H$\alpha$ linewidths) find a normalization that is higher by $\sim 0.3\,$dex in mass (the systematic shift in $M_{\rm 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_{\rm gal}$ dynamically by performing a halo/disc decomposition (the other authors measured $M_{\rm 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_{\rm rot}$ at exactly $r_{\rm opt}=3.2r_{d}$. The H$\alpha$ measurements by Reyes et al. (2011) are intermediate but closer to the results of other H$\alpha$ studies. Therefore, for a same $M_{\rm stars}$ or $M_{\rm gal}$, the rotation speeds from Hi studies tend to be systematically higher than those from H$\alpha$ studies.

This finding has a simple explanation. The H$\alpha$ emission probes star formation, which is mainly concentrated in the inner regions of galaxies ($r<r_{\rm opt}$), while the Hi emission extends beyond the limit $r_{\rm opt}$ of the stellar disc. The rotation curves of discs with M_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{10.5}\,M_{\odot} are usually rising (Persic et al., 1996). Hence, the gas seen in Hi has higher rotation speeds than the gas seen in H$\alpha$. This also explains why the stellar TFR in GalICS 2.0 is intermediate between the Hi determination by Bradford et al. (2016) and the H$\alpha$ determination by Pizagno et al. (2005), though all our models are consistent with the H$\alpha$ 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 measuremts of $M_{\rm gas}/M_{\rm 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, Hi data overestimate the rotation speed at $r_{\rm 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_{\rm rot}$ is equal to the circular velocity $v_{c}$ required to support a particle on a circular orbit against gravity. However, $v_{\rm 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$ (gray 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_{\rm stars}/M_{\rm 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_{\rm vir}$ with $z$ for a fixed $M_{\rm vir}$.

The stellar TFR at $z\sim 1$ can be compared to H$\alpha$ data from the KMOSS 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 $3.68$, which they then used to fit their $z\sim 1$ data and measure the change 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 $3.68$ over the $v_{\rm 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_{\rm rot}$ is very similar to ours. They use the radius that contains $80\%$ of the light. We use the radius that contains $83\%$ 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 of elliptical galaxies will not be as in depth as for spiral galaxies because elliptical galaxies are not the focus of this article 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. 19 compares the FJR predicted by the four models in Table 1 with the data points of Cappellari et al. (2013). The structure coefficient $C$ (Eq. 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); Cappellari et al. (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, section B). Therefore, what our comparison with the FJR really probes is the dependence of $R_{e}$ on $M_{\rm stars}$.

GalICS 2.0 computes bulge radii from energy conservation in major mergers (Eq. 22). This assumption gives results in reasonable agreement with the FJR (Fig. 19). The agreement is better for the default model, which allows the accretion of gas onto the central galaxies of haloes with masses up to $M_{\rm 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. 19 shows that our SAM tends to overestimate the lower envelope of the FJR for $M_{\rm 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 two. At constant mass, shrinking radii translate into higher velocity dispersions.

Dissipation will be more significant for ellipticals with M_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}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_{\rm 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_{\rm 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. 19 show the maximum increase in $\sigma_{e}$ that dissipation could plausibly cause in bulges with M_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{11}\,M_{\odot}. They correspond to a contraction in radius by a factor of two (the maximum value allowed by Covington et al., 2008). Fig. 19 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_{\rm 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 semianalytic 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 article, 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 onto 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 reaccreted. Hence, the shock-heated fraction determines the accretion rate (reionization heating is included but neglible 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_{\rm vir}=10^{10.7}\,M_{\odot}$ and is complete when $M_{\rm vir}=10^{12.7}\,M_{\odot}$. This is the default assumption to compute the accretion rate onto 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 reaccretion 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 ${\rm SFR}=\epsilon_{\rm sf}M_{\rm gas}/t_{\rm 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_{\rm SN}$ of the power output from SNe is used to heat gas to $T_{\rm vir}$ and/or to blow it out of haloes. This fraction depends on the depth of gravitational potential well (measured by $v_{\rm 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_{\rm SN}$ must depend strongly on $v_{\rm vir}$ to reproduce the shallow slope of the local SMF at $10^{9}\,M_{\odot}<M_{\rm stars}<10^{11}\,M_{\odot}$.

Below $M_{\rm stars}\sim 10^{9}\,M_{\odot}$, the slope of the observed SMF steepens again (Baldry et al., 2008; Baldry et al., 2012). This is unlikely to 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_{\rm SN}$ cannot exceed a maximum value. The efficiency of SN feedback increases as $\epsilon_{\rm SN}\propto v_{\rm vir}^{-4}$ to keep the slope of the SMF shallow until it reaches a saturation value $\epsilon_{\rm max}$, after which it the slope starts rising again. For our default choice $\epsilon_{\rm max}=1$, the saturation scale is lower than the N-body resolution. However, lowering the maximum feedback efficiency to $\epsilon_{\rm max}=0.12$ 333Our 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 $\sim 1.5$ if the ejected gas is heated to the virial temperature and mixed with the hot atmosphere. (model 2) changes the slope at M_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}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_{\rm 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 to reproduce 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 ${\rm SFR}>30\,M_{\odot}{\rm\,yr}^{-1}$ (Fig. 9). This excess is also visible on the SFR - $M_{\rm stars}$ diagramme (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_{\rm shutdown}=10^{12.3}\,M_{\odot}$ and $\epsilon_{\rm max}=0.12$ designed to fit the local SMF by Baldry et al. (2008); Baldry et al. (2012). This choice is in better agreement with the number density of galaxies with ${\rm SFR}>30\,M_{\odot}{\rm\,yr}^{-1}$, though it underpredicts the number density of galaxies with ${\rm SFR}=10-20\,M_{\odot}{\rm\,yr}^{-1}$, which the default model reproduces it 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 onto galaxies, the dynamical time, which sets the star formation timescale, 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_{\rm stars}\sim 10^{10.5}\,M_{\odot}$ arises because the slope of the $M_{\rm stars}$ - $M_{\rm vir}$ relation changes from approximately $M_{\rm stars}\propto M_{\rm vir}^{2}$ to approximately $M_{\rm stars}\propto M_{\rm vir}^{1/2}$ around $M_{\rm vir}\sim 7\times 10^{11}\,M_{\odot}$ (Fig. 5). Therefore, assuming $v_{\rm rot}\propto v_{\rm vir}\propto M_{\rm 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_{\rm rot}\propto v_{\rm 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 identifed 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 a 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 $\sim 30\%$ 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 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 article neglect adiabatic contraction but we have made the experiment to run our default model with adiabatic contraction modeled 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_{\rm gal}/M_{\rm 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_{\rm rot}\simeq v_{\rm vir}$ by assuming that $v_{\rm 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_{\rm rot}/v_{\rm vir}$ on $M_{\rm gal}/M_{\rm vir}$, and this is the essential point to reproduce the galaxy SMF and the TFR simultaneously.

Fig. 17 shows that $v_{\rm rot}\simeq v_{\rm vir}$ at $r_{\rm opt}=3.2\lambda r_{\rm vir}/2\sim 0.08r_{\rm vir}$ for dwarf galaxies with $M_{\rm gal}/M_{\rm vir}\sim 10^{-3}-10^{-2}$ (black dashed curve), while $v_{\rm rot}\simeq 1.5v_{\rm vir}$ at the same radius for galaxies with masses comparable to the Milky Way, in which $M_{\rm gal}/M_{\rm vir}\sim 0.04$ (Fig. 5). The relation between $v_{\rm rot}$ and $v_{\rm vir}$ that the dependence of $v_{\rm rot}/v_{\rm vir}$ on $M_{\rm gal}/M_{\rm vir}$ produces (Fig. 18) is in agreement with previous findings from rotation curve studies (Cattaneo et al., 2014) and compensates the relation between $M_{\rm stars}$ and $v_{\rm 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_{\rm 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_{\rm rot}\sim 80{\rm\,km\,s}^{-1}$. In models 2 and 3, which correspond to the steep low-mass slope of Baldry et al. (2008); Baldry et al. (2012), the predicted TFR is slightly convex below $v_{\rm rot}\sim 80{\rm\,km\,s}^{-1}$. Three caveats must be taken into account when discussing these differences. First, observational errors on $M_{\rm stars}$ and $v_{\rm rot}$ may blur these differences. Second, the comparison with observations at $v_{\rm rot}<80{\rm\,km\,s}^{-1}$ should also have in made that the data for dwarf galaxies are all based on Hi measurements. The Hi disc may extend well beyond the optical radius at which we measure $v_{\rm 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_{\rm stars}\sim 10^{9}\,M_{\odot}$ (discussion in Sections 3.1 and 3.2). Hence, any result at M_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}19^{9}\,M_{\odot} (below the horizontal dashed-dotted line in Fig. 14) should be taken with caution to avoid the risk of overinterpretation.

Fig. 18 differs from the relation between $v_{\rm rot}$ and $v_{\rm vir}$ that Sales et al. (2017) measure in the eagle and apostle cosmological hydrodynamic simulations. They find $v_{\rm rot}\simeq 1.15v_{\rm vir}$ over a broad range of virial velocities ($v_{\rm vir}=30-200{\rm\,km\,s}^{-1}$) possibly because their $M_{\rm gal}/M_{\rm vir}$ ratios do not vary with $M_{\rm vir}$ as much as ours (a natural consequence of less extreme feedback). The proportionality of $v_{\rm rot}$ with $v_{\rm 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_{\rm stars}$ and shallower than the observations at high $M_{\rm 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_{\rm opt}=3.2r_{d}$ at which we measure $v_{\rm rot}$ in GalICS 2.0 depends on scale-lengths 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_{\rm rot}$ on $r_{\rm 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_{\rm stars}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}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_{\rm 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_{\rm stars}>10^{9}\,M_{\odot}$.

The FJR is the equivalent of the TFR for elliptical galaxies. The focus of this article 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_{\rm 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.