ABSTRACT

Galaxies display several well-behaved scaling relations between their properties, such as the star formation rate–stellar mass relation (the main sequence, MS) and the stellar mass–halo mass relation (SHMR). In principle, these scaling relations could imply different star formation histories (SFHs) of galaxies and different constraints on galaxy formation physics. In this paper, we derive the SFHs of galaxies by assuming that they always follow the SHMRs at different redshifts and use an empirical model to constrain key processes in their baryon cycle. It is found that, besides cold accretion due to halo growth, outflow of gas produced by stellar feedback has to be recycled to sustain the derived SFHs of galaxies. The recycled fraction is strongly affected by the baryon fraction in accreted low-mass haloes and the mass loading factor that quantifies the ratio between the galactic outflow rate and star formation rate. Our fiducial model predicts that around 20–60 per cent of outflow is recycled in |$\sim 0.5\!-\!4\, \mathrm{Gyr}$|⁠, while simulations predict a slightly higher recycle fraction and a lower recycle time. We argue that strong constraints on the baryon cycle process can be obtained from future observation of the circum-galactic medium (CGM) of galaxies, such as the gas cooling rate of CGM. We also find that the implied SFHs from the SHMRs indicate that galaxies stay on the MS only for part of their lifetimes. Our model reproduces the evolution of the mass–metallicity relation as well.

1 INTRODUCTION

Over the past two decades, huge volume of observation data on properties of galaxies has been obtained from large galaxy surveys, such as the SDSS (York et al. 2000), 2dFGRS (Colless et al. 2001), GOODS (Giavalisco et al. 2004), DEEP (Vogt et al. 2005; Weiner et al. 2005), and COSMOS (Scoville et al. 2007). These data reveal various global scaling relations of galaxies and some of them are closely related to star formation process, such as the star formation rate–stellar mass relation (SFR–M*) for star-forming galaxies (e.g. Oliver et al. 2010; Karim et al. 2011; Leitner 2012; Whitaker et al. 2012, 2014; Speagle et al. 2014; Schreiber et al. 2015), the mean fraction of cold gas (fgasMgas/M*) in the interstellar medium (ISM; e.g. Saintonge et al. 2011, 2017; Boselli et al. 2014; Genzel et al. 2015; Scoville et al. 2017; Tacconi et al. 2018), the mass–metallicity relation (MZR; e.g. Tremonti et al. 2004; Maiolino et al. 2008; Andrews & Martini 2013; Sanders et al. 2018, 2020, 2021; Cullen et al. 2019), and the stellar mass–halo mass relation (SHMR; e.g. Peacock & Smith 2000; Berlind & Weinberg 2002; van den Bosch, Yang & Mo 2003; Yang, Mo & van den Bosch 2003; Mandelbaum et al. 2006; Zheng, Coil & Zehavi 2007; Moster et al. 2010; Yang et al. 2012; Moster, Naab & White 2013; Girelli et al. 2020). These scaling relations are particularly useful to constrain the process of galaxy formation. For example, the SFR–M* relation and the SHMR can separately predict star formation history (SFH) of a galaxy population (e.g. Conroy & Wechsler 2009; Leitner 2012). The gas fraction in ISM and the mass–metallicity relation are strongly affected by the baryon cycle, such as supernovae feedback, gas accretion, outflow by stellar feedback and metal transport. Modelling those scaling relations and understanding their connections remain to be key questions in current theory of galaxy formation.

In the Λ cold dark matter (ΛCDM) framework, it is generally believed that galaxies form at the centre of halos (and subhalos) and their formation is mainly driven by the cooling and condensation of gas in the centre of dark matter halos (White & Rees 1978). Thanks to the advance in numerical methodologies and computing speed, models to follow galaxy evolution in a cosmological context have been developed significantly in the past 20 years, including hydrodynamical simulations (e.g. Hopkins et al. 2014; Vogelsberger et al. 2014a,b; Muratov et al. 2015; Schaye et al. 2015; Wang et al. 2015; Pillepich et al. 2018) and semi-analytic models (SAMs; e.g. Bower et al. 2006; Croton et al. 2006; Kang et al. 2005; Guo et al. 2011; Lagos et al. 2012; Fu et al. 2013; De Lucia et al. 2014; Henriques et al. 2015, 2020; Hirschmann, De Lucia & Fontanot 2016; Lacey et al. 2016).

Undoubtedly, both hydrodynamical simulations and SAMs have achieved great success to match the global properties of galaxies, such as the stellar mass functions (SMF) and some scaling relations mentioned above (for a review, see Somerville & Davé 2015). It is, however, apparent that these models have some limitations and suffer from uncertainties in the physical ingredients of the modelling. Specifically, hydrodynamical simulations rely on assumptions for sub-grid physics, and are limited by numerical resolutions and expensive computational costs, while the SAMs are also becoming enormously complicated, with considerable number of free parameters describing various physical processes. Therefore, the degeneration of model parameters make it difficult to intuitively understand the influence of different physical processes on the evolution of galaxies in a straightforward way. In addition, most models are often tuned to fit the observational scaling relation of galaxy population at a given epoch, usually the SMF at z = 0; it is thus not so apparent to infer how galaxy scaling relations are connected with each other and what are their separate constraints on galaxy SFH.

To investigate if galaxies could evolve along different scaling relations simultaneously and the resulting constraints on their baryon cycle, one first has to obtain SFH of galaxies from their scaling relations. Usually two scaling relations, namely the SHMR and the SFR–M* relation, can be separately used to derive SFHs of galaxies. First, the SHMRs at different redshifts, combined with the average growth of halo mass, can be used to get SFHs of galaxies. This approach was first proposed by Conroy & Wechsler (2009) to derive the star formation rate of galaxies, which are found to be broadly consistent with the data at z < 1. Another one is the observed SFR–M* relation, which can also directly predict the SFHs if galaxies are assumed to follow the SFR–M* relations at different redshifts (Leitner 2012). However, it has been pointed out that the SFHs derived from the SFR–M* relations are inconsistent with the evolution of galaxy SMF, which, in contrast, requires a lower SFR at given stellar mass and an SFR–M* relation that is steeper than observed (Leja et al. 2015; Tomczak et al. 2016; Contini et al. 2017). The resulted SFHs from the SFR–M* relations show that galaxies do not always follow the main sequences (MS) at different redshifts.

In this work, we follow Conroy & Wechsler (2009) to use the SHMRs to track the evolution of galaxies. The main reason for our choice is that the SHMR has been shown as a very fundamental relation for galaxy evolution. Specifically, once the SHMR is fixed, other galaxy properties, such as the SMFs and spatial clusterings, can be well reproduced (e.g. Guo et al. 2010; Moster et al. 2010; Kang et al. 2012). We suggest the readers to refer the review paper by Wechsler & Tinker (2018) for progress on relevant studies of the SHMRs. Our work is different from Conroy & Wechsler (2009) in two aspects. First, we extend our studies to higher redshift as more data at z > 1 is available in recent years, and we aim to study if the resulted SFR–M* relations are consistent with the data across a wide range of cosmic time. Secondly, and most importantly, we use the resulted SFHs to constrain galactic baryon cycle using an empirical model of galaxy formation, which is a simplified but clean version of the comprehensive SAM. Similar approach was adopted by some studies, such as Rodríguez-Puebla et al. (2016). In their work, the SFR for central galaxies can be determined by combining SHMR and the halo mass accretion rate. Then using the empirical model of galaxy formation, with the assumption that the interstellar gas mass is constant for each galaxy, some relations between the net outflow parameter and infall efficiency can be derived. Compared with their work, we pay more attention to the galactic baryon physics, especially the recycle process of the ejected gas. It is well known that galaxy evolution is mainly governed by some key physical processes, such as fresh gas accretion due to halo growth, star formation from cold gas, galactic wind or outflow by stellar feedback, and recycle of gas outflow (e.g. Mo, van den Bosch & White 2010). Numerical studies have independently achieved strong constraints on these processes, such as the gas infall pattern (e.g. Kereš et al. 2005; Brooks et al. 2009; Dekel et al. 2009; Faucher-Giguère, Kereš & Ma 2011; van de Voort et al. 2011), galactic outflow rate (e.g. Muratov et al. 2015; Mitchell et al. 2018, 2020a), and the wind recycle fraction (e.g. Oppenheimer et al. 2010; Christensen et al. 2016; Anglés-Alcázar et al. 2017; Grand et al. 2019; Tollet et al. 2019; Mitchell, Schaye & Bower 2020b). We aim to investigate in this work, by requiring our model to follow the SFHs derived from the observed SHMRs at different redshifts, how the parameters on baryon cycle could be constrained and whether they are consistent with results from simulations.

The structure of this paper is organized as follows: In Section 2, we describe the approach to model galaxy evolution, mainly including two parts: the derivation of galaxy SFH from the SHMRs at different redshifts and an empirical model to describe the evolution of cold gas mass and chemical enrichment. Our main results of the fiducial model are presented in Section 3, including the comparison between predicted star formation rate with the observed MS (Section 3.1), the prediction of gas fraction at high redshift (Section 3.2), the recycle of gas outflow (Section 3.3), and the mass–metallicity relation (MZR) (Section 3.4). We shortly discuss the SFHs of massive galaxies and the constraints on stellar feedback in Section 4. Section 5 summarizes the main results of our work.

Throughout this paper, we assume a standard flat ΛCDM cosmology with |$\Omega _{\rm m}=0.315, \Omega _\Lambda =0.685$|⁠, and |$H_0=67.3 \, \mathrm{km\, s^{-1}\, Mpc}^{-1}$|⁠, and a Chabrier (2003) initial mass function (IMF).

2 METHOD

The details of our method are described in this section. We start with a brief overview of our procedure in Section 2.1. In Sections 2.2 and 2.3, we present the two main ingredients of the method: the SHMR and the average growth of halo/stellar mass for our sample galaxies. A simple empirical model to describe the key physical process of baryon cycle is introduced in Section 2.4.

2.1 Overview

In this part, we briefly introduce the method to derive the SFH of galaxies and how to use them to constrain the physical parameters about the baryon cycle. We select 10 model galaxies with their dark matter halo mass (Mh) evenly distributed between 1011 and |$10^{12}\, \mathrm{M}_\odot$| at z = 0. Note that in Section 4, we will shortly discuss the evolution of massive galaxies, but in most of this work, we focus on galaxies in this mass range for two reasons. The first is that our model neglects merger of dark matter halos (also galaxy mergers) as we try to minimize the uncertainties from modelling galaxy mergers and associated star formation. Previous studies have shown that for galaxies with |$M_{\ast } \lt 10^{11}\, \mathrm{M}_{\odot }$| or |$M_{\rm h} \lt 10^{12}\, \mathrm{M}_{\odot }$|⁠, galaxy mergers contribute very little to the stellar mass growth of galaxies (e.g. Conroy & Wechsler 2009; Moster et al. 2013). Conroy & Wechsler (2009) use a merger model and a no-merger model to derive the SFR of the galaxies, and find that at least for galaxies with |$M_{\ast } \lt 10^{11}\, \mathrm{M}_{\odot }$| and z < 1, the results obtained by the no-merger model are more consistent with the observations (see fig. 7 of their work). Fig. 10 of Moster et al. (2013) also shows that for galaxies with |$M_{\rm h} \lt 10^{12}\, \mathrm{M}_{\odot }$|⁠, the fraction of stellar mass growth in galaxies due to mergers is less than |$5{{\ \rm per\ cent}}$| at z > 0.5, although, at z = 0, this fraction increases to |$\sim 20{{\ \rm per\ cent}}$|⁠. The second reason is that we focus on the predicted evolution of the SFR–M* relation, which is basically held for star-forming galaxies with |$M_* \lt 10^{11}\, \mathrm{M}_{\odot }$|⁠. Throughout this work, the 10 model galaxies are named as G0, G1, G2, ..., G9 with corresponding halo mass at z = 0 as |$M_{\rm h0}=1\times 10^{11}\, \mathrm{M}_\odot , 2\times 10^{11}\, \mathrm{M}_\odot , 3\times 10^{11}\, \mathrm{M}_\odot ,...,1\times 10^{12}\, \mathrm{M}_\odot$|⁠, respectively.

The stellar mass of the model galaxies and their evolution are set using the SHMRs. In this work, the growth of halo mass is determined from a Monte Carlo method, with details to be addressed in Section 2.3. We follow Conroy & Wechsler (2009) to derive the star formation rate of the galaxies based on their halo mass. For example, for the galaxy G0, its progenitor halo mass at a given redshift z is Mh(z), we use the SHMR at redshift z (see equation 1) to get its stellar mass M*(z). By performing the same procedure at a slightly later time-step, say z + dz, we can then obtain a new stellar mass M*(z + dz) based on its new halo mass, Mh(z + dz) at z + dz. We ascribe the stellar mass growth to star formation in the galaxy G0, so it is straight to get the evolution of its star formation rate.

Once the SFH of our model galaxies are obtained, we use a simple empirical model in Section 2.4 to constrain the physical parameters of baryon cycle. As we mentioned before, the SHMRs of galaxies at different redshifts have been shown as best benchmarks for galaxy evolution (e.g. Wechsler & Tinker 2018), it has never been clearly shown whether the implied SFHs can be properly reproduced by a model of galaxy formation. In Section 2.4, we then investigate whether the required baryonic processes, such as gas infall rate, recycled gas fraction from the outflow, are reasonable and consistent with simulation constraints.

2.2 The stellar mass–halo mass relation: SHMR

In this paper, we use the latest SHMR obtained using the abundance matching technique by Girelli et al. (2020). Although there are another two methods to get the SHMR, namely the halo occupation distribution (HOD; e.g. Peacock & Smith 2000; Berlind & Weinberg 2002; Zheng et al. 2007) and the conditional luminosity function (CLF; e.g. van den Bosch et al. 2003; Yang et al. 2003), these two approaches are usually applied to data at low redshift as measurements of galaxy clustering are not reliable at high redshift. The abundance matching technique is straightforward to link galaxy with halo by assuming that a galaxy property monotonically relates to a halo property, such as the stellar mass of the galaxy and the dark matter mass of the halo (e.g. Vale & Ostriker 2004). Introduction on the progress of the SHMR is beyond the scope of this paper, and the readers are referred to the review paper by Wechsler & Tinker (2018). Recently, Girelli et al. (2020) determined the SHMRs up to z = 4 using the COSMOS data, and they adopted the simple double power-law function proposed by Moster et al. (2010):
(1)
where A is the normalization at the characteristic halo mass MA, while β and γ are the slopes of the relation at low- and high-mass ends, respectively. In our work, we use the best-fitting parameters in table 4 of Girelli et al. (2020), and we linearly extrapolate their relations to lower halo/galaxy mass if needed. The scatter in stellar mass at given halo mass is a function of halo mass (e.g. Behroozi, Conroy & Wechsler 2010), and in our work, we follow Girelli et al. (2020) to use a typical constant with σR = 0.2 dex to their best-fitting SHMRs.

2.3 Growth of halo and stellar mass

The formation history of dark matter halo is obtained using the Monte Carlo merger tree code developed by Parkinson et al. (2008). For each model galaxy, from G0 to G9, we produce 100 realizations from the Monte Carlo code and use the average of these realizations as the halo growth history. We then adopt the simple function of Wechsler et al. (2002) to fit the average halo growth history as
(2)
where Mh0 is the halo mass at z = 0 and α is a free parameter used to fit the data. Although any individual halo trajectory may deviate from this one-parameter function, this one-parameter model can indeed provide a pretty good description of the average halo mass accretion history. Our best fitting result is shown in Fig. 1 and the best-fitting parameter ranges from 0.56 to 0.62, which is consistent with the result of Wechsler et al. (2002) for halos in this mass range. In fact, the halo mass growth history can also be obtained by other methods, such as integrating the mean mass accretion rate of dark matter obtained from the simulation, that is, integrating |$\langle {{\rm d}M_{\rm h}}/{{\rm d}t} \rangle$| (e.g. Genel et al. 2008; McBride, Fakhouri & Ma 2009; Fakhouri, Ma & Boylan-Kolchin 2010). After calculation and comparison, we found that there is no obvious difference between our data and those results in our selected halo mass range.
The growth history of halo mass for our 10 model galaxies. The solid lines show the mass evolution predicted from the Monte Carlo merger tree developed by Parkinson, Cole & Helly (2008). The dashed lines are the best fits using equation (2).
Figure 1.

The growth history of halo mass for our 10 model galaxies. The solid lines show the mass evolution predicted from the Monte Carlo merger tree developed by Parkinson, Cole & Helly (2008). The dashed lines are the best fits using equation (2).

With the evolution of both halo mass (equation 2) and the SHMR (equation 1) at hand, it is easy to obtain the stellar mass evolution for each model galaxy. As a check for consistency with previous result, in Fig. 2 we show the evolution of the normalized halo/stellar mass in the left-hand panels and the star formation rate history in the right-hand panel. The left-hand panels are for four model galaxies with increasing mass from the top left- to lower right-hand panel. A general trend seen is that within the mass range of our 10 model galaxies, the halo mass assembles much earlier than the stellar mass, more noticeable for low-mass galaxy. It can also be seen that galaxy assembles earlier with mass increases. The right-hand panel shows the derived SFHs for the 10 model galaxies. It is obvious to find that the peak star formation rate of a high-mass galaxy occurs at high redshift, and for a low-mass galaxy, such as the galaxy G0, its peak star formation has not reached yet at z = 0. The above results obtained from Fig. 2, whether it is about the evolution of halo/stellar mass or about the star formation rate histories, are in broad agreement with previous related studies (e.g. Conroy & Wechsler 2009; Moster et al. 2013; Behroozi, Wechsler & Conroy 2013). In the following sections, we will then use the derived change rate of stellar mass, |$\dot{M}_*(z)$|⁠, to constrain the baryonic process.

Left-hand panels: the halo and stellar mass evolution for four model galaxies. The mass is normalized using the mass at z = 0, which is labelled in each small panel. The red lines are for the stellar mass evolution and the dashed lines for the halo mass. It shows that within the mass range of our 10 model galaxies halo mass grows faster than the stellar mass, and massive galaxies grow more quickly at high redshift. Right-hand panel: the evolution of star formation rates for the 10 model galaxies. Each line is for one galaxy and is normalized by the peak star formation rate of the galaxy during its evolution. The inserted panel shows the absolute star formation rate. This plot shows a strong mass dependence that massive galaxy forms their stars first and has a peak star formation at z ∼ 1. Low-mass galaxy reaches its star formation peak later or is still not reaching its peak until z = 0. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M⊙.
Figure 2.

Left-hand panels: the halo and stellar mass evolution for four model galaxies. The mass is normalized using the mass at z = 0, which is labelled in each small panel. The red lines are for the stellar mass evolution and the dashed lines for the halo mass. It shows that within the mass range of our 10 model galaxies halo mass grows faster than the stellar mass, and massive galaxies grow more quickly at high redshift. Right-hand panel: the evolution of star formation rates for the 10 model galaxies. Each line is for one galaxy and is normalized by the peak star formation rate of the galaxy during its evolution. The inserted panel shows the absolute star formation rate. This plot shows a strong mass dependence that massive galaxy forms their stars first and has a peak star formation at z ∼ 1. Low-mass galaxy reaches its star formation peak later or is still not reaching its peak until z = 0. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M.

2.4 Empirical model for galaxy evolution

Galaxy formation involves very complicated processes with many details yet to be well understood. Never the less, the basic and key processes, such as gas accretion or inflow, cooling of hot gaseous halo, star formation from cold and dense molecular gas, outflow of gas due to stellar wind, and supernova feedback, black hole accretion and associated energetic feedback, etc. are extensively studied and great progress has been made in past years. Analytical descriptions of baryonic processes can be easily found from any typical SAM (e.g. Cole et al. 2000; Kang et al. 2005; Guo et al. 2011).

In this work, we use a simple empirical model to describe the star formation process in our model galaxies. As mentioned before about our selection of the 10 model galaxies, we neglect galaxy mergers and black hole physics in the model. We mainly include the processes about gas cycle, star formation, and its feedback. For star formation rate, we assume a simplified, linear Kennicutt–Schmidt (KS) law (Kennicutt 1998):
(3)
where Mgas is the cold gas content, and SFE is the star formation efficiency, assumed to be a function of galaxy mass and redshift, and in Section 3.2, we present how to set its value. During the evolution of stars, some mass will return to the ISM by stellar wind and supernovae, so the change rate of stellar mass |$\dot{M}_*$| in a galaxy is
(4)
where R is the returned mass fraction from evolving stellar population, which is dependent on the initial SMF (IMF) and can be calculated using the method given by Vincenzo et al. (2016). For our selected Chabrier IMF (Chabrier 2003), we take R = 0.44.

It has been well known that stellar feedback, via both wind and supernova, plays a very important role in galaxy formation and evolution, and all modern cosmological hydrodynamical simulations of galaxy formation have included stellar feedback though the details are quite different (e.g. Vogelsberger et al. 2014a; Schaye et al. 2015; Hopkins et al. 2014; Wang et al. 2015). The major consequence of stellar feedback is the resulted galactic gas outflow that significantly influences the baryon cycle in a galaxy. However, both the relationship between the star formation rate and the outflow mass rate, the so-called mass loading factor, and the fate of the outflow gas are all dependent on the details of modelling stellar feedback in simulations. Even in one particular simulation, the mass loading factor also depends on the properties of the galaxy, and it varies with mass and redshift.

The rate of gas outflow due to stellar feedback is often assumed to be proportional to SFR with a factor depending on galaxy halo mass (or halo velocity) and cosmic time as
(5)
where η is the mass loading factor. Studies using data from cosmological simulations have shown that η depends sensitively on the details of feedback in different models (e.g. Muratov et al. 2015; Christensen et al. 2016; Beckmann et al. 2017; Anglés-Alcázar et al. 2017; Tollet et al. 2019; Nelson et al. 2019b; Mitchell et al. 2020a), and we refer the readers to the paper by Mitchell et al. (2020a) for comparison of η between different cosmological simulations and semi-analytical models. In this work, we use the fitting formula by Muratov et al. (2015) to the FIRE simulation (Hopkins et al. 2014), and it is dependent on halo mass and redshift as
(6)
where Mh60 is the virial mass of a halo with virial velocity (⁠|$V_{\rm vir}=\sqrt{GM_{\rm h}/R_{\rm vir}}$|⁠) equals to 60 km s−1, and Rvir is the virial radius radio of the halo within which the mean density is 200 times the critical density of the Universe.
Now we describe the mass change of the gas reservoir due to gas infall and galactic outflow, which is given by the following equation:
(7)
Here |$\dot{M}_{\rm gas,in}(t)$| is the total gas infall rate, with sources from fresh gas accretion or recycle of outflow. In our model, the gas infall rate is not a free parameter, but is obtained from the above equations. As explained in Section 2.3, once the evolution of stellar mass change rate, |$\dot{M}_*(t)$|⁠, is determined for each model galaxy, the evolution of star formation rate SFR(t) can be obtained by equation (4). Then assuming an SFE, the total cold gas content, Mgas(t), and outflow rate can be calculated from equations (3) and (5), respectively. Thus the gas infall rate can be calculated from equation (7).

We will compare the derived total gas infall rate obtained from equation (7) with model predictions. One main source of gas infall in a galaxy is the fresh gas accretion due to the growth of its dark matter halo. In the classical model (e.g. White & Rees 1978), the newly accreted gas is shocked heated to the virial temperature of the halo, and it then cools down via radiative cooling. This kind of process is called as hot mode accretion. Later studies found that the accretion is not always in hot mode. In low-mass galaxies(⁠|$M_{\rm h} \lt 10^{12}\, \mathrm{M}_{\odot }$|⁠), newly accreted gas from the external halo is not shock heated, but is directly accreted to the halo centre. This cold accretion is rapid and limited by the free-fall time of the halo (e.g. Kereš et al. 2005; Croton et al. 2006; Brooks et al. 2009; Dekel et al. 2009; Faucher-Giguère et al. 2011; van de Voort et al. 2011). The other source of gas infall is the re-accretion/recycle of the outflow. Recent studies using hydrodynamical simulation of galaxy formation have shown that the fraction of recycled gas from outflow is dependent on the properties of the outflow, and it is diverse among different simulations (e.g. Oppenheimer et al. 2010; Christensen et al. 2016; Anglés-Alcázar et al. 2017; Grand et al. 2019; Tollet et al. 2019; Mitchell et al. 2020b). In Section 3.3, we will compare our predictions on gas recycle with simulation results.

We also predict how metals evolve in a galaxy. The production of metal is accompanied by the formation and evolution of stars. Then newly produced metal is returned to the ISM by stellar feedback or stellar wind, and it is assumed to uniformly mix with the ISM. A yield parameter yZ is often used to characterize this process, which is the ratio of the total mass of metals that a stellar population releases into the ISM and the remained stellar mass of that stellar population. This can be intuitively understood as
(8)
where the |$\dot{M}_{{\rm new}\_{\rm metal}}$| is the metal mass growth rate in ISM due to the metal production. Obviously, the value of yZ is also dependent on the selected IMF. For Chabrier IMF (Chabrier 2003), we take yZ as 0.06 from the work by Vincenzo et al. (2016). For more definition or calculation details of R and yZ, we refer readers to Vincenzo et al. (2016). The metallicity of the ISM is defined as
(9)
Here MZ(t) is the total metal content in the ISM. The metal will also be reduced by outflow or diluted by the accretion of cold gas, so it can be calculated as
(10)
Here we assume the metallicity of new infall gas is Zin and it is assumed to be primordial with Zin = 10−4 in our fiducial model, but we will also show predictions of the mass–metallicity relation with different Zin in Section 3.4.

3 PREDICTIONS OF THE FIDUCIAL MODEL

In this section, we present results from our fiducial model. In Section 3.1, we show the trajectories of model galaxies on the MS relations at different redshifts, and we predict the gas fraction at high redshift in Section 3.2. Constraints on the process of gas recycle are shown in Section 3.3. We present the evolution of the mass–metallicity relation (MZR) in Section 3.4.

3.1 Trajectories of model galaxies on the MS relations

Observational studies have found a tight relation between SFR and stellar mass of star-forming galaxies (SFR–M*), often called the ‘Main Sequence (MS)’, and this relation is observed to hold from z = 0 to 6 (e.g. Oliver et al. 2010; Karim et al. 2011; Leitner 2012; Whitaker et al. 2012, 2014; Speagle et al. 2014; Schreiber et al. 2015; Popesso et al. 2023). As the MS relation is a snapshot of star-forming galaxies at one cosmic epoch, one may wonder how a given galaxy evolves along the MS relations at different redshifts. One simple assumption is that ‘A star-forming galaxy is always a star-forming galaxy’, and it then predicts a unique SFH of galaxy using the MS relations at different redshifts. This technique is called as the ‘Main Sequence Integration (MSI)’ (e.g. Leitner 2012; Ciesla, Elbaz & Fensch 2017). Here we show the trajectories of our model galaxies on the MS diagram and investigate whether they agree with the MSI scenario.

Before showing the trajectories of model galaxy in the MS diagram, one should bear in mind two facts. The first is that the MS relation is only for star-forming galaxy, and quenched galaxies are usually not included. The second is that measurements of both star formation rate and stellar mass are tricky, as they depend on star formation rate indicator and techniques to derive the stellar mass of galaxies. To ensure a fair comparison with data, here we use the best-fitting MS relations from three studies, which are Leitner (2012), Speagle et al. (2014), and Popesso et al. (2023). In Leitner (2012), the fitting was done separately to the radio survey data from Karim et al. (2011, hereafter K11) and the IR survey data from Oliver et al. (2010, hereafter O10). Here we use their fitting result for O10 data. Although for O10 data, the fitting is only for galaxies with stellar mass of |$\sim 10^{11}\, \mathrm{M}_\odot$|⁠, compared with the fitting results of K11 data, the fitting result using O10 data is more in line with the latest observational result at the low mass end. Speagle et al. (2014) converted observational data from 25 studies into a common set of calibrations and gave the best-fitting MS relations within the fitted stellar mass range of |$10^{9.7}\!-\!10^{11.1}\, \mathrm{M}_\odot$|⁠. A recent work by Popesso et al. (2023) also did similar analysis by combing data from literature and obtained the best-fitting MS relations in the widest range of redshift (0 < z < 6) and stellar mass (⁠|$10^{8.5}\!-\!10^{11.5}\, \mathrm{M}_\odot$|⁠), and we plot one set of their best fitting formula in the right-hand panel of fig. 3. In the upper row of Fig. 3, we show the evolution of the 10 model galaxies on the MS diagram, and we compare them with data from the aforementioned three observational work in the left-, middle, and right-hand panels, respectively. In each panel, the straight colour lines are the best-fitting MS relations to the data at different redshifts, and the curved colour lines are the trajectories of the 10 model galaxies. Each curved colour line is for one galaxy (from G0 to G9, see labels in the lower row) and the circles along the line mark the positions of the galaxy on the MS diagram at different redshifts. For example, the curved yellow line (the lowest line in each panel) is for the galaxy G0, and the most left circle show its stellar mass and star formation rate at z = 2, and the next circle shows the position at z = 1.5, and so on. It is seen that, compared with data in any of the three panels, our model galaxies are below the MS relations at high redshift and are above the MS relations at low redshift.

The evolution of model galaxies on the MS diagram, or the SFR–M* relation. Upper panels: The curved lines are the trajectories of the 10 model galaxies (each line for one galaxy; see the label in lower right-hand panel), and the circles along each curved line show the star formation rate and stellar mass at several redshifts between z = 0 and 2, indicated by the label in upper panels. The straight colour lines are the best-fitting relations to the observational data: Leitner (2012) , left-hand panel), Speagle et al. (2014, middle panel), and Popesso et al. (2023, right-hand panel) at different redshifts, also indicated by the label in each panel. Lower panel: the offset between the predicted star formation rate and the best-fitting MS relations. It is seen that a galaxy is not always a star-forming galaxy. For galaxies in the mass range ($10^{11}\!-\!10^{12}\, \mathrm{M}_{\odot }$), they have lower star formation rate than the data at high redshift, but are above the MS at z ∼ 0. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M⊙.
Figure 3.

The evolution of model galaxies on the MS diagram, or the SFR–M* relation. Upper panels: The curved lines are the trajectories of the 10 model galaxies (each line for one galaxy; see the label in lower right-hand panel), and the circles along each curved line show the star formation rate and stellar mass at several redshifts between z = 0 and 2, indicated by the label in upper panels. The straight colour lines are the best-fitting relations to the observational data: Leitner (2012) , left-hand panel), Speagle et al. (2014, middle panel), and Popesso et al. (2023, right-hand panel) at different redshifts, also indicated by the label in each panel. Lower panel: the offset between the predicted star formation rate and the best-fitting MS relations. It is seen that a galaxy is not always a star-forming galaxy. For galaxies in the mass range (⁠|$10^{11}\!-\!10^{12}\, \mathrm{M}_{\odot }$|⁠), they have lower star formation rate than the data at high redshift, but are above the MS at z ∼ 0. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M.

The derivation of the model galaxies from the MS relation can be more clearly seen in the lower row of Fig. 3 where we show the offset between our predicted star formation rate and the observed MS, as a function of redshift. As the typical scatter of the observed MS relation is around 0.2 dex with weak dependence on galaxy mass, we plot two horizontal lines to indicate the scatter. It is found that our results agree better with the MS relations from Leitner (2012), and it shows that most of our model galaxies selected at z = 0 were star-forming galaxies between z = 1.5 and 0.5, and their star formation rates were close to the MS relation. Below z = 0.5 the model galaxies are still above the MS relations. Larger deviation is seen in the left-hand panel where it shows that our model galaxies all have star formation rate lower than the data at z > 1. The results show that a star-forming galaxy is not always a star-forming galaxy, and it remains on the MS for only a limited lifetime.

The predicted lower star formation rate than the data at z > 1.5 agrees with previous conclusion using similar methods (e.g. Conroy & Wechsler 2009; Leja et al. 2015). In fact, similar discrepancies are also seen in SAMs or hydrodynamical simulations, that is, at middle and high redshifts, the MS obtained in theoretical models is systematically lower than the observed MS by ∼0.2–0.5 dex (e.g. Kang et al. 2010; Mitchell et al. 2014; Sparre et al. 2015; Furlong et al. 2015; Tomczak et al. 2016; Donnari et al. 2019). It is recently suggested that the tension could be alleviated if similar methods are used to both data and model galaxy to derive the MS relations (Katsianis et al. 2020). Leja et al. (2022) also show that, if a more flexible neural network is adopted to process the data, the normalization of the MS relation can be lowered by ∼0.2–0.5 dex, which is consistent with predictions from theoretical models of galaxy formation.

Before we conclude the discussion on the MS relation, we note that here we only show the trajectories of our model galaxies with halo mass lower than |$10^{12}\, \mathrm{M}_{\odot }$|⁠. It is well realized that low-mass galaxies form later and are now in the process of rapid star formation. This is why we found they are all above the MS at z = 0. In Section 4 we will show the MS relations of more massive galaxies for comparison, and we will see that massive galaxies are basically quenched with star formation rate lower than the MS relations at z ∼ 0.

3.2 Evolution of the cold gas fraction

Using our model, we can also predict the evolution of cold gas fraction, defined as |$\mu _{\rm gas}={M_{\rm gas}}/{M_{*}}$|⁠, for the model galaxies. Using the above derived SFR for each galaxy, the cold gas amount and its evolution can be calculated from equation (3) once the SFE is specified. In principle, SFE could be a function of galaxy mass and redshift, written as SFE(M*, z). In most models, SFE is intrinsically related to the free-fall time or the dynamical time of the galaxy, see the discussion in Somerville & Davé (2015). Here we first use the fitting formula shown in table 3 of Boselli et al. (2014) to get the gas fraction μgas(z = 0) and Mgas(z = 0). Then the value of SFE(M*, z = 0) is determined using our equation (3). It is then found that the SFE has a weak mass dependence as |$SFE\sim M_{*}^{0.5}$|⁠. This mass dependence is the same as that by Schreiber et al. (2016). In addition, observational work have found the redshift dependence of SFE (e.g. Genzel et al. 2015; Scoville et al. 2017; Darvish et al. 2018; Tacconi et al. 2018). Genzel et al. (2015) combined the molecular gas data between z = 0 and 3, and found that the gas depletion time-scale tdepl has a weak redshift dependence as tdepl ∼ (1 + z)−0.3. As SFE is inverse to the gas depletion time, we use a combined |$SFE\sim M^{0.5}_{\ast }(1+z)^{0.3}$| in our model and predict the cold gas content using equation (3).

Fig. 4 shows the evolution of cold gas mass in each model galaxy. The left-hand panel shows the amount of cold gas and the right-hand panel is for its fraction. As one can see from the left-hand panel that for most galaxies, the cold gas content remains almost constant during their evolution, indicating that the gas depleted by star formation and outflow is basically balanced by fresh gas infall. For massive galaxy, there is a steady decrease of cold gas with decreasing redshift, showing that the supply of cold gas is lower than the depletion. This is in general consistent theoretical prediction that quick gas accretion is maintained by cold accretion in low-mass galaxy and massive galaxy starts to quench at low redshift. The right-hand panel shows the gas fraction as a function of redshift. It is found that almost all galaxies were gas-rich at high redshift, and the fraction is higher for low-mass galaxies. For example, for the galaxy G0, its cold gas mass is about 50 times of its stellar mass at z = 2.

The predicted evolution of cold gas mass (left-hand panel) and the cold gas fraction, defined as $\mu _{\rm gas}={M_{\rm gas}}/{M_{*}}$ (the right-hand panel) for the 10 model galaxies from our empirical model. It can be found that for low-mass galaxies, their gas mass remain almost unchanged since z = 2 and they were very gas rich in the past, with gas mass about 10 times higher than stellar mass at z = 2. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M⊙.
Figure 4.

The predicted evolution of cold gas mass (left-hand panel) and the cold gas fraction, defined as |$\mu _{\rm gas}={M_{\rm gas}}/{M_{*}}$| (the right-hand panel) for the 10 model galaxies from our empirical model. It can be found that for low-mass galaxies, their gas mass remain almost unchanged since z = 2 and they were very gas rich in the past, with gas mass about 10 times higher than stellar mass at z = 2. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M.

In Fig. 5, we compare the predicted cold gas fraction with the TNG50 simulation (Nelson et al. 2019a,b; Pillepich et al. 2019) and the observations at z = 0 and 2. The small scatter points are from the TNG50-1 simulation,1 with colour lines to show the median and the error bars to show the 40 percentage from the median. The dash–dotted line is the best-fitting line by Genzel et al. (2015) to the observed molecular gas fraction in galaxies with |$M_{\ast } \gt 10^{9.8}\, \mathrm{M}_{\odot }$| at z = 2. We extrapolated their fit to low-mass galaxies. Our model galaxies at z = 0 and  2 are shown as filled colour circles, connected by black solid and dotted lines, respectively. As mentioned before, our gas fraction at z = 0 is set using the observational results by Boselli et al. (2014). This panel shows that the gas fraction predicted from the TNG50-1 is similar to the data at z = 0, but there is almost no evolution with redshift in the simulation. The gas fraction at z = 2 from our model is slightly higher than the data of Genzel et al. (2015). We note that the data at z = 2 are only for molecular gas, which should be treated as a lower limit of the cold gas fraction. The thick dashed line is the universal gas fraction, which is higher than our model prediction by a factor of 2 at z = 2. Overall, our results are consistent with the data in the sense that most galaxies were gas-rich at early times.

The gas fraction as a function of stellar mass at z = 2 and 0. The model predictions are shown as colour circles connected by the black solid line (z = 0), the black dashed line (z = 2). The small dots are from the TNG50 simulation at z = 2 (blue dots) and 0 (orange dots), with two lines with error bars showing the medium and the scatter. The dash–dotted line is the best fit to the data at z = 2 for the molecular gas fraction (Genzel et al. 2015), and we extrapolate it to low-mass galaxy. The thick dashed line is a theoretical prediction by assuming that total baryon mass (star plus cold gas) in a model galaxy is the universal baryon fraction, $f_{\rm b}={\Omega _{\rm b}}/{\Omega _{\rm m}}$. It is found that the predicted gas fraction at z = 2 is lower than the universal baryon fraction by a factor of 2, and slight higher than the data of Genzel et al. (2015), which should be treated as a lower limit (only molecular gas included).
Figure 5.

The gas fraction as a function of stellar mass at z = 2 and 0. The model predictions are shown as colour circles connected by the black solid line (z = 0), the black dashed line (z = 2). The small dots are from the TNG50 simulation at z = 2 (blue dots) and 0 (orange dots), with two lines with error bars showing the medium and the scatter. The dash–dotted line is the best fit to the data at z = 2 for the molecular gas fraction (Genzel et al. 2015), and we extrapolate it to low-mass galaxy. The thick dashed line is a theoretical prediction by assuming that total baryon mass (star plus cold gas) in a model galaxy is the universal baryon fraction, |$f_{\rm b}={\Omega _{\rm b}}/{\Omega _{\rm m}}$|⁠. It is found that the predicted gas fraction at z = 2 is lower than the universal baryon fraction by a factor of 2, and slight higher than the data of Genzel et al. (2015), which should be treated as a lower limit (only molecular gas included).

3.3 Constraints on the gas cycle

In this section, we present our model prediction on the gas accretion rate, the recycle fraction of gas outflow by stellar feedback, and compare them with constraints from simulations.

3.3.1 The gas infall rate

In Section 2.4, we mentioned that the total gas accretion rate, |$\dot{M}_{\rm gas,in}$|⁠, in a galaxy is calculated using equation (7). This accretion rate is mainly contributed by two sources. One is the gas accretion due to halo growth and the other is the recycle of gas outflow. In our selected galaxy mass range, the gas accretion due to halo growth is in the regime of cold accretion, which is fast and determined by the free-fall time. We approximate this cold accretion rate as
(11)
where fb(Mh) is the baryon fraction in the accreted dark matter halo and tff is the free-fall time of the galaxy. For simple consideration, we set fb as the universal baryon fraction, thus |$f_{\rm b}(M_{\rm h}) = {\Omega _{\rm b}}/{\Omega _{\rm m}} =0.15$|⁠, although the baryon fraction could be influenced by other effects, such as the re-ionization effect that is dependent on the mass of the accreted halo (e.g. Gnedin 2000). We approximate the free fall time as the dynamic time of the halo, which is calculated as tff = tdyn = Rvir/Vvir. The gas accretion rate from the recycle of outflow is then calculated as
(12)
In Fig. 6, we compare the total gas infall rate (solid lines) with the cold accretion rate (dashed lines) for four representative model galaxies, which are galaxy G0, G1, G2, and G9. It is seen that the required total gas infall rate is almost always higher than the cold accretion rate except for the low-mass galaxy, G0, for which the cold accretion is higher than the infall rate at z > 1.1. For most galaxies, the total infall rate is a factor of 2 ∼ 3 higher than the cold accretion rate, showing that the cold accretion is not high enough to sustain the star formation in the galaxy; thus, additional supply is provided by the recycle of ejected gas. In the following section, we study the recycle process in more detail.
The evolution of the total gas infall rate (solid lines) and cold gas accretion rate due to halo growth (dashed lines) for four model galaxies. The difference between the two rates is the accretion rate from the recycle of gas outflow. It is seen that for most galaxies, the total gas infall rate is larger than the cold accretion, indicating that the gas recycle has to be effective. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G9, 1 × 1012) M⊙.
Figure 6.

The evolution of the total gas infall rate (solid lines) and cold gas accretion rate due to halo growth (dashed lines) for four model galaxies. The difference between the two rates is the accretion rate from the recycle of gas outflow. It is seen that for most galaxies, the total gas infall rate is larger than the cold accretion, indicating that the gas recycle has to be effective. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G9, 1 × 1012) M.

3.3.2 Constraints on gas recycle from the outflow

The fate of gas outflow is key for galaxy formation, but it is difficult to model. In early SAMs, it was assumed that the gas outflow is heated to the virial temperature of the halo and mixed with the hot gaseous halo for further cooling (e.g. Springel et al. 2001). A later version of SAMs assumed that part of the outflow remains in the halo, and part leaves the halo and is later re-incorporated in due time, both of which are dependent on halo virial velocity (e.g. Guo et al. 2011; Henriques et al. 2015; Somerville & Davé 2015). However, hydrodynamical simulations have found that the gas outflow is in multi-phase thermal state and the cycle of gas is very complicated with rich diversity among different simulations (e.g. Oppenheimer et al. 2010; Christensen et al. 2016; Anglés-Alcázar et al. 2017; Tollet et al. 2019; Grand et al. 2019; Mitchell et al. 2020b). In this section, we compare our model constraints on the recycle of gas outflow with results from both SAMs and simulations.

We first show the evolution of the outflow rate and the gas recycle rate for our 10 model galaxies in Fig. 7. Each panel is for one galaxy with mass increasing from the upper left to the lower right. In each panel, the vertical lines indicate the peaks for the outflow and the recycle rate. As we can see for all model galaxies, their outflow rates (dashed lines) increase with the look back time and reach a peak roughly at 8 ∼ 10 Gyr ago. This behavior is similar to the evolution of the star formation rate shown in Fig. 2 except for low-mass galaxies. For the low-mass galaxies, their star formation has not reached the peak at z = 0, but due to their low virial velocity at high redshift, the outflow rate was also higher in the past. The solid lines are for the gas recycle rate. We note that here the recycle rate at a given time should not be considered as an instantaneous rate of gas recycled from the outflow produced at the same time, it is the recycle rate of gas from the cumulative outflow in the past.

Evolution of gas outflow rate and recycle rate from the outflow. Each panel is for one galaxy, from G0 to G9. The dashed line is for the gas outflow rate, $\dot{M}_{\rm gas,out}$, and the solid line is for gas recycle rate $\dot{M}_{\rm gas,rec}$. The vertical lines mark the peak of the two lines and their time delay (Δt) can be regarded as the gas recycle time and is labelled in each panel. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M⊙.
Figure 7.

Evolution of gas outflow rate and recycle rate from the outflow. Each panel is for one galaxy, from G0 to G9. The dashed line is for the gas outflow rate, |$\dot{M}_{\rm gas,out}$|⁠, and the solid line is for gas recycle rate |$\dot{M}_{\rm gas,rec}$|⁠. The vertical lines mark the peak of the two lines and their time delay (Δt) can be regarded as the gas recycle time and is labelled in each panel. To make it easier to understand the label, Mh0 = (G0, 1 × 1011), (G1, 2 × 1011), (G2, 3 × 1011), (G3, 4 × 1011), (G4, 5 × 1011), (G5, 6 × 1011), (G6, 7 × 1011), (G7, 8 × 1011), (G8, 9 × 1011), (G9, 1 × 1012) M.

Fig. 7 shows two important features. One is that the recycle rate is always lower than the outflow rate, indicating that only part of the gas from the outflow is recycled to supply the star formation. We define the total recycle fraction as
(13)
In our model, we start the evolution from the initial redshift z = 2, and this equation gives the total recycled fraction of gas outflow by z = 0. In fact, we find that the results are only slightly changed if we start the evolution from z = 4. The other feature is that there is a delay between the peak of outflow rate and the recycle rate. We label the time delay in each panel, and we can see that it is small for massive galaxies and being larger for low-mass galaxies. For the lowest mass galaxy G0, its peak of recycle is not reached yet at z = 0. Although we do not model the recycle process in our model, the time delay between the two curves can be viewed as a typical recycle time of the gas outflow. This can be verified in a simple way. If we shift the two curves in each panel to overlap their peaks, we find that they have similar shapes and the ratio between their maxima is very close to the total recycle fraction calculated from equation (13).

In Fig. 8, we compare the predicted total gas recycle fraction and recycle time from our model with other studies. In the left-hand panel we show the recycle fraction as a function of galaxy halo mass. The filled circles connected with a black solid line are the predictions of our fiducial model, and the other points are from a few simulation results, see the caption for references. It is seen that our results are marginally consistent with the simulation results. Both the NIHAO (Tollet et al. 2019) and the FIRE (Anglés-Alcázar et al. 2017) simulations indicate the recycle fraction has a peak around |$M_{\rm vir} \sim 10^{12}\, \mathrm{M}_{\odot }$|⁠, sightly higher than our peak mass at around |$6 \times 10^{11}\, \mathrm{M}_{\odot }$|⁠. The results from Christensen et al. (2016) do not have similar behavior, and their gas recycle fraction is between 20 and |$70{{\ \rm per\ cent}}$| in the plotted mass range. Overall, our gas recycle fraction is lower than the simulation results. The right-hand panel shows the recycle time for the gas outflow. Here we include results from a few other simulations: the Auriga Simulation Grand et al. (2019), the EAGLE simulation Mitchell et al. (2020b). The red dash–dotted line is the prediction by (Henriques et al. 2013) using the L-Galaxies model. It is found that the model, either our empirical model or the SAM of (Henriques et al. 2013), overpredicts the gas recycle time for low-mass galaxies. The typical recycle times from simulations are between 0.5 and 1.5 Gyr, with a weak dependence on mass than the models. Overall, this panel shows that the recycle time is short in massive galaxy. This has been pointed out by Oppenheimer et al. (2010) that the deeper gravity potential in massive galaxies make the gas recycle faster.

The total recycled mass fraction from the outflow since z = 2 (left-hand panel) and the recycle time (right-hand panel). Both the black solid lines in two panels, connected with colour circles, are our fiducial model predictions with fb = 0.15. The dashed lines with colour circled connected are the model predictions if the baryon fraction in newly accreted halo is half of the universal baryon fraction, fb = 0.075. All other lines or points are results from hydrodynamical simulations and the semi-analytical model; see the label and text for references.
Figure 8.

The total recycled mass fraction from the outflow since z = 2 (left-hand panel) and the recycle time (right-hand panel). Both the black solid lines in two panels, connected with colour circles, are our fiducial model predictions with fb = 0.15. The dashed lines with colour circled connected are the model predictions if the baryon fraction in newly accreted halo is half of the universal baryon fraction, fb = 0.075. All other lines or points are results from hydrodynamical simulations and the semi-analytical model; see the label and text for references.

In principle, hydrodynamical simulation can trace gas dynamics more accurately; thus, it is very likely that both our model and the L-Galaxies model of Henriques et al. (2013) do overestimate the gas recycle time compared to the simulation results. In our model, the high recycle time arises from the low recycle fraction of gas from the outflow. Seen from equation (12) that the gas recycle rate is the difference between the total gas infall rate and the cold accretion rate, and thus our underestimation of the gas recycle fraction is due to the overestimation of cold accretion rate calculated from equation (11). Indeed, it has been well-known that the cosmic re-ionization reduces baryon fraction in low-mass halos (e.g. Gnedin 2000). Other than the cosmic re-ionization effect, galactic supernova feedback also disrupts the incoming cold flow and reduces the gas infall near the halo virial radius (e.g. Mitchell et al. 2018, 2020b). Tollet et al. (2019) used the NIHAO simulation to show that galactic outflow can disrupt the infall gas up to a distance of ∼6Rvir from the host galaxy, and it can reduce the baryon fraction to half of the universal one for halo with mass |$\sim 3 \times 10^{10}\, \mathrm{M}_{\odot }$|⁠. They named this effect as pre-emptivefeedback. We do not attempt to model both the cosmic re-ionization and pre-emptivefeedback effects in this work, but to simply illustrate how these effects will impact our results on the gas recycle process. We test it by setting the baryon fraction in all newly accreted halos as 0.5fb in equation (11). In Fig. 8, we show this test using the black dashed lines in both panels. It is seen that decreasing the incoming baryon fraction from cold accretion will indeed result in a higher recycle fraction of gas from outflow (left-hand panel) and a shorter recycle time (right-hand panel). The agreement with simulations is now better, although the mass dependence in our model is still stronger than the trend in simulations.

3.4 The mass–metallicity relation (MZR)

In this section, we compare the predicted mass–metallicity relation (MZR) with the data at z = 0. As shown by many studies (e.g. Tremonti et al. 2004), the gas-phase metallicity depends critically on the past SFH, the gas outflow, and inflow of the galaxy. Therefore, the MZR can be utilized to constrain these processes in galaxy formation and evolution. In our model, the metal content is calculated using equation (10). We note that, to proceed with that equation, one has to first set the initial metallicity of the gas. As we start the evolution of our model galaxies at z = 2, we use the best-fitting MZR at z = 2 by Maiolino et al. (2008), which is very similar to the recent results of Sanders et al. (2020). The metallicity of infall gas is also a free parameter, and in our model, we test a few assumptions on its value:

  • The metal content in the total gas infall is set as primordial metallicity with Zin = 10−4, no matter where the infall gas is from, either the cold accretion or the recycled gas from the outflow.

  • The metallicity of gas from cold accretion is still primordial, but all metals from the outflow are recycled instantaneously, no matter how much gas is recycled. It means that metal is actually not ejected out, but only the gas by stellar feedback.

  • Only part of the metal in the outflow gas is recycled, and the fraction is the same as the recycled fraction of gas in the outflow.

With the above assumptions, we can calculate the evolution of metallicity in each model galaxy. To compare with the data, we adopt the converting formula in Fu et al. (2013):
(14)
where 8.69 is the solar value in terms of 12 + log10(O/H) (Asplund et al. 2009) and Z = 0.02.

The predicted MZR relations at z = 0 are shown in Fig. 9, in which the three panels are for the above three different assumptions of the metal content in the infall gas, respectively. Different lines are from different observational data (e.g. Tremonti et al. 2004; Maiolino et al. 2008; Andrews & Martini 2013; Sanders et al. 2021). Our model results are shown as filled circles. The upper panel shows that the predicted MZR is below all the data if the infall gas has primordial metal content. Thus, the metal from the gas outflow must be recycled. The middle panel shows that if all metal in the outflow is recycled, the predicted MZR will be higher than the data. This indicates that metal loss must be included in the model. In the lower panel in which the recycle fraction of metal is the same as the recycled gas from the outflow, the predicted MZR is roughly consistent with the data, albeit with a slightly lower normalization. This is linked to our underestimation of gas recycle fraction. We have tested that if a slightly higher recycle fraction about 70–80 per cent is used for the metal recycle, the predicted MZR agrees better with the data, but we omit the plot here.

The mass–metallicity relation at z = 0. The model predictions are shown as colour circles. The three panels are for three different assumptions on the metallicity of the infall gas; see the text for details. All lines and shaded region are from different data: the solid, dashed, and dotted gray lines show the median, $68{{\ \rm per\ cent}}$ contour, and $95{{\ \rm per\ cent}}$ contour, respectively, of the MZR at z ∼ 0.1 (Tremonti et al. 2004). The green dash–dotted line is at z = 0.07 (Maiolino et al. 2008), the purple dotted line is at z = 0.027 ∼ 0.25 (Andrews & Martini 2013), and the brown dash–dotted line is at z ∼ 0 (Sanders et al. 2021).
Figure 9.

The mass–metallicity relation at z = 0. The model predictions are shown as colour circles. The three panels are for three different assumptions on the metallicity of the infall gas; see the text for details. All lines and shaded region are from different data: the solid, dashed, and dotted gray lines show the median, |$68{{\ \rm per\ cent}}$| contour, and |$95{{\ \rm per\ cent}}$| contour, respectively, of the MZR at z ∼ 0.1 (Tremonti et al. 2004). The green dash–dotted line is at z = 0.07 (Maiolino et al. 2008), the purple dotted line is at z = 0.027 ∼ 0.25 (Andrews & Martini 2013), and the brown dash–dotted line is at z ∼ 0 (Sanders et al. 2021).

4 DISCUSSIONS

In this section, we briefly discuss a few cases that we do not include in previous sections. One is the mass range of our model galaxy. In principle, our model can be extended to more massive halos, but as we explained before, galaxy mergers and AGN activities are becoming more effective in massive galaxies, so they hamper us to constrain the baryonic cycle from our simple empirical model. We will shortly show the predicted MS of massive galaxies and whether they agree with the data. Another effect to discuss is the gas outflow by stellar feedback. In our fiducial model, we use the fitting formula of Muratov et al. (2015) for the mass loading factor. As this is a very key parameter governing galaxy formation and most SAMs take it as a free parameter, we will discuss if our empirical model can put some constraints on this factor.

We have seen from Fig. 3 that the star formation rates of our model galaxies (G0–G9) are higher than the MS data at z = 0, but are lower than the data at z > 1.5. This seems to contradict with the common perception that galaxies gradually quench at lower redshifts. In fact, evolution of a galaxy on the MS diagram depends on its mass. To demonstrate this, we select 10 model galaxies with halo mass between |$2\times 10^{11}$| and |$10^{13}\, \mathrm{M}_{\odot }$|⁠, and we show the offset between the predicted star formation rate and the observed MS relations in Fig. 10. The mass of each model galaxy is labelled in the plot. Note that here we only show the comparison with the data of Speagle et al. (2014). The other two fitting formulae used in the paper (Leitner 2012; Popesso et al. 2023) produce slightly different results, but the overall trends are similar.

The offset between the predicted star formation rate and that of the observed MS by Speagle et al. (2014) for the 10 model galaxies with current (z = 0) halo mass ranging from $2\times 10^{11}$ to $10^{13}\, \mathrm{M}_{\odot }$. The dashed lines indicate the typical scatter of 0.2 dex. It is clearly seen low-mass galaxies are recently becoming star forming. Massive galaxies were star forming in the past (z > 2), and are now quenched. More massive galaxy quenched first, in agreement with the ‘downsizing’ scenario. Note that the upturns at z < 0.5 for massive galaxies are artificial due to our neglect of galaxy mergers in massive haloes, see the text for details.
Figure 10.

The offset between the predicted star formation rate and that of the observed MS by Speagle et al. (2014) for the 10 model galaxies with current (z = 0) halo mass ranging from |$2\times 10^{11}$| to |$10^{13}\, \mathrm{M}_{\odot }$|⁠. The dashed lines indicate the typical scatter of 0.2 dex. It is clearly seen low-mass galaxies are recently becoming star forming. Massive galaxies were star forming in the past (z > 2), and are now quenched. More massive galaxy quenched first, in agreement with the ‘downsizing’ scenario. Note that the upturns at z < 0.5 for massive galaxies are artificial due to our neglect of galaxy mergers in massive haloes, see the text for details.

As we can from the plot, at z > 2, low-mass galaxies (⁠|$M_{\rm h} \lt 10^{12}\, \mathrm{M}_{\odot }$|⁠) had star formation rate lower than the MS relation and began to cross the MS at low redshift. This has already been shown in Fig. 3. However, for massive galaxies (⁠|$M_{\rm h} \gt 10^{12}\, \mathrm{M}_{\odot }$|⁠), they were star forming at high redshift and more massive galaxies are more closer to the MS relation. The most massive galaxy (⁠|$M_{\rm h} = 10^{13}\, \mathrm{M}_{\odot }$|⁠) was actually above the MS at z > 3. It is interesting to see that massive galaxies began to depart from the MS and became quenched at low redshift. The larger the mass is, the earlier is the quenching time. This is consistent with the ‘downsizing’ scenario that massive galaxies form first and quench earlier. This is also consistent with the recent studies on the evolution of the MS relation about the ‘bending’ mass, defined as the galaxy mass where the MS relation becomes flat, is also decreasing with redshift and it is around |$M_{\rm h} \sim 10^{12}\, \mathrm{M}_{\odot }$| at z = 0 (e.g. Behroozi et al. 2013; Daddi et al. 2022). This ‘bending’ mass is interpreted as the transitional mass where gas accretion is switched from cold accretion to hot accretion. We note that massive galaxies in our model show upturns of their star formation rates relative to the MS relations at z < 1. This behavior is artificial. We have mentioned that we neglected galaxy mergers in our model, so all stellar mass growth is ascribed to star formation which might be overestimated. Conroy & Wechsler (2009) and Moster et al. (2013) also find that for massive galaxies at low redshifts, the contribution of galaxy mergers to the stellar mass growth of galaxies is non-negligible.

Now we discuss the effect of changing the mass loading factor, η. The mass loading factor we used in equation (6) is from the best-fitting relation by Muratov et al. (2015) to the FIRE simulations. Mitchell et al. (2020a) have compared the mass loading factor from EAGLE simulation with other hydrodynamical simulations and different versions of SAM. They found large discrepancy between simulations and SAMs, mainly due to different definitions of outflow and physics implemented in different simulations. Overall, the mass loading factor from EAGLE (Mitchell et al. 2020a) is larger than the Horizon-AGN simulation (Beckmann et al. 2017), but lower than FIRE simulation (Hopkins et al. 2014), TNG50 simulation (Nelson et al. 2019b), and the NIHAO simulation (Tollet et al. 2019). Among these hydrodynamical simulations, the mass loading factor in NIHAO is the strongest. We refer the readers to the paper by Mitchell et al. (2020a) for more details. In Fig. 11, we show two additional cases where we increase/decrease the mass loading factor, η, by a factor of 2. It is interesting to find that increasing the stellar feedback efficiency, or the mass loading factor, results in a higher recycle fraction from the outflow (dashed line in the left-hand panel) and a shorter recycle time (dashed line with circles in the right-hand panel). Reducing the stellar feedback brings an opposite effect (dotted lines with filled circles). Using our empirical model introduced in Section 2.4, it is straightforward to show that the recycle fraction is related to the outflow as
(15)
In our model, the cold accretion rate (equation 11), cold gas, and stellar mass are all independent of the parameter η; thus, it is understandable from equation (15) that a larger η will result in a larger fREC and a lower recycle time.
The same as Fig. 8, but here we increase/decrease the fiducial mass loading factor from Muratov et al. (2015) by a factor of 2. It is found that a stronger mass loading factor, or higher stellar feedback efficiency, agrees better with the simulation results on both gas recycle fraction (dashed line with circle connected in the left-hand panel) and recycle time (the same line style in the right-hand panel). Our results suggest that the mass loading factor from Muratov et al. (2015) could be a lower limit, and a boost by a factor of 2 is more favoured by the simulation data.
Figure 11.

The same as Fig. 8, but here we increase/decrease the fiducial mass loading factor from Muratov et al. (2015) by a factor of 2. It is found that a stronger mass loading factor, or higher stellar feedback efficiency, agrees better with the simulation results on both gas recycle fraction (dashed line with circle connected in the left-hand panel) and recycle time (the same line style in the right-hand panel). Our results suggest that the mass loading factor from Muratov et al. (2015) could be a lower limit, and a boost by a factor of 2 is more favoured by the simulation data.

In general, most hydrodynamical simulations have implemented different stellar feedback physics and tuned their parameter to best fit different set of observational data, it is naturally to see there is a large diversity among their predictions on gas recycle. As we mentioned before, although the mass loading factor is largely diverse in different simulations (Mitchell et al. 2020a), the predicted gas recycle fraction and the gas recycle time from the outflow are not significantly different (points in Fig. 11). One possible reason is that gas recycle fraction does not depend strongly on the details of the stellar feedback, but the dynamical properties of the halo where a galaxy stays. Our results in Fig. 11 indicate that the mass loading factor we used (Muratov et al. 2015) could be a lower limit. In fact, our results favour a stronger mass loading factor and a recycle time, which are similar to the constraints from the L-Galaxies model (Henriques et al. 2013).

5 SUMMARY

In this paper, we use the evolution of the galaxy stellar-to-halo mass relation to obtain SFH of a few model galaxies with halo mass ranging from 1011 to |$10^{13}\, \mathrm{M}_{\odot }$| at z = 0. We combine their SFH with a simple empirical model for galaxy evolution to constrain the gas cycle. Our results can be summarized as the following:

  • Galaxies are not always evolving along the MS throughout their life times. Low-mass galaxies (⁠|$M_{\rm h} \lt 10^{12}\, \mathrm{M}_{\odot }$| at z = 0) had star formation rate lower than the MS at z > 1, and they gradually crossed the MS at z < 0.5 and stay above the MS at z = 0. High-mass galaxies (⁠|$M_{\rm h} \gt 10^{12}\, \mathrm{M}_{\odot }$| at z = 0) were star-forming galaxies at z > 2 with star formation rate close to the MS relations. With decreasing redshift, massive galaxies began to depart from the MS and gradually quenched. More massive galaxies quenched earlier. The evolution of galaxy SFH generally agrees with the ’downsizing’ scenario that massive galaxy forms first and is quenched early.

  • By incorporating a star formation efficiency with dependence on both stellar mass and redshift as |$SFE \sim M_{\ast }^{0.5}(1+z)^{0.3}$|⁠, we find that the cold gas fraction in all model galaxies at z = 2 is lower than the universal baryon fraction by a factor of 2, and it is consistent with the observations. The predicted cold gas fraction from TNG50-1 simulation does not evolve with redshift and is lower than the data at z = 2.

  • It is found that, in addition to the cold gas accretion due to halo growth, gas recycle from outflow is required to sustain the star formation in our model galaxy. The recycle fraction is between 20 and 60 per cent for galaxy with current halo mass between 1011 and |$10^{12}\, \mathrm{M}_{\odot }$|⁠. The peak recycle fraction is at halo mass around |$6\times 10^{11}\, \mathrm{M}_{\odot }$| and low-mass galaxy has lower recycle fraction. This is quantitatively consistent with simulation results. The gas recycle time for massive galaxy is around 0.5 Gyr, in consistent with simulation results, but it is about 4 Gyr for low-mass galaxy, higher than the simulation results.

  • We predicted the mass–metallicity relation using our model. In the model, the metal in ISM is enriched by newly formed stars. The metal could be replenished by newly accreted gas or diminished by gas outflow. It is found that if all metal in the outflow is lost and infall gas has primordial gas metallicity, Zgas ∼ 10−4, the predicted metallicity is lower than the data at z = 0 even we have tuned the initial mass–metallicity relation to match the data at z = 2. However, if all metal in the outflow is recycled to the ISM, the predicted MZR is higher than the data at z = 0 and can be excluded. An equal fraction of recycle for both gas outflow and metal produces a MZR that agrees with the data slightly better.

We also find that both the gas recycle fraction and recycle time are affected by the processes of cold gas accretion and stellar feedback. From the discussion of Section 3.3.2 and Fig. 8, we can conclude that if the gas fraction in newly accreted halo is lower than the universal baryon fraction (e.g. fb = 0.075), as indicated by previous simulation results (e.g. Mitchell et al. 2018; Tollet et al. 2019), the required gas recycle fraction from the outflow is higher and the recycle time is shorter, which are more consistent with the simulation results shown in Fig. 8. From the discussion on the effect of stellar feedback in Section 4 and Fig. 11, we find that a stronger feedback efficiency or a higher mass loading factor than that from the FIRE simulation (Muratov et al. 2015) will lead to the same result as above. So in our model, the cold accretion rate and mass loading factor can be constrained by the cooling of the recycled gas.

It has been long recognized that the circum-galactic medium (CGM) of a galaxy is composed of gas from both incoming accretion and ejected gas from star formation, and CGM plays an important role in governing galaxy formation and evolution (e.g. Tumlinson, Peeples & Werk 2017). Great efforts have been paid to observe the CGM and characterize its fundamental properties (e.g. total mass, cooling rate). For example, Zhang et al. (2021) recently measured the mass cooling rates for CGM gas as a function of galactic stellar mass using the observations of H α emission line. Simply considering the galactic gas outflow forms the major part of the CGM of a galaxy (Tumlinson et al. 2017), the mass cooling rate of CGM gas can be greatly simplified as the mass recycle rate of outflow gas (⁠|$\dot{M}_{\rm gas,rec}$|⁠). Then using our equation (15), with the recycle fraction (fREC) from simulations, the mass loading factor η will be uniquely determined. Therefore, the cold accretion rate for galaxies (⁠|$\dot{M}_{\rm gas,acc}$|⁠) will also be uniquely determined. In addition, with the mass cooling rate of CGM gas and the total mass of CGM gas, the cooling time can be roughly estimated. Of course, the above intuitive derivation is based on a few simple assumptions, and more accurate and complex models are needed to realize and test it. Never the less, by this simple derivation, we can see the important role that future observation of CGM observations can greatly constrain the baryon cycle process in galaxy formation.

On the other hand, it would also be interesting to incorporate the constraints on gas recycle fraction and recycle time obtained in this study into a semi-analytical model to see if it can reproduce the global properties of galaxy, such as the evolution of galaxy SMF and star formation MS. We will investigate them in a future work.

ACKNOWLEDGEMENTS

This work is partly supported by the National Key Research and Development Program (No. 2022FYA1602903), NSFC (No. 11825303, 11861131006), the Fundamental Research Fund for Chinese Central Universities (No.226-2022-00216), the science research grants from the China Manned Space project with NO.CMS-CSST-2021-A03, CMS-CSST-2021-B01, and the cosmology simulation database (CSD) in the National Basic Science Data Center (NBSDC-DB-10). We thank Miao Li, Yu Luo, Qingxin Wen, Yisheng Qiu, and Liping Wei for helpful discussions.

DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

Footnotes

1

The cold gas data here refer to the sum of molecular hydrogen and atomic hydrogen, from the calculations in Diemer et al. (2018, 2019).

REFERENCES

Andrews
B. H.
,
Martini
P.
,
2013
,
ApJ
,
765
,
140

Anglés-Alcázar
D.
,
Faucher-Giguère
C.-A.
,
Kereš
D.
,
Hopkins
P. F.
,
Quataert
E.
,
Murray
N.
,
2017
,
MNRAS
,
470
,
4698

Asplund
M.
,
Grevesse
N.
,
Sauval
A. J.
,
Scott
P.
,
2009
,
ARA&A
,
47
,
481

Beckmann
R. S.
et al. ,
2017
,
MNRAS
,
472
,
949

Behroozi
P. S.
,
Conroy
C.
,
Wechsler
R. H.
,
2010
,
ApJ
,
717
,
379

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

Berlind
A. A.
,
Weinberg
D. H.
,
2002
,
ApJ
,
575
,
587

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

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

Brooks
A. M.
,
Governato
F.
,
Quinn
T.
,
Brook
C. B.
,
Wadsley
J.
,
2009
,
ApJ
,
694
,
396

Chabrier
G.
,
2003
,
PASP
,
115
,
763

Christensen
C. R.
,
Davé
R.
,
Governato
F.
,
Pontzen
A.
,
Brooks
A.
,
Munshi
F.
,
Quinn
T.
,
Wadsley
J.
,
2016
,
ApJ
,
824
,
57

Ciesla
L.
,
Elbaz
D.
,
Fensch
J.
,
2017
,
A&A
,
608
,
A41

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

Colless
M.
et al. ,
2001
,
MNRAS
,
328
,
1039

Conroy
C.
,
Wechsler
R. H.
,
2009
,
ApJ
,
696
,
620

Contini
E.
,
Kang
X.
,
Romeo
A. D.
,
Xia
Q.
,
2017
,
ApJ
,
837
,
27

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

Cullen
F.
et al. ,
2019
,
MNRAS
,
487
,
2038

Daddi
E.
et al. ,
2022
,
A&A
,
661
,
A7

Darvish
B.
,
Scoville
N. Z.
,
Martin
C.
,
Mobasher
B.
,
Diaz-Santos
T.
,
Shen
L.
,
2018
,
ApJ
,
860
,
111

De Lucia
G.
,
Tornatore
L.
,
Frenk
C. S.
,
Helmi
A.
,
Navarro
J. F.
,
White
S. D. M.
,
2014
,
MNRAS
,
445
,
970

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

Diemer
B.
et al. ,
2018
,
ApJS
,
238
,
33

Diemer
B.
et al. ,
2019
,
MNRAS
,
487
,
1529

Donnari
M.
et al. ,
2019
,
MNRAS
,
485
,
4817

Fakhouri
O.
,
Ma
C.-P.
,
Boylan-Kolchin
M.
,
2010
,
MNRAS
,
406
,
2267

Faucher-Giguère
C.-A.
,
Kereš
D.
,
Ma
C.-P.
,
2011
,
MNRAS
,
417
,
2982

Fu
J.
et al. ,
2013
,
MNRAS
,
434
,
1531

Furlong
M.
et al. ,
2015
,
MNRAS
,
450
,
4486

Genel
S.
et al. ,
2008
,
ApJ
,
688
,
789

Genzel
R.
et al. ,
2015
,
ApJ
,
800
,
20

Giavalisco
M.
et al. ,
2004
,
ApJ
,
600
,
L93

Girelli
G.
,
Pozzetti
L.
,
Bolzonella
M.
,
Giocoli
C.
,
Marulli
F.
,
Baldi
M.
,
2020
,
A&A
,
634
,
A135

Gnedin
N. Y.
,
2000
,
ApJ
,
542
,
535

Grand
R. J. J.
et al. ,
2019
,
MNRAS
,
490
,
4786

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

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

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

Henriques
B. M. B.
,
White
S. D. M.
,
Thomas
P. A.
,
Angulo
R.
,
Guo
Q.
,
Lemson
G.
,
Springel
V.
,
Overzier
R.
,
2015
,
MNRAS
,
451
,
2663

Henriques
B. M. B.
,
Yates
R. M.
,
Fu
J.
,
Guo
Q.
,
Kauffmann
G.
,
Srisawat
C.
,
Thomas
P. A.
,
White
S. D. M.
,
2020
,
MNRAS
,
491
,
5795

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

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

Kang
X.
,
Jing
Y. P.
,
Mo
H. J.
,
Börner
G.
,
2005
,
ApJ
,
631
,
21

Kang
X.
,
Lin
W. P.
,
Skibba
R.
,
Chen
D. N.
,
2010
,
ApJ
,
713
,
1301

Kang
X.
,
Li
M.
,
Lin
W. P.
,
Elahi
P. J.
,
2012
,
MNRAS
,
422
,
804

Karim
A.
et al. ,
2011
,
ApJ
,
730
,
61
(K11)

Katsianis
A.
et al. ,
2020
,
MNRAS
,
492
,
5592

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

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

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

Lagos
C. d. P.
,
Bayet
E.
,
Baugh
C. M.
,
Lacey
C. G.
,
Bell
T. A.
,
Fanidakis
N.
,
Geach
J. E.
,
2012
,
MNRAS
,
426
,
2142

Leitner
S. N.
,
2012
,
ApJ
,
745
,
149

Leja
J.
,
van Dokkum
P. G.
,
Franx
M.
,
Whitaker
K. E.
,
2015
,
ApJ
,
798
,
115

Leja
J.
et al. ,
2022
,
ApJ
,
936
,
165

McBride
J.
,
Fakhouri
O.
,
Ma
C.-P.
,
2009
,
MNRAS
,
398
,
1858

Maiolino
R.
et al. ,
2008
,
A&A
,
488
,
463

Mandelbaum
R.
,
Seljak
U.
,
Kauffmann
G.
,
Hirata
C. M.
,
Brinkmann
J.
,
2006
,
MNRAS
,
368
,
715

Mitchell
P. D.
,
Lacey
C. G.
,
Cole
S.
,
Baugh
C. M.
,
2014
,
MNRAS
,
444
,
2637

Mitchell
P. D.
,
Blaizot
J.
,
Devriendt
J.
,
Kimm
T.
,
Michel-Dansac
L.
,
Rosdahl
J.
,
Slyz
A.
,
2018
,
MNRAS
,
474
,
4279

Mitchell
P. D.
,
Schaye
J.
,
Bower
R. G.
,
Crain
R. A.
,
2020a
,
MNRAS
,
494
,
3971

Mitchell
P. D.
,
Schaye
J.
,
Bower
R. G.
,
2020b
,
MNRAS
,
497
,
4495

Mo
H.
,
van den Bosch
F. C.
,
White
S.
,
2010
,
Galaxy Formation and Evolution
.
Cambridge Univ. Press
,
Cambridge

Moster
B. P.
,
Somerville
R. S.
,
Maulbetsch
C.
,
van den Bosch
F. C.
,
Macciò
A. V.
,
Naab
T.
,
Oser
L.
,
2010
,
ApJ
,
710
,
903

Moster
B. P.
,
Naab
T.
,
White
S. D. M.
,
2013
,
MNRAS
,
428
,
3121

Muratov
A. L.
,
Kereš
D.
,
Faucher-Giguère
C.-A.
,
Hopkins
P. F.
,
Quataert
E.
,
Murray
N.
,
2015
,
MNRAS
,
454
,
2691

Nelson
D.
et al. ,
2019a
,
Comput. Astrophys. Cosmol.
,
6
,
2

Nelson
D.
et al. ,
2019b
,
MNRAS
,
490
,
3234

Oliver
S.
et al. ,
2010
,
MNRAS
,
405
,
2279
(O10)

Oppenheimer
B. D.
,
Davé
R.
,
Kereš
D.
,
Fardal
M.
,
Katz
N.
,
Kollmeier
J. A.
,
Weinberg
D. H.
,
2010
,
MNRAS
,
406
,
2325

Parkinson
H.
,
Cole
S.
,
Helly
J.
,
2008
,
MNRAS
,
383
,
557

Peacock
J. A.
,
Smith
R. E.
,
2000
,
MNRAS
,
318
,
1144

Pillepich
A.
et al. ,
2018
,
MNRAS
,
473
,
4077

Pillepich
A.
et al. ,
2019
,
MNRAS
,
490
,
3196

Rodríguez-Puebla
A.
,
Primack
J. R.
,
Behroozi
P.
,
Faber
S. M.
,
2016
,
MNRAS
,
455
,
2592

Saintonge
A.
et al. ,
2011
,
MNRAS
,
415
,
32

Saintonge
A.
et al. ,
2017
,
ApJS
,
233
,
22

Sanders
R. L.
et al. ,
2018
,
ApJ
,
858
,
99

Sanders
R. L.
et al. ,
2020
,
MNRAS
,
491
,
1427

Sanders
R. L.
et al. ,
2021
,
ApJ
,
914
,
19

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

Schreiber
C.
et al. ,
2015
,
A&A
,
575
,
A74

Schreiber
C.
,
Elbaz
D.
,
Pannella
M.
,
Ciesla
L.
,
Wang
T.
,
Koekemoer
A.
,
Rafelski
M.
,
Daddi
E.
,
2016
,
A&A
,
589
,
A35

Scoville
N.
et al. ,
2007
,
ApJS
,
172
,
38

Scoville
N.
et al. ,
2017
,
ApJ
,
837
,
150

Somerville
R. S.
,
Davé
R.
,
2015
,
ARA&A
,
53
,
51

Sparre
M.
et al. ,
2015
,
MNRAS
,
447
,
3548

Speagle
J. S.
,
Steinhardt
C. L.
,
Capak
P. L.
,
Silverman
J. D.
,
2014
,
ApJS
,
214
,
15

Springel
V.
,
White
S. D. M.
,
Tormen
G.
,
Kauffmann
G.
,
2001
,
MNRAS
,
328
,
726

Tacconi
L. J.
et al. ,
2018
,
ApJ
,
853
,
179

Tollet
É.
,
Cattaneo
A.
,
Macciò
A. V.
,
Dutton
A. A.
,
Kang
X.
,
2019
,
MNRAS
,
485
,
2511

Tomczak
A. R.
et al. ,
2016
,
ApJ
,
817
,
118

Tremonti
C. A.
et al. ,
2004
,
ApJ
,
613
,
898

Tumlinson
J.
,
Peeples
M. S.
,
Werk
J. K.
,
2017
,
ARA&A
,
55
,
389

Vale
A.
,
Ostriker
J. P.
,
2004
,
MNRAS
,
353
,
189

van de Voort
F.
,
Schaye
J.
,
Booth
C. M.
,
Haas
M. R.
,
Dalla Vecchia
C.
,
2011
,
MNRAS
,
414
,
2458

van den Bosch
F. C.
,
Yang
X.
,
Mo
H. J.
,
2003
,
MNRAS
,
340
,
771

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

Vogelsberger
M.
et al. ,
2014a
,
MNRAS
,
444
,
1518

Vogelsberger
M.
et al. ,
2014b
,
Nature
,
509
,
177

Vogt
N. P.
et al. ,
2005
,
ApJS
,
159
,
41

Wang
L.
,
Dutton
A. A.
,
Stinson
G. S.
,
Macciò
A. V.
,
Penzo
C.
,
Kang
X.
,
Keller
B. W.
,
Wadsley
J.
,
2015
,
MNRAS
,
454
,
83

Wechsler
R. H.
,
Tinker
J. L.
,
2018
,
ARA&A
,
56
,
435

Wechsler
R. H.
,
Bullock
J. S.
,
Primack
J. R.
,
Kravtsov
A. V.
,
Dekel
A.
,
2002
,
ApJ
,
568
,
52

Weiner
B. J.
et al. ,
2005
,
ApJ
,
620
,
595

Whitaker
K. E.
,
van Dokkum
P. G.
,
Brammer
G.
,
Franx
M.
,
2012
,
ApJ
,
754
,
L29

Whitaker
K. E.
et al. ,
2014
,
ApJ
,
795
,
104

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

Yang
X.
,
Mo
H. J.
,
van den Bosch
F. C.
,
2003
,
MNRAS
,
339
,
1057

Yang
X.
,
Mo
H. J.
,
van den Bosch
F. C.
,
Zhang
Y.
,
Han
J.
,
2012
,
ApJ
,
752
,
41

York
D. G.
et al. ,
2000
,
AJ
,
120
,
1579

Zhang
H.
et al. ,
2021
,
ApJ
,
916
,
101

Zheng
Z.
,
Coil
A. L.
,
Zehavi
I.
,
2007
,
ApJ
,
667
,
760

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)