G.A.S. II: Dust extinction in galaxies; Luminosity functions and InfraRed eXcess

Dust is a crucial component of the interstellar medium of galaxies. The presence of dust strongly affects the light produced by stars within a galaxy. As these photons are our main information vector to explore the stellar mass assembly and therefore understand a galaxy's evolution, modeling the luminous properties of galaxies and taking into account the impact of the dust is a fundamental challenge for semi-analytical models.We present the complete prescription of dust attenuation implemented in the new semi-analytical model: G.A.S. This model is based on a two-phase medium originating from a physically motivated turbulent model of gas structuring (G.A.S. I paper). Dust impact is treated by taking into account three dust components: Polycyclic Aromatic Hydrocarbons, Very Small Grains, and Big Grains. All three components evolve in both a diffuse and a fragmented/dense gas phase. Each phase has its own stars, dust content and geometry. Dust content evolves according to the metallicity of it associated phase.The G.A.S. model is used to predict both the UV and the IR luminosity functions from $z=9.0$ to $z=0.1$. Our two-phase ISM prescription catches very well the evolution of UV and IR luminosity functions. We note a small overproduction of the IR luminosity at low redshift ($z<0.5$). We also focus on the Infrared-Excess (IRX) and explore its dependency with the stellar mass, UV slope, stellar age, metallicity and slope of the attenuation curves. Our model predicts large scatters for relations based on IRX, especially for the IRX-$\beta$ relation. Our analysis reveals that the slope of the attenuation curve is more driven by absolute attenuation in the FUV band than by disk inclination. We confirm that the age of the stellar population and the slope of the attenuation curve can both shift galaxies below the fiducial star-birth relation in the IRX-$\beta$ diagram.


Introduction
Galaxies are defined by their stellar populations -the '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 characterize in galaxy evolution models, essentially because the formation of stars involves many non-linear processes that occur over a large range of temporal and spatial scales in e.g., density, velocity, and magnetic field strength and regularity Krumholz & McKee 2005). Observations show that star formation in galaxies is a very inefficient process, with typically 0.1%-10% of the available gas being 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 the star-formation efficiency is highly dependent on how ionization and kinetic feedback is injected into the interstellar medium (ISM; e.g., Krumholz & Thompson 2012;Gatto et al. 2015;Hennebelle & Iffrig 2014;Dale 2015;Geen et al. 2017;Gavagnin et al. 2017).
The difficulty in simulating feedback-regulated star formation, as well as 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, make the global regulation of star formation one of the greatest challenges in modelling galaxy formation and evolution. 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 supernovaedriven shocks. High resolution hydrodynamic simulations, with implementation of sub-grid turbulence models (Schmidt et al. 2013;Semenov et al. 2016), show complex, multi-phase, turbu-lent structures within the ISM with a realistic global Kennicutt-Schmidt relation on kpc scale, 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 (τ depl = M g /Ṁ ≈ 1 − 10 Gyr, where M g andṀ are respectively the gas mass and the star-formation rate), 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 re-cycling 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 the mass of the galaxy, different prescriptions for various feedback mechanisms have been invoked to regulate star formation. For low stellar masses, M < 10 9 M , feedback is used to either regulate gas accretion, mostly by heating the gas through photo-ionisation (e.g., Doroshkevich et al. 1967;Couchman & Rees 1986;Ikeuchi 1986;Rees 1986), or ejecting the gas using the mechanical energy generated by supernovae (SN, e.g., White & Rees 1978). For high stellar masses, M > 10 10 M , models rely on the action of Super-Massive Black-Holes (SMBH) to inhibit gas accretion onto galaxies. A significant fraction of the power generated by AGN is used to limit the cooling of the hot gas phase surrounding the galaxy.
Despite including 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 and dynamical coupling between gas phases are still lacking. As explained in Cousin et al. (2015a), current semi-analytical model 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 11 M observed at redshifts greater than 3, are not robustly reproduced in current models.
We present here a new semi-analytical model G.A.S.the Galaxy Assembler from dark-matter Simulation -which is based, in part, on previous version described in Cousin et al. (2015b) and Cousin et al. (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: G.A.S. II: in which we describe and model the mechanisms that leads to dust attenuation of galaxian light and G.A.S. III in which we explore the panchromatic emission of galaxies from the FUV to the sub-millimetre 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. 2015a). Here 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 following. 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. Sect. 6 focuses 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 large velocity dispersions. This allows us to define an effective cooling rate. In section 7 we present and discuss our results, mainly focusing 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 halos.

From dark-matter to baryons
2.1. Dark-matter G.A.S. is built upon a set of dark-matter merger trees extracted from a pure N-body simulation. The current simulation uses WMAP-5yr cosmology (Ω m = 0.28, Ω Λ = 0.72, f b = 0.16, h = 0.70, Komatsu et al. 2009) with a volume of [100/h] 3 M pc in which 1024 3 particles evolve. Each particle has a mass of m p = 1.025 10 8 M . Halos and sub-structures (satellites) are identified by using the HaloMaker code (Tweed et al. 2009). In our merger trees, we only consider halos with at least 20 darkmatter particles leading to a minimal dark-matter halo mass of 2.050 × 10 9 M . Dark-matter halos grow from smooth accretion. The dark-matter accretion rate,Ṁ dm , only includes particles that are newly detected in the halo, and that have never been previously identified in another halo.

Baryonic accretion
Based on the dark-matter accretion rateṀ dm , the accretion rate of baryons is defined by,  (M h , z) is the effective baryonic fraction. The baryon accretion rate depends on the gas ionisation state and we adopt here the photo-ionisation model based on the Gnedin (2000) and Kravtsov et al. (2004) prescription, but using the effective filtering mass given by Okamoto et al. (2008). We assume that re-ionisation occurs at z = 7.0, we therefore set Fig. 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 ( f ph−ion b / f b < 0.5) of the photo-ionisation prescription occurs at low redshift (z < 0.9).
As initially proposed by Khochfar & Silk (2009), in G.A.S. we define two different modes of accretion, a cold 1 mode, and a hot mode. Depending on the dark-matter halo mass, the fraction of accreted, hot gas is computed as 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 in two parts that feed two different reservoirs: M cold for the cold mode and M hot for the hot mode. Both the cold and the 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 1 The temperature of the accreted gas in this mode is close to 10 4 K. 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.

The 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:Ṁ 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 .

Hot accretion: Radiative cooling
As in the previous versions of this model (Cousin et al. 2015b(Cousin et al. , 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 atmosphere a constant temperature T hot , and a average gas metallicity Z hot . For more details, please refer to Cousin et al. (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, 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 3 K and 10 8 K, and for gas metallicities over the range 10 −4 Z and 2Z . We assume the solar metal mass fraction of Z = 0.02. Fig. 2 shows the cooling efficiency as a function of both gas temperature and average gas metallicity.
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: In this equation, T hot cool is a "cooling clock", which is a measure of the effective cooling time of the hot gas. After each timestep ∆t, the mass-weighted cooling time is updated via: where M is the total mass in the hot phase after the latest timestep, ∆t. ∆M is the net mass variation of hot gas during this timestep (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 2 , 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. Knowing the cooling radius r cool (deduced from Eq.5), we can write the condensation rate of the gas: 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 . Fig. 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 halos with M h = 10 10.5 M , the average accretion due to the cold mode decreases from 20M yr −1 at z 9.0 to 0.2M yr −1 at z = 0.3; for the hot mode, accretion rate decreases from 30M yr −1 to less than 0.3M yr −1 between z 9.0 and z = 0.3.  The mean values are computed by selecting only halos which are actively accreting gas. Overall, as expected, the both the cold and hot mode accretion rates decline with decreasing redshift but the relative fraction of hot mode accretion increases with decreasing redshift and for the most massive halos (M h > 10 10.5 M ). This trend is the result of the cooling regulation related to thermal instabilities. Colour code indicates redshift. The grey solid vertical line marks the dark-matter halo mass where contributions of cold and hot modes are equivalent.

The 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 development of multi-phase 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 small 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 phenomenon 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).

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. They link the mass surface density µ and the 1D-velocity dispersion σ, respectively, to the wave number k = 1/l: 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): 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: b = 1 3 (2 − a). Following Romeo et al. (2010), we assume (a, b) = (1/5, 3/5). These values Fig. 4. Diagram illustrating how the gas fragments during the turbulent inertial cascade. When a structure which 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,Ṁ f rag , as discussed in Sect. 4.1 and given by Eq. 14.
reproduce observations of the dynamics of clouds at the atomicmolecular transition. The energy transfer rate during the inertial cascade is now proportional to: The normalisation parameters used in Eq. 8 are listed in Table 2, and correspond to standard values for dense, supersonic, compressible gases (Jog & Solomon 1984b,a;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 pre-stellar cores (e.g., 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 −1 (e.g., Orkisz et al. 2017), a value slightly higher than the speed of sound for a molecular gas at 10 K (c s (10K) = 0.2 km s −1 ), 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 −2 (Gao & Solomon 2004;Heiderman et al. 2010;Lada et al. 2012).

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 kpc 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 Fig. 5. The gas fragmentation rate as a function of both the maximum energy injection scale, l max , which we assume is equal to the disk scaleheight, 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 bottom right corner of each box. The mass bins range from M = 10 9.0 M to 10 11.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 indicated the value of the median of the distribution. The regions shown are for modelled galaxies at z = 2.1.
(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. The Bonnor-Ebert mass at a given scale 1/k, M BE 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 : We derive the mass transfer rateṀ k between wave numbers k and 2k, using the conservation of energy: where V k = π 6 k −3 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).

Gas fragmentation timescale
We define the gas fragmentation timescale, T f rag , 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Ṁ f rag . During the process of fragmentation, we assume that the disk scale-height is constant. . Diagram illustrating the flow of mass between the three gas reservoirs we considered in our multi-phase 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,Ṁ s f g (Eq. 18). The gas reservoir fragments at a rate,Ṁ f rag (Eq. 14). The diffuse gas reservoir is fed through the rates of two mechanisms, (1) gas accretion rate, (1 − f in f rag )Ṁ b (Eq. 1), and, (2) the rate at which the fragmented and star-forming gas is disrupted via the energy injected by SN and AGN,Ṁ disrupt (Eqs. 27 and 28). Each of these three gas reservoirs contributes to the outflow rate,Ṁ wind (Eqs. 24 and 26).
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. Fig. 5 shows our estimated gas fragmentation timescale T f rag [h d ,Ṁ f rag ] as a function of the instantaneous disk scaleheight h d and mass flow rate at which the largest structure is fed, M f rag (see Sect. 4.1 and Eq. 14). Both the average disk scaleheight and the GMC 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Ṁ f rag . To take this into account, we will define an effective disk fragmentation timescale T f rag , which depends on the history of the disk (see Sect. 4.2 and Eq. 19).

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 multi-phase and we set the mass fraction of fragmented gas to f in f rag = 1/3 (van de Voort & Schaye 2012). We further discuss the physical origin and consequences of that assumption in Sect. 8.

The diffuse gas reservoir
The first main reservoir stores a diffuse ( 1 cm 3 ) and warm (10 4 K) gas which is traced by emission lines of the warm, ionised medium such as [Oiii], [Nii], etc. (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 due to supernovae and/or an actively accreting super-massive black hole (see Sect. 5.4.3). The initial temperature of the warm, diffuse gas is assumed to be 10 4 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 T unstr,n cool is computed after each time step ∆t, as the sum of the mass-weighted cooling times of the halo and the newly incoming gas: where M is the total mass of diffuse gas after the time-step ∆t. ∆M is the mass increase of diffuse gas, coming from both 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Ṁ f rag is computed using the cooling timescale,Ṁ where M di f f is the mass of diffuse gas. The cooling timescale is computed as a function of the average metallicity Z di f f (Eq. 4).
We assume a temperature of 10 4 K and an average volume of the warm, diffuse gas component 3 V di f f = 22πr 2 d h d . The other parameters in Eq. 14 are computed as follows: φ m , the mass fraction of diffuse gas that condenses in an effective cooling time T di f f cool , depends on the characteristic cooling timescale of this reservoir t warm cool . φ m was calculated in Cornuault et al. (2018), and for computational purposes we fitted their computation by an error function given by, where x = T di f f cool /t warm cool . The best fit gives x t = 0.55 and s = 0.13.
where M di f f , M f rag and M s f g 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, 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.

The fragmented gas reservoir
The clumpy gas phase has low filling factor (typically less than ≈ 10%) with densities n H = 1 − 10 3 cm −3 and temperatures 10 3 -10 K. It is traced by atomic and molecular gas lines such as CO, 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 M 1/h BE 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:Ṁ T f rag is the effective disk fragmentation timescale. This timescale 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Ṁ f rag ). If M is the total mass of fragmented gas after the last time-step ∆t, and ∆M = ∆tṀ f rag is the gas mass incorporated in the fragmented gas reservoir, thus the disk fragmentation timescale is updated as: This timescale is a function of T f rag [h d ,Ṁ f rag ], which depends on values of h d andṀ f rag , the disk scale-height and the growth rate of GMCs, respectively (see Sect. 3.3).

The 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 s f = l /σ = 0.1Myr. The star-formation rate is then simply given by: In our prescription, the rate of star formation is mainly limited by the rate at which gas becomes clumpy.

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 starformation rate surface density Σ S FR is computed in the same way. Using these two definitions, our predictions of the galaxyscale Schmidt-Kennicutt law is shown in Fig. 7.
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 10 Σ 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 over the redshift range z=1.5-2.1, we estimate a powerlaw index, N = 1.232 ± 0.002. This slope is slightly shallower than the standard Kennicutt (1998) slope. In addition, we find a clear trend for 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: 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 follow: 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 (Eqs. 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.

Model Feedback
Massive stars and active galactic nuclei (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. 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.

Morphology and the efficiency of outflows
Galaxies in the early universe frequently have a clumpy morphology that suggest there are interacting regions of dense gas and stars (e.g., Elmegreen et al. 2009;Elmegreen 2009  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). Blue stars and points are from the starburst and normal disk galaxy samples from Kennicutt (1998), respectively. Central panel: The solid black line is the median trend of our star-forming galaxy sample. 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: The distribution of galaxies in the Σ gas − Σ S FR 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.
clumps in these distant galaxies are not typically observed over 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 −1 , and V w = 200 km s −1 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 25 th , 50 th and 75 th percentiles for all starforming galaxies, which values areṀ wind /SFR= 6.5, 11.5, and 26.4, respectively.

Supernovae feedback
In Cousin et al. (2016) the instantaneous SN event rate (η S N Gyr −1 ) is proportional to the star-formation rate, which is a function of time (related to the star-formation history; see Sects. 3.5.2). Based on the rate of SNe, the instantaneous power generated by SN is simply given by Q sn = η sn E sn . We assume the standard value E sn = 10 44 erg s −1 (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. (2016). Please see that paper for details.

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 Super Massive Black Holes (SMBHs) are created during a major merger if the remnant galaxy have a bulge of at least 10 6 M . The seed of the SMBH is given an initial mass, M init • = 300M . We associate a gas torus to 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: where: µ g is the mass fraction of the gas that is diffuse in the remnant galaxy disk; -µ m = MIN(M 1 ,M 2 ) MAX(M 1 ,M 2 ) 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 di f f (r < 3r torus ) is the mass of diffuse gas enclosed in r < 3r torus . The torus radius, r torus = 10pc. 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 − f ME )Ṁ acc of the infalling gas goes to increasing the mass of the SMBH, while the remaining gas Q agn = f MEṀacc is converted into energy (see Sect. 5.3 for details). The radiative and mechanical energy of the AGN is then distributed in the various ways as illustrated in Fig. 9.
From Eq. 22 it is clear that major mergers with gas-rich progenitors will accrete the largest amount of gas onto the torus of the remnant, while minor mergers of gas-poor progenitors will 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: 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). 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 on to 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 SMBḢ M e j . This yields the total infall rateṀ in f all =Ṁ e j +Ṁ acc . In our model, we assume that the relative fraction of ejected to accreted mass flow rates is η • =Ṁ e j M acc = 1.0. This division results in equal shares of the infalling gas to be: (i) accreted onto the SMBH, driving the increase the SMBH mass; and (ii) generates 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Ṁacc c 2 .

Distribution of feedback power
We use the power output from AGN and SN 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 S N -are each divided in two parts: kinetic power fraction f k,AGN and f k,S N , 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,S N respectively, are used to heat the ejected gas (Sect. 5.4.4).

Large scale ejecta
We derive the instantaneous ejection rate using the conservation of the kinetic energy released by SNs (e.g., Dekel & Silk 1986;Kauffmann et al. 1993;Efstathiou 2000): where v w is the average velocity of large scale wind (Sect. 5.1). 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 (Ṁ e j ; Sect. 5.3). We derive the velocity of the jet using conservation of energy: 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: The total instantaneous ejection rate of the combined action of SN and AGN isṀ wind =Ṁ wind,S N +Ṁ wind,AGN . Each gas reservoir, M di f f , M f rag and M s f , contributes to the instantaneous ejection rate in proportion of its mass fraction (Fig. 6).

Overestimating the gas escape fraction: Re-accretion timescale
When the gas is ejected from the disc, a fraction of the gas will remain in the hot circumgalactic medium. The remaining fraction is ejected from the dark matter potential. As in Cousin et al. (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, e.g., 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. Due to 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) adopt 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 .

Disruption rate
The power generated by SNs 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 6 yr; Cousin et al. 2016). The fraction of the SN energy which is injected into the fragmented and star-forming gas f f rag disrupts this phase, feeding the diffuse gas reservoir (Eq. 16 and Fig. 6). The disruption rate is defined as: 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 SNs power 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 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: As for SNs, the residual power 1 − f f rag 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,S N +Ṁ disrupt,AGN .

Radiative heating and bolometric luminosity
In the previous two sections we discussed our prescriptions for the kinetic power of SN and AGN. The non-kinetic fraction of the total power is also divided in two different parts. A fraction f th,S N ( 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 6 and 10 7 K. The residual power (1 − f k,S N )(1 − f th,S N ) for SNs 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).

Average velocity dispersion
The gas disruption rates, as given in Eqs. 27, 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, SNs, and AGN.
Our prescription is based on the evolution of the "turbulent" kinetic energy budget E σ v : σ v being 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 2∆E acc σ v = f incr ∆Mv 2 acc , 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 di f f , of diffuse gas. 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 SNs and AGN, Considering all of these energy terms, after a time step ∆t, we have, 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).
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, where M i is the mass (baryon + dark-matter) included in the galaxy half mass radius, and r i d 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 programs 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 to do reproduce some of the most extreme dispersions observed (both small and large values; Fig. 10). At z = 1.1, the median value of the velocity distribution of the diffuse gas is close to 35 km s −1 .

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 parameters of the turbulent cascade. The disk scale height h d is deduced from σ v = σ(h −1 d ). 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).

Thermal instabilities
For galaxies with large masses, previous semi-analytical models invoked powerful AGN feedback to greatly reduce or even completely quench accretion of cooling gas in their halos (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 circum-galactic gas, which results in powerful outflows expelling the gas and heating the halo gas. However, observations are not entirely clear on the precise impact of outflows on distant galaxies(see, e.g., Mullaney et al. 2015;Netzer et al. 2014Netzer et al. , 2016Scholtz et al. 2018;Falkendal et al. 2018). In Cousin et al. (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 efficiently reduce the gas accretion onto massive galaxies, M > 10 11 M . This new mechanism is based on the growth of thermal instabilities in the hot halo phase surrounding the galaxy.

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 circum-galactic medium. In the standard model of cooling presented previously and also used in e.g., Croton et al. (2006); Baugh (2006); 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 4 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, find that the cloud-cloud velocity dispersion can be comparable to the characteristic virial velocity of the darkmatter 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.

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, T T I . This timescales evolves concomitantly with the cooling timescale. After each time step, ∆t, the thermal instability clock is updated following a mass-weighted prescription, similar to the one applied to the cooling clock, namely, where M is the total mass stored in the hot phase after the last time step, ∆t. ∆M is the mass accreted by the hot gas phase and/or ejected from the galaxy. For each time step, ∆t, the time increment ∆t T I is calculated as, where F(Λ) is the thermal instability function (Sharma et al. 2012), depending on the logarithmic derivative of the cooling efficiency function Λ, 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 4 . 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 T I is set to 0. According to the preexisting/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 3 K and 10 8 K and gas metallicity between 10 −4 Z and 2Z . 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.

The 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 center of the halo. We define the instantaneous size of the mixing zone as, 4 Following a scheme where gas is continuously added into the hot gas reservoir.
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 ε T I 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 will 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 circumgalactic medium. 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 T I = r cool − ∆r T I .

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, Mean value of the volume filling factor, φ, of gas that is thermally unstable (Eq. 36). φ v = 1 implies that all the region within r = r cool is thermally unstable in our model. Lower panel: The fraction of galaxies that are accreting radiatively cooling halo gas as a function of halo virial mass, M vir . 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 gray area spans the range of the mean fraction of galaxies that are accreting hot gas for each halo mass summed over redshifts 0.1 to 10 (Eq. 2). For example, for a halo mass of 10 10.5 M , the fraction of galaxies that have hot gas in the halo which is radiatively cooling can range between 30% to ∼55% depending on redshift.
φ 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. φ v = 1 indicates that all the 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. Fig. 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 halos, M vir ¿10 9 M , have φ v dis-tributed between 0.3 and 0.8. At z ≥ 9.0, the average φ v of the first halos 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 halos 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 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 halos, M vir 10 9 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 halos in the low mass bins at very high redshift, z 9, is low and thus these bins are greatly affected by statistical noise.
Over the redshift range, z = 4.0 to z = 1.5, the volume filling factor of the most massive halos reaches very high φ v , 0.95. At z < 1.0, the average volume filling factor of thermally unstable gas in massive halos starts to decrease. In these massive halos, some accretion still occurs. This accretion produces a new star formation and therefore some large scale ejection events. In these massive halos 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. Fig. 12 also shows the change 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 it strongly decreases and reaches values smaller 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 halos, the mass of hot halo gas is small, M hot ≤ 10 7 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 1M /yr (Fig. 3). The fraction of dark matter halos hosting actively radiatively cooling hot gas increases (or remains constant) as a function of the dark-matter virial mass. The largest population of actively cooling hot halos is reached between 10 10.25 M and 10 12 M . Between z = 4.0 and z = 3.0, in the most massive halos > 10 12 M the fraction reaches 1.0. However by z < 3.0 for the most massive halos, we find a strong decrease in the fraction of galaxies with actively cooling halos. This decrease is directly linked to the high φ v reached in these halos and this decrease is responsible for the reduction in the effective cooling rate estimated for the most massive halos at z < 2 (Fig. 3).

Results and analysis
7.1. The stellar mass assembly and co-moving density of galaxies We compare the predictions of the G.A.S. model for the comoving 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 are in agreement with observations. Predictions are in reasonable agreement with estimates even at the highest redshifts, z > 3, although 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 to 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. (2015a). These comparisons allow us to gauge the impact of including the starformation and gas-dissipation prescriptions that are in G.A.S. with models that do not include such prescriptions (Fig. 15).
Between z 1.0 and z 6.0, our model agrees very well with observational measurements of both the cosmic comoving 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 comoving density of intermediate mass galaxies and over-predicts the comoving density of low-mass galaxies. We find that in the low and intermediate mass range, 10 8.5 ≥ M ≤ 10 11.5 M , the comoving 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 the 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 due to 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, ε T I , could limit the growth of these very massive model galaxies.
Compared to the ad-hoc recipe (∝ M 3 h ) presented previously in Cousin et al. (2015a), the new regulation prescription based on the turbulent cascade is in better agreement with the comoving 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 comoving density of low/intermediate 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/intermediate mass galaxies, hence producing an excess. Furthermore, a lack of gas in more massive galaxies means that their models under-estimate 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 photo-ionisation 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 3 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 3 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 the number of lower mass galaxies, those with M < ∼ 10 11 M essentially unchanged. Photo-ionisation 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 photo-ionisation prescription (Sect. 2.2, Fig. 1). The process with the largest impact on the co-moving number density of galaxies appears to be the reaccretion 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 comoving density of almost all stellar masses but especially so for intermediate-mass galaxies.

Evolution of disk properties
We discuss with Fig. 17 the evolution with both the stellar mass and the redshift of five main properties of star-forming galaxies. We present the evolution of: -The gas mass fraction, f gas = M gas M gas +M . 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 50 .
-The diffuse gas density, ρ di f f . 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, σ v V . 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 implies 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 11 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 larger fraction 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 larger than 10 10 M stabilise their fragmented gas mass fraction at around 42%. For less massive galaxies, < 10 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 disc half mass radius is a decreasing function of the redshift. At redshift z > 3.0, the disc half mass radius is a decreasing function of the stellar mass. Massive galaxies appear more compact than low-massive 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 3 . The less massive galaxies (≤ 10 10 M ) see their diffuse gas density strongly decrease and reach a value of < 0.1atoms/cm 3 . This evolution is mainly linked to the decrease of the fragmented gasmass fraction discussed previously.
This decrease also impacts the dispersion to orbital velocities ratio, σ v V . This ratio appears to be a decreasing function of both the redshift and the stellar mass. At z ≤ 3.0, we measure a strong increase of σ v V in the two lowest stellar mass bins. This 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, indian-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.  Fig. 16. The impact of our main prescriptions for regulating star formation and gas accretion, photo-ionisation, 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 presented here. The dot-dashed gold lines indicate our model without photo-ionisation. In the three right-most panels, the dotted brown lines and the dashed orange lines are 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 are for four different redshift bins, z = 0.3, 1.5, 2.8, and 4 (from the left-most to the right-most panel respectively).
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 σ v V 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 σ v V ratio slightly decreases through the evolution.

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 time-scale, t cool , given by Eq. 4; the fragmentation time-scale, t f rag = M f rag +M s f ġ M f rag , 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 SN and/or AGN; the gas depletion time of the fragmented gas, t s f g = M f raġ M s f g , defined as the time required to convert all the fragmented gas and star-forming gas into stars. As the star-formation timescale is less than all other timescales, t s f g , 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 given in four different stellar mass bins, 10 9 , 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 Tab. 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 Myrs and 300 Myrs. 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 disc 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 helps 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 1 2 3 4 5 6 7 redshift (z) values of the ratio, T 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 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 Myrs. Those values correspond to 10% of the disk orbital timescale. As for the cooling timescale, the fragmentation and the disruption timescales of our starforming 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 SNs 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 simu-lations 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.

Impact of fragmentation and disruption on the life cycle of GMCs
The small differences between the condensation and disruption timescales implies that the 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, those formed in galaxies with short or intermediate disk-scale heights, h d ≤50pc, 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 SN. It is the competition between accretion and disruption which regulates the mass of GMCs (Vázquez-Semadeni et al. 2010). Due to 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.
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).

Gas depletion timescale
Depending on the redshift and the stellar mass of a model galaxy, the gas depletion timescale, t s f g , ranges from 80 Myr to 2 Gyr (Fig. 18). The upper end of this range is consistent with with 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 the stellar mass. Above M ≥ 10 10.5 M and at z ¿ 3.0, the average gas depletion timescale reaches a lower-limit close to 80 Myr.
The t s f g /t f rag 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 reach 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 -5%. Such values are consistent with observational estimates (e.g., Wong & Blitz 2002;Leroy et al. 2008;Murray et al. 2010;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 Myrs of these exchange rates is only a small fraction of the orbital time-scale. 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.

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 pre-stellar 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 multi-phase ISM where 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 time-scale, we measure a large scatter in both the depletion time-scale and the disk orbital time-scale. This proves that prescriptions that rely solely on the disk orbital timescale do not capture properly the complexity of gas cycles in galaxies. Smaller characteristic timescales, such as the cooling and/or energy injection timescales due to supernovae 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 SNs/AGN kinetic energy injection. A fraction of the gas can also be ejected from disks by SN explosions but the characteristic time-scale of ejection is in average 10 times larger than the fragmenting/disruption cycle time-scale. While the SN/AGN kinetic energy injection in the ISM regulates the star-formation rate in low mass galaxies, our model galaxies retain a large fraction of the gas. This "local" gas is then used with a progressively increasing efficiency necessary to build massive, M star > 10 11 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 FUV fluxes ...) predicted by our model have been post-processed to successfully predict the [CII] luminosity functions at high redshifts and 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 photo-ionization 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 darkmatter halos, 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 will could contribute in quenching the radiative cooling and dissipation in the most massive halos.