© ESO, 2015

1. Introduction

After more than twenty years, (hybrid)-semi analytical models (SAMs) are the best tool for the physical interpretation of large surveys. Under the hypothesis that baryonic processes cannot strongly influence the structuration of the dark-matter, the dark-matter evolution is decoupled from the computation of the baryonic component. Even if this decoupling is a strong assumption (see van Daalen et al. 2011), it allows a wide range of physical prescriptions to be explored in a realistic cosmological context and in a short computational time.

Originally proposed by White & Frenk (1991), SAMs are still being developed (e.g. Cole 1991; Cole et al. 2000; Hatton et al. 2003; Baugh 2006; Croton et al. 2006; Cattaneo et al. 2006; Somerville et al. 2008; Guo et al. 2011; Henriques et al. 2013). After a first version described in a dedicated paper, the updates of these models are often fragmented into a large number of publications. In this context, we found it necessary to clarify all the steps of the modelling and build an up-to-date model based on the latest results obtained by hydrodynamic simulations and the comparison of empirical models with observations. Inspired by the previous GalICS model (Hatton et al. 2003), we revisit all the standard mechanisms acting on the formation and the evolution of galaxies in our new SAM, step by step, and in a single paper.

We mainly focus on the gas cycle, which is the main actor in galaxy stellar mass assembly. It is thus described from the cosmological smooth accretion to the complex feedback mechanisms. Obviously, all physical processes act in the gas cycle with different time scales. Indeed, gas-accretion, star formation, and clump migration time scales are all different. For this reason, we propose an adaptive time-step scheme to solve the gas evolution equations.

Like its predecessor (Hatton et al. 2003), the new model is based on a dark-matter N-body simulation computed in the Λ cold dark-matter paradigm (Λ-CDM). Indeed, even if the growth of dark-matter haloes can be followed using the Press & Schechter (1974) formalism and their extensions (e.g. Cole 1991; Kauffmann et al. 1993; Cole et al. 1994; Somerville & Kolatt 1999), the evolution of the dark-matter component is now commonly extracted from huge dark-matter simulations (e.g. Kauffmann et al. 1999; Helly et al. 2003; Hatton et al. 2003). In such high-resolution simulations, the dark-matter halo population covers a wide range of mass (M_h ∈ [ 10⁷:10¹³ ] M_⊙) and scale, from dwarf haloes to groups. To follow its time evolution, the hierarchical growth of dark-matter structures is commonly represented by merger trees in which each branch represents the growth of a given halo. The analysis of dark-matter simulations have set the current paradigm in which the dark-matter halo growth follows two modes (e.g. Lacey & Cole 1994; Tormen et al. 1997; Wechsler et al. 2002; van den Bosch 2002; Zhao et al. 2009; Fakhouri & Ma 2010; Fakhouri et al. 2010; Genel et al. 2010):

• A smooth accretion that is difficult to define with accuracy. Indeed, it is closely linked to the mass resolution. For example, a 10⁶ M_⊙ dark-matter particle can be seen (or defined) as an unresolved dark-matter halo. However Genel et al. (2010) show that the smooth accretion process converges to non-zero values when resolution increases. More precisely, they demonstrate that the smooth accretion can account for up to 40% of the total mass of a dark-matter halo. These results prove that the naive expectation that all accretion would come through mergers when the resolution increases is wrong. The smooth dark-matter accretion is therefore a major vector of the growth of dark-matter structures. It is subject of a detailed description in this paper.

• The second mode is obviously linked to the merging of pre-existing structures.

Applied to the evolution of the dark-matter component, the baryonic post-treatment is articulated around three new prescriptions:

• A two-phase smooth baryonic accretion, with a cold and a hot mode that are built on the smooth dark-matter accretion.

• A complete monitoring of the hot phase. The evolution of the mean temperature and density profile parameters are followed using explicit energy conservation laws, under the hydrostatic equilibrium hypothesis.

• A new approach for describing of the disc gravitational instabilities. Based on the formation, migration and disruption of giant-clumps, this new prescription adapted for SAMs relies on the work by Dekel et al. (2009b).

Beyond these three new prescriptions, the model examines the standard mechanisms of star formation regulation. For the low-mass structures two processes are in play:

• Supernovae (SNe) feedback. Originally proposed by Larson (1974) and White & Rees (1978), this process is based on the injection of SN energy in the ISM, which generates a wind that expels a fraction of the gas from the disc structure. This mechanism has been widely used in SAM during the past 15 years (e.g. Kauffmann et al. 1993; Cole et al. 1994, 2000; Silk 2003; Benson et al. 2003; Hatton et al. 2003; Croton et al. 2006; Cattaneo et al. 2006; Somerville et al. 2008; Guo et al. 2011; Henriques et al. 2013).

• Photoionisation. Originally proposed by Doroshkevich et al. (1967), photoionisation has been introduced in the context of the CDM paradigm by Couchman & Rees (1986), Ikeuchi (1986), and Rees (1986). The ultraviolet background, created by the first generations of stars, heats the gas surrounding the halo. In the smallest structures (M_h< 10¹⁰ M_⊙), the temperature reached by the gas is high enough to prevent it from collapsing. The gas accretion on galaxies, hence the star formation, hosted by these small structures, are thus reduced. In the context of SAM, the impact of this mechanism has been previously explored by Benson et al. (2002, 2003); Somerville (2002); Croton et al. (2006); and Guo et al. (2011), among others.

The galaxy/dark-matter potential well is too deep for high-mass structures (M_h> 10¹² M_⊙). In this case, SNe and photoionisation processes are not powerful enough to eject some gas or to prevent the infall of baryons. SNe still produce winds but ejecta velocities are not sufficient to expel the gas. It is thus rapidly re-accreted in the disc. A popular interpretation that can explain the low star formation efficiency in M_h> 10¹² M_⊙ haloes is based on the activity of a super massive black hole (SMBH) hosted by the galaxy (e.g. Croton et al. 2006; Cattaneo et al. 2006; Bower et al. 2006; Malbon et al. 2007; Somerville et al. 2008; Ostriker et al. 2010; Bower et al. 2012). If a fraction of the energy produced by the accretion on the SMBH is used to heat the hot gas phase, the efficiency of the cooling process is decreased. The accretion of the gas in the galaxy, and thus the star formation, are then reduced.

With these standard recipes some profound disagreements between model predictions and observations are still present, even if observed local (z = 0) galaxy properties are generally reproduced well. For example, one of the most important problems is the strong excess of low-mass galaxies predicted at high z. Indeed the stellar mass functions measured at 1 <z< 4, for example by Ilbert et al. (2010, 2013), or Caputi et al. (2011), demonstrate that SAMs predict a number of low-mass galaxies roughly ten times larger than what observed. Models cannot also easily explain the cosmic star formation history in detail and the bulk of star formation activity observed at 1 <z< 3. In addition, recent studies show that the link between dark-matter haloes and their host galaxies is not well understood (e.g. Leauthaud et al. 2012). Some empirical models (e.g. Béthermin et al. 2012; Behroozi et al. 2013a,b) suggest a non-monotonically growth of the stellar mass with the dark-matter halo mass. Indeed, the star formation activity derived from observations seems to be very inefficient both in low-mass (M_h< 10¹⁰ M_⊙) and high-mass haloes (M_h> 10¹³ M_⊙). The star formation activity seems to peak for galaxies evolving at intermediate halo mass (M_h ≃ 10¹² M_⊙) (Eke et al. 2004, 2005; Guo et al. 2010).

While the current paper describes all the standard prescriptions used in the new model in detail, the reasons for the tensions between models and observations are discussed in a companion paper where we compare the predictions with fundamental observations, such as the stellar mass function, M_h versus M_⋆ relation, IR luminosity function, and specific star formation rate (Cousin et al. 2015). The companion paper describes, in particular, the impact of SNe feedback and photoionisation process on the stellar mass assembly. It emphasises the difficulty matching observations by using only the standard star formation regulation processes and proposes, in this context, an ad hoc modification of the standard star formation cycle.

The paper is organised as followed. In Sect. 2 we describe the dark-matter simulation used to process the dark-matter evolution. We then present an explicit description of the smooth dark-matter process based on background particles. In Sect. 3, we describe the bimodal baryonic accretion prescription adopted in the model. It is based on a photoionisation model and on two gas reservoirs dedicated to cold and hot phases. Section 4 presents a complete modelling of the monitoring of the hot gas phase. Ejecta processes coming from SNe and/or SMBH activity are also presented, together with the prescription for computating the gas escape fraction and the cooling rate. In Sect. 5, we present all the properties and evolution laws of discs and bulges. The new clumpy component that derives from disc instabilities is also fully described in this section. In Sect. 8 we describe the effect of merger events onto galaxy evolution. Section 9 closes the paper by giving an overview and some research perspectives. In addition, we have added in Appendix A.4, a description of the adaptive time-step scheme used to follow the different dynamical time scales of the gas cycle.

2. dark-matter

In this section, we first describe the dark-matter merger-tree structure and then define the smooth accretion process. We also give the useful dark-matter halo properties. All symbols and their definition are listed in Table 1.

2.1. Dark-matter evolution

Like its predecessor, the new version eGalICS uses a pure N-body simulation to follow the dark-matter evolution. We use an N-body cosmological simulation based on the WMAP-3yr cosmology (Ω_m = 0.24, Ω_Λ = 0.76, f_b = 0.16, h = 0.73). In a volume of (100h^-1)³ ≃ 150 Mpc³, 1024³ particles evolve with an elementary mass of m_p = 8.536 × 10⁷ M_⊙. Dark-matter merger trees (see Fig. 1) are extracted with the halo-finder program AdaptaHOP proposed by Aubert et al. (2004) and extended by Tweed et al. (2009). The dark-matter assembly is therefore pre-computed for the different time steps. Its evolution cannot be modified by baryonic processes.

The halo-finder version used here (Tweed et al. 2009) allows us to extract and separate different levels of structures: main haloes (M-H) and subhaloes (S-H or satellites). In the two cases, we only consider dark-matter structures containing at least 20 dark-matter particles. This limit gives a minimum dark-matter halo mass . Subhaloes are extracted and identified as overdensities in a pre-existing structure that can be a M-H or another S-H. Currently, it is possible to extract up to six different levels of S-H. The halo finder generates all the links between all haloes in time and in structuration level.

2.2. Dark-matter background accretion

To follow the dark-matter halo growth and, later on, the baryonic accretion process, we identify at a given time step n, all particles that are detected in a halo ℋ, and that have never been identified in an other halo ℋ’ before. These background particles () play an important role. Indeed, we consider them as the only source of new baryonic mass. In this framework, the accretion rate on a given structure at a given time step is defined as (1)Figure 1 shows the evolution of the dark-matter accretion rate Ṁ_acc. We can see that before merging with a larger structure (M-H), a halo is systematically identified as a S-H of this more massive halo. During this pre-merger period, it is clear that the dark-matter accretion rate (Ṁ_acc) decreases strongly and may even stop. In these conditions, the diffuse accretion is then supported by the most massive structure. The decrease in background accretion is also visible during fly-by events, as identified in the 18th branch of the merger tree shown in Fig. 1.

By following with time the background accretion process for a given structure i, we can define an integrated accretion mass: (2)In this equation:

• I is the integrated background mass accreted by the halo i up to the previous time step n − 1.

• II is the background mass accreted during the time step [n − 1, n]. This mass is computed using all particles identified in the halo i at the time step n and that have never been identified in an other halo before.

• III corresponds to the integrated background mass accreted by all the progenitors j of the halo i.

2.3. Remark on the dark-matter background accretion

For M-H, the integrated accretion mass (M_acc) is closed to the well known halo finder mass M_fof, which is defined as the mass of all particles instantaneously detected in a structure. In contrast, the integrated background mass M_acc for S-H is often higher than M_fof. Indeed, S-H are sensitive to the stripping process, and a fraction of the halo mass (M_fof) can be ripped by gravitational interactions acting before a merger event. In the context of our dark-matter background accretion process, one of the key points of the model resides in the fact that the stripped particles that have already been identified in a structure cannot be considered in any other structure.

The mass estimator M_acc summarises the accretion history of a given halo. When the instantaneous mass is needed, we use the virial mass M_vir, as in other SAMs. This mass is computed using all particles that follow the virial theorem for a given halo.

2.4. Dark-matter halo properties

2.4.1. Density profile

The spatial distribution of a dark-matter halo mass is described by a Navarro Frenk and White (NFW) density profile (Navarro et al. 1995, 1996, 1997) defined as: (3)where r_dm is the characteristic scale radius of the dark-matter halo and ρ_dm/ 4 is the density at r = r_dm. With this density profile, the dark-matter halo mass M_h( <r) contained in the radius r is given by (e.g. Suto et al. 1998; Makino et al. 1998; Komatsu & Seljak 2001) (4)where φ(x) is a geometrical function defined as (5)Obviously, if r_vir is the virial radius of the dark-matter halo extracted from the dark-matter simulation (see Sect. 2.3 of Hatton et al. 2003 for more information): (6)where c_dm is the dark-matter concentration parameter defined as c_dm = r_vir/r_dm. With the dark-matter virial mass and the virial radius, it is possible to define the dynamical time of the dark-matter structure: (7)This time will be used in the next sections to define the baryonic accretion rate.

2.4.2. Virial temperature and escape velocity

In the hot atmosphere evolution monitoring, we use the virial temperature T_vir linked to the dark-matter halo: (8)We also use the escape velocity of the dark-matter halo defined as (9)and the circular velocity of the dark-matter halo defined as (10)

3. Baryonic accretion

The evolution of the dark-matter component is therefore known using a large set of merger trees. The dark-matter halo properties are listed and defined in the previous section. We now have to add the baryonic content on the hierarchical structure of the dark-matter. This section describes the link between the dark-matter smooth accretion and the baryons that will feed the galaxy. All parameters used in this section and their definition are listed in Table 2.

As explained previously, we consider that the dark-matter background-accreted mass (Ṁ_acc) is the only vector of baryonic accretion. The link between the dark-matter and baryonic accretion is defined through the effective accreted baryonic fraction, f_b (Eq. (11)).

3.1. Photoionisation

The value taken by this parameter f_b follows two distinct regimes. Before reionisation (z_reion), we consider that the baryons feed the dark-matter structures following the universal baryonic fraction . After reionisation, the heating of the gas generated by the first generation of stars limits the baryonic gas accretion in the smallest structures (Kauffmann et al. 1993). Originally proposed by Doroshkevich et al. (1967), the photoionisation mechanism has been developed in the CDM paradigm by Couchman & Rees (1986), Ikeuchi (1986), and Rees (1986). This process is widely used in the context of galaxy evolution (Efstathiou 1992; Babul & Rees 1992; Shapiro et al. 1994; Quinn et al. 1996; Thoul & Weinberg 1996; Bullock et al. 2000; Gnedin 2000; Benson et al. 2002; Somerville 2002; Somerville et al. 2008, 2012; Croton et al. 2006; Hoeft et al. 2006; Okamoto et al. 2008; Guo et al. 2011; Henriques et al. 2013).

The effective baryonic fraction f_b therefore depends on both the redshift z and dark-matter halo mass M_vir. The most commonly used formulation is that proposed by Gnedin (2000) and Kravtsov et al. (2004): (11)Following Okamoto et al. (2008), we use α = 2. The characteristic halo mass follows¹: (12)According to this filtering mass, haloes with M_vir = M_c really collect only half of the baryons coupled with the accreted dark-matter.

Figure 2 shows the evolution of the effective accreted baryonic fraction as a function of both the redshift and the dark-matter halo mass. The mass resolution limit is 20 dm-particules. At this mass scale, the effective baryonic fraction is strongly reduced (f_b< ⟨f_b⟩ / 2) only in the smallest structures evolving at z< 2. As demonstrated in Cousin et al. (2015), the impact of photoionisation is limited. Even if the amount of accreted gas is reduced, it is not sufficient to significantly limit the star formation activity.

3.2. Bimodal accretion

Following the definitions summarised by Eqs. (1) and (11), the baryonic accretion rate on a dark-matter halo is given by (13)As proposed by Khochfar & Silk (2009) we adopt an explicit bimodal accretion composed of a cold and a hot gas flow. These two parallel accretion modes have been discussed in Kereš et al. (2005); Dekel et al. (2009a,b); Khochfar & Silk (2009); van de Voort et al. (2011), and Faucher-Giguère et al. (2011).

On the first hand, we distinguish a cold mode, for which the gas is accreted through filamentary streams. This mode is very efficient. Even if the gas is previously shock-heated (Nelson et al. 2013), it cools with a very high efficiency (t_cool ≪ t_dyn) and falls in the centre of the halo in an approximately free-fall time (t_dyn, see Eq. (7)). The first galaxy discs are formed quickly with this accretion mode. This mode dominates the growth of galaxies at high redshifts (z> 3) and the growth of galaxies in the lowest mass haloes (M_vir< 10¹¹ M_⊙) at any times.

On the other hand, in more massive haloes (M_vir> 10¹² M_⊙) and at late times, the accretion is dominated by a hot mode. For these haloes, the cosmological accreted gas falls into the dark-matter structure and is shock-heated to temperatures close to the virial temperature (T_vir). This gas participates in the development of a hot stable atmosphere (T> 10⁶ K) around the central host galaxy.

Consequently, and as shown in Fig. 3, the baryonic cosmological accretion is divided into two phases (M_cold and M_hot). At a given redshift and for a given halo, the hot mass fraction (f_sh(M_vir,z), shock heated) is computed following Lu et al. (2011) prescription (their Eqs. (24) and (25)): (14)Figure 4 gives the evolution of f_sh with redshift and mass. As explained previously, the cold mode is dominated in the smallest structures (M_vir< 10¹¹ M_⊙), and thus f_sh(M_vir,z) ≪ 1. Unlike for these small objects, the cosmological accretion onto the largest haloes is dominated by the hot mode, and f_sh(M_vir,z) ≃ 1. It is interesting to note that, at high redshift (z> 4), even for the largest haloes (M_vir> 10¹² M_⊙), a small fraction (≃10%) of the accreted gas is always attributed to the cold mode. This residual cold fraction decreases with time and tends to zero in the local Universe.

3.2.1. Cold gas accretion, the filamentary part

Before feeding the galaxy, the cold gas falls in the dark-matter potential well and forms some filamentary structures. This cold gas follows a very collimated distribution. The formation rate of these filamentary structures (Ṁ_flt) is deduced from the dark-matter accretion rate (Ṁ_acc, Eq. (1)) and from the fraction of shock-heated gas (f_sh, Eq. (14)). We define (15)The cold gas is represented by the cold reservoir M_cold (Fig. 3). Once collimated in the filamentary structure, the cold gas falls into the centre of the dark-matter halo, and feeds a galaxy disc with a rate Ṁ_ff (see Fig. 3) close to the free-fall rate (ff). We define (16)where ε_ff = 1.0 is an efficiency parameter.

The evolution terms of the cold gas reservoir are thus completely known. They follow (17)

3.2.2. Hot gas accretion, the shock-heated part

The other part of the baryonic cosmological accretion on a given dark-matter halo is an isotropically distributed gas that is shock-heated when its falls in the potential well. The shock is localised around the virial radius (r_vir) of the dark-matter structure. The accretion rate (Ṁ_sh) of this hot isotropic gas phase is given by (18)At any given time, the total hot gas mass (M_hot in Fig. 3) located around the galaxy is the sum of this cosmological shock-heated gas and the ejected gas produced by the galaxy wind (Ṁ_wind in Fig. 3 and Sect. 4.1).

In summary, the time evolution of the hot atmosphere follows the differential equation (19)where Ṁ_sh is given by Eq. (18), Ṁ_wind is the ejecta rate (wind) produced by the galaxy (see Sect. 4.1), and Ṁ_cool is the cooling rate computed by taking the radiative cooling of the hot gas into account (see Sect. 4.6). Finally, Ṁ_IGM is an output rate corresponding to the mass fraction (wind + pre-existing hot gas phase), which leaves the hot atmosphere due to its high velocity (see Sect. 4.3 for more details). Each term of this equation is shown in Fig. 3.

3.2.3. Total baryonic accretion on the galaxy

Figure 5 shows the total baryonic accretion on the galaxy (Ṁ_cool + Ṁ_ff) for a set of redshifts (z ∈ [ 0:6 ]). The grey vertical arrow indicates the transition mass (10^11.4 M_⊙) used in Lu et al. (2011) (Eq. (14)) for the accretion mode repartition. For haloes with mass lower (higher) than this threshold the accretion is dominated by cold mode (cooling process). The colour histogram is built with our reference model. This reference model uses Okamoto et al. (2008) photoionisation prescription (see Sect. 3.1). The decrease in accretion efficiency for massive haloes is clearly visible. For comparison, we plot for haloes with mass lower than M_vir = 10^11.4 M_⊙, the average accretion on the galaxy computed in the case of Gnedin (2000) photoionisation prescription. As explained in Cousin et al. (2015) the Gnedin (2000) photoionisation prescription is more efficient, and in this case, the average accretion is lower than with our reference model. We also show in the figure the mean accretion on the galaxy in a model without any photoionisation. At high redshift (z> 2), the difference between our reference model and the model without any photoionisation is very small. For z< 2, the differences are larger, especially for the smallest structures (by a factor 10 at z = 1). Again, this confirms that the Okamoto et al. (2008) photoionisation prescription only has a real impact on the smallest structures (M_vir< 10¹⁰ M_⊙) and at low redshift. Finally, we plot the mean baryonic accretion computed in hydrodynamic simulations by Ceverino et al. (2010, their Eq. (7), Ṁ) and Fakhouri et al. (2010, their Eq. (2), ) for comparison. We see that our model is in good agreement with these two simulations.

3.2.4. Some remarks concerning the bimodal accretion prescription

Our accretion model is based on an explicit bimodal baryonic accretion. This prescription allows us to follow the mass of cold and/or hot gas evolving around the galaxy. Indeed, even if the original calculation proposed by White & Frenk (1991) included the cooling and the free-fall timescales (Benson & Bower 2011), it was difficult to follow the mass of gas associated with each of these modes. In the context of our explicit two-phase accretion, it is easier to identify the different galaxy populations associated with a given accretion mode. The total fresh gas accreted on the galaxy () in our model has been compared with the prediction from the original calculation proposed by White & Frenk (1991). It differs by less than 10%.

Some recent works based on a new computational hydrodynamic technique (e.g. Nelson et al. 2013), are questioning the existence of cold streams. It seems that the different computational methods (SPH, AMR, or moving mesh) do not give the same results. In parallel, the impact of the numerical resolution is also very important. At the time of writing, the debate is still open.

4. The hot halo phase

We have presented the hot gas accretion mode that is the first input term of the hot gas phase. In a classical galaxy evolution model, SNe and the SMBH activity heat a fraction of the gas. This gas can be ejected from the galaxy (Ṁ_wind) and thus has to be added to the hot gas phase. It is considered as the second input term of the hot gas phase. This section describes the two distinct feedback processes, SN and SMBH. All symbols used in the equations of this section and their definition are summarised in Table 3.

We present in the following (Sects. 4.1 to 4.6) a self-consistent model of the hot atmosphere. We begin (Sect. 4.1) by a complete description of the ejecta rate due to SNe and/or AGN. We continue (Sect. 4.2) by a description of the kinetic and thermal energy carried by the ejected mass. We add to this energy calculation a model of energy transfer from the ejecta to the pre-existing hot gas phase. These energy considerations allow us to compute (Sect. 4.3) an evolving hot halo temperature and the fraction of mass that can stay in equilibrium in the dark-matter halo potential well. The other part of the mass, which leaves the hot atmosphere, is removed from the hot gas reservoir M_hot and is transferred to (see Fig. 3):

• A passive IGM reservoir if the structure is identified as a main halo,

• the hot gas reservoir of the host M-H if the structure is a S-H (satellite).

We give the formalism for this hot atmosphere model for two different equations of state. In our reference model, we use the ideal gas case (Sect. 4.5.1). We also give the equations for the polytropic gas case in Appendix B.

SNe and SMBH (see Figs. 8 and 11) generate outflows (Ṁ_wind in Fig. 3). The total galaxy ejecta rate is given by the following relations: (20)We develop these terms in the following section.

4.1. Ejecta processes

4.1.1. SNe ejecta

The instantaneous ejecta rate produced by SNe (Ṁ_wind,sn) is deduced with the principle of energy conservation. This has already been used (see Dekel & Silk 1986; Kauffmann et al. 1993) and Efstathiou (2000) for a review. We link the kinetic energy per time unit to the power produced by SNe following (21)where

• (1 − ε_sn)is the fraction of SNe kinetic energy used to feed the disc velocity dispersion (see Sect. 5.3.1);

• f_Kin,snE_snis the SNe energy converted in kinetic energy. We use E_sn = 10⁴⁴ Joules and f_Kin,sn = 0.3 (Kahn 1975; Aguirre et al. 2001);

• η_snis the number of SNe produced per unit of stellar mass (). This parameter is linked to the initial mass function IMF²,

• Ṁ_⋆is the star formation rate.

This formulation of the kinetic energy per unit of time leads to a degeneracy between the ejected mass rate (Ṁ_wind,sn) and the wind velocity (V_wind). To break the degeneracy, we compute the wind velocity using a model extracted from Bertone et al. (2005) and based on Efstathiou (2000) and Shu et al. (2005). The wind velocity can be linked to the star formation rate (Martin 1999) and seems to be independent on the galaxy morphology (Heckman et al. 2000; Frye et al. 2002). We therefore use the following relation (Bertone et al. 2005, their Eq. (9)) (22)

4.1.2. Super-massive-black hole activity

An accreting SMBH produces wind. The outflow rate (Ṁ_•,out) is linked to the effective accretion rate on the SMBH (Ṁ_•,acc) by the relation (Ostriker et al. 2010) (23)This parameter η_• is currently distributed in the range [ 0.1:1 ] (Ostriker et al. 2010). In our reference model we use η_• = 0.6. If Ṁ_•,inf is the total infall rate (see Sect. 6 and Eq. (88) for explicit prescriptions), the conservation of the mass gives obviously (24)and therefore (25)Using the same energy conservation criterion as that used for SNe ejecta (see Sect. 4.1.1 and Eq. (21)), we obtain (26)where

• f_Kin,•is the fraction of energy produced by the accretion on the SMBH that is converted into kinetic energy. This parameter is not well known but observations and numerical simulations give a range of plausible values, f_Kin,• ∈ [ 10^-4:10^-3 ] (Ostriker et al. 2010; Proga et al. 2000; Stoll et al. 2009). In our reference model we use f_Kin,• = 10^-3. This value is compatible with observed outflow velocities, as given for example by Emonts et al. (2005) and Morganti et al. 2005a,b.

• V_jetis the outflow velocity.

Observations of ionised and neutral gas outflows in radio galaxies show that the SMBH activities and, more precisely, the jet have an impact on the gas content (e.g. Morganti et al. 2005a,b; Emonts et al. 2005; Nesvadba et al. 2006, 2008; Lehnert et al. 2011; Guillard et al. 2012). However, it is still unclear how the gas is affected. To compute the galaxy ejecta rate linked to the SMBH activity (Ṁ_wind,•), we use the momentum conservation between the SMBH jet and the gas in the galaxy. We assume that a fraction of the gas can reach a velocity equal to the galaxy escape velocity () thanks to momentum transfer. The total ejecta coming from the SMBH activity is therefore (27)where

• ε_• = 0.6is the coupling factor between the jet momentum and the gas reservoir, This factor is adjusted to produce an outflow rate of a few solar masses per year during the secular evolution of galaxies and a few tens of solar masses per year during the merger induced activity (e.g. Morganti et al. 2005a,b; Emonts et al. 2005).

• is the escape velocity of the galaxy (Eq. (56)).

4.2. Energy of the ejecta, transfers to the hot gas phase

As developed in subsequent sections, we introduce in our model an explicit description of the energetic content of the hot gas surrounding the galaxy. Indeed, accretion and ejecta processes generate energy. A fraction of this energy (kinetic, thermal, and luminous) leaves the galaxy and is distributed in the hot halo gas phase. This energy contributes to the heating of the hot gas.

We follow the internal energy of

• the hot gas evolving around the galaxy: E_Th,atm (28);

• the hot accreted gas: E_Th,acc (29);

• the wind generated by the galaxy: E_Th,wind (32).

With these three terms, it is possible to follow the evolution of the mean temperature of the hot gas phase and estimate the mean temperature of the wind produced by the galaxy.

In this section, we separately discuss the three terms that are considered as three different internal energy sources. The properties of the hot gas are then computed using the following steps:

• 1.

  The internal energy of the accreted gas E_Th,acc is added to the internalenergy of the pre-existing hot gas E_Th,atm.

• 2.

  We then deduce a mean temperature for the pre-existing hot gas phase.

• 3.

  We deduce a mean temperature for the wind phase from its internal energy (E_Th,wind).

• 4.

  We deduce a mean temperature for the wind phase, by taking these two mean temperatures ( and ) into account in a distribution function, f^†. We finally compute the escape mass fraction f_esc and the effective mean temperature .

4.2.1. Thermal energy of the pre-existing hot phase

At a given time t, the gas mass in the hot pre-existing atmosphere is M_hot. This gas is considered in thermal equilibrium in the dark-matter potential well. If the mean temperature of the gas³ is , and removing the mass that has condensed during the time step (Ṁ_coolΔt, Eq. (54)), the total internal energy of the pre-existing hot gas phase at the end of the time step (t + Δt) is given by (28)

4.2.2. Thermal energy in the gas accretion

As explained in Sect. 3.2.2, the hot isotropic atmosphere is also fed by cosmological hot (shock-heated) accretion. The thermal energy of the gas accreted during the previous time step Δt is (29)where T_vir is the virial temperature of the dark-matter halo (Eq. (8)) and Ṁ_shΔt (Eq. (18)) is the cosmological hot gas mass accreted during the time step Δt.

4.2.3. Thermal energy in the ejecta

During a SNe explosion, a fraction f_Kin,sn of the total energy (E_sn) is converted into kinetic energy. The remaining fraction is distributed in luminous energy (photons) and in thermal energy in the ejected gas. We can write the instantaneous non-kinetic power produced by SNe: (30)The instantaneous non-kinetic power produced by SMBH activity follows (31)Equations (30) and (31) are used to compute the instantaneous power. We then compute the (time-)average on the evolution time step Δt of

• the non-kinetic power, ;

• the SMBH non-kinetic power, .

With these two quantities, we can derive the thermal energy carried by the galactic winds during Δt following (32)where we define f_Th,sn and f_Th,• as two free parameters describing the fraction of the non-kinetic energy distributed in thermal energy by SNe and SMBH feedback, respectively. The reference values are set to f_Th,sn = 0.1 and f_Th,• = 10^-3. These values have been chosen so as to produce a temperature distribution centred on:

• 10⁶K for SN wind. Indeed, typical SN remnant emission shows temperatures close to k_BT ∈ [ 0.1:0.7 ] keV, which correspond to T ∈ [ 1:8 ] 10⁶ K (Koo et al. 2002; Sasaki et al. 2014);

• 10⁷K for SMBH wind.

4.2.4. Temperatures

During a time step Δt, the hot gas phase is fed by accretion, cooling, and ejecta processes. In this context, the mean temperature of the hot atmosphere after a time step Δt is given by (33)In parallel we consider that the temperature of the wind is constant during the time step Δt. It is given by (34)

4.3. Escape fraction and mean temperature

4.3.1. Average velocity and mass of the galaxy wind

On the galaxy scale, the instantaneous total wind velocity (SN and SMBH) is computed using velocity winds of each component weighted by their ejecta rate. We define (35)This instantaneous wind velocity is only valid on the galaxy scale. As described in Eq. (A.4) the evolution of the galaxy components (disc and bulge) is performed with an adaptative time step δt shorter than the hot atmosphere time step Δt = ∑ _iδt_i. To take the impact of the hot gas wind on the hot gas atmosphere into account, such as the permanently lost fraction (see Sect. 4.3 and Eq. (41)), it is necessary to estimate the average wind velocity for the time scale used for the hot phase. The wind velocity is therefore computed following (36)The total mass contained in the wind is therefore given by (37)

4.3.2. Velocity distributions

We assume that the hot atmosphere in which the wind extends is in hydrostatic equilibrium in the dark-matter halo potential well. The hot gas phase has a mean temperature . In these conditions, the velocity probability distribution of a particle in this hot gas phase is given by the well-known Maxwell-Boltzman distribution, (see Eq. (C.1)).

Concerning the wind, in the gas frame of reference, the probability velocity distribution is also given by a Maxwell-Boltzman distribution, . In the fixed halo (or galaxy) referential, if we take the mean velocity of the wind into account, the velocity probability distribution is also given by the conditional following relation (38)with v′ = v − ⟨V_wind⟩.

At any given time, we can consider that the hot gas phase is composed of a pre-existing hot atmosphere in which the hot gas produced by feedback processes is injected. Therefore, the global velocity probability distribution is given by, (39)with (40)

4.3.3. The effective escape fraction

At a given time some particles in the pre-existing hot gas phase and even more in the wind phase can have higher velocities than the escape velocity of the dark-matter structure V_esc,dm (Eq. (9)). We assume that only the dark-matter mass and therefore the dark-matter gravitational well have a direct impact on the hot gas phase confinement. Suto et al. (1998) demonstrate that the self-gravity of the hot gas phase has an impact on the hot gas phase density profile but only at large radius (r>r_dm). The density profile of the hot gas is an important quantity in the gas cooling modelling, but if we consider as Suto et al. (1998) that only the core part of the hot atmosphere can cool in a short time, the impact of not taking the self-gravity on cooling processes into account is small.

The mass fraction f_esc that can leave the dark-matter potential well is given by (41)Figure 6 shows the distribution of (Eq. (39)) for two different dark-matter halo masses. In the two cases the wind velocity is close to V_wind = 250 km s^-1 but for M_vir = 10¹⁰ M_⊙, the temperature of the wind is close to K, while for M_vir = 10¹³ M_⊙, it is close to K. As expected, the fraction of the wind mass that can definitively leave the structure is larger for the lowest dark-matter halo mass. Indeed, for such small structures, all the mass contained in the wind leaves the potential well.

After ejection and evaporation, the distribution and the average temperature are recomputed. In Fig. 6 the new distribution is shown in red. For the halo with M_vir = 10¹³ M_⊙, the temperature has increased because the hot wind has participated in heating the hot isotropic atmosphere.

4.3.4. Mass transfer and mass loss

The mass M_IGM = f_escM_hot leaves the hot halo owing to its high velocity. If the galaxy that generates the outflows is hosted by an S-H (satellite), we recompute the ejection processes, following the previous description, in the host M-H referential. Indeed, the halo finder algorithm creates links between S-H and their host M-H. We compute the escape velocity of M-H, and deduce the fraction of ejected mass that can leave its dark-matter potential well. The difference between the mass lost by S-H and M-H is added to the hot reservoir of M-H. The fraction of the mass that can leave M-H is definitively removed from the hot reservoir and added to a passive IGM reservoir (see Fig. 3).

After mass ejection, the average temperature of the gas is modified. This new mean temperature is then given by (42)

4.4. Some remarks concerning the SMBH feedback efficiency

SMBH feedback has been implemented in different models (e.g. Croton et al. 2006; Cattaneo et al. 2006; Bower et al. 2006; Malbon et al. 2007; Somerville et al. 2008; Ostriker et al. 2010; Bower et al. 2012). For example, in Croton et al. (2006; their Eqs. (10) and (11)) and/or Somerville et al. (2008; their Eq. (21)) SMBH activity directly affects the cooling of the hot gas phase. Indeed in these two models, a fraction of the power generated by the accretion is considered as a heating gas rate Ṁ_heat. This heating rate is then used to compute an effective cooling rate, Ṁ_cool,eff = MAX(0,Ṁ_cool,eff − Ṁ_heat). This approach leads to a very efficient SMBH feedback in which a large fraction of the power generated by the accretion is directly used to reduce the cooling process.

In our approach, SMBH feedback contributes to gas ejection (Eq. (27)), but the impact on the cooling rate is only given by a possible increase in the hot gas phase temperature due to wind (thermal energy in the wind, Eqs. (31) and (32)). In this context, even if the SMBH activity produces high-temperature wind (≃10⁷ K), our SMBH feedback implementation leads to less efficiency and therefore has a weaker impact on the cooling process.

4.5. The hot gas density profile

To describe the hot isotropic atmosphere located around the galaxy and to take the impact of the cooling process into account, we need to define a density profile. This section is dedicated to describing of i) the density profile and ii) the cooling process. All the parameters and their definitions are listed in Table 4.

We first assume that the hot stable mass M_hot (in hydrostatic equilibrium) is enclosed in the virial radius r_vir of the dark-matter structure. To model the mass distribution of the hot gas we assume a spherical geometry and thus define (43)The density distribution ρ_g(r) of the hot atmosphere is deduced from the hydrostatic equilibrium condition (HEC) computed in a potential dominated by the dark-matter, (44)where G is the gravitational constant, P(r) the pressure radial profile and M_h( <r) the dark-matter mass enclosed in the radius r. In our model it is possible to use two different density profiles deduced from an ideal or a polytropic gas. In the reference model, we use the isothermal gas case. For completeness we describe the polytropic gas case in Appendix B.

4.5.1. The isothermal ideal gas case

For an ideal (id) gas, the equation of state links the pressure, the density, and the temperature by the well-known following relation (45)where k_B is the Boltzman constant, μ the mean molecular weight of the gas, m_p the proton mass, P^id(r) the radial pressure profile, the density profile, and T₀ the central temperature of the gas. In this simple case of an ideal gas we assume, that the radial profile of the temperature T(r) is set to a constant value and that this constant value is equal to the central value T₀. Here, T₀ is equal to the average value . The HEC in the dark-matter halo potential (Suto et al. 1998; Makino et al. 1998; Komatsu & Seljak 2001; Capelo et al. 2012) gives (46)where ρ_dm and r_dm are the core density and the characteristic radius of the dark-matter halo, respectively and φ(x) is the geometrical function given in Eq. (5).

For an NFW profile (Navarro et al. 1995, 1996, 1997), the hot gas profile that follows the HEC is Suto et al. (1998); Capelo et al. (2012)(47)with ℱ^id(x), a geometrical function defined as (48)where Φ(x) is another geometrical function based on the previous one φ(x) (Eq. (5)) and defined by (49)In the case of an ideal gas, the geometrical distribution of the hot gas is completely set by the dark-matter properties. The last parameter that remains to be determined is the gas core density ρ_g,0. It is extracted from the total mass contained in the hot halo (Eq. (43)) as (50)

4.6. Cooling process

The cooling and condensation of the hot atmosphere is computed using the classical model initially proposed by White & Frenk (1991). We recall the cooling time function here, which gives the cooling time of the gas mass enclosed in a radius r and which follows a given density ρ_g(r) and temperature profile T(r), (51)The cooling time function also depends on the mean molecular weight of the gas μ, the proton mass m_p, and on the cooling function Λ that depends on the temperature T(r) and metallicity Z_g of the hot gas. We use the cooling function Λ from Sutherland & Dopita (1993), for temperature between 10⁴ K and 10^8.5 K and metallicity between 10^-30Z_⊙ and 10^0.5Z_⊙.

4.6.1. Cooling time, radius, and rate

To compute the cooling rate, and therefore the condensed mass of gas that can fall into the galaxy, we need to define the cooling radius, r_cool which the radial extension of the hot atmosphere that can condense during a given time t_cool. The cooling radius r_cool is the solution of the following equation (52)where t(r) is the cooling time function (Eq. (51)).

Different definitions for t_cool exist in the literature (e.g. Cole et al. 2000; Hatton et al. 2003), and thus for r_cool (see De Lucia et al. 2010, for a discussion). In our model, we link t_cool to a cooling clock:

• The cooling clock starts when the hot reservoir M_hot receives some mass. The clock runs as long as the hot phase contains gas. Therefore, t_cool increases with the evolving time of the hot atmosphere.

• The cooling clock is set to 0 when a major merger event occurs. We consider a major merger event when the mass ratio of the two hot gas reservoirs ( and ) (53)

• After a minor merger, the cooling clock is set to the cooling time t_cool of the most massive progenitor⁴.

The cooling rate of the hot atmosphere is then computed following (54)The mass enclosed in the cooling radius r_cool falls into the galaxy with a rate close to the free-fall rate. We assume that the hot atmosphere extends up to the virial radius, so that the cooling radius cannot be larger than this virial radius. The free efficiency parameter is set to ε_cool = 1.0.

5. Galaxy components

We present the relations and definitions linked to galaxy formation and evolution in the next sections. An eGalICS galaxy (Fig. 7) is composed of two distinct but interacting components: a disc (Fig. 8) and a bulge or bulge-like component ℬ (Fig. 11).

5.1. Accretion on the galaxy

At the first steps of galaxy formation, the cold gas (from filamentary collimated structures and/or the isotropic cooling flow) falls into the center of the dark-matter halo. We assume that this cold gas initially forms a thin exponential disc. The gas acquires angular momentum during the mass transfer (Peebles 1969). After its formation, the disc is supported by its angular momentum. This paradigm is based on the prescription given by Blumenthal et al. (1986) and Mo et al. (1998) and has been frequently used in SAMs (e.g. Cole 1991; Cole et al. 2000; Hatton et al. 2003; Somerville et al. 2008). During the secular evolution of the galaxy, some mass may be transferred to a pseudo-bulge component by disc instabilities (see Sects. 5.4.3 and 5.4.4 for more information). After a major merger event (see Sect. 8.3), the remnant galaxy has a spheroidal morphology (see Sect. 5.5). The gas contained in the remnant galaxy and the freshly accreted gas, which will fall into the centre of the halo in the next time step, will then form a new disc.

5.2. The circular and escape velocity of the galaxy

The radial profile of the circular velocity of the galaxy is deduced from the dynamic of all its components, the dark-matter halo , disc , and bulge ℬ. We obtain (55)All the terms will be defined in their corresponding description section (Eqs. (10), (60), and (83)).

The escape velocity of the galaxy is computed following (56)where and are the galaxy and dark-matter mass enclosed in the galaxy half-mass radius , respectively.

5.3. The disc

Table 6 summarised the parameters of the disc. In the classical scenario adopted here, the disc component is the first structure formed in a new galaxy. We use the standard approach as proposed by Cole et al. (2000); Hatton et al. (2003), among others.

We assume that the disc component is infinitely thin. Its mass radial distribution is given by (57)where is the dimensionless radius .

As explained previously, we assume that the gas accretion is only supported by the disc. At a given time n, the disc is defined by its total mass and its exponential radius . During a time step Δt = t_{n + 1} − t_n, we assume that the new accreted mass acquires angular momentum and therefore forms a very thin disc structure with a mass (Ṁ_ff + Ṁ_cool)Δt and an exponential radius deduced from the dark-matter halo virial radius r_vir using the well-know spin parameter λ, (Blumenthal et al. 1986; Sellwood & McGaugh 2005; Macciò et al. 2008; Cole et al. 2008; Antonuccio-Delogu et al. 2010; Muñoz-Cuartas et al. 2011).

The assembly of the pre-existing disc with the accreted mass leads to a secular evolution of the disc exponential radius that follows (58)In the model, the disc is defined as an object that contains two different parts: (59)The radial profile of the disc circular velocity is given by (Freeman 1970) (60)where I_r and K_r are modified Bessel functions of rank r.

5.3.1. Disc velocity dispersion and dynamical time

In this section, we use a new method based on energy conservation that is used to monitor the average velocity dispersion of the gas in the disc component. At a given time n, we assume that the energy dispersion of the gas is given by (61)where is the total mass of the gas in the disc and the average velocity dispersion in the disc component. At time n the kinetic energy contained in the disc rotation is (62)where is the circular velocity of the galaxy (Eq. (55)). Here the integral is computed from 0 to , which encloses 99.99% of the disc mass.

As proposed by Khochfar & Silk (2009) or Ocvirk et al. (2008), we assume that the velocity dispersion is generated by gas infall and SNe energy injection. Since the infall mass is structured through a very thin disc in our model, assuming that the dynamics of the gas is given by the total potential well gives the following definition for the kinetic energy transported by gas fuelling (63)where Ṁ_cool and Ṁ_ff are the cooling rate and the free-fall rate fuelling the galaxy, respectively.

In addition we take the SNe kinetic energy injection in the gas into account, (64)where

• is the fraction of SNe kinetic energy used to feed disc velocity dispersion;

• η_snis the number of SNe produced per unit of stellar mass formed (). This parameter is linked to the initial mass function IMF;

• f_Kin,snE_snis the SNe energy converted in kinetic energy. We use E_sn = 10⁴⁴ Joules and f_Kin,sn = 0.3 (Kahn 1975; Aguirre et al. 2001);

• is the star formation rate in the disc.

To compute the evolution of the average velocity dispersion, we assume that the variation of the rotational energy between times n and n + 1, () is supported by gas infall energy (E_inf ≥ ΔE_V) and that the variation of the energy dispersion () is then given by (65)where f_disp = 0.05 is a dissipation factor. The mean velocity dispersion in the disc is therefore given by (66)Figure 9 shows the evolution with redshift of the normalized distribution of . We have assumed an isotropic three-dimensional velocity dispersion, . While high-redshift discs present an average ratio close to 0.6, local discs have ratios close to 0.1. The disturbed disc population disappears progressively with time and a more stable population is generated.

The disc dynamical time is set to the minimal time between the disc’s full rotational time and local velocity dispersion. We use (67)

5.4. Disc instability, bulge-like component formation

For more than one decade, observations (e.g. Cowie et al. 1995; van den Bergh 1996; Elmegreen & Elmegreen 2005; Genzel et al. 2008; Bournaud et al. 2008) and hydrodynamic simulations (e.g. Bournaud et al. 2007; Ceverino et al. 2010, 2012) have shown that there are of gas-rich, turbulent discs. These kinds of discs are unstable, and undergo gravitational fragmentation that forms giant clumps. Such clumps interact and migrate to the centre of the galaxy where they form a bulge-like component (e.g. Elmegreen 2009; Dekel et al. 2009b).

Based on these works, we assume in our model that mass overdensities and low velocity dispersion (Eq. (66)) lead to the formation and the migration of giant clumps. We present an analytic self-consistent model of disc instabilities here. All parameters used in this section are described in Table 7.

5.4.1. Origin of instabilities

Following Toomre (1963, 1964), a disc becomes unstable if its gas surface density () is high and its circular velocity and/or mean velocity dispersion is low. The disc stability is controlled by the following parameter (68)where (69)is the mass-weighted epicyclic frequency with Ω the angular velocity and the gas mass surface density.

In Eq. (68), σ_r is the mean radial velocity dispersion. As for the disc dynamical time computation (Eq. (67)), we assume an isotropic three-dimensional velocity dispersion, . The gaseous component of the disc is considered as stable or marginally stable if Q>Q_crit = 1.0. We use this stability criterium to evaluate the mass that will form clumps by gravitational instabilities.

5.4.2. The clump component,

We assume that disc instabilities are linked to the formation of giant clumps. These giant clumps are formed into the disc and migrate to the centre of the disc to form a pseudo-bulge component (see Fig. (8)). As previously, all parameters used in the description of disc instabilities linked to clumps are listed in the dedicated Table 7.

In our model, the clumpy phase () is not considered as a new separated component but is considered as an inherent part of the disc. The amount of gas in clumps is therefore given by a fraction of the disc mass. As for the disc component, clumps contain a cold star-forming gas and a stellar population . We monitor the stellar and the gaseous mass in clumps by taking their formation (Sect. 5.4.3), transfer (Sect. 5.4.4), and disruption (Sect. 5.4.5) into account. Star formation in the disc is divided into a clumpy and a homogeneous component.

5.4.3. Clump formation

We consider that disc instabilities (and therefore clumps formation) are only due to an excess of mass. When the homogeneous (non-clumpy) gas fraction of the disc ( is unstable (Q<Q_crit), we divide the gaseous mass of the disc in two parts, a stable and an unstable part . We have (70)with Q′ is a reduced-Toomre parameter that takes into account only the fraction of the disc homogeneous gas mass (the gaseous mass of the clumps being subtracted in the gas surface density computation). Within this framework, at a given time, a new amount of unstable gas feeds the clumpy part, (see Fig. 8). The unstable gas mass is given by (71)Knowing the new unstable mass in the disc, we can compute the instantaneous formation rate of the clump component (, see Fig. 8) following (72)Clumps can form stars. The dynamical time for star formation in clumps (73)where and are the characteristic radius and mass of one clump, which we extract from Dekel et al. 2009b (their Eqs. (8) and (9)), (74)and (75)In these relations, δ is the disc mass fraction defined as (76)where , , and are the mass, enclosed in (radius for 90% of the disc mass), of the dark-matter halo, disc, and bulge component, respectively. As explained in the following section, the individual mass is used as a threshold mass that allows us to compute the clump transfer to the pseudo-bulge component.

5.4.4. Clumps transfer

To take the migrations into account, we define as the mass of clumps that is transferred to the pseudo-bulge. The transfer rate, , writes as (77)where is the total mass in clumps, and the mass of the clumps that is already in a transfer process. We add as the typical clump transfer time scale. Like other clump properties, we compute this characteristic time following Dekel et al. (2009b; their Eq. (19)) (78)When the transferred mass reservoir becomes larger than the reference mass , the stellar and gas mass contained in this reservoir () is subtracted from the disc and added to the pseudo-bulge component. The transfer is considered as a micro-merger event. When the micro-merger is finished, we reset the transferred mass, and recompute the characteristic individual mass of clumps . This new characteristic mass will be the new threshold mass for the next clump migration process.

5.4.5. Clumps disruption

At a given time, a disc may have giant clumps (). If, at the same time Q>Q_crit then we compute the disruption rate (79)where is the disruption time computed as proposed by Dekel et al. (2009b; their Eq. (14)): (80)In Fig. 10 we show the evolution of the Toomre stability parameter for one disc chosen arbitrarily for z< 11. We also show the evolution of the number and individual mass of clumps for the unstable period. We see that for 0.8 <z< 1.5 the disc instability generates numerous clumps that are destroyed after each merger event and systematically formed again but with higher individual mass. During this period of instabilities, the decrease in the number of clumps has to be linked to the migration process. After z = 0.8 only one giant clumps can be formed, but with a higher mass. For z< 0.5 the disc is stable (Q>Q_crit = 1), the clumps disruption process is therefore turned on and the clumps mass decrease. After the last merger event, the disc does cannot form any clumps owing to the lack of gas.

5.4.6. Summary of mass transfers

The evolution of the mass transfers of the different components of the disc is governed by the following equations: (81)

5.5. Pseudo-bulge or spheroidal galaxy ℬ

In our model, we assume the same properties for a post-major merger galaxy (with spheroidal morphology) and for a pseudo-bulge formed by giant-clumps migration. The mass transfers of these components ℬ are illustrated in Fig. 11, and the parameters are listed in Table 8.

As for the disc, we treat the (pseudo-)bulge as a numerical object composed of a cold gas reservoir (M_ℬ,g^⋆) and a stellar population (M_ℬ,⋆). In a pseudo bulge the cold gas reservoir is fed by clump migration. This gas can be used to form new stars. In a spheroidal galaxy, formed during a merger event, the gas reservoir is only fed by the stellar cycle (see Fig. 13). In this context, the amount of gas is very low and therefore no star formation is possible. At each merger, the gas contained in the bulge component is added to the disc component of the remnant galaxy.

The mass distribution of the spheroidal galaxy or the pseudo-bulge is described by a Hernquist (1990) model, (82)where r_ℬ is the characteristic radius of the bulge.

5.5.1. Bulge circular and escape velocity

The circular velocity of the bulge or pseudo-bulge is given by (Hernquist 1990) (83)with (84)where M_ℬ is the total mass (gas + stars) of the bulge, M_• the mass of the SMBH, and r_ℬ the characteristic radius. We give also here the escape velocity of the bulge, (85)

5.5.2. Bulge dynamical time and velocity dispersion

As commonly used in SAMs, we define the dynamical time of the bulge using its velocity dispersion (σ_v,ℬ) instead of its circular velocity, (86)where r_ℬ,50 is the radius that enclosed 50% of the bulge mass, and where we used the definition of the velocity dispersion proposed by Hernquist (1990; their Eq. (10)), computed at the half-mass radius.

5.5.3. Summary of the mass transfers

The physical processes acting in the bulges are linked to clump migration and gas ejection (Ṁ_wind,agn due to SMBH activity). As in the disc, the gas in the bulge is converted into stars and these stars feed the gas reservoir (Ṁ_wind,⋆, Sect. 7 and Fig. 13). We summarize the mass transfer rates with the following equations: (87)

6. The SMBH component

An SMBH can be formed and evolve in the centre of the galaxy. This section describes the hypothesis behind the formation of an SMBH and the accretion mechanisms that govern its evolution. The parameters used in this section are listed in a dedicated table (Table 9).

We decided to form the SMBH by converting a fraction of the bulge mass into an SMBH when the bulge mass (M_ℬ) becomes greater than a threshold M_ℬ,• = 10⁶ M_⊙. We assume that the formation process of an SMBH is only possible during a merger. Indeed, we assume that the SMBH formation needs powerful processes. Therefore, after a merger, if the total mass of the bulge is higher than M_ℬ,• an SMBH with a mass M_•,min = 10³M_sun is created. We subtract the equivalent mass in their respective proportions from the gas phase and from the stellar population. As can be expected, the hypothesis of a constant seed mass as a threshold to create an SMBH generates a scatter in the M_• versus M_ℬ relation. This scatter is similar to what is obtained by Croton et al. (2006).

The SMBH evolves in the centre of the bulge and grows by accretion processes. The total accretion rate is given by (88)where is an accretion rate linked to the evolution of a gaseous torus formed around the black hole and defined as (89)In this last equation, is the total mass of gas in the torus and t_dyn,ℬ the dynamical time of the bulge. During the galaxy evolution processes the torus is fed by

• the gas coming from clumps migration (),

• a fraction of the bulge gas ()⁵, instantaneously added during a merger event.

In Eq. (88), Ṁ_Bondi is the classical Bondi (1952) accretion rate. We use the same calculation as proposed by Somerville et al. (2008), (90)where T_• ≃ 10⁷ K is the mean temperature in the accretion torus, Λ the cooling function, and M_• the SMBH mass. At any time, the accretion onto the SMBH is possible only if the torus contains gas () and only in the Eddington accretion rate limit, (91)We recall that the complete description of ejecta processes driven by the SMBH is given in Sects. 4.1.2 and 4.2.3.

The instantaneous bolometric luminosity produced by the SMBH is deduced from Eq. (31) and is set to (92)We show in Fig. 12 the evolution of the fraction of dark-matter haloes that hosts an SMBH, as a function of the halo mass. The fraction increases with M_h. In the model presented here, all structures with dark-matter halo mass larger than ≃4 × 10¹¹ M_⊙ host an SMBH, regardless of the redshift.

We also show the fraction of active SMBH in Fig. 12. An SMBH is considered as active if L_•,bol> 0.05L_Edd. For redshift z> 2 this fraction is large then >40% but it decreases for structures evolving at lower redshift (z≤1). This decrease is mainly visible at the highest and lowest mass. Haloes with intermediate mass (10¹¹<M_vir< 10¹² M_⊙) host most of the active SMBH.

7. The stellar component

In the disc, bulge, and clumps, the gas (, M_g^⋆,ℬ, , respectively) is converted into stars (, M_⋆,ℬ, , respectively). We adopt the same definition for the star formation rate in all these components, (93)In this formulation, ε_⋆,i is a free efficiency parameter, depending on the component. We adopt ε_⋆,i = 0.02 in discs and bulge and 0.2 in clumps. These values allow to be reproduced the Kennicutt (1998) law on both galactic and clump scales. Here M_g^⋆,i is the total mass of star-forming, gas and t_dyn,i is the dynamical time (Eqs. (67), (86), and (73) for the disc, bulge and, clumps, respectively).

If the disc contains some clumps, the star formation is divided into two parts: i) the homogeneous gas part in the disc (, Sect. 5.4.3 and Eq. (70)) and ii) the clumpy part (). The overall disc star formation rate is computed by considering the fraction of mass in each of the two components (94)As shown in Fig. 13 the evolution of stellar population is closely related to the gas reservoirs. At the end of their life, stars re-injects gas into the galaxy’s interstellar medium. This gas is enriched in metals. As in its predecessor (Hatton et al. 2003), the metals enrichment is followed using the Devriendt et al. (1999) prescription.

8. Mergers

8.1. Major or minor merger?

We define the merger type (major or minor) using the following mass ratio, (95)where is the total mass (galaxy and dark-matter halo), enclosed in the galaxy half-mass radius . We consider a minor merger event when η_merger< 0.25. Table 10 gives symbols and their definition associated with mergers.

8.2. Minor mergers

During a minor merger, we keep the disc and the bulge. The stellar populations of each disc (bulge), are added to form the remnant disc (bulge). In contrast, cold gas reservoirs are added in the remnant disc. The size of the remnant disc (bulge) is set to the larger disc (bulge) progenitor size. The velocities (dispersion and circular) are recomputed with the properties of the remnant galaxy and its dark-matter halo.

We consider that the minor merger of two discs has a strong impact on the evolution and formation of clumps. We decided to destroy all clumps in the remnant disc. New clumps have to be formed in the new disc after the merger, by taking the post-merger dark-matter and disc properties into account.

8.3. Major mergers

In the case of a major merger (η_merger> 0.25), the galaxy structure and dynamic are completely modified. The stellar population of the remnant galaxy presents an elliptical form modelled by the same Hernquist (1990) mass distribution than the bulge component. In parallel the gas components of the two progenitors are added and a new disc is formed. The characteristics of this new disc are computed in the context of the post-merger dark-matter halo (r_vir,λ).

At the merger time, the half-mass bulge radius r_ℬ,50 is computed using (Hatton et al. 2003), (96)where M_{1 / 2,i} has been introduced in Eq. (95). The characteristic bulge radius (r_ℬ) is then deduced following (97)In the case of a micro-merger event, when a clump migrates to the centre, we use the same framework as in Eq. (95), M_{1 / 2,1} being in this case the bulge plus dark-matter halo mass enclosed in the bulge half mass radius, and the transferred clump mass.

8.4. SMBH: evolution and/or formation

For the SMBH component, if the two progenitors contain an SMBH, their masses are added, and a fraction of the gas contained in the remnant galaxy is added to the SMBH accretion torus (). If only one of the two bulge components contains an SMBH, the remnant galaxy keeps this SMBH component. If the progenitors do not contain any SMBH, but if the remnant bulge mass becomes higher than the mass threshold, an SMBH is created with a mass M_•,min = 10³ M_⊙.

8.5. Merger boost factor

Mergers are violent events. During a short time-scale galaxy evolution processes are strongly modified and secular evolution laws (efficiencies) are no longer valid. To take these modifications into account we define a merger boost factor following (98)In this definition, η_merger is the progenitor mass ratio (Eq. (95)), η_gas the gas fraction in the post-merger structure, τ the time elapsed since the last merger event, and τ_merger = 0.05 Gyr is the characteristic merger time scale. This factor is used to modify the efficiencies just after a merger event. We use MAX[ 1,ε_boost(τ) ] as a correction factor to the star formation law (Eq. (93)), for each component of the galaxy (disc and bulge) and for the SMBH torus accretion rate (Eq. (89)). The formulation used here allows simulating a time-dependent gas compression (decreasing with time), and takes the gas content of the two progenitors into account. The more gas they contain, the higher the gas compression.

8.6. The hot gas phase

During a merger, the hot gas masses of the two haloes are added, and the equilibrium temperature is computed with the thermal energy of each hot atmosphere, (99)which is added to the gravitational interaction energy, (100)where r_hot,i is the half-mass radius. The equilibrium temperature is then given by (101)

9. Conclusions

We have presented in this paper the physical mechanisms acting on the formation and evolution of a galaxy in a dark-matter halo in a semi-analytical framework. The model follows the dark-matter structuration through merger trees. The growth of dark-matter haloes was followed using a set of properties, including the virial mass and radius and density profile parameters. The dark-matter accretion rate is computed using background particles, which were identified as particles that have never been assigned to any other halo at a previous time step. These particles allowed following the baryonic accretion. Using a photoionisation model, the baryonic accretion rate was deduced from the smooth dark-matter accretion rate, built with the background particles. The baryonic accretion is then divided into two parts:

• A cold mode, that dominates the growth of small structures. It relies on the fact that the cooling time of the gas accreted by low-mass haloes (M_vir< 5 × 10¹⁰ M_⊙) is shorter than their dynamical time.

• A hot mode: for massive haloes (M_vir> 5 × 10¹⁰ M_⊙), the accreted hot gas has to cool before collapsing.

In the model, the evolution of the hot atmosphere, in equilibrium in the dark-matter potential well, is based on a detailed monitoring of its average temperature and density profile. These two quantities were followed by using energy conservation laws based on SN feedback and hot-gas accretion. In addition, the coupling between the pre-existing hot-gas phase and the wind phase generated by the galaxy, allowed us to compute the gas escape fraction. The cooling rate, computed using the standard White & Frenk (1991) prescription, relies on this explicit monitoring of the hot-gas phase properties.

Coming from the cooling process or directly from the cold accretion mode, the cold gas forms a disc at the centre of the dark-matter halo. The gas in the disc can form stars. In the fiducial approach, we used the Kennicutt law based on the total gas mass and the dynamical time of the disc.

A feature of this new model is the treatment of instabilities in the disc. Based on Dekel et al. (2009b), the model we developed accounts for the formation, migration, and disruption of giant clumps. The dynamic of the disc is followed using the rotational profile that is deduced from the entire mass distribution (dark-matter halo, disc, and bulge). To this rotational profile is added a model to follow the average gas velocity dispersion. The evolution through cosmic time of the ratio indicates that the population of highly disturbed discs (), formed at high z, decreases progressively with time, and gives birth to a more quiet disc population (). The mean velocity dispersion, added to the gas surface density computation, allows the stability behaviour of the disc to be derived. Using a method based on the Toomre parameter, we compute at each time-step the unstable mass, which is stored in clumps. The clumpy mass then follows a migration process. The mass that is transferred participates to the formation of a pseudo-bulge component.

This pseudo-bulge, or bulge formed during mergers, can host the formation and the evolution of an SMBH. The accretion of mass onto the SMBH generates energy that is at the origin of mass ejection through wind. Together with the mass ejected due to SNe, the SMBH activity participates in the heating of the hot atmosphere around the galaxy. SN and SMBH are the two main process that allow the star formation to be regulated. Indeed gas ejection induces gas depletion in the galaxy and the heating of the hot atmosphere reduces the cooling rate.

In our model all the up-to-date mechanisms used to regulate the star formation process were considered. Some variations of the fiducial recipes are discussed in detail in Cousin et al. (2015). We showed in this paper that even using strong efficiencies for the photoionisation and feedback processes, the star formation activity remains too efficient in low-mass galaxies. We thus introduce an ad hoc modification of the standard paradigm, based on the presence of a no-star-forming gas component, and an artificial concentration of the star-forming gas in galaxy discs. The no-star-forming gas reservoir generates a delay between the accretion of the gas and star formation. This model agrees much better with the observations of the stellar mass function in the low-mass range than the reference model, and agrees quite well with a large set of observations, including the redshift evolution of the specific star formation rate.

SAMs are currently the best models for following a large number of galaxies in their dark-matter halo. Their predictions, deduced from large scale dark-matter simulations, allow interpreting the recent surveys of galaxy population. However, even if the decoupling between dark-matter evolution (merger tree) and the baryonic physics allows exploring a large set of prescriptions in a short computation time, SAMs cannot follow all the complexity of the galaxy evolution. This applies particularly to the computation of the star formation rate. SAMs generally assume that the gas distribution follows an exponential profile, in a infinitely thin disc. All the gas content in the disc is available to form stars. In contrast, real galaxies host a complex gas structuration; high-z structures are generally gas-rich discs and show a finite scale-height. Galaxy discs are a multi-fluid system, in which multi-scale interactions (2D, 3D Shi et al. 2011) generate instabilities that participate in the ISM structuration (e.g. Jog & Solomon 1984; Hoffmann & Romeo 2012). A better description of these ISM structuration mechanisms is essential for better understanding the star formation regulation processes. In a future work, we plan to follow the gas structuration in the disc from the large (λ ≫ h) to the small scales (λ ≪ h), assuming a finite scale-height (h) and using the turbulent cascade description. Assuming that star formation occurs only on the smallest scales, we could estimate the star-forming gas fraction and give a better description of the star formation activity.

---

Acknowledgments

We thank Nicole Nesvabda, Pierre Guillard, and Nicolas Bouché for useful discussions of the SNe and SMBH feedback. We also thank Jérémy Blaizot and Andrea Cattaneo for very insightful discussions.

References

Appendix A: An adaptive time-step scheme

In a galaxy, different processes act at the same time but on different time scales. It is not efficient, for example, to compute the evolution of the cold filamentary phase that evolves with a typical dark-matter dynamical time (t_dyn = 10⁶ yr) and with a time step dedicated to the ejecta rate of the galaxy (5 × 10³ yr). If we use the same time step to follow these two components, we generate numerical errors on the component with the shorter dynamical time.

In this context, each part of the gas cycle (in the baryonic halo phase or in the galaxy) must evolve with a time step that is as close as possible to its dynamical time. Accordingly, we have developed an adaptive time-step scheme. Each component evolves passively by receiving and transfers mass until the evolution of an other component significantly affects its own evolution. We consider that a component is modified if the variation ΔM of its mass (gas or stars) is greater than p% = 10%. This value is low enough to accurately follow star formation and ejection processes in gas-rich galaxies up to z ≃ 4. Naturally, for passive galaxy (like spheroidal at low z), a higherr value could be acceptable.

Appendix A.1: Description

At any given time step each component with a mass M must evolve with an input and output rate (Ḟ_in and Ḟ_out). With these two gas flows, the optimal time step, δt^MAX, for a component with a mass M is therefore computed as (A.1)For example, for the cold gas component M_cold (Fig. (3)), Ḟ_in = Ṁ_flt (Eq. (15)) and Ḟ_out = Ṁ_ff (Eq. (16)); therefore, (A.2)At the first time step, the mass is not defined and the reference is then fixed to (A.3)

Appendix A.2: The hot-gas phase

The case of the hot atmosphere is more complicated. Indeed, the hot phase:

• receives gas from the intergalactic medium: Ṁ_sh (Eq. (18));

• follows cooling process: Ṁ_cool (Eq. (54));

• receives gas coming from the galaxy: Ṁ_wind (Eq. (20)).

This last term is not linked to the mass of the hot phase M_hot but depends on the star formation activity of the galaxy. To follow the evolution of the hot atmosphere with a high accuracy, and therefore compute the optimal time step , we must take possible variations in the ejecta rate into account.

In this context, it is important to check whether the evolution of the hot phase is dominated by input (Ṁ_sh + Ṁ_wind) or output (Ṁ_cool) gas flows. In the first case Ṁ_in = Ṁ_sh + Ṁ_wind>Ṁ_out = Ṁ_cool if we take instantaneous values of all input and output rates to compute , an increase, in the next time step, in the galaxy ejection rate would lead to a violation, by excess, of the “stability criterium”: ΔM_hot/M_hot ≤ p%. On the contrary, if the hot phase is dominated by the cooling process Ṁ_in = Ṁ_sh + Ṁ_wind<Ṁ_out = Ṁ_cool, a decrease in the galaxy ejection rate would lead to a larger decrease in the hot gas mass than p%.

We thus have (A.4)where we authorize, depending on the case, an increase or a decrease of p% of the instantaneous galaxy ejection rate Ṁ_wind. These two cases are imprinted in two limits (a lower limit: and a upper limit: ) beyond which the criterion of p% variation is no longer valid.

For the hot-gas phase, with this optimal time step, we can deduce the two lower and upper limits for the galaxy ejection rate. The lower limit reads as (A.5)and the upper limit as (A.6)Between these two limits, the stability criterion is followed. Figure A.1 illustrates this optimal time step computation. In the case chosen for the figure, the evolution of the hot gas phase is dominated by the input flux (Ṁ_sh + Ṁ_wind>Ṁ_cool). The optimal time step for the hot gas phase evolution is therefore computed by assuming an increase in the galaxy ejecta rate.

Appendix A.3: The stop and restart mechanism

The two galaxy ejection rate limits are and . Above and below , the evolution of the hot gas phase violates the p% criterion. These two limits are transferred to the galaxy evolution procedure. During the evolution of the galaxy components, if the average⁶ of the ejecta rate is above or below these two limits, the evolution process is stopped. The properties of the hot gas phase and of the galaxy are computed again at this time, and a new time step starts from this point.

Appendix A.4: Cold and hot phase coupling

In the new eGalICS algorithm, it is therefore possible to compute the optimal time step for the cold () and the hot phases (). The same computation is performed for all galaxy components, such as for the disc and the bulge components, which composed a galaxy. We present in this section the coupling that exists between the hot and cold phases.

We assume that Ṁ_ff and Ṁ_cool are constant during and , respectively. As shown in Fig. A.2 and as explained previously, the computation of and allows some barriers to be established: and . These barriers mark the evolution of the cold and the hot phase. If the galaxy has already evolved until t_gal = t_cold = t_hot, the next galaxy time step is computed following (A.7)During the first part of the evolution presented in Fig. A.2 the time step of the galaxy is limited by the hot gas phase evolution. Indeed, the cold time step is longer (). The hot gas phase and the galaxy therefore evolve on the same time grid. Later on, the galaxy evolution is limited by the cold-phase evolution (). Indeed, Ṁ_ff has to be recomputed after . Owing to the strong coupling between the galaxy and the hot gas phase (ejection processes Ṁ_wind), the evolution of these two components are synchronized at . In these conditions, even if the properties of the hot-gas phase should allow an evolution over a longer time step (until ), its evolution is stopped. The three components are updated, their time evolution are then similar (t_gal = t_cold = t_hot), and the evolution process is started again.

Appendix A.5: Effective time step duration

Figure A.3 shows the distribution of the time step of the dark-matter component, the cold or/and hot phase and of the galaxy component (disc and bulge). The histograms are built following the analysis of a dark-matter halo and its galaxy between z = 11 and z = 0. The time step durations used most for the hot or cold phases are identical to those used for the dark-matter. Indeed, the bulk of the blue histogram coincides with the green histogram. Minimum values coincides with the synchronisation process ( in Fig. A.2). Concerning the distribution of time steps used for the galaxy components (disc or bulge), the maximum value is linked to the stellar evolution process (δt_⋆ = 10^-3 Gyr, Devriendt et al. 1999). The second bump of the grey histogram corresponds to merger time-scale events (τ_merger = 5 × 10^-4 Gyr). Lower values are again linked to the synchronisation process (clump transfers) between the disc and the bulge.

Appendix B: The polytropic hot gas case

For a polytropic gas (poly), the temperature distribution is non-isothermal. The pressure, density, and temperature profiles are linked with the following relations (B.1)where T₀, P₀, and ρ_g,0 are the central temperature, pressure, and density, respectively, and Γ is the polytropic index of the gas. The HEC in the dark-matter halo potential (Suto et al. 1998; Makino et al. 1998; Komatsu & Seljak 2001; Capelo et al. 2012) gives (B.2)with ℱ^poly(x,T₀), a geometrical function defined as (B.3)where Φ(x) is the geometrical function defined by Eq. (49) (and see Table B.1 for parameter definitions).

In the polytropic gas case, T₀ is not equal to the mean temperature of the gas phase but to the mass-average temperature (B.4)Following Eq. (50), the core density ρ_g,0 is given by (B.5)Using Eqs. (B.4) and (B.5), we can compute an implicit relation between and T₀(B.6)where ℱ^poly(x,T₀) (Eq. (B.3)) is a function of the central temperature T₀. Since is a strictly increasing function of T₀ for a set of dark-matter properties (r_dm and ρ_dm), the implicit relation given above can be inverted using a dichotomy-research algorithm. We can thus compute the central temperature T₀ of the polytropic gas connected to the mean temperature deduced from energy exchanges (see Sect. (4.2)).

The geometrical profile defined by Eq. (B.3) must be strictly positive to ensure density (Eq. (B.2)) and temperature (Eq. (B.1)) consistency. This condition leads to a minimum threshold value of the central temperature . Indeed all other parameters are linked to dark-matter halo properties and are already fixed. This critical central temperature depends on the geometrical function Φ(x) (Eq. (49)), which is a strictly increasing function of x. As (B.7)then (B.8)If the central temperature T₀ of the polytropic gas derived from the average temperature is lower than this minimum threshold value , the hot atmosphere is not considered stable enough to activate the cooling processes. The gas contained in this hot phase again falls on the galaxy with a feeding rate close to the free-fall rate (as defined for the filamentary structure). We set in this case (B.9)where M_hot is the total mass of gas contained in the unstable hot phase, t₀ the dynamical time of the dark-matter structure and ε_ff = 1.0 a free parameter (the same as used in Eq. (16)).

Appendix C: Maxwell-Boltzman distribution

We use the following formulation of the well-known Maxwell-Boltzman velocity distribution (C.1)

Appendix D: Data repository

Outputs from the model are available at the CDS under the reference m₁. Data are distributed under *.fits format and are therefore compatible with the TOPCAT software (http://www.star.bris.ac.uk/~mbt/topcat/).

All Tables

All Figures