© M.Cousin et al. 2019, A&A, 627, A132

1. Introduction

Galaxies are defined by their stellar populations, that is “when and where” their stars formed. Therefore, if models are to capture accurately the process of galaxy formation and evolution, researchers must determine how star formation is regulated locally and globally in galaxies. However, star formation is one of the most challenging processes to characterise in galaxy evolution models, essentially because the formation of stars involves many non-linear processes that occur over a wide range of temporal and spatial scales in for example density, velocity, and magnetic field strength and regularity (Kritsuk et al. 2013; Krumholz & McKee 2005). Observations show that star formation in galaxies is a very inefficient process. Only 0.1%–10% of the available gas is converted into stars per local free-fall time (e.g. Lada et al. 2010; Hennebelle & Falgarone 2012; Agertz & Kravtsov 2015; Lee et al. 2016). On the other hand, numerical simulations of molecular clouds indicate that star formation efficiency is highly dependent on how ionisation and kinetic feedback is injected into the interstellar (ISM; e.g. Krumholz & Thompson 2012; Gatto et al. 2015; Hennebelle & Iffrig 2014; Dale 2015; Geen et al. 2017; Gavagnin et al. 2017).

The global regulation of star formation is one of the greatest challenges in modelling galaxy formation and evolution due to the difficulty in simulating feedback-regulated star formation, and the absence of a detailed physical description of processes responsible for large-scale feedback [whether driven by active galactic nuclei (AGN) or intense star formation or both] and gas accretion onto galaxies. Recently, there is growing observational and theoretical evidence that turbulent pressure injected by young stars is comparable to gravitational pressure in distant disks, enabling self-regulated star formation with low efficiency (Lehnert et al. 2013). Warm diffuse and cold dense gas evolve under the influence of compression by passing spiral arms, thermal and gravitational instabilities, and supernovae-driven shocks. High resolution hydrodynamic simulations, using sub-grid turbulence models (Schmidt et al. 2013; Semenov et al. 2016), show complex, multiphase, turbulent structures within the ISM with a realistic global Kennicutt–Schmidt relation on kiloparsec scales and gas depletion times in star-forming regions over scales of 10–50 pc consistent with observations. In these models, the global gas depletion time is long, i.e. τ_depl = M_g/Ṁ_⋆ ≈ 1−10 Gyr, where M_g and Ṁ_⋆ are the gas mass and star formation rate, respectively. This lengthy depletion time is because the gas spends most of the time in a state that does not rapidly lead to star formation. The gas is recycled many times since the lifetime of star-forming clouds or the local gas depletion time in star-forming regions is typically 1 − 500 Myr. This wide range of local depletion times, which lead to significant gas cycling and recycling is due to dynamical disruption, dispersal by feedback (Raskutti et al. 2016; Semenov et al. 2017), and supersonic turbulence (Guillard et al. 2009). This complexity must be captured in some way in galaxy evolution models to generate realistic galaxies.

In high resolution hydrodynamical simulations (e.g. Hopkins et al. 2014; Schaye et al. 2015; Kimm et al. 2017; Mitchell et al. 2018) or in semi-analytical models (SAMs; 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), depending on the model and mass of the galaxy, different prescriptions for various feedback mechanisms have been invoked to regulate star formation. For low stellar masses, M_⋆ <  10⁹ M_⊙, feedback is used to either regulate gas accretion, mostly by heating the gas through photoionisation (e.g. Doroshkevich et al. 1967; Couchman & Rees 1986; Ikeuchi 1986; Rees 1986), or ejecting the gas using the mechanical energy generated by supernovae (SNe; e.g. White & Rees 1978). For high stellar masses, M_⋆ >  10¹⁰ M_⊙, models rely on the action of supermassive black holes (SMBH) to inhibit gas accretion onto galaxies. A significant part of the power generated by AGN is used to limit the cooling of the hot gas phase surrounding the galaxy.

Despite the inclusion of some of these processes to regulate the gas content of galaxies, galaxy evolution models fail to reproduce the star formation histories and physical properties of galaxies, mainly because a robust theory of star formation including a dynamical coupling between gas phases is still lacking. As explained in Cousin et al. 2019, A&A, 627, A132 (2015a), current semi-analytical models overestimate the number of low-mass galaxies, especially at high redshift (z >  2.0), where the gap between models and observations is roughly an order of magnitude or more. For high-mass galaxies, AGN feedback, initially used to limit the growth of massive galaxies at low redshift, has also an impact on the star formation rate and history at higher redshift. Consequently, even the massive galaxies, those with M_⋆ ≃ 10¹¹ M_⊙ observed at redshifts greater than 3, are not robustly reproduced in current models.

We present a new semi-analytical model G.A.S. – the Galaxy Assembler from dark-matter Simulation – which is based, in part, on a previous version described in Cousin et al. 2019, A&A, 627, A132 (2015b, 2016). This paper (Paper I) provides an overview of the physical processes considered in G.A.S. and how they are implemented as phenomenological rate equations. In addition, we have two complementary companion papers: Cousin et al. 2019, A&A, 627, A132 (2019) G.A.S. II, in which we describe and model the mechanisms of dust attenuation and Cousin (2019) G.A.S. III, in which we explore the panchromatic emission of galaxies from the far-ultraviolet to the submillimetre bands.

In this first paper we focus mainly on the regulation of star formation. In previous models, we have adopted an ad hoc recipe to generate a delay between accretion and star formation (Cousin et al. 2019, A&A, 627, A132 2015a). We implement a physical prescription based on the inertial cascade of turbulent energy from large to small scales. Accreted gas onto galaxies is initially considered as mainly diffuse. We compute the mass fraction of the gas subject to phase separation and fragmentation following Sharma et al. (2012) and Cornuault et al. (2018) and references therein, which also depends on the disk properties. Star formation occurs in the fragmented gas at a scale of 0.1 pc. In a large set of semi-analytical models (e.g. Croton et al. 2006; Cattaneo et al. 2006; Somerville et al. 2008), star formation in massive galaxies is regulated by a strong reduction, or even complete suppression, of gas accretion by AGN feedback. In our new model, we do not need efficient AGN feedback, but instead our regulation process is a natural outcome of both the growth of thermal instabilities in the hot halo phase and the dissipation of turbulent energy within the denser, fragmented gas reservoir. Those two processes delay gas accretion onto galactic disks and star formation.

The paper is organised as follows. In Sect. 2, we provide a brief description of the dark-matter simulation we use, as well as the prescription we adopt to implement baryonic accretion rates. In Sect. 3, we focus on the turbulent inertial cascade. We describe how we compute the energy and mass transfer rates between physical scales; we define and compute the gas fragmentation timescale and show that it is a key parameter for the regulation of star formation in low-mass galaxies. In Sect. 4, we describe our model for the gas cycle as it accretes onto galactic disks, from diffuse accreted gas to potentially star-forming gas. In Sect. 5, we describe our implementation of SN and AGN feedback. Section 6 focusses on the thermal instabilities arising in the hot gas phase in the halo. We assume that gas accretion onto galaxies is limited by turbulent mixing in the range of radii, where thermal instabilities develop because gas acquires high velocity dispersions. This allows us to define an effective cooling rate. In Sect. 7 we present and discuss our results, mainly focussing on the evolution of the galaxy stellar mass function with redshift, and the evolution of relevant timescales of physical processes (e.g. gas cooling, fragmentation, and orbital timescales). We also discuss the impact of our implementation of thermal instabilities on the quenching of massive galaxies in massive haloes.

2. From dark matter to baryons

2.1. Dark matter

The G.A.S. model is built upon a set of dark-matter merger trees extracted from a pure N-body simulation. The current simulation uses WMAP-5 yr cosmology (Ω_m = 0.28, Ω_Λ = 0.72, f_b = 0.16, h = 0.70, Komatsu et al. 2009) with a volume of [100 h]³ Mpc in which 1024³ particles evolve. Each particle has a mass of m_p = 1.025 × 10⁸ M_⊙. Haloes and substructures (satellites) are identified via the HaloMaker code (Tweed et al. 2009). In our merger trees, we only consider haloes with at least 20 dark-matter particles leading to a minimal dark-matter halo mass of 2.050 × 10⁹ M_⊙. Dark-matter haloes grow from smooth accretion. The dark-matter accretion rate, Ṁ_dm, only includes particles that are newly detected in the halo and have never been previously identified in another halo.

2.2. Baryonic accretion

Based on the dark-matter accretion rate Ṁ_dm, the accretion rate of baryons is defined by

(1)

where is the effective baryonic fraction. The baryon accretion rate depends on the gas ionisation state and we adopt in this work the photoionisation model based on the Gnedin (2000) and Kravtsov et al. (2004) prescription, but we use the effective filtering mass given by Okamoto et al. (2008). We assume that re-ionisation occurs at z = 7.0, we therefore set .

Figure 1 shows the normalised mass baryonic fraction that is associated with the dark-matter smooth accretion; we assume a universal baryonic fraction f_b = 0.16. Following the minimal dark-matter halo mass used in our model, the main impact () of the photoionisation prescription occurs at low redshift (z <  0.9).

Fig. 1. Relative mass fraction of baryons which are accreted smoothly. The fraction is plotted as a function of both the dark-matter virial mass and redshift. The redshift of each curve is colour coded by the colour bar on the right. Each curve, starting at z = 0.0, is separated by Δz = 0.1. For z >  zreion = 7.0, we set . The grey vertical bar indicates the minimum halo mass in our dark-matter model. The grey horizontal line indicates f = 0.5.

As initially proposed by Khochfar & Silk (2009), in G.A.S. we define two different modes of accretion: a cold¹ mode, and a hot mode. Depending on the dark-matter halo mass, the fraction of accreted, hot gas is computed as

(2)

where M_mix is the transition mass when the cold and hot gas mass accretion rates are equal and M_vir is the halo virial mass. The evolution of the hot gas fraction is inspired by the study of Lu et al. (2011; see their Eqs. (24) and (25)), but we do not account for evolution with redshift since it is very weak and does not strongly impact how the gas is accreted in our model.

The baryonic accretion is divided into two parts that feed two different reservoirs: M_cold for the cold mode and M_hot for the hot mode. Both the cold and hot reservoirs are fed by metal-free gas. During the evolution of any galaxy, metal-rich ejecta coming from the galaxy are added to the hot reservoir. The metal content of the hot gas therefore depends directly on the rates and timescales over which galaxies create metals. The metallicity of the hot reservoir evolves with time. The chemo-dynamical model included in G.A.S. tracks the abundance of the main elements in the gas phase. The production and the re-injection of these metals are taken into account for stars with initial masses between 0.1 M_⊙ and 100 M_⊙ over metallicities from zero to super-solar.

2.2.1. Cold accretion mode

For the cold accretion mode, we assume that gas falls directly onto the galaxy (Dekel & Birnboim 2006), and we compute the cold accretion rate via

(3)

where M_cold is the mass stored in the cold reservoir and t_dyn is the dynamical time of the dark-matter halo: t_dyn = r_vir/v_vir.

2.2.2. Hot accretion: Radiative cooling

As in the previous versions of this model (Cousin et al. 2019, A&A, 627, A132 2015b, 2016), we assume that the hot gas surrounding the galaxy is confined in the potential well of the dark-matter halo and in hydrostatic equilibrium. The hot gas density profile ρ_hot(r) is computed following the prescriptions in Suto et al. (1998), Makino et al. (1998), Komatsu & Seljak (2001) and Capelo et al. (2012). We assume for such hot atmospheres a constant temperature T_hot, and a average gas metallicity Z_hot. Details are provided in Cousin et al. 2019, A&A, 627, A132 (2015b, Sect. 4). Radiative cooling and the associated gas condensation is computed using the prescription in White & Frenk (1991). The cooling time of the hot gas phase is defined – as a function of the radius – as

(4)

In previous versions of this model, we adopted the Sutherland & Dopita (1993) cooling efficiencies Λ(T, Z). In the present version, we use those computed by De Rijcke et al. (2013), tabulated as a function of gas temperature and metallicity, which we interpolate between T = 10³ K and 10⁸ K, and for gas metallicities between the range 10⁻⁴ Z_⊙ and 2 Z_⊙. We assume the solar metal mass fraction of Z_⊙ = 0.02. Figure 2 shows the cooling efficiency as a function of both gas temperature and average gas metallicity.

Fig. 2. Cooling efficiency, Λ(Thot, Zhot), as a function of both gas temperature and gas metallicity (see De Rijcke et al. 2013). The values of Λ(Thot, Zhot) are colour coded using the colour bar on the right (in units of erg s−1 cm−3).

At a given time, the mass of warm gas that can condensate and feed the galaxy is enclosed within the cooling radius, r_cool. This radius is calculated using the cooling time equation

(5)

In this equation, is a “cooling clock”, which is a measure of the effective cooling time of the hot gas. After each time step, Δt, the mass-weighted cooling time is updated via

(6)

where M is the total mass in the hot phase after the latest time step, Δt. The quantity ΔM is the net mass variation of hot gas during this time step (accretion–ejection). We assume that cooling takes place during the overall last time step for the gas already in the hot atmosphere. However, accretion is continuous², and incoming gas also starts to cool. Therefore, taking into account the time for the incoming gas to enter the hot phase, this new gas cools during only half a time step on average. Therefore, the effective cooling time of the hot gas halo can increase or decrease between two time steps depending on the relative fraction of halo gas to incoming gas. During mergers, the cooling clock of the remnant hot phase is set to the value of the most massive progenitor at the time of the merger (see Table 1).

Knowing the cooling radius r_cool (deduced from Eq. (5)), we can write the condensation rate of the gas,

(7)

The mass within r_cool decreases with a timescale r_cool/v(r_cool), where v(r_cool) is the circular velocity of the dark-matter halo measured at r = r_cool. We assume that the hot atmosphere extends up to the virial radius r_vir of the dark-matter halo, thus the cooling radius cannot be larger than r_vir.

Figure 3 shows the average galaxy gas accretion rates produced by the two different modes as a function of both the dark-matter virial mass and the redshift. At low dark-matter halo masses, M_vir <  10^10.5 M_⊙, accretion on the galaxy is dominated by the cold mode. The contribution to the accreted mass of the hot mode increases progressively with halo mass. Around M_vir ∼ 10^10.5 M_⊙, the contribution of the two modes are equal. This transition occurs at approximately the same halo mass for all redshifts we considered. For both the cold mode and the hot mode accretion, the average accretion rate decreases with the redshift. For dark-matter haloes with M_h = 10^10.5 M_⊙, the average accretion due to the cold mode decreases from 20 M_⊙ yr⁻¹ at z ≃ 9.0 to 0.2 M_⊙ yr⁻¹ at z = 0.3; for the hot mode, the accretion rate decreases from 30 M_⊙ yr⁻¹ to less than 0.3 M_⊙ yr⁻¹ between z ≃ 9.0 and z = 0.3.

Fig. 3. Mean accretion rates into the galaxy halo galaxy from the cold (upper panel) and hot mode accretion (lower panel) as a function of virial mass. The mean values are computed by selecting only haloes which are actively accreting gas. Colour code indicates redshift. The grey solid vertical line inidcates the dark-matter halo mass for which contributions of cold and hot modes are equivalent.

3. Turbulent inertial cascade

Both observations and numerical simulations of the ISM at high spatial resolutions show that turbulence plays a fundamental role in star formation (e.g. Bergin & Tafalla 2007; Miville-Deschênes et al. 2010; Renaud et al. 2014). Turbulence controls the rate at which kinetic energy is dissipated (e.g. Krumholz & McKee 2005; Padoan & Nordlund 2011; Federrath & Klessen 2012) and leads to the development of multiphase morphological structures in the gas (e.g. André 2013; Levrier et al. 2018). Recent numerical models of galaxy formation have adopted a gravo-turbulent sub-grid model for star formation (e.g. Hopkins 2012; Kimm et al. 2017), but those time-consuming simulations are limited to low cosmological volumes. In this section, we develop our modelling of the mass and energy transport from the large scale of injection to the small dissipative scale, through a hierarchy of structures, using a phenomenological prescription for the cascade of turbulent energy. This allows us to compute the mass of gas that may form stars and to test the impact of these phenomena on the properties of galaxies. In our prescription, the mass of gas that is able to form stars is the mass of gas reaching the dissipation scale l_⋆ (see Sect. 3.1), calculated from the mass flow rate between scales, under the assumption of a constant energy transfer (Sect. 3.2).

3.1. Self-similar scaling relations and energy transfer rate

We assume that the two Larson (1981) self-similar scaling relations are satisfied to compute the mass and energy transfer rates during the inertial turbulent cascade. The two scaling relations link the mass surface density μ and the 1D-velocity dispersion σ, respectively, to the wave number k = 1/l as follows:

(8)

where k_⋆ is the wave number associated with the star formation scale, l_⋆ = 1/k_⋆. μ_⋆ and σ_⋆ are two normalisation parameters and a and b are the two slope indices; both depend on the gas phase considered (Fleck 1996; Hoffmann & Romeo 2012). The second self-similar scaling relation gives the energy transfer rate per unit volume, as in Kritsuk et al. (2013),

(9)

We set the energy transfer rate ė to be constant (e.g. Hennebelle & Falgarone 2012), which gives a relation between the two slopes of Larson’s relations as follows: . Following Romeo et al. (2010), we assume (a, b)=(1/5, 3/5). These values reproduce observations of the dynamics of clouds at the atomic-molecular transition. The energy transfer rate during the inertial cascade is now proportional to

(10)

The normalisation parameters used in Eq. (8) are listed in Table 2 and correspond to standard values for dense, supersonic, and compressible gases (Jog & Solomon 1984a,b; Romeo et al. 2010; Hoffmann & Romeo 2012). We assume l_⋆ = 0.1 pc for the star formation scale, which corresponds to the characteristic width of interstellar filaments hosting prestellar cores (e.g. André 2013; Palmeirim et al. 2013). Although the detailed physical interpretation of this width is still debated, this scale is of the order of the scale below which the turbulence becomes subsonic in star-forming filaments (e.g. Padoan et al. 2001). At this scale, we adopt a typical velocity dispersion σ_⋆ = 0.3 km s⁻¹ (e.g. Orkisz et al. 2017), a value slightly higher than the speed of sound for a molecular gas at 10 K (c_s(10 K)=0.2 km s⁻¹), which corresponds to the observed transition between bound and unbound filaments (Arzoumanian et al. 2013). The critical mass surface density above which the gas is gravitationally unstable and converted into stars is set to μ_⋆ ≃ 150 M_⊙ pc⁻² (Gao & Solomon 2004; Heiderman et al. 2010; Lada et al. 2012).

3.2. Mass transfer rate during the inertial cascade

In star-forming galactic disks, the structure of the gas is observed to be self-similar over a wide range of scales, from a few kiloparsec to the length scale l_⋆ over which star formation occurs (Dickman et al. 1990; Miville-Deschênes et al. 2010). In this section, we compute the mass transfer rate from the energy transfer rate (Eq. (10)). Starting from the scale of the disk scale height, h_d, which is the largest possible injection scale, the gas mass is progressively distributed over smaller scales. Following the inertial cascade, when going from a scale 1/k to 1/2k, large structures break into smaller ones. A diagram of this progressive fragmentation of the gas in the ISM of our modelled galaxies is shown in Fig. 4.

Fig. 4. Diagram illustrating how the gas fragments during the turbulent inertial cascade. When a structure formed at the scale k reaches the Bonnor–Ebert mass, it breaks into smaller structures at the next level, 2k. The mass transfer rate, Ṁk, between scales is given by Eq. (12). The largest structures are fed at a rate, Ṁfrag, as discussed in Sect. 4.1 and given by Eq. (14).

The Bonnor–Ebert mass at a given scale 1/k, , sets the critical mass above which the structure becomes unstable and collapses into smaller structures (Bonnor 1956), until we reach the dissipation scale l_⋆, i.e.

(11)

We derive the mass transfer rate Ṁ_k between wave numbers k and 2k, using the conservation of energy

(12)

where is the volume of the cloud at scale k, M^k is the total mass stored at scale k, and ė is the constant kinetic turbulent energy transfer rate per unit volume (Eq. (10)).

3.3. Gas fragmentation timescale

We define the gas fragmentation timescale, T_frag, as the time needed to transfer all the gas from the disk scale height, h_d, to the star formation length scale, l_⋆, assuming that the fragmented gas reservoir is fed at a constant rate Ṁ_frag. During the process of fragmentation, we assume that the disk scale height is constant. In practice, we track the mass of the initial Bonnor-Ebert sphere along the cascade given by Eq. (12), until a steady state is reached.

Figure 5 shows our estimated gas fragmentation timescale T_frag[h_d, Ṁ_frag] as a function of the instantaneous disk scale height h_d and mass flow rate at which the largest structure is fed, Ṁ_frag (see Sect. 4.1 and Eq. (14)). Both the average disk scale height and the giant molecular clouds (GMCs) growth rate increase with the stellar mass, and the average fragmentation timescale decreases. Obviously, gas accretion history, star formation, gas ejection modify those two parameters, h_d and Ṁ_frag. To take this into account, we define an effective disk fragmentation timescale 𝒯_frag, which depends on the history of the disk (see Sect. 4.2 and Eq. (19)).

Fig. 5. Gas fragmentation rate as a function of both the maximum energy injection scale, lmax, which we assume is equal to the disk scale-height, and the growth rate of GMCs (Eq. (14)). Each coloured box shown in the diagram represents the gas fragmentation rate for a specific stellar mass bin as indicated by the labels at the lower right corner of each box. The mass bins range from M⋆ = 109.0 M⊙ to 1011.5 M⊙. The limits of each box encompass the 15% and 85% percentiles of the distribution for each mass bin. Coloured points within each box indicate the value of the median of the distribution. The regions shown are for modelled galaxies at z = 2.1.

4. G.A.S. cycle: Evolution of the gas reservoirs

To model the gas cycle in galactic disks, we follow the mass content of three gas reservoirs: (i) a diffuse gas reservoir; (ii) a fragmented gas reservoir; and (iii) a star-forming gas reservoir. In the following, we describe each of these gas reservoirs and their mass flow rates, of which we give a schematic view in Fig. 6. We consider that the accreted gas onto the galaxy is multiphase and we set the mass fraction of fragmented gas to (van de Voort & Schaye 2012). We further discuss the physical origin and consequences of that assumption in Sect. 8.

Fig. 6. Diagram illustrating the flow of mass between the three gas reservoirs considered in our multiphase G.A.S. model. The star-forming reservoir contains the mass that is immediately available for star formation. It is fed by a turbulent cascade and gas fragmentation at a mass flow rate, Ṁsfg (Eq. (18)). The gas reservoir fragments at a rate, Ṁfrag (Eq. (14)). The diffuse gas reservoir is fed through the rates of two mechanisms: (1) gas accretion rate (Eq. (1)); and (2) the rate at which the fragmented and star-forming gas is disrupted via the energy injected by SNe and AGN, Ṁdisrupt (Eqs. (27) and (28)). Each of these three gas reservoirs contributes to the outflow rate, Ṁwind (Eqs. (24) and (26)).

4.1. Diffuse gas reservoir

The first main reservoir stores a diffuse (≃1 cm³) and warm (10⁴ K) gas which is traced by emission lines of the warm, ionised medium such as [O III] and [N II] (e.g. Béthermin et al. 2016; Laporte et al. 2017; Suzuki et al. 2017). As illustrated in Fig. 6, the diffuse gas reservoir is fed by two different sources: (i) the accretion of warm gas coming from both the cold and the hot mode, and (ii) the disruption of fragmented gas by the injection of energy from a SN or an actively accreting SMBH (see Sect. 5.4.3).

The initial temperature of the warm, diffuse gas is assumed to be 10⁴ K, and we compute its isobaric cooling using the same approach as with the hot gas phase (Sect. 2.2.2). The effective cooling time of this phase is computed after each time step Δt, as the sum of the mass-weighted cooling times of the halo and the newly incoming gas, i.e.

(13)

where M is the total mass of diffuse gas after the time step Δt. The quantity ΔM is the mass increase of diffuse gas coming from accretion and disrupted gas coming from the fragmented phase. We assume that radiative cooling acts on the warm diffuse gas during all the previous time step (Eq. (13)). However, for the newly incoming gas, we assume that the radiative cooling occurs only during half of the previous time step (as for the hot gas phase). Therefore, the fraction of halo and freshly acquired gas can increase or decrease after each step.

The mass transfer rate between the diffuse gas phase to the fragmented gas phase Ṁ_frag is computed using the cooling timescale

(14)

where M_diff is the mass of diffuse gas. The cooling timescale is computed as a function of the average metallicity Z_diff (Eq. (4)). We assume a temperature of 10⁴ K and an average volume of the warm, diffuse gas component³ . The other parameters in Eq. (14) are computed as follows:

• ϕ_m, the mass fraction of diffuse gas that condenses in an effective cooling time, , depends on the characteristic cooling timescale of this reservoir . The parameter ϕ_m was calculated in Cornuault et al. (2018) and for computational purposes we fitted their computation by an error function given by,

  (15)

  where . The best fit gives x_t = 0.55 and s = 0.13.

• f_frag, the fraction of the gas that is fragmented, is defined as

  (16)

  where M_diff, M_frag and M_sfg are the gas masses stored in the three different reservoirs (Fig. 6). This factor accounts for the fact that the more the gas is fragmented, the lower the mass flow to bound structures.

• f_Q, the Toomre disk instability criterion, is calculated as,

  (17)

  This factor accounts for the fact that the mass flow rate to bound structures increases as the diffuse gas becomes gravitationally unstable. We adopt a standard value of Q_crit = 1.0.

4.2. Fragmented gas reservoir

The clumpy gas phase has low filling factor (typically less than ≈10%) with densities n_H = 1 − 10³ cm⁻³ and temperatures 10³–10 K. It is traced by atomic and molecular gas lines such as CO, [C I], and [C II] (e.g. Aravena et al. 2016; Bothwell et al. 2016, 2017; Decarli et al. 2016; Popping et al. 2017). We assume that the fragmented gas is contained within spherical and bound structures with initial radii r = h_d/2, and must contain at least a mass to be unstable. At each time step, the maximum number of bound structures formed in the disk is .

The gas fragmentation process feeds the star-forming gas reservoir at the rate

(18)

The effective disk fragmentation timescale 𝒯_frag is updated at each time step and allows us to follow the full history of the gas cycle in relation with the current disk properties (h_d and Ṁ_frag). If M is the total mass of fragmented gas after the last time step Δt, and ΔM = ΔtṀ_frag is the gas mass incorporated in the fragmented gas reservoir, thus the disk fragmentation timescale is updated as

(19)

This timescale is a function of T_frag[h_d, Ṁ_frag], which depends on values of h_d and Ṁ_frag, the disk scale height and the growth rate of GMCs, respectively (see Sect. 3.3).

4.3. Star-forming gas reservoir

Progressively, the fragmented gas is converted into star-forming gas. The gas contained in the star-forming reservoir is very dense and cold, typically traced by molecular lines as HCN (e.g. Oteo et al. 2017). The star-forming gas reservoir is characterised by its very short timescale before it forms stars. We assume a constant star formation timescale which is linearly dependent on the length scale over which star formation occurs, t_sf = l_⋆/σ_⋆ = 0.1 Myr. The star formation rate is then simply given by

(20)

In our prescription, the rate of star formation is mainly limited by the rate at which gas becomes clumpy.

4.4. Schmidt–Kennicutt laws

The G.A.S. model follows the evolution of three interacting gas reservoirs. Star formation is assumed to occur in the reservoir containing the densest gas after a progressive structuring starting at disk scale height in some GMCs. We can estimate the Schmidt–Kennicutt laws at two different scales within this context. At the galaxy scale, the gas surface density Σ_gas is computed assuming that half of the mass of the fragmented gas is enclosed in the half mass radius (1.68r_d) of the disk. The star formation rate surface density Σ_SFR is computed in the same way. Using these two definitions, our predictions of the galaxy-scale Schmidt–Kennicutt law is shown in Fig. 7.

Fig. 7. Schmidt–Kennicutt relations estimated from our model galaxies at both the galaxy scale and GMC scale. In the three panels the solid blue line indicates the Kennicutt (1998) relation. The grey shaded area is the distribution of a full sample of star-forming galaxies in the redshift range z = 1.5–2.1. Each coloured point indicates the position of an individual star-forming galaxy, uniformly selected within the galaxy-scale Σgas − ΣSFR plan. The colour bar indicates the fragmented gas fraction in each representative galaxy. Left panel: the solid orange line is the best fit computed for our entire sample of star-forming galaxies within this redshift range. The dashed contour indicates measurements of nearby galaxies (Bigiel et al. 2008). The blue stars and points are from the starburst and normal disk galaxy samples from Kennicutt (1998), respectively. Central panel: the solid black line indicates the median trend of our star-forming galaxy sample. The dashed black lines indicate the 15% and 85% percentiles ranges of our model galaxies. We estimate an average scatter in the distribution of 0.46 dex. Right panel: distribution of galaxies in the Σgas − ΣSFR estimated over the GMC-scale of star formation (see text for details). We show this estimate for the same galaxy sample used in the galaxy-scale Schmidt–Kennicutt relation in the two leftmost plots in this figure.

The median trend of our star-forming galaxy sample at z = 2.1 is in good agreement with the Kennicutt (1998) relation over four orders of magnitude (Fig. 7). Below log₁₀Σ_gas = 1.0, our star formation rates are slightly higher than those deduced from Kennicutt (1998). We measure an average scatter of ≃0.46 dex. Our star-forming galaxy distribution is fully consistent with individual measurements galaxies measurements (Kennicutt 1998; Bigiel et al. 2008). Using the whole star-forming sample of galaxies in the redshift range z = 1.5–2.1, we estimate a power-law index, N = 1.232 ± 0.002. This slope is slightly shallower than the standard Kennicutt (1998) slope. In addition, we find a clear trend of an increasing gas surface density as the fragmented gas fraction increases (Fig. 7).

In our model, we assume that star formation is initiated in GMCs by the progressive fragmentation of the gas. GMCs are scaled to the disk scale height h_d. By assuming that the fragmented gas mass and the star formation rate is homogeneously distributed in all GMCs formed into the disk, we define the gas surface density and the star formation rate surface density as follows:

(21)

The clear correlation we found at the galaxy scale disappears at the GMC scale (Fig. 7). Compared to galaxy scale, the fragmented gas fraction shows the opposite trend: the higher the gas surface density in a GMC, the lower its fragmented gas fraction. This trend can be translated as follows: galaxies with a relatively higher (lower) fraction of fragmented gas host relatively more (fewer) GMCs. The mass of gas and star formation is homogeneously distributed in all GMCs (Eq. (21)). The average gas surface density is therefore lower in galaxies with highly fragmented gas, which also happen to host more GMCs than galaxies with relatively low levels of fragmentation.

5. Model feedback

Massive stars and AGN in galaxies inject significant amounts of energy into the ISM and the circumgalactic medium (CGM) of galaxies. One of main challenges in galaxy evolution models is to distribute this power into the various gas phases in these media (efficiency parameters are listed in Table 1). We now describe how we distribute the SN and AGN power within the galaxy and its surroundings in the G.A.S. model, and how this power is used to regulate star formation.

5.1. Morphology and the efficiency of outflows

Galaxies in the early universe frequently have a clumpy morphology, which suggests there are interacting regions of dense gas and stars (e.g. Elmegreen et al. 2009; Elmegreen 2009). The clumps in these distant galaxies are not typically observed in the same range of mass and size as star-forming regions in local disk galaxies. This morphological evolution provides clues as to how star formation in distant galaxies may have proceeded and was regulated, but overall observations are consistent with mass and energy injected into the ISM playing an important role in regulating star formation in galaxies (self-regulation; e.g. Lehnert et al. 2015, and references therein).

To account for this morphological evolution, we assign to each galaxy a specific morphological type. Galaxies that formed recently are assumed to be clumpy. The morphological type can then change, from clumpy to a smooth(er) disk, only during a merger. If two clumpy galaxies merge, the remnant galaxy can be converted to a smooth disk with a probability of 25%. If two different morphologies merge, the remnant galaxy is assigned the morphology of the most massive progenitor (clumpy or smooth). Following these rules, clumpy galaxies become progressively regular, smoother disk galaxies. Furthermore, galaxies with these two different morphologies also behave differently as they evolve. For low-mass, clumpy galaxies, we assume that gas is more compact and therefore that ejection of gas through outflows is more efficient. The terminal velocity of a wind is proportional to the square root of the energy injection rate divided by the total mass flow rate of the wind Ṁ_wind. Thus, we can translate this efficiency in terms of the average terminal wind velocity. For clumpy galaxies, we assume the terminal velocity of the wind V_w = 100 km s⁻¹, and V_w = 200 km s⁻¹ for disk galaxies. The mass loading factor Ṁ_wind/SFR is higher in clumpy galaxies than in rotating smooth disks. At z = 2.1, distributions of the mass loading factor are characterised by the probabilities that they lie above the 25th, 50th, and 75th percentiles for all star-forming galaxies, which values are Ṁ_wind/SFR = 6.5, 11.5, and 26.4, respectively.

5.2. Supernovae feedback

In Cousin et al. 2019, A&A, 627, A132 (2016) the instantaneous SN event rate (η_SN Gyr⁻¹) is proportional to the star formation rate, which is a function of time (related to the star formation history; see Sect. 3.5.2). Based on the SN rate, the instantaneous power generated by SNe is simply given by Q_sn = η_snE_sn. We assume the standard value E_sn = 10⁴⁴ erg s⁻¹ (e.g. Aguirre et al. 2001). In G.A.S., outside of the morphological division in the efficiency of outflows, we use the same prescriptions as in Cousin et al. 2019, A&A, 627, A132 (2016). Please see that paper for details.

5.3. Active nuclei feedback

In addition to the energy input from SNe, the different gas phases are also impacted by the energy produced by AGN. We assume that SMBHs are created during a major merger if the remnant galaxy have a bulge of at least 10⁶ M_⊙. The seed of the SMBH is given an initial mass, . We associate a gas torus with each SMBH formed. The torus is the gas reservoir that feeds the growth and energy output of the SMBH. This torus is fed by diffuse gas during merger events using the following prescription:

(22)

where μ_g is the mass fraction of the gas that is diffuse in the remnant galaxy disk; is the merger mass ratio and M_i is the total mass (dark-matter + galaxy) enclosed in the half mass radius of the halo; M_diff(r <  3r_torus) is the mass of diffuse gas enclosed in r <  3r_torus. The torus radius, r_torus = 10 pc.

From Eq. (22) it is clear that major mergers with gas-rich progenitors accrete the largest amount of gas onto the torus of the remnant, while minor mergers of gas-poor progenitors lead to very little accretion. The infall rate is chosen to be the maximum value of the Bondi infall rate (Bondi 1952) and the free-fall rate given by

(23)

where the gas temperature of the torus is fixed to T_torus = 10^6.5 K. Z_torus is the gas-phase metallicity of the torus. The metallicity and temperature determine the cooling efficiency Λ(T, Z) (see Sect. 2.2.2). The values M_• and M_torus are the SMBH mass and the torus mass, respectively; t_• is the orbital time of gas at the radius of the torus r_torus.

Our prescription for the infall of gas onto the SMBH follows closely that presented in Ostriker et al. (2010). The total mass infall flux into the region of influence of the SMBH is the sum of the mass flux that is accreted onto the SMBH Ṁ_acc, and the fraction that is driven out of the region of influence of the SMBH Ṁ_ej. This yields the total infall rate Ṁ_infall = Ṁ_ej + Ṁ_acc. In our model, we assume that the relative fraction of ejected to accreted mass flow rates is . This division results in equal shares of the infalling gas to (i) be accreted onto the SMBH, driving the increase the SMBH mass; and (ii) generate power by converting a fraction of the accreted mass f_ME = 0.1 (Fig. 8). The power produced by the AGN is then Q_AGN = f_MEṀ_accc².

Fig. 8. Illustration showing the distribution of accretion and energy output of a model AGN. Gas infalling from the torus to the AGN is divided between accretion onto the SMBH Ṁacc and ejection Ṁjet. We assume the diffuse gas in the disk can be driven outwards by the AGN to give the overall ejection rate of gas Ṁwind,agn. The accretion onto the SMBH Ṁacc is itself divided in two parts: a fraction (1−fME)Ṁacc of the infalling gas goes to increasing the mass of the SMBH, while the remaining gas Qagn = fMEṀacc is converted into energy (see Sect. 5.3 for details). The radiative and mechanical energy of the AGN is then distributed in various ways as illustrated in Fig. 9.

5.4. Distribution of feedback power

We use the power output from AGN and SNe to regulate the gas cycle in galaxies in three different ways in our model. The two different sources of power – the AGN power Q_AGN and the mechanical energy of SN Q_SN – are each divided in two parts: kinetic power fraction f_k, AGN and f_k, SN, and the bolometric power (Fig. 9). A fraction of the kinetic power f_w is used to drive a large-scale wind (Sect. 5.4.1). The residual fraction 1 − f_w is used to disrupt the fragmented gas of the disk (Sect. 5.4.3) and power the turbulence of the diffuse gas. A fraction of both the AGN and starburst bolometric luminosity, f_th, AGN and f_th, SN, respectively, are used to heat the ejected gas (Sect. 5.4.4).

Fig. 9. Illustration of how the energy output from AGN and SNe are distributed. The SNe and AGN power QSN and QAGN are shown in red and yellow, respectively. Specifically, the power is distributed through four different channels: (1) as a large-scale wind, (2) in generating turbulence, (3) in heating the gas, and (4) output bolometric luminosity.

5.4.1. Large-scale ejecta

We derive the instantaneous ejection rate using the conservation of the kinetic energy released by SNe (e.g. Dekel & Silk 1986; Kauffmann et al. 1993; Efstathiou 2000), i.e.

(24)

where v_w is the average velocity of large-scale wind (Sect. 5.1) and f_esc is the fraction of mass that escapes the disk.

The ejection rate due to the energy output of the AGN is determined by the fraction of the mass infalling towards the SMBH driven out (Ṁ_ej; Sect. 5.3). We derive the velocity of the jet using conservation of energy,

(25)

The outflowing jet is coupled to the diffuse gas in the disk which reduces its velocity. To obtain the coupling efficiency between the jet and the diffuse gas, we assume that jet velocity is equal to the escape velocity of the galaxy V_esc. Thus, the kinetic energy of the jet is simply used to drive the diffuse gas out of the galaxy. This leads to an mass outflow rate given by

(26)

The total instantaneous ejection rate of the combined action of SNe and AGN is Ṁ_wind = Ṁ_wind,SN + Ṁ_wind,AGN. Each gas reservoir, M_diff, M_frag, and M_sf, contributes to the instantaneous ejection rate in proportion of its mass fraction (Fig. 6).

5.4.2. Overestimating the gas escape fraction: Re-accretion timescale

When the gas is ejected from the disk, a fraction of the gas remains in the hot CGM. The remaining fraction is ejected from the dark-matter potential. As in Cousin et al. 2019, A&A, 627, A132 (2015b, see their Sect. 4.3 and Fig. 6) to compute the escape fraction of the gas, we adopt a “ballistic” approach based on the comparison of the dark-matter escape fraction to a shifted Maxwell–Boltzmann velocity distributions. However, as shown by for example Oppenheimer & Davé (2008), the ejected gas is not only affected by gravitational forces (V_esc) but also by ram pressure of gas in the CGM. Because of this affect, our ballistic approach overestimates the escape fraction. To correct this bias, we store all the hot gas ejected from a halo in a extended circumgalactic reservoir M_ex-cir. The mass stored in this reservoir is then progressively re-accreted and added to the hot gas trapped into the dark-matter potential well.

The gas that is re-accreted out of this reservoir has been implemented in various semi-analytical models (e.g. De Lucia et al. 2004; Somerville et al. 2008; Guo et al. 2011; Henriques et al. 2013) and potentially plays a major role in contributing to the overall gas supply of galaxies. The formulation in the various semi-analytical models can vary. For example, Guo et al. (2011) adopted a prescription depending on both of the dark-matter halo mass and redshift. In Henriques et al. (2013) the re-accretion of gas is inversely proportional to the dark-matter virial mass without any dependence on redshift. In G.A.S., we adopt a prescription similar to De Lucia et al. (2004) and Somerville et al. (2008). We assume that the gas is re-accreted on a timescale which is twice the halo crossing time t_reacc = 2r_vir/v_vir.

5.4.3. Disruption rate

The power generated by SNe that does not contribute to driving a wind is assumed to be injected directly into the ISM. This remaining power is distributed between the different ISM gas phases in proportion to their mass. Massive stars are assumed to remain mostly embedded in their dense birth clouds during their short lifetimes (3 × 10⁶ yr; Cousin et al. 2019, A&A, 627, A132 2016). The fraction of the SN energy which is injected into the fragmented and star-forming gas f_frag disrupts this phase, feeding the diffuse gas reservoir (Eq. (16) and Fig. 6). The disruption rate is defined as

(27)

where σ_v is the average velocity dispersion of the diffuse gas (Sect. 5.5.1). We assume that the velocity dispersion of the disrupted gas is equal to that of the diffuse gas. Assuming this implies that if the diffuse gas is highly turbulent (high σ_v) then it takes more energy to disrupt the fragmented gas. The disruption rate depends on two terms specifically related to the gas phase. The first term, 1 − f_w, corresponds to the minimum fraction of SN power which disrupts the fragmented gas. The second term, f_w(1 − f_esc), corresponds to the fraction of power that remains in the gas because of the limit imposed by the galaxy escape fraction. The fraction of power which does not contribute to the wind is therefore re-injected to disrupt the fragmented gas. This fraction increases with the galaxy mass. The fraction of the power of SNe injected into the diffuse gas phase maintains or increases the level of turbulence of the diffuse gas (i.e. σ_v). During a time step Δt, the possible increase of the turbulent energy of the diffuse gas is given by (see Sect. 5.5.1 for all other contributions).

Simultaneously, we also include the energy output from AGN in disrupting the gas. The contribution from any AGN to the disruption rate is given by

(28)

As for SNe, the residual power 1 − f_frag is injected into the diffuse gas as turbulent energy, . The total mass disruption rate is the sum of the SN and AGN contributions, i.e. Ṁ_disrupt = Ṁ_disrupt,SN + Ṁ_disrupt,AGN.

5.4.4. Radiative heating and bolometric luminosity

In the previous two sections we discussed our prescriptions for the kinetic power of SNe and AGN. The non-kinetic fraction of the total power is also divided in two different parts. A fraction f_th, SN (f_th, AGN) of the non-kinetic SN (AGN) power is used to heat the ejected gas. These fractions are adjusted so that the average temperature of the ejected gas is between 10⁶ and 10⁷ K. The residual power (1 − f_k, SN)(1 − f_th, SN) for SNe is then assumed to be emitted as bolometric luminosity in each galaxy (i.e. the part of the total radiative power that does not go directly into heating the gas).

5.5. Average velocity dispersion and disk scale height

5.5.1. Average velocity dispersion

The gas disruption rates, as given in Eqs. (27) and (28), depend on the average velocity dispersion of the diffuse gas, σ_v. We compute and continuously update the velocity dispersion by taking into account simultaneously the kinetic energy injected by gas accretion, SNe, and AGN.

Our prescription is based on the evolution of the “turbulent” kinetic energy budget E_σv as follows:

(29)

where σ_v is the 3D gas velocity dispersion. Between two time steps, we assume that

• A fraction, f_disp = 1/2, of the turbulent energy is dissipated per orbital time.

• The turbulent energy is increased by , corresponding to a fractional increase, f_incr = 1/3, per dynamical time, of the kinetic rotational energy of the newly accreted diffuse gas, ΔM. We assume that the freshly accreted diffuse gas (ΔM) forms a thin rotating disk (e.g. Mo et al. 1998). This thin gas disk merges with the pre-existing disk hosting a mass, M_diff, of diffuse gas. The quantity v_acc is the orbital velocity of the freshly accreted diffuse gas computed at the half mass radius.

• The total budget of the turbulent energy is finally increased by the energy injected in the diffuse gas by SNe and AGN, and (see Sect. 5.4.3).

Considering all of these energy terms, after a time step Δt, we have

(30)

Based on this updated total turbulent kinetic energy, we calculate the average velocity dispersion of the diffuse gas σ_v (Eq. (29)). The two free parameters, f_incr and f_disp, have been adjusted to reproduced observed values of velocity dispersion (see Fig. 10).

Fig. 10. Probability distribution of the velocity dispersion of the diffuse gas. This distribution represents model galaxies with M⋆ >  107 M⊙ and at z = 1.1. We compare this distribution with observational measurements from Swinbank et al. (2017, dark blue histrogram) and Pelliccia et al. (2017, light blue histrogram). The solid and dashed vertical lines represent the median and mean value of the distribution of our modelled galaxies, respectively.

During a merger, the total turbulent kinetic energy of the two progenitors are added. We also add a fraction of the gravitational energy due to the interaction between the two galaxies,

(31)

where M_i is the mass (baryon + dark matter) included in the galaxy half mass radius and is the disk exponential radius (Hatton et al. 2003).

We compare our predictions of the gas velocity dispersion with recent observational measurements of the velocity dispersion of the warm ionised media from Swinbank et al. (2017) and Pelliccia et al. (2017). Objects targeted by these two programmes are distributed in redshift between z ≃ 0.3 and z ≃ 1.7. Our predictions are in good agreement with these observational results but our results reproduce some of the most extreme dispersions observed (both low and high values; Fig. 10). At z = 1.1, the median value of the velocity distribution of the diffuse gas is close to 35 km s⁻¹.

5.5.2. Disk scale height

The disk fragmenting timescale (Eq. (19)) depends on the average disk scale height. We use the disk scale height to define the initial energy injection scale of the inertial turbulent cascade. Our gas fragmenting scenario is mainly based on scaling relations and the disk scale height is crucial in calculating this parameter of the turbulent cascade. The disk scale height h_d is deduced from . At z = 1.1, the median value of the velocity dispersion leads to a median value of the disk scale height is ≈160 pc (Fig. 10).

6. Thermal instabilities

For galaxies with high masses, previous semi-analytical models invoked powerful AGN feedback to reduce greatly or even completely quench accretion of cooling gas in their haloes (e.g. Cattaneo et al. 2006; Somerville et al. 2008; Guo et al. 2011; Benson 2012; Henriques et al. 2013). Such prescriptions assume a strong coupling between AGN power, the ISM, and the hot circumgalactic gas, which results in powerful outflows expelling the gas and heating the halo gas. However, observations are not entirely clear concerning the precise impact of outflows on distant galaxies (see e.g. Mullaney et al. 2015; Netzer et al. 2014, 2016; Scholtz et al. 2018; Falkendal et al. 2019). In Cousin et al. 2019, A&A, 627, A132 (2015a,b), we also assume that AGN can have an impact on the hot halo phase but we have limited the effect to simply heating the halo gas.

In the following we present a new prescription to reduce efficiently the gas accretion onto massive galaxies, M_⋆ >  10¹¹ M_⊙. This new mechanism is based on the growth of thermal instabilities in the hot halo phase surrounding the galaxy.

6.1. Thermal instabilities: Description

Our new mechanism assumes that the condensation of gas and therefore the gas accretion onto the galaxy is progressively and strongly limited by the growth of thermal instabilities in the hot gaseous CGM. In the standard model of cooling presented previously and also used in, for example Croton et al. (2006), Baugh (2006), and Somerville et al. (2008), there is a strong dichotomy between the gas stored within and beyond the cooling radius, r_cool. At radii greater than the cooling radius, the gas is assumed to remain hot, while below r_cool, the gas is assumed to be warm, 10⁴ K, and dense enough to condense and feed the galaxy. Of course, the reality is more complex, and obviously hot and warm gas co-exist around the cooling radius. In this transition region, the gas can be thermally unstable. These thermal instabilities can generate warm clouds that are orbiting in the hotter gas. Cornuault et al. (2018), in a phenomenological model of accreting gas, have found that the cloud–cloud velocity dispersion can be comparable to the characteristic virial velocity of the dark-matter halo in which they formed. The mixing and dynamics of warm clouds into the hot atmosphere can moderate the effective accretion rate of gas onto a galaxy.

6.1.1. Thermal instability clock

As shown by Sharma et al. (2012), the timescale over which gas becomes thermally unstable is related to the cooling time. Following their approach, and similar to gas cooling generally, we define a thermal instability timescale, 𝒯_TI. This timescale evolves concomitantly with the cooling timescale. After each time step, Δt, the thermal instability clock is updated following a mass-weighted prescription, similar to that applied to the cooling clock, namely,

(32)

where M is the total mass stored in the hot phase after the last time step, Δt. The quantity ΔM is the mass accreted by the hot gas phase or ejected from the galaxy. For each time step, Δt, the time increment Δt_TI is calculated as

(33)

where F(Λ) is the thermal instability function (Sharma et al. 2012), depending on the logarithmic derivative of the cooling efficiency function Λ,

(34)

We assume that the hot gas becomes thermally unstable during all the previous time steps (Eq. (32)). However, as for the effective cooling clock, we consider that thermal instabilities can evolve with the addition of freshly accreted gas only during half a time step⁴. The thermal instability function, F(Λ), determines if the hot gas is stable or unstable. Gas is assumed to be unstable if F(Λ) ≥ 0 and stable if F(Λ) < 0. When hot gas is stable, the time increment Δt_TI is set to 0. According to the pre-existing/incoming mass ratio, the thermal instability clock can increase or decrease during a given time step.

As for the cooling efficiency, the thermal instability function has been interpolated and tabulated. Figure 11 shows the thermal instability function dependence on gas temperature between 10³ K and 10⁸ K and gas metallicity between 10⁻⁴ Z_⊙ and 2 Z_⊙. As for the effective cooling clock, the thermal instability clock of the hot phase is assumed to have the value that the most massive progenitor had prior to merging.

Fig. 11. Thermal instability efficiency, F(Λ), as a function of both gas temperature and gas metallicity. The amplitudes of thermal instabilities are computed as given in Sharma et al. (2012) and the values are indicated in the colour bar on the right side of the plot. When F(Λ) < 0, the gas is thermally unstable. The solid and dashed blue contours indicate where 2 − dlnΛ/dlnT = 0 and −3.0, respectively.

6.1.2. Mixing zone

The thermal instability clock provides at any time the effective timescale of the growth of thermal instabilities. Starting at the cooling radius, we assume that thermal instabilities are propagating through the gas and ultimately reach the centre of the halo. We define the instantaneous size of the mixing zone as

(35)

We assume an adiabatic index γ = 1.4 and a mean molecular mass μ = 0.62. R is the specific ideal gas constant and T_hot is the average temperature of the hot gas. We assume that the propagation speed of these instabilities is close to the sound speed of the gas. Following this hypothesis, the value of the free parameter ε_TI is set to 0.63.

In the mixing zone, we assume that thermal instabilities lead to the formation of warm gas clouds orbiting around the galaxy. As already stated, the cloud–cloud velocity dispersion is assumed to be close to the virial velocity (Cornuault et al. 2018). In addition, even if we do not take into account explicitly such a mechanism, the warm gas clouds probably interact strongly with the large-scale wind escaping the galaxy. Through this interaction, we can reasonably assume that kinetic energy of the outflowing gas is injected into a dynamic and cloudy CGM. Thus we assumed that the condensation of warm clouds is suppressed in the mixing zone. The condensation of clouds is only efficient in the inner region of the hot atmosphere and the effective cooling rate is then computed using Eq. (7). For this effective cooling rate, the initial cooling radius calculated is substituted with the expression, r_TI = r_cool − Δr_TI.

6.2. Impact of thermal instabilities on the gas cooling

To understand and quantify the progression of thermal instabilities in the hot gas phase, we define the thermal instability volume filling factor as

(36)

The parameter ϕ_v measures the volume fraction in which thermal instabilities are developed; ϕ_v = 0 indicates that there is no impact on the gas cooling. As in the original model, all the gas enclosed in the cooling radius can condense and therefore feed the galaxy. The case ϕ_v = 1 indicates that the entire region with r <  r_cool is thermally unstable. When ϕ_v = 1, since the clouds are orbiting at the virial velocity, is the point at which the cooling flow would otherwise occur is effectively stopped.

Figure 12 shows the mean trend of ϕ_v (Eq. (36)) as a function of both the dark-matter halo mass and redshift. Depending on the redshift, the first resolved haloes, M_vir >  10⁹ M_⊙, have ϕ_v distributed between 0.3 and 0.8. At z ≥ 9.0, the average ϕ_v of the first haloes formed is close to 0.6. This value progressively decreases with the redshift. At z ≃ 4.0, the average ϕ_v reaches a minimum around 0.3. Then the average ϕ_v of first haloes detected increases as the redshift decreases. The maximum value, 0.8, is reached at low redshift. In this low-mass regime, hot atmospheres are only formed via large-scale winds coming from the galaxy. The average ϕ_v is closely linked to galaxies driving winds. Indeed, as mentioned in Sect. 6.1.1, the effective thermal instability (TI) clock runs in proportion to the freshly added to pre-existing gas mass ratio in the hot halo phase (Eq. (32)). The higher the proportion of freshly added ejecta, the slower the TI clock advances and thus the slower ϕ_v increases. Between z ≃ 9.0 and z = 4.0, the intensity of ejecta affecting the first resolved haloes, M_vir ≃ 10⁹ M_⊙, increases continuously (related to both the gas accretion and the SFR). Then the intensity progressively decreases until the lowest redshifts. We caution that the number of haloes in the low-mass bins at very high redshift, z ≃ 9, is low and thus these bins are greatly affected by statistical noise.

Fig. 12. Upper panel: mean value of the volume filling factor, ϕ, of gas that is thermally unstable (Eq. (36)). The case ϕv = 1 implies that all the region within r = rcool is thermally unstable in our model. Lower panel: fraction of galaxies that are accreting radiatively cooling halo gas as a function of halo virial mass, Mvir. Different colours indicate this fraction as a function of redshift. The value of the redshift for each line is indicated in the colour bar on the right side of the upper panel. The grey area spans the range of the mean fraction of galaxies that are accreting hot gas for each halo mass summed over redshifts 0.1–10 (Eq. (2)). For example, for a halo mass of 1010.5 M⊙, the fraction of galaxies that have hot gas in the halo which is radiatively cooling can range between 30% and ∼55% depending on redshift.

Between the redshift range, z = 4.0 to z = 1.5, the volume filling factor of the most massive haloes reaches very high ϕ_v, ≃0.95. At z <  1.0, the average volume filling factor of thermally unstable gas in massive haloes starts to decrease. In these massive haloes, some accretion still occurs. This accretion produces a new star formation and therefore some large-scale ejection events. In these massive haloes we therefore observe an increase of freshly acquired hot halo gas. This increase in freshly acquired hot halo gas results in a progressive decrease of the volume filling factor of thermally unstable gas.

Figure 12 also shows the change in the dark-matter halo mass and the evolution with redshift of the fraction of galaxies in which condensation of the hot halo gas is occurring, i.e. Ṁ_cool > 0. (Eq. (7)). At the halo detection mass threshold of our model, radiatively cooling gas is accreting over the whole mass range. The higher the redshift, the higher the fraction of galaxies that accrete radiatively cooled hot halo gas. For z >  5.0, the fraction remains high, > 0.8. Then this fraction sharply decreases and reaches values lower than 0.1 at z <  3.0. At low virial halo masses, the hot gas phase is only fed by large-scale ejecta from galactic outflows (f_hot ≤ 0.18). In low virial-mass haloes, the mass of hot halo gas is small, M_hot ≤ 10⁷ M_⊙. The cold mode is the dominant accretion mode at low masses. Even if at z >  5.0 radiative cooling produces an effective accretion in a large fraction of the galaxies, the average radiative cooling rate is less than 1 M_⊙ yr⁻¹ (Fig. 3). The fraction of dark-matter haloes hosting actively and radiatively cooling hot gas increases (or remains constant) as a function of the dark-matter virial mass. The largest population of actively cooling hot haloes is reached between ≃10^10.25 M_⊙ and ≃10¹² M_⊙. Between z = 4.0 and z = 3.0, in the most massive haloes > 10¹² M_⊙ the fraction reaches 1.0. However by z <  3.0 for the most massive haloes, we find a strong decrease in the fraction of galaxies with actively cooling haloes. This decrease is directly linked to the high ϕ_v reached in these haloes and this decrease is responsible for the reduction in the effective cooling rate estimated for the most massive haloes at z <  2 (Fig. 3).

7. Results and analysis

7.1. Stellar mass assembly and co-moving density of galaxies

We compare the predictions of the G.A.S. model for the co-moving stellar mass density as a function of redshift with a compilation of observations (Fig. 13). Overall, the stellar mass growth of the Universe predicted in our model agree with observations. These predictions are in reasonable agreement with estimates even at the highest redshifts, z >  3, although they are generally consistent with the highest total co-moving stellar mass density estimates. Furthermore, this general agreement between G.A.S. and observations is also reflected in a comparison of the co-moving number density of galaxies binned by redshift (Fig. 14). We compare our results with a variety of observed distributions covering the redshift range z ≃ 0.1–z ≃ 6 (Baldry et al. 2008; Yang et al. 2009; Caputi et al. 2011; Ilbert et al. 2013; Grazian et al. 2015; Duncan & Conselice 2015; Song et al. 2016; Davidzon et al. 2016). In addition, we compare our best G.A.S. predictions for the co-moving number density with those of Guo et al. (2011) and also those obtained with the previous ad hoc regulation prescriptions of Cousin et al. 2019, A&A, 627, A132 (2015a). These comparisons allow us to gauge the impact of including the star formation and gas-dissipation prescriptions that are in G.A.S. with models that do not include such prescriptions (Fig. 15).

Fig. 13. Redshift evolution of the co-moving stellar mass volume density. The solid orange line is the prediction of the G.A.S. model. We compare our predictions with the compilation of observational results (see Madau & Dickinson 2014, and references therein).

Fig. 14. Comparison between the predictions of G.A.S. with observational estimates of the co-moving number density of galaxies for 9 redshift bins spanning the range, z ≃ 0.1–z ≃ 6 (from upper left to lower right and as indicated at the top of each panel). We compare the predictions of the G.A.S. model (solid grey lines) with various studies: Baldry et al. (2008, red squares), Yang et al. (2009, orange circles), Caputi et al. (2011, brown diamonds), Grazian et al. (2015, chocolate squares), Duncan & Conselice (2015, chestnut-red triangles), Song et al. (2016, salmon circles), Ilbert et al. (2013, yellow dashed-line), and Davidzon et al. (2016, red dot-dashed line). For this comparison with observations, we indicate binned measurements using points and best-fit stellar mass functions as dashed lines (Schechter and double-Schechter). All stellar mass functions are for a Chabrier IMF. The best fits of Ilbert et al. (2013) and Davidzon et al. (2016) are corrected for the Eddington bias which affects the high-mass galaxy bins in particular. The grey dotted-line indicates the total gas plus stars, baryonic mass function predicted by G.A.S.The grey dashed-line indicates the stellar mass function of the previous panel for the next highest redshift bin. The tick marks placed on the left side of each panel indicate the number of galaxies used in the stellar mass function predicted by the G.A.S. model.

Fig. 15. Comparison of the stellar mass functions predicted by our best fiducial model with those predicted by Guo et al. (2011). Each panel shows a different redshift as indicated at the upper left in each panel; solid red, dashed orange, and dashed gold lines indicate the predictions of G.A.S., Guo et al. (2011), and Cousin et al. 2019, A&A, 627, A132 (2015a), respectively (as indicated in the legend of the left-most panel). We compare these models with observations (dashed grey lines; Davidzon et al. 2016). At z = 4.0, we also plot additional observational measurements (salmon circles, red upward pointing triangles, and orange squares; Song et al. 2016; Duncan & Conselice 2015; Grazian et al. 2015, respectively, and as indicated in the legend of the right-most panel)

Between z ≃ 1.0 and z ≃ 6.0, our model agrees very well with observational measurements of both the cosmic co-moving mass densities and shape of the co-moving number density of galaxies for a wide range of redshifts (cf. Figs. 13 and 14). At z <  1.0, our model systematically under-predicts the co-moving density of intermediate-mass galaxies and over-predicts the co-moving density of low-mass galaxies. We find that in the low- and intermediate-mass range, 10^8.5 ≤ M_⋆ ≤ 10^11.5 M_⊙, the co-moving density functions predicted by our model show a double power-law shape, with a steeper slope at the low-mass end of the model distribution. At the highest redshifts, z ≃ 6.0, the density function keeps this shape even for the most massive model galaxies, M_⋆ ≥ 10^11.5 M_⊙. This power-law shape is the signature of the continuous competition between the disruption and progressive fragmentation of the gas through the turbulent cascade. Some galaxies, those with stellar masses above 10^10.5 M_⊙, have already substantially formed at z ≃ 6.0. The turbulent cascade regulates the star formation in intermediate- and low-mass objects and keeps a sufficient amount of diffuse and fragmented non-star-forming gas that actively feeds star formation until the formation of the most massive observed galaxies.

As the redshift decreases, the power-law shape is broken at the high-mass end of the distribution and an exponentially declining function is progressively formed. In agreement with observational estimates, the knee of our predicted co-moving number density functions develops around 10^10.75 M_⊙ at z ≃ 5.0 and progressively shifts to 10^11.0 M_⊙ by z ≃ 0.1. The break appears when the star formation is strongly reduced or quenched, i.e. when all the fresh gas contained in the disk is consumed. At the high-mass end of the galaxy distribution, gas is accreted via the hot mode (Sect. 2.2 and Fig. 3). In, for example, Guo et al. (2011), Croton et al. (2006), or Somerville et al. (2008), the reduction of the gas accretion is because of the power of the AGN heating the halo gas and halting accretion. To obtain the necessary heating, such prescriptions need an increasing and constant production of power from AGN. In G.A.S., much of the reduction or halting of gas accretion is a result of the increase in the thermally unstable hot gas phase surrounding galaxies. As shown in Fig. 16, the impact of thermal instabilities appears to be progressive and its impact starts to become significant at z ∼ 4.0 for the most massive galaxies. The effect is clearly visible for the galaxy mass function at z = 2.1. Without such regulation, the continuous accretion of gas leads to the formation of very massive, 10^12.0 M_⊙, galaxies at z ≃ 1.5. At very high masses, the reduction of the gas accretion onto galaxies is not sufficient to quench completely the star formation. Indeed, some galaxies have stellar masses more than 10^11.5 M_⊙ are predicted in our model but are not observed (e.g. Davidzon et al. 2016; Ilbert et al. 2013; Yang et al. 2009; Baldry et al. 2008). This is due to a progressive decrease of the volume filling factor of thermal instabilities in the most massive galaxies (Fig. 12). A slight increase of the thermal instabilities propagation efficiency, ε_TI, could limit the growth of these very massive model galaxies.

Fig. 16. Impact of our main prescriptions for regulating star formation and gas accretion, photoionisation, re-accretion, and thermal instabilities on the stellar mass assembly of model galaxies (i.e. the co-moving density of model galaxies as a function of redshift). The solid red line shows our fiducial model using the full suite of prescriptions. The dot-dashed gold lines indicate our model without photoionisation. In the three right-most panels, the dotted brown lines and dashed orange lines indicate predictions of the stellar mass functions without including the prescriptions for re-accretion of gas over two halo crossing times and thermal instabilities, respectively (as indicated in the legends of the two left-most panels). These are compared to the observational estimates of the co-moving density of galaxies from Davidzon et al. (2016). The four different panels represent four different redshift bins, z = 0.3, 1.5, 2.8, and 4 (from the left-most to the right-most panel, respectively).

Compared to the ad hoc recipe () presented previously in Cousin et al. 2019, A&A, 627, A132 (2015a), the new regulation prescription based on the turbulent cascade is in better agreement with the co-moving densities of low-mass galaxies. The ad hoc recipe was clearly not appropriate for capturing the necessary ingredients to produce sufficient numbers of low-mass galaxies and especially those at high redshifts (Fig. 15). On the contrary, the stellar mass functions predicted by Guo et al. (2011) clearly over-produce the co-moving density of low- or intermediate-mass galaxies, M_⋆ <  10^10.5 M_⊙. Their prescriptions for regulating star formation do not have the correct dependencies as a function of the stellar mass. Gas is converted to stars too efficiently in low- or intermediate-mass galaxies, hence producing an excess. Furthermore, a lack of gas in more massive galaxies means that their models underestimate the number of massive galaxies.

In fact, to elucidate the role played by various physical prescriptions we now include in G.A.S., we compare four different model configurations: our best complete model and three alternative versions in which: (i) the gas re-accretion prescription is disabled; (ii) thermal instabilities are not considered; and, (iii) the photoionisation prescription is switched off. These three alternative versions allow us to determine the effective impact of those prescriptions in regulating the mass growth of the ensemble of model galaxies (Fig. 16). At the highest redshift for which we made this comparison, z = 4, we see that turning off all three of these prescriptions has little impact on the modelled stellar mass distribution function. As we decrease in redshift, we see progressively greater impact for these three processes on the stellar mass distribution. The lack of thermal instabilities in the hot gas in the halo results in overly massive galaxies by z ≈ 2 − 3, but leaves essentially unchanged the number of lower mass galaxies, i.e. those with M_⋆ ≲ 10¹¹ M_⊙. Photoionisation appears to only impact the co-moving number density of galaxies at low redshift and for galaxies with low stellar masses. At z <  1, the power-law shape in the co-moving density of low-mass galaxies, i.e. those with M_⋆ <  10^8.5 M_⊙, is broken and becomes flatter as the redshift decreases. This trend is a result of a progressive reduction in the effective gas accretion with decreasing redshift and is a direct consequence of the photoionisation prescription (Sect. 2.2, Fig. 1). The process with the largest impact on the co-moving number density of galaxies appears to be the re-accretion of gas (Sect. 5.4.2). Its impact on the galaxy co-moving density distribution starts to become evident at z ≈ 3 (Fig. 16). Without re-accretion, the model strongly under-predicts the co-moving density of almost all stellar masses but especially so for intermediate-mass galaxies.

7.2. Evolution of disk properties

We discuss with Fig. 17 the evolution with both the stellar mass and redshift of five main properties of star-forming galaxies. We present the evolution of the following disk properties:

Fig. 17. Evolution with redshift of 5 main disk properties of disk galaxies which lie along the star-forming galaxy main sequence shown for a range of mass bins. Legend in the lower right gives the stellar mass values for each curve in the panels. Upper left: gas mass fraction, fgas, as a function of redshift for 5 stellar mass bins. Upper centre: fraction of the gas which is fragmented, ffrag. Upper right: average half mass radius, r50, of the ensemble of star-forming model galaxies. Lower left: diffuse gas density, ρdiff. Lower right: ratio of the gas velocity dispersion to orbital velocity, .

• The gas mass fraction, . The value M_gas is the total mass of the three gas reservoirs combined for a given disk galaxy.

• The fragmented gas fraction, .

• The galaxy half mass radius, r₅₀.

• The diffuse gas density, ρ_diff. This density is computed based on the mass of diffuse gas, a disk with an external radius of 22 × r_d and a scale height, h_d. We also assume a mean molecular mass, μ = 0.62.

• The ratio of the gas velocity dispersion to orbital velocity, . The dispersion velocity is that of the diffuse gas. The orbital velocity is computed at 2.2r_d (Pelliccia et al. 2017).

The following analysis is based on model star-forming galaxies that have been extracted at a variety of redshifts. We define star forming as those galaxies which lie above 25% of the mean of the relation between the star formation rate and stellar mass (the “main sequence” of star-forming galaxies Schreiber et al. 2015, their Eq. (9)). Within this sample of model galaxies, we calculate median of the ensemble binned in redshift and stellar mass for each of the quantities listed above. We first focus on the gas mass fraction, f_gas. At a given stellar mass, the highest redshift galaxies have the highest gas fractions while galaxies with larger stellar masses have lower gas fractions. If we consider both the increase in mass with decreasing redshift, these tracks imply that the gas mass fraction globally decreases as galaxies grow and evolve. At low redshift, z <  0.5, the gas mass fractions of star-forming galaxies lie between 15%–50% depending on the stellar mass. In two samples of star-forming galaxies, evolving around z = 1.5, Daddi et al. (2010) measured a gas mass fraction higher than our predictions. At z = 1.5, for an average stellar mass of M_⋆ = 10^10.5 M_⊙, the gas mass fraction is distributed as [p15, p50, p85] = [0.18, 0.25, 0.36] in comparison to f_gas ≃ 0.6 measured by Daddi et al. (2010). In the same redshift slice, for an average stellar mass of M_⋆ = 10¹¹ M_⊙ we measured [p15, p50, p85] = [0.11, 0.16, 0.26] in comparison to f_gas ≃ 0.5 measured by Daddi et al. (2010).

In parallel to the progressive decrease of the gas-mass fraction, the distribution of this gas between the diffuse and the fragmented/dense gas also evolves. We first note that at a given stellar mass, the fragmented gas mass fraction is an increasing function of the redshift. Galaxies formed at higher redshift contain a higher percentage of fragmented gas than galaxies evolving at lower redshift with a similar mass. At high redshift z >  2.5, the fragmented gas mass fraction is a clear decreasing function of the stellar mass. The hierarchy in mass is less clear at lower redshifts. We note that galaxies hosting a stellar mass higher than 10¹⁰ M_⊙ stabilise their fragmented gas mass fraction at around 42%. For less massive galaxies, < 10¹⁰ M_⊙, the fragmented gas fraction strongly decreases and reaches low values distributed between 25% and 33% at z = 0.1.

The gas content clearly evolves with time: gas-rich and strongly fragmented galaxies evolve through structures dominated by stars and in which the gas is mainly diffuse. This evolution in gaseous and stellar contents is following the evolution of galaxy size: At a given stellar mass, the disk half mass radius is a decreasing function of the redshift. At redshift z >  3.0, the disk half mass radius is a decreasing function of the stellar mass. Massive galaxies appear more compact than low-mass galaxies. For all stellar mass bin analysed, the average half mass radius is strictly lower than 1 kpc. At z <  3.0, we note the opposite trend: the average half mass radius is strictly larger than 1 kpc and is an increasing function of the stellar mass.

The anti-correlated evolution of the gas content and disk size directly impacts the density of the diffuse gas, which is the main driver of the GMC feeding rate (Sect. 4.1, Eq. (14)). At a given stellar mass, the diffuse gas density appears to be an increasing function of the redshift. In galaxies formed at higher redshift, the gas is denser than in galaxies of a same stellar mass formed at lower redshift. At a z >  3, the diffuse gas density is an increasing function of the stellar mass. At z <  3, the most massive galaxies see their diffuse gas density converge to 0.2–0.3 atoms cm⁻³. The less massive galaxies (≤10¹⁰ M_⊙) see their diffuse gas density strongly decrease and reach a value of < 0.1 atoms cm⁻³. This evolution is mainly linked to the decrease of the fragmented gas-mass fraction discussed previously.

This decrease also impacts the dispersion to orbital velocities ratio, . This ratio appears to be a decreasing function of both the redshift and stellar mass. At z ≤ 3.0, we measure a strong increase of in the two lowest stellar mass bins. This strong evolution observed at low redshift is linked to the injection of kinetic energy by the latest-formed SN populations in a gas that becomes rare and diffuse. We recall that the velocity dispersion is computed at the disk scale height and takes into account only the diffuse gas mass. Less gas means a turbulent energy budget (Sect. 5.4.3) divided into a smaller number of gas atoms and therefore a higher average velocity per mass unit. In parallel, the decrease with the redshift of the ratio is also due to the increase with the redshift of the orbital velocity traced by the decrease of the average disk orbital timescale (see upper left panel of Fig. 18). By tracing average evolutions from low mass at high redshift to high mass at low redshift, the average galaxy trends indicate that the ratio slightly decreases through the evolution.

Fig. 18. Evolution with redshift of 5 main timescales for regulating the growth of disk galaxies which lie along the star-forming galaxy main sequence shown for a range of redshift bins. The legend in the lower right panel gives the stellar mass value for each curve in the panels. Upper left: orbital timescale, torb. Upper centre: characteristic cooling timescale, tcool. Upper right: fragmentation timescale of the gas, tfrag. Lower left: gas disruption timescale, tdisrupt. Lower left: gas depletion timescale, tsfg.

7.3. Timescales

The gas cycle presented in Sect. 4 is based on three gas reservoirs and different exchange rates between these reservoirs. From these transfer rates and knowing the mass in each of these gas reservoirs, we can deduce the various timescales. Evolution with both the redshift and the stellar mass of five main timescales are shown in Fig. 18. We focus on

• the orbital time t_orb estimated at a radius, 2.2r_d (e.g. Pelliccia et al. 2017);

• the cooling timescale, t_cool, given by Eq. (4);

• the fragmentation timescale, , defined as the time required to double the fragmented gas mass;

• the disruption timescale, , defined as the time required to deplete the fragmented gas via the kinetic energy provided by SNe and/or AGN; and

• the gas depletion time of the fragmented gas, , defined as the time required to convert all the fragmented gas and star-forming gas into stars. As the star formation timescale is shorter than all other timescales, t_sfg, is essentially the time required to convert all of the fragmented gas into stars.

We present details of the model output by providing the medians and the 15% and 85% percentiles of these characteristic timescales (Table 3). We compile the statistics for both the full and star-forming samples of model galaxies at z = 2.1. These estimates are given in the four different stellar mass bins 10⁹, 10¹⁰, 10^10.5, and 10^11.25 M_⊙. Even if the median trends presented in Fig. 18 suggest some regularity in galaxy properties as a function of mass, all timescales at constant mass have significant scatter (see Table 3). This large scatter indicates that the dynamics for any one galaxy or ensemble of galaxies is not regular or simple, but is a complex interaction of the various processes included in the model.

First, we first focus on the disk orbital timescale which plays a key role in controlling processes that occur over long timescales. Depending on redshift and stellar mass, the orbital timescale ranges between ≃4 Myr and ≃300 Myr. At fixed stellar mass, the orbital timescale increases with the redshift. At fixed redshift, it is a decreasing function of the stellar mass. The range of orbital timescales calculated using the G.A.S. model are in good agreement with orbital timescales estimated for local galaxies (e.g. Kennicutt 1998; Leroy et al. 2008; Colombo et al. 2018).

The fragmentation of the ISM is driven by energy injection on scales of and larger than the disk height via the condensation of the diffuse gas. The rate of fragmentation depends on the radiative cooling rate of the diffuse warm gas. The characteristic cooling time is dependent on redshift and galaxy mass and ranges over five orders of magnitude (Eq. (4)). At a given stellar mass, the characteristic cooling timescale is a decreasing function of the redshift (Fig. 18). This trend is driven by the increase with the redshift of the diffuse gas density. At a given redshift, the characteristic cooling time is a decreasing function of the mass. These two trends together means that the characteristic cooling timescale of a median galaxy track decreases with both redshift and stellar mass (Fig. 18). At z = 2.1, the characteristic cooling timescale of our star-forming galaxies sample is distributed between ≃0.8 Myr and ≃6.2 Myr. These values are slightly smaller than in the fully resolved galaxy sample of ≃1.3 Myr and ≃8.0 Myr (Table 3). At all redshifts and for all stellar masses, the characteristic cooling timescale is significantly shorter than the orbital timescale (< 6%). If no additional mass were added to the reservoir of diffuse gas, its condensation, which is mainly governed by the characteristic cooling timescale, would deplete completely in less than an orbital time. The short cooling and condensation timescales help to drive the overall complex dynamics of the gas cycle in model galaxies as alluded to earlier.

To provide an indication of how much of the diffuse gas is condensing at any given time, we tabulate the mass fraction of the diffuse gas (Table 3). The mass fraction of the diffuse gas that is condensing can be estimated by comparing the effective cooling time (Eq. (6)) to the characteristic cooling timescale (Eq. (4)). At z = 2.1, for the five stellar mass bins tabulated, the median values of the ratio, 𝒯_cool/t_cool, range between ≃2.7 and ≃3.5 for both model galaxy samples. Low-mass galaxies have the lowest ratios. The ratio is approximately constant, ≃3.5, for galaxies with M_⋆ ≥ 10¹⁰ M_⊙.

The fragmentation and disruption timescales have very similar dependencies on redshift and stellar mass (Fig. 18). These two timescales show a decreasing trend with both redshift and stellar mass. At z = 2.1, for the star-forming galaxy sample, the fragmenting timescale and the disrupting timescale are distributed between 4 and 15 Myr. Those values correspond to ≃10% of the disk orbital timescale. As for the cooling timescale, the fragmentation and disruption timescales of our star-forming galaxy sample are marginally smaller than those of the full galaxy sample. The disruption timescale is always (slightly) larger than the fragmentation timescale (Table 3) and thus gas fragmentation is always faster than gas disruption in G.A.S.The short fragmentation timescales indicate that the growth of GMCs is a continuous and efficient process. Timescales measured in our model are fully compatible with GMC gas accretion rates found via numerical simulations (e.g. Vázquez-Semadeni et al. 2010; Colin et al. 2013).

Working against the gas cooling and fragmentation is the injection of energy by SNe and AGN. This energy injection is used in G.A.S. to disrupt the gas and transfer some of the fragmented gas back into the reservoir of diffuse gas. Hydrodynamical simulations find that disruption timescales of about 10–20 Myr (Colin et al. 2013; Dobbs 2015). This range of values is fully consistent with our predictions.

7.3.1. Impact of fragmentation and disruption on the life cycle of GMCs

The small differences between the condensation and disruption timescales implies that the total mass of fragmented gas in GMCs stays approximately constant and that significant fraction of the fragmented mass is continuously regenerated. From this equilibrium, there are two evolutionary paths for GMCs that depend on both the scale and mass of GMCs.

• Low-mass GMCs, i.e. those formed in galaxies with short or intermediate disk-scale heights, h_d ≤ 50 pc; the mass disrupted by SN kinetic energy injection is close to the total of low-mass GMCs. In such circumstances, after only few SN cycles, GMCs are fully disrupted. These disrupted GMCs are constantly replenished through the formation of new GMCs through the condensation of diffuse gas. Thus the disruption/fragmentation timescales are approximately the lifetime of the GMCs in our model. Our results are in reasonable agreement with the sizes and short GMC lifetimes estimated in local galaxies, 25–70 pc and 17 ± 4 Myr (Murray 2011; Miura et al. 2012; Meidt et al. 2015).

• More massive GMCs resist being disrupted by SNe. The competition between accretion and disruption regulates the mass of GMCs (Vázquez-Semadeni et al. 2010). Because of this competition, GMCs have to be therefore considered as dynamical structures where the gas is just passing through from the diffuse-gas state to the star-forming gas state.

The SN/AGN-driven gas disruption slows significantly the rate at which gas fragments. In small GMCs, the cloud structure and therefore the sites of star formation can be completely destroyed after only a few SN cycles. In more massive GMCs, the fragmentation is also significantly slowed. Indeed, gas recently added to the reservoir of diffuse gas starts to fragment on larger scale than the gas within the GMCs. In our prescription, this constant renewal and exchange of the gas in these reservoirs and the impact of the added gas to the overall cascade of fragmenting gas is fully accounted for via the mass-weighted fragmentation clock (Eq. (19)).

7.3.2. Gas depletion timescale

Depending on the redshift and stellar mass of a model galaxy, the gas depletion timescale, t_sfg, ranges from 80 Myr to 2 Gyr (Fig. 18). The upper end of this range is consistent with the gas depletion timescale measured in local spiral galaxies (e.g. Leroy et al. 2013; Colombo et al. 2018) and the values predicted at high redshift, z ≈ 2 − 3, are also generally consistent with the values estimated for distant star-forming galaxies (e.g. Daddi et al. 2010; Genzel et al. 2010). Gas depletion timescales are a decreasing function of both the redshift and stellar mass. Above M_⋆ ≥ 10^10.5 M_⊙ and at z >  3.0, the average gas depletion timescale reaches a lower limit of close to 80 Myr.

The t_sfg/t_frag ratio provides an estimate of the efficiency of star formation in GMCs. The number of condensation/disruption cycles and the depletion timescale of fragmented gas imply that faster cycles lead to shorter depletion timescales. At z = 2.1, median values extracted from our star-forming sample indicate that between ≃30 and ≃65 cycles are required to convert fragmented gas into star-forming gas (and therefore into stars). This number of gas cycles is consistent with gas dynamics and number of cycles estimated in hydrodynamic simulations (e.g. Semenov et al. 2017). In the star-forming sample, the average number of condensation/disruption-cycles is strictly a decreasing function of the stellar mass. However, in the full galaxy sample, a minimum is reached for models galaxies with mass of 10^10.25 M_⊙. The fragmenting gas in galaxies less and more massive than 10^10.25 M_⊙ need more cycles to convert the fragmented gas into stars. During each cycle only a small fraction, 1%–5%, of the fragmented gas stops fragmenting and is available for the star formation. In fact, the majority of the gas that forms GMCs is recycled in of order one to a few dynamical times or simply remains in a state that cannot form stars. We obtain star formation efficiencies between ≃1% and ≃5%. Such values are consistent with observational estimates (e.g. Wong & Blitz 2002; Leroy et al. 2008; Murray et al. 2010; André 2013; Heiderman et al. 2010) and efficiencies measured in hydrodynamic simulations (e.g. Semenov et al. 2017; Kimm et al. 2017).

In summary, the gas cycle implemented in the new G.A.S. model attempts to capture the complex dynamics of the gas in galaxies. In our prescriptions, gas is continuously exchanged between the diffuse and the fragmented non-star-forming phases. A large number of condensation and disruption cycles, 30–70, are needed to progressively convert diffuse gas into very dense star-forming gas and then stars. The characteristic timescales of ≃15 Myr of these exchange rates is only a small fraction of the orbital timescale. Under such conditions, the fragmented gas is converted into star-forming gas and stars continuously and the gas spends the majority of its time in a non-star-forming gas phase.

8. Discussion and conclusions

We present a new semi-analytical model, G.A.S., in which we implemented a more realistic gas cycle than has been previously implemented in a semi-analytical model. We introduced a prescription for delaying star formation which is underpinned by progressively fragmenting the gas. The formation of giant molecular clouds, filaments, and cold cores takes time and therefore at a given instant only a small fraction of the total gas is in the form of prestellar cores. The majority of the gas is not in a form that is immediately available to form stars. Within this framework, we account for the continuous dissipation of the turbulent kinetic energy through the different ISM scales and phases. We implemented an approximate multiphase ISM in which the gas cycles between a warm diffuse phase, a cold fragmented non-star-forming phase, and a very dense star-forming gas phase. Diffuse warm gas is progressively fragmented following the effective cooling time. Even if the overall depletion timescale of the fragmented gas is, on average, proportional to the disk orbital timescale, we measure a large scatter in both the depletion timescale and the disk orbital timescale. This proves that prescriptions that rely solely on the disk orbital timescale do not properly capture the complexity of gas cycles in galaxies. Smaller characteristic timescales, such as the cooling and/or energy injection timescales due to SNe and AGN are also important. This energy input on shorter timescales, efficiently disrupts the star-forming gas. The large-scale velocity dispersion of the diffuse gas is also maintained by SNe/AGN kinetic energy injection. A fraction of the gas can also be ejected from disks by SN explosions but the characteristic timescale of ejection is on average ten times larger than the fragmenting/disruption cycle timescale. While the SN/AGN kinetic energy injection in the ISM regulates the star formation rate in low-mass galaxies, our model galaxies retain a high percentage of the gas. This “local” gas is then used with a progressively increasing efficiency necessary to build massive, M_⋆ >  10¹¹ M_⊙, galaxies in the early Universe, z ∈ [4; 6].

Star formation occurs only in the very dense gas, which in our model is a product of continuous fragmentation over a wide range of scales. This new complete gas cycle strongly regulates the star formation in our modelled galaxies. Only a few percent of the available fragmented gas is converted into stars during a GMC life cycle. By taking into account the fragmented gas on a galaxy scale, our model is able to reproduce the standard Schmidt–Kennicutt law. Our estimated gas depletion timescales, which are directly related to the time it takes for gas to fragment, are in good agreement with observational estimates.

This new gas regulation cycle leads to very good agreement with the observed stellar mass function over a wide range of redshifts, z ∈ [0.8; 6]. The ability of our model to catch the stellar mass assembly in high redshift galaxies as been already used with success in Lagache et al. (2018). Galaxy properties (gas/stars contents, metallicities, and far-ultraviolet fluxes ...) predicted by our model have been post-processed to successfully predict the [CII] luminosity functions at high redshifts and have allowed us to explore the main characteristics of this emission.

At z <  0.8, we find some discrepancies between the predicted and the observed stellar mass functions. Specifically, an over-density of low-mass galaxies and the under-density of intermediate-mass galaxies modelled at z <  0.8 could probably be solved by slightly increasing the efficiency of photoionisation in our prescription.

At z <  4.0, to reduce the efficiency of or to stop the growth of stellar mass in massive galaxies, the accretion of gas and its transformation into stars has to be reduced or stopped. Previous semi-analytical models implemented strong AGN feedback to quench gas accretion onto massive galaxies. In those models, a significant fraction of the power produced by the AGN is directly used to reduce the cooling of the hot halo gas. This implies a constant AGN power production, which may not be consistent with our understanding of AGN variability (see e.g. Hickox et al. 2014; Stanley et al. 2015; Volonteri et al. 2015, and references therein). For massive galaxies, we propose an alternative model which regulates the radiative cooling and gas accretion onto galaxies. In parallel to radiative cooling, we implemented the development of thermal instabilities creating warm gas surrounding the galaxy. The growth of these thermal instabilities in the central region of the halo progressively reduces or halts the cooling and dissipation of kinetic energy in the halo and therefore reduces or stops the accretion onto galaxies. This process of regulation is a natural outcome of the growth of thermal instabilities in the hot halo gas and does not depend on the power produced by AGN.

However, our actual efficiency parameters do not fully stop gas accretion onto very massive galaxies. In some massive dark-matter haloes, radiative cooling restarts (VFF < 1.0) even if its hot gaseous halo has been fully quenched previously (VFF = 1.0). The recovery of the gas accretion leads to the formation of some (< 10) unobserved overly massive galaxies at z ≃ 0.3. Our prescription of thermal instability growth in the hot gas phase is governed by a set of two parameters which can be further refined to overcome this problem.

In our model of the gas cycle, we assume that the gas initially fragments at the disk scale height with the progressive and continuous formation of over-densities which are akin to observed GMCs. The fragmentation of the gas will progress down to the scale at which stars form or ∼0.1 pc. However, this fragmentation can start at larger scales in the cold streams or in the hot(warm) gas phase surrounding the galaxy. Some clumps of warm gas could already be formed around the galaxy. We assume in our model that a fraction of the newly accreted gas is already fragmented but we do not include any interaction between this already fragmented gas in the halo and the large-scale wind. This kind of interaction could reduce, perhaps substantially, the gas-accretion efficiency (Cornuault et al. 2018). Indeed such coupling could increase the cloud-cloud velocity dispersion and maintain the turbulence in the hot and warm gas contained in the CGM of the most massive galaxy. Such a hypothesised mechanism will be complementary to the development of thermal instabilities, could in fact act as a catalyst to the formation of more clouds in the CGM, and could contribute in quenching the radiative cooling and dissipation in the most massive haloes.

Acknowledgments

MC acknowledges Thomas Fenouillet for his much appreciated help in the use of the LAM computation clusters and wishes to thank Mathieu Génois for his help in the proofreading of this article and for all physical and technical discussions. MC thanks Olivier Ilbert for the numerous and very helpful discussions about stellar mass function (estimations, errors, limits) and Benoit Epinat for a useful discussion about galaxy disk dynamics. MC also wishes to express his appreciation to Alexandre Beelen and Yanick Roehlly for very useful physical and technical discussions. The authors thank the Centre National d’Etudes Spatiales (CNES), Aix Marseille Université, Sorbonne Université, the Programme National de Cosmologie and Galaxies (PNCG) and Programme de Physico-chimie du Milieu Interstellaire (PCMI) of CNRS/INSU for their financial support. PG gratefully acknowledges the support of the Institut Universitaire de France (IUF).

References

Appendix A: The GALAKSIENN library

The GALAKSIENN library stores the main results produced by our new G.A.S. semi-analytical model, especially MOCK galaxy catalogues and sky maps. It is available online through the ZENODO platform: https://zenodo.org/, DOI: 10.5281/zenodo.1451229. A complete description of the GALAKSIENN library is given in Paper III. In association with this Paper I, we distribute the ASCII tables of the stellar mass functions (Fig. 14).

All Tables

All Figures

Fig. 1. Relative mass fraction of baryons which are accreted smoothly. The fraction is plotted as a function of both the dark-matter virial mass and redshift. The redshift of each curve is colour coded by the colour bar on the right. Each curve, starting at z = 0.0, is separated by Δz = 0.1. For z >  zreion = 7.0, we set . The grey vertical bar indicates the minimum halo mass in our dark-matter model. The grey horizontal line indicates f = 0.5.

In the text

	Fig. 2. Cooling efficiency, Λ(Thot, Zhot), as a function of both gas temperature and gas metallicity (see De Rijcke et al. 2013). The values of Λ(Thot, Zhot) are colour coded using the colour bar on the right (in units of erg s−1 cm−3).
In the text

	Fig. 3. Mean accretion rates into the galaxy halo galaxy from the cold (upper panel) and hot mode accretion (lower panel) as a function of virial mass. The mean values are computed by selecting only haloes which are actively accreting gas. Colour code indicates redshift. The grey solid vertical line inidcates the dark-matter halo mass for which contributions of cold and hot modes are equivalent.
In the text

	Fig. 4. Diagram illustrating how the gas fragments during the turbulent inertial cascade. When a structure formed at the scale k reaches the Bonnor–Ebert mass, it breaks into smaller structures at the next level, 2k. The mass transfer rate, Ṁk, between scales is given by Eq. (12). The largest structures are fed at a rate, Ṁfrag, as discussed in Sect. 4.1 and given by Eq. (14).
In the text

Fig. 5. Gas fragmentation rate as a function of both the maximum energy injection scale, lmax, which we assume is equal to the disk scale-height, and the growth rate of GMCs (Eq. (14)). Each coloured box shown in the diagram represents the gas fragmentation rate for a specific stellar mass bin as indicated by the labels at the lower right corner of each box. The mass bins range from M⋆ = 109.0 M⊙ to 1011.5 M⊙. The limits of each box encompass the 15% and 85% percentiles of the distribution for each mass bin. Coloured points within each box indicate the value of the median of the distribution. The regions shown are for modelled galaxies at z = 2.1.

In the text

Fig. 6. Diagram illustrating the flow of mass between the three gas reservoirs considered in our multiphase G.A.S. model. The star-forming reservoir contains the mass that is immediately available for star formation. It is fed by a turbulent cascade and gas fragmentation at a mass flow rate, Ṁsfg (Eq. (18)). The gas reservoir fragments at a rate, Ṁfrag (Eq. (14)). The diffuse gas reservoir is fed through the rates of two mechanisms: (1) gas accretion rate (Eq. (1)); and (2) the rate at which the fragmented and star-forming gas is disrupted via the energy injected by SNe and AGN, Ṁdisrupt (Eqs. (27) and (28)). Each of these three gas reservoirs contributes to the outflow rate, Ṁwind (Eqs. (24) and (26)).

In the text

Fig. 7. Schmidt–Kennicutt relations estimated from our model galaxies at both the galaxy scale and GMC scale. In the three panels the solid blue line indicates the Kennicutt (1998) relation. The grey shaded area is the distribution of a full sample of star-forming galaxies in the redshift range z = 1.5–2.1. Each coloured point indicates the position of an individual star-forming galaxy, uniformly selected within the galaxy-scale Σgas − ΣSFR plan. The colour bar indicates the fragmented gas fraction in each representative galaxy. Left panel: the solid orange line is the best fit computed for our entire sample of star-forming galaxies within this redshift range. The dashed contour indicates measurements of nearby galaxies (Bigiel et al. 2008). The blue stars and points are from the starburst and normal disk galaxy samples from Kennicutt (1998), respectively. Central panel: the solid black line indicates the median trend of our star-forming galaxy sample. The dashed black lines indicate the 15% and 85% percentiles ranges of our model galaxies. We estimate an average scatter in the distribution of 0.46 dex. Right panel: distribution of galaxies in the Σgas − ΣSFR estimated over the GMC-scale of star formation (see text for details). We show this estimate for the same galaxy sample used in the galaxy-scale Schmidt–Kennicutt relation in the two leftmost plots in this figure.

In the text

Fig. 8. Illustration showing the distribution of accretion and energy output of a model AGN. Gas infalling from the torus to the AGN is divided between accretion onto the SMBH Ṁacc and ejection Ṁjet. We assume the diffuse gas in the disk can be driven outwards by the AGN to give the overall ejection rate of gas Ṁwind,agn. The accretion onto the SMBH Ṁacc is itself divided in two parts: a fraction (1−fME)Ṁacc of the infalling gas goes to increasing the mass of the SMBH, while the remaining gas Qagn = fMEṀacc is converted into energy (see Sect. 5.3 for details). The radiative and mechanical energy of the AGN is then distributed in various ways as illustrated in Fig. 9.

In the text

	Fig. 9. Illustration of how the energy output from AGN and SNe are distributed. The SNe and AGN power QSN and QAGN are shown in red and yellow, respectively. Specifically, the power is distributed through four different channels: (1) as a large-scale wind, (2) in generating turbulence, (3) in heating the gas, and (4) output bolometric luminosity.
In the text

Fig. 10. Probability distribution of the velocity dispersion of the diffuse gas. This distribution represents model galaxies with M⋆ >  107 M⊙ and at z = 1.1. We compare this distribution with observational measurements from Swinbank et al. (2017, dark blue histrogram) and Pelliccia et al. (2017, light blue histrogram). The solid and dashed vertical lines represent the median and mean value of the distribution of our modelled galaxies, respectively.

In the text

Fig. 11. Thermal instability efficiency, F(Λ), as a function of both gas temperature and gas metallicity. The amplitudes of thermal instabilities are computed as given in Sharma et al. (2012) and the values are indicated in the colour bar on the right side of the plot. When F(Λ) < 0, the gas is thermally unstable. The solid and dashed blue contours indicate where 2 − dlnΛ/dlnT = 0 and −3.0, respectively.

In the text

Fig. 12. Upper panel: mean value of the volume filling factor, ϕ, of gas that is thermally unstable (Eq. (36)). The case ϕv = 1 implies that all the region within r = rcool is thermally unstable in our model. Lower panel: fraction of galaxies that are accreting radiatively cooling halo gas as a function of halo virial mass, Mvir. Different colours indicate this fraction as a function of redshift. The value of the redshift for each line is indicated in the colour bar on the right side of the upper panel. The grey area spans the range of the mean fraction of galaxies that are accreting hot gas for each halo mass summed over redshifts 0.1–10 (Eq. (2)). For example, for a halo mass of 1010.5 M⊙, the fraction of galaxies that have hot gas in the halo which is radiatively cooling can range between 30% and ∼55% depending on redshift.

In the text

	Fig. 13. Redshift evolution of the co-moving stellar mass volume density. The solid orange line is the prediction of the G.A.S. model. We compare our predictions with the compilation of observational results (see Madau & Dickinson 2014, and references therein).
In the text

Fig. 14. Comparison between the predictions of G.A.S. with observational estimates of the co-moving number density of galaxies for 9 redshift bins spanning the range, z ≃ 0.1–z ≃ 6 (from upper left to lower right and as indicated at the top of each panel). We compare the predictions of the G.A.S. model (solid grey lines) with various studies: Baldry et al. (2008, red squares), Yang et al. (2009, orange circles), Caputi et al. (2011, brown diamonds), Grazian et al. (2015, chocolate squares), Duncan & Conselice (2015, chestnut-red triangles), Song et al. (2016, salmon circles), Ilbert et al. (2013, yellow dashed-line), and Davidzon et al. (2016, red dot-dashed line). For this comparison with observations, we indicate binned measurements using points and best-fit stellar mass functions as dashed lines (Schechter and double-Schechter). All stellar mass functions are for a Chabrier IMF. The best fits of Ilbert et al. (2013) and Davidzon et al. (2016) are corrected for the Eddington bias which affects the high-mass galaxy bins in particular. The grey dotted-line indicates the total gas plus stars, baryonic mass function predicted by G.A.S.The grey dashed-line indicates the stellar mass function of the previous panel for the next highest redshift bin. The tick marks placed on the left side of each panel indicate the number of galaxies used in the stellar mass function predicted by the G.A.S. model.

In the text

Fig. 15. Comparison of the stellar mass functions predicted by our best fiducial model with those predicted by Guo et al. (2011). Each panel shows a different redshift as indicated at the upper left in each panel; solid red, dashed orange, and dashed gold lines indicate the predictions of G.A.S., Guo et al. (2011), and Cousin et al. 2019, A&A, 627, A132 (2015a), respectively (as indicated in the legend of the left-most panel). We compare these models with observations (dashed grey lines; Davidzon et al. 2016). At z = 4.0, we also plot additional observational measurements (salmon circles, red upward pointing triangles, and orange squares; Song et al. 2016; Duncan & Conselice 2015; Grazian et al. 2015, respectively, and as indicated in the legend of the right-most panel)

In the text

Fig. 16. Impact of our main prescriptions for regulating star formation and gas accretion, photoionisation, re-accretion, and thermal instabilities on the stellar mass assembly of model galaxies (i.e. the co-moving density of model galaxies as a function of redshift). The solid red line shows our fiducial model using the full suite of prescriptions. The dot-dashed gold lines indicate our model without photoionisation. In the three right-most panels, the dotted brown lines and dashed orange lines indicate predictions of the stellar mass functions without including the prescriptions for re-accretion of gas over two halo crossing times and thermal instabilities, respectively (as indicated in the legends of the two left-most panels). These are compared to the observational estimates of the co-moving density of galaxies from Davidzon et al. (2016). The four different panels represent four different redshift bins, z = 0.3, 1.5, 2.8, and 4 (from the left-most to the right-most panel, respectively).

In the text

Fig. 17. Evolution with redshift of 5 main disk properties of disk galaxies which lie along the star-forming galaxy main sequence shown for a range of mass bins. Legend in the lower right gives the stellar mass values for each curve in the panels. Upper left: gas mass fraction, fgas, as a function of redshift for 5 stellar mass bins. Upper centre: fraction of the gas which is fragmented, ffrag. Upper right: average half mass radius, r50, of the ensemble of star-forming model galaxies. Lower left: diffuse gas density, ρdiff. Lower right: ratio of the gas velocity dispersion to orbital velocity, .

In the text

Fig. 18. Evolution with redshift of 5 main timescales for regulating the growth of disk galaxies which lie along the star-forming galaxy main sequence shown for a range of redshift bins. The legend in the lower right panel gives the stellar mass value for each curve in the panels. Upper left: orbital timescale, torb. Upper centre: characteristic cooling timescale, tcool. Upper right: fragmentation timescale of the gas, tfrag. Lower left: gas disruption timescale, tdisrupt. Lower left: gas depletion timescale, tsfg.

In the text