Anatomy of galactic star formation history: Roles of different modes of gas accretion, feedback, and recycling

We take the parametrized form for the mass growth of the dark matter halo given in equations (H1)-(H6) in Behroozi, Wechsler & Conroy (2013), which is tested against Bolshoi, Consuelo, and Multidark cosmological simulations.

Gas accretion on to the galactic discs and subsequent star formation from the interstellar gas is the most important driver of the galaxy formation and evolution. In recent development of the hydrodynamical cosmological simulations, we have seen a significant revision of the traditional view about the physics of cosmic gas accretion.

The widely-used scenario for gas accretion on to forming galaxies, namely, the shock heating paradigm, assumes that the gas entering a dark matter halo is heated by shock waves to roughly virial temperature of the halo and subsequent radiative cooling induces the accretion of this gas to the inner parts where the galaxy grows (e.g. Rees & Ostriker, 1977). In this picture, the gas accretes with the free-fall time if it cools rapidly while the accretion proceeds with the radiative cooling time when the cooling is not efficient. The cooling argument underlying this scenario succeeded in explaining the segregation of galaxies and galaxy clusters on the density-temperature plane (Blumenthal et al., 1984) and has constituted a backbone for many theoretical studies.

However, realistic simulations for hydrodynamical evolution of the primordial gas in the $\Lambda$CDM Universe (e.g. Fardal et al., 2001; Kereš et al., 2005; Dekel & Birnboim, 2006; Ocvirk, Pichon & Teyssier, 2008; van de Voort & Schaye, 2012; Nelson et al., 2013) found that the assumption of the stable shock underlying this picture is not always correct. The shock wave is unstable and no shock heating occurs in low mass halos. The unheated cold gas directly accretes to the inner parts in free-fall fashion in this case (e.g. Fardal et al., 2001; Birnboim & Dekel, 2003; Kereš et al., 2005; Dekel & Birnboim, 2006). In addition to this cold flow in low mass halos, high mass halos at high redshifts develop filaments of cold gas inflowing through the surrounding hot gas created by shocks (e.g. Dekel & Birnboim, 2006; Ocvirk, Pichon & Teyssier, 2008; van de Voort & Schaye, 2012; Nelson et al., 2013). The new picture, sometimes called the cold accretion scenario, has a large impact on our understanding of galaxy evolution, For example, it helps interpret the high abundance of actively star-forming galaxies at high redshifts (e.g. Cattaneo et al., 2006; Dekel et al., 2009). Lyman limit systems may be observational manifestations of the predicted cold streams (Fumagalli et al., 2011). The Ly$\alpha$ emitting filaments at high redshift observed by Martin et al. (2016) may provide a more direct example. Extending this picture to larger spatial scales, Overzier (2016) proposes that the transition from the cold accretion regime to the regime dominated by the hot intracluster gas in most massive halos marks the division between protoclusters and clusters of galaxies.

Noguchi (2018) demonstrated that the cold accretion scenario, combined with a chemical evolution model, reproduces the bimodality in the stellar abundance distributions of the Milky Way disc stars (for new observations, see Queiroz et al. (2020)). The transition from the early cold-mode to the late hot-mode gas accretion sets the stage for two distinctive star formation episodes (see also Birnboim, Dekel & Neistein (2007)). The short period active star formation in early times produces metal-poor (i.e., deficient in Fe) stars with relatively enriched $\alpha$-elements by the dominant action of type II supernovae. The late-time star formation proceeds slowly for a longer period allowing much iron enrichment by type Ia supernovae, leading to a metal-rich $\alpha$-element-poor (low $[\alpha/{\rm Fe}]$) population. Several cosmological simulations also show signs of such two-stage star formation in Milky Way mass galaxies (e.g. Hopkins et al., 2014; Grand et al., 2018), although details are often masked by other baryonic processes included. For galaxy morphology, the cold flow scenario was found to reproduce the relative dominance of thin discs, thick discs, and bulges in the present-day disc galaxies as a function of galaxy mass by associating the hot mode with thin disc formation and the cold mode with bulge and thick disc formation, respectively (Noguchi, 2020). The observed galaxy-mass dependence of the mean age and the age spread of the bulge (Breda & Papaderos, 2018; Breda et al., 2020) is also naturally explained by this model. We adopt this cold flow picture in our fiducial scheme described later and examine its essential role by comparing it with the accretion scheme based on the shock-heating picture.

The accreted halo gas turns into interstellar gas in the galactic disc and serves as raw material for star formation. To calculate the star formation rate (SFR), we adopt the formulation by Blitz & Rosolowsky (2006) in which the SFR is essentially determined by the molecular gas content in the interstellar medium (ISM), with the molecular fraction in turn controlled by the disc midplane pressure. This treatment enables us to trace the evolution of the atomic and molecular contents separately. Those results are compared with the observation. Mass-loss associated with the evolution of massive stars is evaluated by assuming Chabrier IMF with lower and upper initial mass limits set to $1{\rm M}_{\odot}$ and $100{\rm M}_{\odot}$, respectively. The resulting mass fraction, $R=0.44$ is returned to the ISM instantaneously (see Appendix A3).

It is well established that the evolution of low-mass galaxies is strongly affected by the ejection of the ISM by supernovae on both observational and theoretical grounds (e.g. Dekel & Silk, 1986; Mac Low & Ferrara, 1999; Brooks et al., 2007; McQuinn, van Zee & Skillman, 2019). It not only helps explain observations such as the mass-metallicity relation but also solves the over-production problem in the low galaxy mass domain inherent in the $\Lambda$CDM paradigm (see Bullock & Boylan-Kolchin (2017) for a review and Weinmann et al. (2012) for detailed discussion). We take the energy-driven feedback as the standard choice although the model with the momentum-driven feedback is also discussed.

Recent studies (Torrey et al., 2014; Somerville, Popping & Trager, 2015, e.g.) suggest that supernova feedback alone cannot produce the masses of low-mass galaxies consistent with their metallicity. In order to realize the observed metallicity, those galaxies should be more massive than observed. Another type of feedback, called preventive feedback, was found to resolve this tension. Preventive feedback refers to the suppression of gas accretion on to halos or galactic discs by some mechanism(s). It can thus reduce the stellar mass of these galaxies without affecting metal enrichment driven by ejective feedback. Lu et al. (2017) point out the importance of preventive feedback in reproducing the galaxy stellar mass function and the mass-metallicity relation of the Milky Way satellites simultaneously. Several cosmological simulations (Angléz-Alcázar et al., 2017; Christensen et al., 2016; Oppenheimer et al., 2010) show that this type of preventive feedback indeed occurs in galaxy formation although its detailed nature is not completely clarified. Candidate mechanisms include radio-mode feedback by active galactic nuclei (AGNs) (Kereš et al., 2009), heating from outflowing winds (Oppenheimer et al., 2010), and heating by ionizing radiation field (Efstathiou, 1992; Gnedin, 2000; Hoeft et al., 2006; Okamoto, Gao & Theuns, 2008; Christensen et al., 2016). Preventive feedback can take two different forms. One is the reduction of the gas entering dark matter halos and another is the suppression of the gas supply from parent halos to galactic discs. Although this should be discriminated in the strict sense (e.g. Christensen et al., 2016), we include this process by reducing the accretion of the halo gas to the disc by a specified factor. We confirmed that the reduction of the gas entering the halo by the same factor gives indistinguishable results.

Outflow of the interstellar gas is ubiquitously observed in star-forming galaxies with a large mass range and taken as a signature of ejective feedback. Because the outflow velocity attained by supernova feedback is typically less than several hundred kilometre per second (e.g Chisholm et al., 2015), the outflowing gas may not leave the host galaxy entirely but a part of it may eventually fall back to the galaxy. Many numerical simulations demonstrate that such a fall-back actually occurs in actively star-forming galaxies (e.g. Christensen et al., 2016; Grand et al., 2019). This type of gas recycling has also been incorporated in the semi-analytic models (SAMs) (e.g. Croton et al., 2006; Henriques et al., 2013; White, Somerville & Ferguson, 2015) and its effect on galaxy evolution has been intensively studied. These studies show that recycling affects not only global properties of galaxies but also their internal structures (e.g. Brook et al., 2012; Marasco, Fraternali & Binney, 2012; Grand et al., 2019). On the other hand, observational confirmation of gas recycling turned out to be inconclusive (see Fraternali & Pezzulli (2018) for a review). Although the infall of cold gas like galactic fountains is observed, several possibilities can be considered for its origin, including the cooling flow of the hot halo gas or the inflowing streams of unheated gas (cold-mode accretion). It is theoretically anticipated that the circulation of the ejected gas over much larger scales than the observed local fountains also takes place (e.g. Angléz-Alcázar et al., 2017; Mitchell, Schaye & Bower, 2020). Recent studies (e.g. Ford et al., 2014; Hasan et al., 2022) suggest that the recycling gas may be observed as quasar absorption lines (especially ${\rm C}_{\rm IV}$) but more work is required to firmly establish their connection. Considering these situations, we include this process simply by returning the ejected gas to the disc component with a specified timescale.

The fiducial scheme assumes the existence of cold-mode accretion as revealed in recent cosmological simulations. Cosmological simulations and analytical consideration predict three different regimes for the properties of the gas distributed in the dark matter halos depending on the virial mass and redshift (e.g. Dekel & Birnboim, 2006; Ocvirk, Pichon & Teyssier, 2008). They are demarcated by two characteristic mass scales as shown in Figure 1 (left panel). One is the shock mass, $M_{\rm shock}$, above which the halo gas develops a stable shock that heats the gas nearly to the virial temperature and the other is the stream mass, $M_{\rm stream}$, below which part of the halo gas remains cold and is confined into narrow filaments that thread the smoothly distributed shock-heated gas. The latter mass scale is valid only for high redshifts.

The halo gas in Domain F is unheated and expected to accrete in free-fall to the inner region (the disc plane). In domain G, the gas heated to the virial temperature attains near hydrostatic equilibrium in the halo gravitational field and the radiative cooling induces cooling flow to the centre with the cooling timescale. The gas behaviour is not so clear in Domain H, where cold gas streams coexist with the surrounding shock-heated hot gas. They may behave independently and accrete with their own timescales or they may interact with each other leading to modification of accretion timescales. Because no detailed information is available, we assume that the cold and hot gases in Domain H accrete with the free-fall time and the radiative cooling time, respectively.

It is difficult to pin down the transition borders, namely $M_{\rm shock}$ and $M_{\rm stream}$, which depend on various physical quantities of accreting gas with the gas metallicity being the most influential parameter. Cosmological simulations using different gas dynamical codes (e.g., SPH and AMR) show different results, with the Eulerian mesh codes (e.g., AMR) tending to suppress the cold fraction compared with the Lagrangian codes (e.g., SPH) (e.g. Nelson et al., 2013). We use the average of the results from several studies summarized in Figure 6 of Ocvirk, Pichon & Teyssier (2008), keeping in mind that the adopted prescription may suffer from ambiguity and affect quantitative detail of the result. Trial with several variations for the transition borders, however, revealed that the qualitative nature of model evolution is conserved within plausible parameter ranges.

Existence of Domain H is a distinguished characteristic of the cold mode accretion. In the previous work (Noguchi, 2020), it was related to the formation of galactic bulges. The transition from the cold to hot modes corresponds to the quenching of star formation because the gas accretion slows down abruptly at this point. The reconstruction of galaxy assembly and SFH for the redshift range $z=1-10$ by Behroozi et al. (2019) indicates increase of the quenching mass with redshift, though less steeply than the theoretical prediction of Dekel & Birnboim (2006). To elucidate the effect of Domain H on the galactic SFH, we run an artificial model family in which this domain is removed (middle panel in Figure 1).

The shock heating picture has been routinely assumed in many theoretical studies of galaxy formation. It is interesting to see the difference between the star formation history derived under the fiducial scheme and this picture. There are two domains in this scheme separated by the equality of the radiative cooling time and the free-fall time (Figure 1, right panel). When the gas cools with timescale shorter than the free-fall time, it is made to accrete with the free-fall time. Otherwise, the gas is assumed to accrete with the cooling timescale.

For each accretion scheme, we run a series of galaxy models with the present halo virial masses varied in the range $10^{10.36}{\rm M}_{\odot}\leq M_{\rm vir,0}\leq 10^{14}{\rm M}_{\odot}$. This upper mass limit roughly corresponds to the halo mass inferred for the most massive (disc) galaxies with the stellar mass $\sim 10^{11.5}{\rm M}_{\odot}$ on the basis of the abundance matching analysis (e.g. Moster, Naab & White, 2013; Rodríguez-Puebla et al., 2015; Behroozi et al., 2019). The lower mass limit corresponds to the stellar mass of $\sim 10^{7}{\rm M}_{\odot}$, below which the observational material becomes sparse. Grey lines in Figure 1 show the evolution of each halo, illustrating how it goes through different regimes of gas accretion as it grows in mass. Because no galaxy merger is considered in this study, the present models are regarded as mainly describing the evolution of galaxies whose growth rates are inflow dominated.

It is remarked here that the schemes illustrated in Figure 1 may need some modification. As seen later, preventive feedback which acts to reduce gas accretion to the galactic disc, is found to be important especially for low mass halos. The red line in each panel of Figure 1 indicates the critical halo mass assumed in this study for which this effect starts to be active. Although the specific mechanism responsible for this process is not clear, previous studies suggest that it affects a wider range of mass than the ionization by background radiative fields. For comparison, we also plot the ’filtering mass’ by Gnedin (2000) and the 90$\%$ suppression line by Okamoto, Gao & Theuns (2008) as a typical estimate for the radiative effect. Preventive feedback, either on to discs or on to halos, may thus significantly alter the picture solely based on the cooling argument (and assumption of a universal cosmic baryon fraction). We revisit this issue later.

Figure 2 shows the stellar-to-halo mass ratio as a function of the halo virial mass for different redshifts. The halo mass here denotes the mass at each redshift. The empirical analysis by Behroozi et al. (2019) using the abundance matching method shows no significant difference between the star forming galaxies and all galaxies including quiescent galaxies. At the present epoch ($z\sim 0$), all schemes give similar results approximately agreeing with the result by Behroozi et al. (2019) and other observational studies. Moving to higher redshifts ($z>2$), the shock-heating scheme progressively underproduces stars especially for more massive halos. This is due to earlier transition of massive halos to the hot-mode gas accretion and associated weakening of star formation. The predicted redshift dependence seems to be significantly stronger than the observationally inferred one even if we make allowance for the possible systematic deviation as suggested by the small ($\sim 0.3$ dex) overproduction in the low mass domain ($M_{\rm vir}<10^{12}{\rm M}_{\odot}$) common to all three schemes.

Figure 3 shows the star formation history (SFH) for different values of the present halo masses. It is seen that the fiducial scheme captures several key features of the observed SFHs. Namely, the SFR in halos with larger present masses attains the peak at earlier epochs and declines more drastically after the peak. The halos in the lowest mass bin (the blue line) show monotonically increasing SFR. Notable dicrepancy is seen at z>2 where the models underpredict SFR compared with the observation. This may indicate that the adopted feedback is too strong. However, the discrepancy for most massive halos appears also in absence of feedback as seen later in section 5. We discuss possible reasons for this discrepancy there.

The flat scheme is unable to reproduce this mass dependent behaviour. This indicates that the existence of the domain H defined by the existence of cold gas streams contributes much in producing the observed mass dependent SFH, especially the early drastic ’quenching’ in massive halos. The shock-heating scheme also fails due to the same reason. It is interesting that Behroozi et al. (2019) compare their results for the quenching fraction (their Figure 30) to the prediction of Dekel & Birnboim (2006) cold accretion model and find rough agreement, i.e., the domain for star forming galaxies extends to higher halo masses at higher redshifts. This is in line with the present result showing the importance of the cold-mode accretion of gas filaments penetrating through the surrounding shock-heated halo gas although some other mass-dependent process (e.g. mergers) might bring about a similar effect.

Differences in SFHs in the three schemes displayed in Figure 3 appear in high-mass domain with ${\rm log}M_{\rm vir}[{\rm M}_{\odot}]>12.5$, where accretion processes are different. Specifically, massive galaxies in the fiducial scheme experience peak SFR in earlier epochs and larger drop in SFR after the peak than those in the flat scheme while massive galaxies in the shock heating scheme show no peak. Figure 4 is meant to clarify the origin of these various SFHs going back to the difference in the prescribed accretion recipes. Top panels show the gas accretion rates (red solid lines) and SFRs (red dotted lines) for ${\rm log}M_{\rm vir}[{\rm M}_{\odot}]=14$ in the three accretion schemes. For aid in analysis, middle panels show masses of the halo gas and bottom panels plot the radiative cooling time (blue) and the free-fall time (green) for the halo gas. We here turn off supernova feedbacks, gas recycling and preventive feedback to explore more directly how the halo gas mass and the two timescales eventually control the SFR under accretion-driven situation.

Top panels indicate that SFRs essentially follow gas accretion rates with small time delay. Gas accretion rates in turn can be approximated by the amount of halo gas divided by the accretion timescales. Middle panels show that the mass of halo gas (solid) is almost the same in the three schemes. Therefore, the differences in SFHs should originate mostly in those in accretion timescales. Green and blue lines in upper panels are the ’accretion rates’ calculated by dividing the true halo gas mass respectively by the free-fall time and radiative cooling times. These are again similar in the three schemes reflecting the similarity in the halo gas mass.

The key for interpreting different SFHs lies in how the transition from the cold-mode accretion to the hot-mode accretion occurs. At the transition, the accretion timescale is switched from the free-fall time to the radiative cooling time. Figure 5 compares this transition for three accretion schemes. In the shock-heating scheme, the transition occurs exactly when the cooling time overshoots the free-fall time, but in the other two schemes it occurs when the initial cold-accretion phase ends, which is later than the equality of two timescales. Instability of shocks in the flat scheme delays the transition compared with that in the shock-heating scheme. Existence of cold streams in high-mass halos further retards the transition in the fiducial scheme. Upper panels in Figure 4 overplot the transition epoch by red circles on the accretion rate curves. It is seen that the actual accretion rates (red lines) follow closely the green lines before this transition epoch but they deviate and approach to blue lines after that. Comparison of three accretion schemes thus reveals that the difference in the accretion history comes from the difference in the transition epoch. Moving from the shock-heating to flat to fiducial schemes, the transition epoch moves progressively to more recent epochs and therefore the drop in accretion rate gets larger because of larger difference in timescales at transition. This explains why the peak SFR occurs latest and the decrement in SFR after the peak is the largest for the fiducial scheme as shown in Figure 3.

The SFHs such as plotted in Figure 3 will contain all information about the star formation activity in galaxies if accurately derived over the whole mass range. However, halo masses in the empirical approach are not obtained for each sample galaxy but derived statistically from abundance matching. On the other hand, it is not easy either to measure the halo masses for individual galaxies more directly, for example, from satellite kinematics (e.g. More et al., 2011) or weak gravitational lensing (e.g. Mandelbaum et al., 2006). Any derived SFHs as a function of halo mass could thus be subject to uncertainty even if SFRs are obtained precisely. We examine in the next section the dependence of star formation activity on the stellar mass, which is more directly assessed than the halo mass in general.

The correlations between star formation properties and galaxy stellar masses have been intensively examined in previous studies. Among various scaling relations, the star formation main sequence (MS) of galaxies, i.e., the SFR-$M_{\star}$ relation, provides useful information on how galaxies with different masses formed stars over the cosmic time (e.g. Whitaker et al., 2014; Johnston et al., 2015; Schreiber et al., 2015; Kurczynski et al., 2016; Tomczak et al., 2016; Schreiber et al., 2017; Popesso et al., 2019; Corcho-Caballero, Ascasibar & López-Sánchez, 2020; Khusanova et al., 2021) and serves as an independent constraint on theoretical modelling (e.g. Feldmann, 2020).

Figure 6 compares the SFR and sSFR predicted by the models with the observational data. It is seen that the fiducial scheme reproduces several important features at least qualitatively. Regarding the ${\rm SFR}-M_{\star}$ relation (top panels), it gives the SFR that increases almost linearly with the stellar mass for the low-mass range, with normalization increasing with redshift. The relation at higher masses bends in progressively larger degree for lower redshifts. This behaviour agrees with the observational results about the redshift evolution of the normalization of the star formation MS and the turn-over in the high mass range (e.g. Whitaker et al., 2014; Tomczak et al., 2016). Both the flat and shock-heating schemes give a reverse trend at high masses with smaller bending at lower redshifts. The middle panels plot the sSFR as a function of the stellar mass. The fiducial scheme shows nearly horizontal relation for low masses and increasingly steeper decline at high masses toward lower redshifts, agreeing with the observational analysis. The other two schemes fail to reproduce the observed trend. Although we take here the observation by Whitaker et al. (2014) as a typical example, other observational studies also reveal similar trends (e.g. Schreiber et al., 2015; Tomczak et al., 2016; Rodríguez-Puebla et al., 2020).

Difference in the turn-over at high masses in different accretion schemes is directly related to difference in the SFHs. In order to clarify this point, Figure 7 overplots the trajectory of evolving model galaxies on the ${\rm SFR}-M_{\star}$ relation. It is evident that large flattening at high mass end in the fiducial scheme is caused by strong drop in the SFR when the gas accretion is switched from the cold mode (solid curve) to the hot mode (dashed curve) in massive halos with $M_{\rm vir}>10^{13}{\rm M}_{\odot}$. During the hot-mode phase, $M_{\star}$ increases by only 0.2 dex while SFR decreases by almost 1 dex in the most massive halo. The decrement in SFR is larger for more massive halos that move down on more right sides of the ${\rm SFR}-M_{\star}$ diagram, leading to the fan shape pattern of redshift sequence in the fiducial scheme. High-mass end of the main-sequence is related to the hot mode accretion also in the flat and sock-heating schemes. In these cases, however, the decrease in SFR is milder and the stellar mass increases with almost constant SFR. Namely, the galaxy moves rightward nearly horizontally in the diagram. This means that the main sequence itself moves rightward. Combination of this movement and the weak turn-over at high stellar masses at each redshift makes a configuration of the redshift sequence converging to higher stellar masses.

There has been much debate about the existence of the turn-over at high masses (e.g. Santini et al., 2017; Popesso et al., 2019). Speagle et al. (2014) argue that the MS is nearly linear and its slope and normalization decrease toward recent epochs, dismissing the existence of turn-over. The slope at the high mass domain likely depends on the selection criteria for star-forming galaxies (Johnston et al., 2015). Galaxy morphology is also suggested to affect the high-mass behaviour of the MS. For example, including galaxies with dominant bulge components, usually characterized by smaller sSFR compared with the bulgeless galaxies of similar masses, was found to create strongly downward bending at high masses as reported by Cook et al. (2020). Teimoorinia, Bluck & Ellison (2016) also suggest close empirical relationship between the bulge mass (or the bulge-to-total mass ratio) and the star formation activity. There is also a possibility that the star formation properties and the main sequence characteristics depend on the galactic environments (e.g. Elbaz et al., 2007; Gnedin & Kravtsov, 2011; Zeimann et al., 2013; Ji et al., 2018; Noirot et al., 2018; Nantais et al., 2020; Old et al., 2020). We return to these issues later in Discussion.

The bottom panels of Figure 6 show the redshift evolution of the sSFR for different stellar masses. The fiducial scheme produces steeper decline of sSFR with time for more massive galaxies in qualitative agreement with data, whereas other two schemes show opposite trends in contrary to observations. This is a reflection of different turn-over behaviours seen in top panels. Absolute value of sSFR is underpredicted commonly in all accretion schemes for ${\rm logM}_{\star}[{\rm M}_{\odot}]$ < 10.4. For 10.4<${\rm logM}_{\star}[{\rm M}_{\odot}]$<11.0, the fiducial scheme gives closer match to the observational data than other schemes. Finally, for 11.0<${\rm logM}_{\star}[{\rm M}_{\odot}]$, the fiducial scheme exhibits smaller sSFR than other schemes.

The present result highlights the importance of gas accretion process in determining the shape and evolution of the star formation MS. Especially, the shock-heating scheme seems to be incapable of explaining the qualitative behaviour of turn-over in the star formation main sequence claimed by several observational works. This result agrees with that from several theoretical studies. The semi-analytic models adopting the shock-heating paradigm for the halo gas (e.g. Mitchell et al., 2014; Brennan et al., 2015; Henriques et al., 2015; Somerville, Popping & Trager, 2015; White, Somerville & Ferguson, 2015; Guo et al., 2016) tend to produce the sSFR-$M_{\star}$ relation that stays almost flat or even rises more steeply with the stellar mass toward more recent epoch, being completely opposite to the observation. The disc galaxy evolution models by Dutton, van den Bosch & Dekel (2010) and Lapi et al. (2020) assuming the shock heating also produce the SFR-$M_{\star}$ relation similar to these models. On the other hand, the model by Bouché et al. (2010) taking into account the cold-mode gas accretion produces the SFR-$M_{\star}$ relation that evolves only in the the normalization. Cosmological simulations (e.g. Hopkins et al., 2014; Vogelsberger et al., 2014; Furlong et al., 2015; Sparre et al., 2015; Schaye et al., 2015; Donnari et al., 2019; Scholz-Díaz, Sánchez Almeida & Dalla Vecchia, 2021) tend to produce turn-over at high masses and show agreement with the observation to some extent but its reason remains unclear. One possibility is that these simulations, treating gas processes more directly, automatically contain the cold-mode phase of gas accretion.

It should be noted that the main sequence turn-over may arise from other reason than specific feature of gas accretion. Generally, any feedback process operating preferentially on massive halos would produce a qualitatively similar effect. For example, the feedback from AGNs (when included) may contribute to the main sequence turn-over as argued in (Donnari et al., 2019). In this specific case, AGN feedback should closely follow star formation. Quasar-mode feedback may be thus promising, though the coupling of feedback energy and the halo gas should be clarified to confirm its feasibility.

We focused so far on the star formation activity represented by the SFR and the sSFR. We here examine the growth of stellar mass for different halo masses. Figure 8 illustrates the stellar mass build-up in terms of the absolute mass (top row) and the fractional mass (middle and bottom rows). It is clear that the fiducial scheme explains the so-called downsizing that more massive galaxies grow faster. Indeed, downsizing behaviour is already hinted in SFHs shown in Figure 3. This trend is more clearly shown by the fractional mass growth plotted in the middle and bottom rows. In the flat and shock-heating schemes, the galaxies with different masses grow almost in concert. The grey shading in the bottom panels indicates the observed range derived by Zhou et al. (2021) for massive red spirals. The fiducial scheme gives better reproduction than the other two schemes for the high-mass range. The reason is that the prolonged period of cold-accretion in this scheme allows more halo gas to accrete to the disc at early times (z>2) and raises SFR above those in other two schemes as illustrated in Figure 4. It nevertheless fails in matching the observation for $t>8$ Gyr. This arises from the upturn in SFR for massive halos at z<0.5 seen in Figure 3 because it brings about significant stellar mass growth in recent epochs and makes early quenching impossible. This discrepancy could be erased by including, for example, AGN radio-mode feedback which likely operates for low gas accretion rate in the hot-mode phase.

Star formation rates in galaxies are controlled by the amount of cold gas that serves as raw material for star formation and the timescale (depletion time) for gas consumption. The amount of gas in turn could be determined by several factors including the rate of gas accretion, ejection or removal rate of ISM, and recycling efficiency of ejected ISM. Confrontation of the model prediction for the evolution of gas contents with observational data therefore helps judge importance of each factor in galactic star formation. Gaseous content of the cold interstellar medium is classified into atomic hydrogens (HI) and molecular hydrogens (${\rm H}_{2}$). While the molecular gas is probed to high redshifts up to $z=2\sim 3$ thanks to the detection of CO lines, the observation of atomic component is currently restricted to the local Universe, $z<\sim 0.2$, due to the bandshift effect of 21 cm emission lines. Recent large surveys (e.g., BIMA SONG, HERACLES, COLD GASS) have revealed that the molecular masses relative to the stellar masses generally increase toward earlier cosmological epochs (e.g. Helfer et al., 2003; Leroy et al., 2008; Genzel et al., 2010; Saintonge et al., 2011; Tacconi et al., 2013). Census for the atomic gas in the local Universe (e.g., ALFALFA, THINGS, GASS) consistently indicates a decreasing HI mass fraction toward increasing stellar masses (e.g. Giovanelli et al., 2005; Walter et al., 2008; Catinella et al., 2010).

Figure 9 compares the model result (thick solid lines) with the currently available observational data (symbols). Redshift is colour-coded as indicated. The top panels show the mass ratio of the atomic hydrogen and the stellar content as a function of the stellar mass. At $z=0$, all the schemes give similar trends in broad agreement (a factor of $\sim$ 2) with the observation. As we move to higher redshifts, the difference between the three schemes becomes clearer, with the fiducial scheme showing the largest decrease. The thin lines indicate the empirical result obtained by Popping, Behroozi & Peeples (2015), which shows very little evolution with redshift. The simulation by Davé et al. (2013) also suggests little evolution in HI content. It is noted that the result by Popping, Behroozi & Peeples (2015) was obtained indirectly from the observed SFR by using the abundance matching technique combined with the locally established star formation law and the result by Davé et al. (2013) is also based on several assumptions in calculating HI fraction. If nevertheless these theoretical studies turn out to give better match to future observations, the shock-heating scheme might seem the most promising. There is a caveat here, however. Stronger decrease in HI content to higher redshifts for the fiducial scheme than in other schemes may have partly resulted from more active cold accretion that produced discs with richer gas content. Gas richness (i.e., higher surface densities) in discs means a relatively higher fraction of ${\rm H}_{2}$ in the disc and therefore poorer HI content. Another possible factor influencing gas contents is their spatial distribution. In actual galaxies, ${\rm H}_{2}$ and HI are spatially decoupled, with the former more concentrated to galactic centers (e.g. Young & Scoville, 1991). The model assumes that two components have the same scalelength corresponding to the halo spin parameter $\lambda$ of 0.06, which gives disc radii more relevant to the molecular component. Increasing $\lambda$ to 0.08 as a compromise between two components significantly lifts up HI content relatively to ${\rm H}_{2}$ due to resulting lower gas densities. It also diminishes difference in redshift dependence between different schemes. weakening apparent advantage of the shock-heating scheme. Argument relying on the spin parameter is admittedly rough but this result suggests a promissing avenue to resolve the difficulty shown here.

Middle panels of Figure 9 show the mass ratio of the molecular hydrogen. Compared with HI, all the schemes show much weaker mass-dependence with only the normalization increasing with redshift. Weak mass dependence agrees with the observed trend. The observational data and the empirical model by Popping, Behroozi & Peeples (2015) (thin lines) indicate the mass ratio at $z=2-3$ larger than the local value by more than $\sim 1.5$ dex. All schemes give good match to this result but tend to overpredict $M_{\rm H2}/M_{\star}$ up to $\sim 1$ dex at $z=0$. This likely originates in overabundance in ISM gas in model galaxies (especially massive ones) at recent epoch. Increasing gas ejection artificially by 1 dex for z<0.5 leads to lower molecular content and better agreement at z=0. This modification also helps erase turn-up of model SFRs at $z<0.5$ in massive halos seen in Figure 3 for the fiducial scheme. and improve agreement with the observation. Although we boosted here the ejection rate by supernova feedback artificially, more promising alternative mechanism to operate is radio-mode AGN feedback that likely acts preferentially for massive halos (Domain G of the fiducial scheme) at low redshifts.

Bottom panels of Figure 9 show the mass fraction of the total gas including atomic and molecular hydrogen. The three schemes all show an increasing fraction toward smaller stellar masses, with the value at fixed mass at $z=2$ larger than in the local Universe. These trends agree qualitatively with the observation but the latter suggests larger ratios and steeper mass dependence at $z=2$ whereas the ratios are smaller at $z=0$. The dashed lines indicate the result of the analysis performed by Popping, Behroozi & Peeples (2015). Their result shows larger redshift dependence than our models, though the slope of mass dependence is similar. The discrepancy of the fiducial scheme and observation can be partly relieved by using combination of the large spin parameter and boosting of gas ejection discussed above. Larger spins contribute to boost HI content at $z=2$ without affecting the molecular content, leading to larger gas-to-stellar ratio. Applied boosting reduces ${\rm H}_{2}$ content at $z=0$ reducing gas fraction. Two effect combine to produce larger redshift evolution in closer agreement with the observational data. Exploring evolution of two gas components consistently requires knowledge of density profiles of those comonents and reliable information about angular momentum distribution of infalling gas, which control disc sizes. Radio mode AGN feedback is believed to operate in low efficiency accretion to massive black holes provided by the hot-mode accretion (Fabian, 2012, e.g.). Including growth of black holes is necessary to treat this process properly. We defer such investigation to future study.

Recent observations of the molecular gas can reach $z>5$. Additional analysis was carried out of the evolution of the molecular gas over large redshift range. The result is shown in Figure 10 for different stellar mass bins. This figure gives the essentially same information as the middle panels of Figure 10 but now shows the redshift dependence to earlier cosmological epochs. The observations suggest the ratio of the molecular mass to the stellar mass increasing with redshift up to $z>5$ for all mass bins and larger ratios for lower mass galaxies at the fixed redshift. The three schemes considered here show a qualitatively similar behaviour to the observation although the discrepancy can reach $\sim 1$ dex at the lowest mass bin, where the three schemes give nearly identical evolution. Although the fiducial scheme cannot reproduce every detail of the observational result, it predicts strong evolution of the molecular content in agreement with the observation. These analyses also suggest that the improvement in the observational data for galactic gas contents (especially the redshift evolution of the atomic gas) is needed to further constrain the models examined here. The data expected to be obtained by the newest instruments such as Atacama Large Millimeter Array (ALMA) and Square Kilometre Array (SKA) will serve as a powerful constraint on theoretical models.

The importance of preventive feedback was suggested by the results discussed above. We introduced this process by reducing the gas accretion rate to the disc. Ejective feedback also helps suppress star formation by removing the interstellar gas. These two processes operate effectively in similar domains in the $M_{\rm vir}-z$ plane (Figure 16), suggesting a possibility that these processes are degenerate. Namely, a question arises whether ejective feedback with properly adjusted parameters can serve as a substitute for the required preventive feedback.

To explore this possibility, we first turned off preventive feedback and increased the efficiency of ejective feedback up to unity (i.e., all the energy liberated by supernovae is used in ejecting the ISM). For any choice in this range ($0.2<\varepsilon_{\rm FB}<1$), the replacement by ejective feedback exacerbates the poor fit to the observed stellar-to-halo mass ratio. Matching the ratio for low mass halos inevitably leads to unacceptable underproduction of the stellar content for high mass halos. This may be due to the steep $V_{\rm vir}^{-2}$ dependence of the energy-driven feedback adopted here ($V_{\rm vir}$ is the virial velocity of the halo). We next tried the momentum feedback with the $V_{\rm vir}^{-1}$ dependence. This choice did not give a satisfactory agreement either. This result means that raising the ejection efficiency above the fiducial value is not a promising replacement for preventive feedback.

This in turn casts a doubt that the efficiency of 20% adopted for the fiducial ejective feedback may already be too large. In order to check this point, we examined the mass loading factor (i.e., the ratio of the mass ejection rate and the SFR, hereafter MLF) in the present model. Figure 19 compares our fiducial-scheme model with available observations and other theoretical models. The MLFs for the efficiency $\varepsilon_{\rm FB}=0.2$ roughly border the upper envelope for the observational data. In the low mass domain where feedback has dominant effect, the semi-analytic model by Lagos, Lacey & Baugh (2013) and Illustris simulation by Vogelsberger et al. (2013) give higher MLFs than our model while other cosmological simulations and the steady state models by Hopkins, Quataert & Murray (2012) show lower values.

Considering this result, we lowered the efficiency to $\varepsilon_{\rm FB}=0.05$. This brings the model in better agreement with the observation and the theoretical results. For this value, massive halos suffer from negligible effect of ejective feedback. In this case, however, we need much stronger action of preventive feedback that reduces the gas accretion rate to $\sim 5\%$ instead of the adopted 30% in order to match the stellar mass in low-mass galaxies. This leads to the contribution in the stellar mass reduction being comparable between two types of feedback while the contribution of ejective feedback is about twice that by preventive feedback in the $\varepsilon_{\rm FB}=0.2$ case. It is not clear if such a strong preventive feedback is possible. Because the observationally derived MLFs likely give only lower limits (but also see the discussion by McQuinn, van Zee & Skillman (2019)), the ejective efficiency of $20\%$ may not necessarily be unrealistic. Incidentally, Dutton & van den Bosch (2009) found that $\sim 25\%$ of the supernova energy is required to reproduce the scaling relations of disc galaxies although this conclusion was derived for the shock-heating picture. In any case, it seems robust that we need a significant contribution from preventive feedback to explain the stellar mass budget over the whole mass range.

The importance of prevented feedback has been stressed also from the chemical point of view (e.g. Mitra, Davé & Finlator, 2015; Christensen et al., 2016). It is well known that there is a tension between low efficiency of making stars and retaining enough metal in low-mass galaxies. The model that relies solely on ejective feedback in reducing stellar masses produces a too steep mass-metallicity relation compared with the observation (e.g. Torrey et al., 2014; Somerville, Popping & Trager, 2015). Preventive feedback can reduce the stellar mass in low mass galaxies with their metallicity nearly untouched (e.g. Torrey et al., 2014). For example, the cosmological simulation by Christensen et al. (2016) that includes preventive feedback produces the mass-metallicity relation at $z=0$ with the slope and normalization agreeing with the observation.

We have seen above that preventive feedback is required also to reproduce the stellar-to-halo mass ratio for low mass galaxies. Nevertheless, preventive feedback is here introduced in an ad-hoc manner tuned for better reproduction of the observational data. It is important to check the plausibility of the introduced prescription to describe this process. Because there is no strong constraint from the observation at the moment, we compare the outcome of our prescription to the cosmological simulations addressing this process. The red line in the left panel of Figure 20 indicates the critical halo mass for preventive feedback adopted in this study. Two asterisks indicate the halo mass below which the inclusion of stellar feedback reduces the gas accretion rate on to halos in the simulation by van de Voort et al. (2011). Simulations by Correa et al. (2018) and Wright et al. (2020) give similar results. Two triangles denote the mass below which the prevented gas starts to dominate the total baryon mass budget in EAGLE simulation performed by Mitchell & Schaye (2022). Our prescription agrees with these cosmological simulations at least up to $z\sim 2$. The FIRE-2 cosmological simulation by Hopkins et al. (2018) and the SAM by Pandya et al. (2020) also suggest the critical mass for preventive feedback caused by supernovae that decreases with redshift. Their critical masses lie in the range $10^{11}-10^{12}{\rm M}_{\odot}$ in agreement with our prescription and other cosmological simulations. Analytical prescription considered in Davé, Finlator & Oppenheimer (2012) adopts systematically lower efficiency than our choice over the redshift range $0<z<2$ and applies to lower halo masses. Nevertheless, their recipe dictates stronger feedback for higher redshifts and less massive halos, which is qualitatively similar to our recipe.

Right panel of Figure 20 shows that the decrease in the gas accretion due to preventive feedback is around $30\sim 40\%$ for $M_{\rm vir}<10^{12}{\rm M}_{\odot}$ in our fiducial-scheme model. This diminishment is comparable to the one obtained in the cosmological simulation by Christensen et al. (2016). For the stellar mass, our model gives much smaller values but in agreement with the empirical estimates for the stellar-to-halo mass ratios. This fact may lend some support to the recipe adopted in the present work.

Lu et al. (2017) parametrized the strength of preventive feedback as a function of halo mass and redshift by fitting their semi-analytic model to the stellar mass function and the mass-metallicity relation of the satellite galaxies of Milky Way. Dashed and solid purple lines in the left panel indicate the halo mass for which the accretion to the halo is reduced to respectively $99\%$ and $30\%$ with respect to $f_{\rm bar}\dot{M}_{\rm vir}$. These masses are much lower than our recipe and other studies mentioned above, rather in agreement with the boundary for UV background suppression given by Gnedin (2000) and Okamoto, Gao & Theuns (2008) (see Figure 1). Their model thus seems to imply that the background radiation field is sufficient in preventing gas accretion in the Milky Way system. By contrast, other studies suggest a close link between ejective and preventive feedback. For example, Oppenheimer et al. (2010), based on the numerical simulations, argue that the suppression of gas accretion by supernova wind energy and momentum has a strong effect for low mass galaxies. van de Voort et al. (2011) and Muratov et al. (2015) also appear to support this picture.

The suggested dominance of preventive feedback for low mass halos may demand significant revision of the cold accretion picture for the low-mass regime. Possible interactions between the outgoing gas ejecta and the inflowing cold gas may change the physical state of the halo gas drastically. In other words, preventive feedback by definition invalidates the assumption of the universal baryon fraction in parent halos and/or the gas accretion rate (on to discs) calculated from either the free-fall time or the radiative cooling time. Nevertheless, we could gauge the required suppression of gas accretion on to galaxies relative to the accretion rate calculated from the universal baryon fraction and appropriate inflow timescales (also see Appendix, section A6).

The present study demonstrates the ability of simple isolated models in reproducing the observed diversity in galactic star formation history. This success is attained by assuming the cold accretion regime of halo gas aquisition and cooperation of ejective feedback, preventive feedback and gas recycling. This study reveals at the same time several shortcomings of these idealized models. The most notable disagreement with the observational data is that the models underpredict SFR at high redshifts. At low redshifts, massive galaxies show higher SFRs than observed. Here we discuss several processes which are neglected in our models but may resolve or alleviate these problems. Specifically, we touch on here AGN feedback, galaxy mergers and wind transfer.

The most promising process to reduce high SFR in massive galaxies in recent epochs is feedback from AGNs. It has been considered to be one of important processes that potentially affect the galaxy evolution and has been taken into account in many theoretical studies (e.g. Sijacki et al., 2007; Vogelsberger et al., 2013; Weinberger et al., 2018). These studies (both the SAMs and cosmological simulations) succeeded in explaining the quenching of massive galaxies by introducing AGN feedback (e.g. Kauffmann & Haehnelt, 2000; Croton et al., 2006; de Lucia et al., 2010; Lu et al., 2014). Quenching mechanisms can act in various forms. For example, radio mode feedback can act as preventive feedback (e.g, Kereš et al., 2009; Zinger et al., 2020) and may lead to ’halo quenching’ as discussed by Lu et al. (2014). The accumulation of AGNs on way from the blue cloud to the red sequence for late type galaxies (e.g. Schawinski et al., 2010; Powell et al., 2017) may suggest AGN-quenching for star forming galaxies from the observational viewpoint. McPartland et al. (2019) report that the inclusion of AGNs (LINERs and Seyfert 2 galaxies) into observational analysis leads to the turn-over of the MS at high masses. On the other hand, de Lucia et al. (2010), using a semi-analytic model, argue that the necessity of AGN feedback is changed by various assumptions adopted for gas cooling. Recent cosmological simulations (e.g. Correa et al., 2018; Nelson et al., 2019; Wright et al., 2020) give more quantitative results for the effect of AGN feedback on gas accretion rates. For example, Wright et al. (2020) indicate that stellar feedback has much stronger preventative effect than AGN feedback (factor $\sim 4$ versus $20\sim 30\%$). They infer that this difference is due partially to the former acting preferentially in low mass halos that have shallow gravitational potentials while the latter operates in massive halos having deep potential wells. Our result also suggests that preventive feedback is more important in low mass halos in line with these simulation result (e.g., Figure 16). Probably, one of the most relevant results to the current issue is provided by the EAGLE simulation (Katsianis et al., 2017) which shows that AGN feedback at low redshifts improves the agreement of the simulation with the observed evolution of the cosmic SFR density (CSFRD) (see their Figure 10). Although CSFRD is the integrated quantity over all galaxy populations, this result hints at important role of AGN feedback in quenching massive galaxies in recent epochs. This result suggests that inclusion of AGN feedback could help cure unrealistic upturn in SFR in low-redshift massive galaxies seen in the present fiducial scheme.

Regrettably, most observational works conducted to date (e.g. Schawinski et al., 2010; McQuinn, van Zee & Skillman, 2019) give only circumstantial evidence for AGN quenching. More direct information could be provided by measuring the mass loading factor for AGN-driven outflows. For example, Qiu et al. (2021) calculated the ratio of outflow rate and black hole accretion rate for massive elliptical galaxies. If some scaling relation between this ratio and the BH accretion rate or BH mass is established by extending such analysis to larger samples including late-type galaxies, it could be combined with the BH growth model to provide better theoretical insight into the importance of AGN feedback for different galaxy masses and cosmological epochs.

On the other hand, overly inactive star formation in model galaxies at high redshift could be alleviated by galaxy mergers which are expected to be frequent at early cosmological epochs. Along with AGN feedback, galaxy mergers have also been considered as an important process affecting mainly massive galaxies. They can provide a promising route for morphological transformation and star formation quenching, both processes possibly leading to the formation of early-type (elliptical) galaxies. One concern with this scenario is that mergers may not quench star formation completely. For example, the simulation by Springel, Di Matteo & Hernquist (2005) indicates that the aid from AGN feedback is necessary to expell the residual ISM completely outside merging galaxies. Despite this caveat, it seems evident that galaxy mergers (especially major ones) have non-negligible effects on the evolution of galaxy populations. Including mergers is difficult in our modelling by construction. Therefore, the present model should be, in the strict sense, regarded as mainly describing evolution of disc-dominated galaxies that have probably suffered little from the effects of major galaxy mergers in their lives. Although the merger-driven starbursts in disc galaxies are thus not handled in our current model, it seems reasonable to expect that frequent mergers in early times would enhance star formation activity by perturbing young gas-rich galactic discs.

Recent cosmological simulations (e.g. Angléz-Alcázar et al., 2017; Grand et al., 2019; Mitchell, Schaye & Bower, 2020) point out the importance of intergalactic transfer of material such as ’wind transfer’ in addition to mergers of galaxies themselves. These simulations suggest comparable contribution of this process to other processes such as gas recycling (fountain flow) and fresh gas accretion. Our model is unable to predict how these extragalactic ingredients affect galaxy evolution. What is naturally expected is that the idealized concept of hot and cold modes of gas accretion depicted in Figure 1 may be significantly modified (in addition to possible modification posed by preventive feedback mentioned in section 6.2). Currently, the definition of these newly introduced processes itself varies among different authors and mutual comparison of simulation results is not straightforward. It is challenging but worthwhile to forge various processes occuring on halo scales into a unified picture as the next target.