How feedback shapes galaxies: an analytic model

TL;DR

This paper presents a simple analytic model linking dark matter halo growth to galaxy stellar mass build-up, successfully reproducing key galaxy evolution observables and explaining the galaxy stellar mass function's shape.

Contribution

It introduces an analytic, time-independent star formation efficiency model that captures galaxy-halo co-evolution and explains the double Schechter function.

Findings

Reproduces cosmic star formation rate density evolution.

Explains the double Schechter function in galaxy stellar mass distribution.

Identifies two efficiency regimes: feedback and black hole regulation.

Abstract

We introduce a simple analytic model of galaxy formation that links the growth of dark matter haloes in a cosmological background to the build-up of stellar mass within them. The model aims to identify the physical processes that drive the galaxy-halo co-evolution through cosmic time. The model restricts the role of baryonic astrophysics to setting the relation between galaxies and their haloes. Using this approach, galaxy properties can be directly predicted from the growth of their host dark matter haloes. We explore models in which the effective star formation efficiency within haloes is a function of mass (or virial temperature) and independent of time. Despite its simplicity, the model reproduces self-consistently the shape and evolution of the cosmic star formation rate density, the specific star formation rate of galaxies, and the galaxy stellar mass function, both at the present time and at high redshifts. By systematically varying the effective star formation efficiency in the model, we explore the emergence of the characteristic shape of the galaxy stellar mass function. The origin of the observed double Schechter function at low redshifts is naturally explained by two efficiency regimes in the stellar to halo mass relation, namely, a stellar feedback regulated stage, and a supermassive black hole regulated stage. By providing a set of analytic differential equations, the model can be easily extended and inverted, allowing the roles and impact of astrophysics and cosmology to be explored and understood.

Tables

Table 1. Table 1: Effective star formation efficiency parameters for the six idealised models. To agree with observational data, the values of α 𝛼 \alpha and β 𝛽 \beta are typically positive. The “Fiducial” model captures this behaviour while the five variations systematically explore the effect of the effective star formation efficiency on the physics of galaxy formation.

	${ϵ}_N$	$M_{\mathrm{crit}}$	$\alpha$	$\beta$
Fiducial	0.125	${10}^{12}$	0.75	0.75
Constant	0.250	${10}^{12}$	0	0
No AGN	0.125	${10}^{12}$	0.75	0
No SN	0.125	${10}^{12}$	0	0.75
$M_{\mathrm{crit}}={10}^{10}$	0.125	${10}^{10}$	0.75	0.75
High Efficiency	0.320	${10}^{12}$	0.75	0.75

Table 2. Table 2: Best fit parameters for the halo mass-dependent ( Model I , ϵ ∗ ​ ( M h ) subscript italic-ϵ subscript 𝑀 ℎ \epsilon_{*}({M_{h}}) ), and virial temperature-dependent ( Model II , ϵ ∗ ​ ( T vir ) subscript italic-ϵ subscript 𝑇 vir \epsilon_{*}({T_{\mathrm{vir}}}) ) star formation efficiency models. For the virial temperature-dependent model, M crit subscript 𝑀 crit M_{\mathrm{crit}} is given at redshift z = 0 𝑧 0 z=0 , which corresponds to a critical virial temperature T crit = 10 5.3 ​ K subscript 𝑇 crit superscript 10 5.3 K T_{\mathrm{crit}}=10^{5.3}\mathrm{K} . As T crit subscript 𝑇 crit T_{\mathrm{crit}} is kept constant, M crit ∝ a ​ ( t ) 3 / 2 proportional-to subscript 𝑀 crit 𝑎 superscript 𝑡 3 2 M_{\mathrm{crit}}\propto a(t)^{3/2} (see • ‣ Section   2.1 ). χ ν 2 superscript subscript 𝜒 𝜈 2 \chi_{\nu}^{2} is the reduced chi-squared statistic used for goodness of fit testing.

Model	${ϵ}_N$	$M_{\mathrm{crit}}\left[M_{\odot }\right]$	$\alpha$	$\beta$	${\chi }_{\nu }^2$
Mass-dependent	0.178	${10}^{11.68}$	1.537	0.656	1.5
Temp-dependent	0.140	${10}^{12}$	2.377	0.834	1.6

Equations

$${\log }_{10}\left(\frac{M_{\ast }}{10^{12}{M}_{\odot }}\right)=\epsilon (M_h,t){\log }_{10}\left(\frac{M_h}{10^{12}{M}_{\odot }}\right)+{\log }_{10}\mathcal{N}(t),$$

$$\varphi (M_{\ast })\equiv \frac{\text{differential}{n}_{\mathrm{g}al}}{\text{differential}{\log }_{10}{M}_{\ast }}={\epsilon }^{-1}\frac{\text{differential}{n}_{\mathrm{h}}}{\text{differential}{\log }_{10}{M}_h},$$

$${ϵ}_{\ast }=\frac{{\dot{M}}_{\ast }}{f_b{\dot{M}}_h}$$

$${ϵ}_{\ast }(M_h)=2{ϵ}_{\mathrm{N}}{\left[{\left(\frac{M_h}{M_{\mathrm{crit}}}\right)}^{-\alpha }+{\left(\frac{M_h}{M_{\mathrm{crit}}}\right)}^{\beta }\right]}^{-1},$$

$$M_{\ast }=$$

$$=$$

$$F(\eta ,z)=$$

$$=$$

$$\epsilon =\frac{\mathrm{dlog}{M}_{\ast }}{\mathrm{dlog}{M}_h}=\frac{1+\alpha }{(1+z)F(\eta ,z)}.$$

$${ϵ}_{\ast }(T_{\mathrm{vir}})=2{ϵ}_{\mathrm{N}}{\left[{\left(\frac{T_{\mathrm{vir}}}{T_{\mathrm{crit}}}\right)}^{-\alpha }+{\left(\frac{T_{\mathrm{vir}}}{T_{\mathrm{crit}}}\right)}^{\beta }\right]}^{-1},$$

$$M_h=200\frac{4\pi }{3}{R}_{200\mathrm{m}}^3\bar{\rho },$$

$$T_{\mathrm{vir}}=\frac{\mu m_{p}G{M}_h}{5k_{\mathrm{B}}{R}_{200\mathrm{m}}},$$

$$f_{\Lambda }(t,A,B)=1+A{\left(\frac{t}{t_{\Lambda }}\right)}^2+B{\left(\frac{t}{t_{\Lambda }}\right)}^4,$$

$$\frac{M_h(t)}{10^{12}{\mathrm{M}}_{\odot }}=$$

$${\left\{{\left(\frac{M_0}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{-\gamma /2}+0.31\gamma \left[{\left(\frac{t}{t_m}\right)}^{-2/3}{f}_{\Lambda }(t,0.16,-0.01)-1.67\right]\right\}}^{-2/\gamma }$$

$$\frac{\text{differential}n(M_h,t)}{{\mathrm{dlog}}_{10}{M}_h}$$

$${\left(\frac{t}{t_m}\right)}^{-2/3}{f}_{\Lambda }(t,0.16,-0.01)\times$$

$$\exp \left[-0.13{\left(\frac{M_h}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{\gamma }{\left(\frac{t}{t_m}\right)}^{-4/3}{f}_{\Lambda }(t,0.32,0)\right].$$

$$\varphi (M)={\varphi }^{\ast }\text{quantity}(\frac{M}{M^{\ast }}{)}^{\alpha }{e}^{-M/M^{\ast }},$$

$$\varphi (M)=\text{quantity}[{\varphi }_1^{\ast }\text{quantity}(\frac{M}{M^{\ast }}{)}^{{\alpha }_1}+{\varphi }_2^{\ast }\text{quantity}(\frac{M}{M^{\ast }}{)}^{{\alpha }_2}]e^{-M/M^{\ast }}.$$

$$\frac{\text{differential}{M}_{\ast }}{\text{differential}t}$$

$$={ϵ}_{\ast }{f}_{b}1.6\times 10^{10}{\mathrm{M}}_{\odot }{Gyr}^{-1}{\left(\frac{t}{t_m}\right)}^{-5/3}\times$$

$$f_{\Lambda }(t,-0.32,0.06){\left(\frac{M_h}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{\left(1+\frac{\gamma }{2}\right)}.$$

$$\mathrm{SFR}$$

$$f_{\mathrm{SFR}}=\{\begin{array}{ll}1& \text{ for }{M}_h\le M_{\mathrm{crit}}\\ \text{quantity}(\frac{M_h}{M_{\mathrm{crit}}}{)}^{\eta }& \text{ for }{M}_h>M_{\mathrm{crit}},\end{array}$$

$$\eta ={\log }_{10}(0.5)/(13-12)=-\mathrm{0.3.}$$

$$\mathrm{SFR}$$

$${\dot{\rho }}_{{}_{\mathrm{SFR}}}(t)$$

$$=\int {ϵ}_{\ast }{f}_b{\dot{M}}_h\frac{f_{\mathrm{SFR}}}{(1-R)}\frac{\text{differential}n(M_h,t)}{{\mathrm{dlog}}_{10}{M}_h}{\mathrm{dlog}}_{10}{M}_h.$$

$$\frac{\text{differential}{\dot{\rho }}_{{}_{\mathrm{SFR}}}}{{\mathrm{dlog}}_{10}{M}_h}$$

$${\left(\frac{M_h}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{\gamma }{\left(\frac{t}{t_m}\right)}^{-7/3}{f}_{\Lambda }(t,-0.16,0)\times$$

$$\exp \left[-0.13{\left(\frac{M_h}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{\gamma }{\left(\frac{t}{t_m}\right)}^{-4/3}{f}_{\Lambda }(t,0.32,0)\right],$$

$$a^{{\color[rgb]{1,0,0}c^{*}}}(t)\propto\left[\frac{3}{2}\frac{t}{t_{m}}\right]^{2/3}\left(1+\frac{1}{4}\left(\frac{t}{t_{\Lambda}}\right)^{2}+\frac{1}{80}\left(\frac{t}{t_{\Lambda}}\right)^{4}+...\right),$$

$$t^{{\color[rgb]{1,0,0}c^{*}}}_{\mathrm{m}}=\sqrt{\frac{3}{8\pi G{\rho}_{0}}}=\frac{1}{{H}_{0}\sqrt{{\Omega}_{\mathrm{m},0}}},$$

$$t^{{\color[rgb]{1,0,0}c^{*}}}_{\Lambda}=\sqrt{\frac{3}{\Lambda c^{2}}}=\frac{1}{H_{0}\sqrt{\Omega_{\Lambda,0}}}.$$

Full text

††pubyear: 2019††pagerange: How feedback shapes galaxies: an analytic model–B

1 Introduction

The co-evolution between galaxies and their haloes is perhaps one of the most fundamental aspects of every galaxy formation model. In the current paradigm of galaxy formation, every galaxy forms within a dark matter halo. However, understanding the relationship between a dark matter halo and the galaxies it hosts is not a trivial exercise due to our lack of detailed understanding of the complex baryonic process involved.

In a standard Lambda Cold Dark Matter ($\Lambda$CDM) cosmology, gravitationally bound dark matter structures build up hierarchically, by a combination of the smooth accretion of surrounding matter and continuous merging with smaller structures (White & Rees, 1978; Qu et al., 2017). The formation and evolution of galaxies within these haloes is thought to be a highly self-regulated process, in which galaxies tend to evolve towards a quasi-equilibrium state where the gas outflow rate balances the difference between the gas inflow rate and the rate at which gas is locked up in stars and black holes (BHs) (e.g. White & Frenk, 1991; Finlator & Davé, 2008; Bouché et al., 2010; Schaye et al., 2010; Davé et al., 2012a; Bower et al., 2017). Consequently, galaxy formation is thought to be determined on the one hand by the formation and growth of dark matter haloes, which depends solely on the cosmological background, and on the other hand, by the regulation of the gas content in these haloes, which in turn depends on complex baryonic processes such as radiative cooling, stellar mass loss, and feedback from stars and accreting BHs. This co-evolution process results in a tight correlation between the properties of galaxies and their dark matter haloes (see e.g. Wechsler & Tinker 2018 for a review).

A fundamental requirement for a successful galaxy formation model, is to reproduce the relation between stellar mass and halo mass inferred from observations. However probing the dark matter distribution and its evolution represents an observational challenge. Direct observational probes include galaxy-galaxy lensing (e.g. Brainerd & Specian, 2003; Hoekstra et al., 2004; Hudson et al., 2015) and the kinematics of satellite galaxies (e.g. Zaritsky et al., 1993; van den Bosch et al., 2004; Norberg et al., 2008). However, direct observation techniques are limited to low redshifts ($z<1$), due to the difficulty of resolving individual distant galaxies. Indirect methods include, for example, comparing the abundance and clustering properties of galaxy samples with predictions from phenomenological halo models (e.g. Neyman & Scott, 1952; Berlind & Weinberg, 2002; Cooray & Sheth, 2002; Cowley et al., 2018). This method however, depends heavily on the underlying modelling and assumptions, for example, the bias with which haloes trace the underlying matter distribution.

From the theoretical point of view, the formation and evolution of dark matter haloes is largely considered a “solved problem” (see however, van den Bosch et al. 2018). Using extremely accurate measurements of the density perturbations imprinted onto the cosmic microwave background radiation fluctuations as initial conditions (e.g. Planck Collaboration et al., 2016), many different groups have produced convergent results using large cosmological N-body simulations (e.g Springel et al., 2005; Klypin et al., 2011; Trujillo-Gomez et al., 2011; Angulo et al., 2012; Fosalba et al., 2015).

On the other hand, the complex physics of galaxy formation still has many open questions. Different approaches have been used to model the intricate baryonic physics of galaxy formation. The most widely used technique combines the evolution of dark matter with either a semi-analytical (e.g. Cole et al., 1994; Somerville et al., 2008; Henriques et al., 2015; Lacey et al., 2016) or hydrodynamical (e.g. Vogelsberger et al., 2014; Schaye et al., 2015; Davé et al., 2016; Dubois et al., 2016; Pillepich et al., 2018a) treatment of the baryonic processes involved. A key ingredient in both methods that has led us to a better understating of the physics of galaxy formation is the use of physically motivated models for feedback processes (see Somerville & Davé 2015; Naab & Ostriker 2017 for a comprehensive review).

An alternative approach known as empirical modelling takes the advantage of the vast number of observational data sets from large galaxy surveys and relate statistical galaxy scaling relations to the evolution of dark matter haloes without assuming strong physical priors (e.g. Behroozi et al., 2013; Moster et al., 2013; Rodríguez-Puebla et al., 2016; Behroozi et al., 2019; Moster et al., 2018; Grylls et al., 2019).

While all of these approaches have been very productive, the increasing complexity of the models and simulations make it difficult to pinpoint and understand the fundamental physics driving the results. For instance, in cosmological simulations, “sub-grid” physics are implemented as micro phenomena that depend only on local gas properties from which macroscopic patterns emerge. However, it is hard to track down the link between what emerges and why (e.g. Bower et al., 2017). In this paper, we examine this issue in detail by adopting the opposite approach. We develop a fully analytic model of galaxy formation derived from a simple relation between the star formation rate and halo growth rate that disentangles the role of cosmology from the role of astrophysics in the galaxy formation process. Our model restricts the role of baryonic astrophysics to setting the relation between galaxies and their haloes. With this simple relation, we can use an analytic approximation to the growth of dark matter haloes to predict galaxy properties. By providing a set of analytic equations, the model can be easily “inverted” and allows for rapid experiments to be conducted, providing a powerful tool to explore the differential effects of baryonic physics, averaged over galaxy scales. Despite its simplicity, the model reproduces self-consistently the shape and evolution of the cosmic star formation rate (SFR) density, the specific star formation rate (sSFR) of galaxies, and the galaxy stellar mass function (GSMF), both at the present time and at high redshift.

We validate our results by comparing to numerical hydrodynamic simulations from the eagle project. The eagle simulation suite111The galaxy and halo catalogues of the simulation suite, as well as the particle data, are publicly available at http://www.eaglesim.org/database.php (McAlpine et al., 2016). (Schaye et al., 2015; Crain et al., 2015) consists of a large number of cosmological hydrodynamical simulations that include different resolutions, simulated volumes and physical models. These simulations use advanced smoothed particle hydrodynamics (SPH) and state-of-the-art subgrid models to capture the unresolved physics. A complete description of the code and physical parameters used can be found in Schaye et al. (2015). Here we compare to the eagle reference simulations that used a flat, $\Lambda$CDM cosmology with parameters ($\Omega_{m}=0.307$, $\Omega_{\Lambda}=0.693$, $h=0.6777$, $\sigma_{8}=0.8288$, $n_{s}=0.9611$) consistent with the Planck Collaboration et al. (2014) results. The calibration strategy of the eagle simulations is described in detail by Crain et al. (2015), who also presented additional simulations to demonstrate the effect of parameter variations.

The layout of this paper is as follows: In Section 2 we introduce the analytic model of galaxy formation. We present two models of the effective star formation efficiency: A time-independent efficiency which depends only on halo mass, and an efficiency that depends on the virial temperature of the halo. In Section 3, we explore the effect of the different efficiency parameters in the galaxy formation outputs. Namely, the cosmic star formation rate density, the specific star formation rate of galaxies, and the galaxy stellar mass function. In Section 4 we compare the results from our model to different observational datasets. We also discuss the need for a time-evolving efficiency in order to reproduce the rapid evolution of the GSMF. We discuss the limitations of our model, and summarise our conclusions in Section 5.

2 An analytic model of galaxy formation

2.1 The effective star formation efficiency

The formation, evolution and abundance of dark matter haloes can be predicted accurately when the cosmology and dark matter model (i.e. cold, warm, self-interacting, etc.) is known. Although these processes are highly non-linear, the underlying physics is well understood. However, the gas and stellar content of haloes is much less well understood because of the intrinsic complexity of the baryonic processes, such as cooling, star formation and feedback, that drive it. An empirical approach to populating dark matter haloes with galaxies is to focus on the relation between stellar mass and halo mass inferred from observations. We write this relation as

[TABLE]

where $M_{*}$ and $M_{h}$ are the central galaxy stellar mass and host halo mass respectively, $\mathcal{N}(t)$ is a normalisation factor, and $\varepsilon(M_{h},t)$ is the logarithmic slope of the stellar to halo mass relation (SHMR). Allowing $\mathcal{N}(t)$ to be a random variable would account for the scatter in the relation, but in this paper we will focus on the mean relation and replace $\mathcal{N}$ by its expectation value. $\varepsilon$ is closely related to the galaxy formation efficiency of haloes, and we will explore this connection further below.

Probing the dark matter distribution and its evolution directly represents an observational challenge. Perhaps the simplest, and most commonly used alternative technique is to use (sub-) halo abundance matching to determine the typical SHMR (e.g. Behroozi et al., 2013; Moster et al., 2013). In essence, the SHMR is derived by mapping the theoretical halo mass function and the observed abundance of galaxies given by the GSMF,

[TABLE]

where $n_{\mathrm{g}al}$ and $n_{\mathrm{h}}$ are the co-moving abundances of galaxies and haloes respectively. A subtlety here is that $M_{*}$ refers to the stellar mass of the central object in the halo. More complex formulations of the abundance matching method allow for the contribution of satellite galaxies to the mass function, but we will keep to the simple approach. This is adequate if the stellar mass function is dominated by central galaxies (e.g. Yang et al., 2009; Lan et al., 2016).

Abundance matching studies have consistently shown that a simple picture of the galaxy population is consistent with much of the observational data. The SHMR is a strong function of halo mass but a weak function of cosmic time; it can be approximated well by two power laws that connect at a stellar mass that corresponds roughly to the knee of the stellar mass function (e.g. Moster et al., 2010; Yang et al., 2012; Mitchell et al., 2016). We use this as the basis of a simple model that couples the build-up of dark matter haloes and the build-up of galaxy stellar mass.

We must be careful, however, to distinguish the instantaneous efficiency with which infalling baryons are converted into stars,

[TABLE]

(where $f_{b}={\Omega_{b}}/{\Omega_{m}}$ is the cosmic baryon fraction) and the integral of this growth over the history of the halo, Eq. 1. Note that $\dot{M}_{*}$ includes both star formation within the central object, and the accretion of infalling stars. We will need to distinguish between the two in order to relate the stellar mass growth to the observed star formation rate density.

In Appendix A, we show that the build up of dark matter haloes can be described analytically using recent developments of the linear theory original described by Press & Schechter (1974). This allows the abundance and growth rates of haloes to be derived from the power spectrum of density fluctuations in the early universe. Equation 2 provides a promising approach to connect the growth of haloes to the formation of galaxies, and the observational results suggest that a good starting point is to consider a time-independent SHMR that depends only on halo mass. In addition, we will examine a model in which the efficiency of star formation depends on the halo’s virial velocity, as it reflects the evolution of the gravitational potential of the halo (e.g. Sharma & Theuns, 2019).

• •

Halo mass-dependent efficiency (Model I)

In order to model $\epsilon_{*}$ in Eq. 3, we begin by assuming that the efficiency of conversion of infalling baryons into stars depends only on halo mass (e.g. Rodríguez-Puebla et al., 2016; Salcido et al., 2018; Tacchella et al., 2018). Motivated by the results from abundance matching techniques (e.g. Behroozi et al., 2013; Rodríguez-Puebla et al., 2016), that estimate a galaxy formation efficiency that peaks at masses similar to Milky-Way sized haloes (${\sim}10^{12}\,{\hbox{M}_{\odot}}$) and falls steeply for higher and lower masses, we model $\epsilon_{*}$ as a double power law, a similar parametrisation as that proposed by Moster et al. (2010):

[TABLE]

where $\epsilon_{\mathrm{N}}$ is a normalisation parameter, and $\alpha$ and $-\beta$ are the power-law slopes at low and high mass respectively. The maximum efficiency occurs at halo mass $M_{\mathrm{crit}}$. To agree with observational data, the values of $\alpha$ and $\beta$ are typically positive, i.e. at low masses, star formation is suppressed because of the efficiency of feedback from star formation, and at higher masses, the cooling of the inflowing gas is suppressed by heating from BHs (e.g. White & Frenk, 1991; Bower et al., 2006; Benson, 2012; Haas et al., 2013).

Because of our assumption that $\epsilon_{*}$ depends only on $M_{h}$, we can integrate to determine $M_{*}$ without needing to know the time evolution of the halo mass.

[TABLE]

where $\eta=(1+\alpha)/(\alpha+\beta)$, $\mathbb{z}=\left({M_{h}}/{M_{\mathrm{crit}}}\right)^{\alpha+\beta}$ and

[TABLE]

where222We have used the symbol $\mathbb{z}$ to differentiate from redshift $z$. ${}_{2}F_{1}(a,b;c;\mathbb{z})$, is the Gaussian hypergeometric function. For values of $\alpha>0$ and $0<\beta<1$, in the limit $M_{h}{\ll}M_{\mathrm{crit}}$, $\lim_{\mathbb{z}\to 0}F(\eta,\mathbb{z})=1$, while for $M_{h}{\gg}M_{\mathrm{crit}}$, $F(\eta,\mathbb{z})\approx 1/((\eta-1)\,\mathbb{z})$. Differentiating Eq. 5, the logarithmic slope of the stellar mass halo mass relation, $\varepsilon$, can be written in a simple analytic form,

[TABLE]

Equation 7 describes a smooth transition in slope from $(1+\alpha)$ for $M_{h}{\ll}M_{\mathrm{crit}}$, to $(1-\beta)$ for $M_{h}{\gg}M_{\mathrm{crit}}$.

Figure 1 shows an illustration of the effective star formation efficiency $\epsilon_{*}$ as a function of $M_{h}$. The parametrisation provides a smooth a transition from the $\alpha$ dominated regime (for low halo masses), to the $\beta$ dominated regime (for high halo masses). The figure shows that $\epsilon_{*}$ has the same slopes as $M_{*}/M_{h}$, i.e. $\alpha$ and $\beta$, but the normalisation of $M_{*}/M_{h}$ is different by a factor of $1/(1+\alpha)$, and $1/(1-\beta)$ for low and high mass haloes respectively. We note that for our chosen parametrisation, $M_{*}/M_{h}$ is closely approximated by a double power law.

• •

Virial temperature-dependent efficiency (Model II)

As we will discuss in Section 3, an efficiency dependent only on halo mass turns out to be a very good approximation of the stellar mass build up of galaxies because most of the stellar mass builds up when the mass of the halo has roughly its current value. However, a time-independent efficiency model significantly under predicts the abundance galaxies at high redshifts ($z>4$), which hints at the need for a time-evolving efficiency model. A purely empirical approach (e.g Moster et al., 2018; Behroozi et al., 2019) would relax the physical priors and let, in this case, the four efficiency parameters in the model to evolve freely in time. Instead, we consider an alternative model in which $\epsilon_{*}$ depends only on the virial temperature of the halo $T_{\mathrm{vir}}$ (and hence the gravitational potential of the halo).

Considering the energetics of galaxy winds suggests that, wind that marginally escape the gravitational binding energy of the galaxy’s halo can carry a higher mass loading in lower mass haloes (Dekel & Silk, 1986; White & Frenk, 1991; Davé et al., 2012b; Sharma & Theuns, 2019). Since the energy required for escape depends on the halo virial temperature, $T_{\mathrm{vir}}$, this leads to an inverse scaling of the mass loading, and hence, star formation efficiency. At sufficiently high mass, the energy associated with individual supernovae becomes smaller than the halo binding energy. This may lead to the accumulation of gas around the central BH, and consequently a wind driven by BH accretion, rather than star formation (Dubois et al., 2015). A related argument can also be made based on the buoyancy of gas heated by star formation. Bower et al. (2017) discusses a possible physical origin of a transition from where star formation driven outflows get hotter than the virial temperature of the halo and can escape (i.e. supernovae energy, or entropy, is much greater than the halo binding energy), to where outflows stall inside the halo triggering star formation and BH growth.

In order to explore these effects, we model the effective star formation efficiency as a function of the halo’s virial temperature using the same double power law parametrisation as in Eq. 4,

[TABLE]

Note that the relation between halo mass and virial temperature depends on redshift (see Eqs. 9 and 10), as the density of collapsed haloes decrease as the universe expands. As a result, it is not possible to analytically determine $M_{*}(M_{h},t)$, but the required integrals can easily be evaluated numerically. We discuss this model further in Section 4.1.

2.2 Halo definition

Dark matter haloes are typically identified by growing a sphere outwards from the potential minimum of the dark matter halo out to a radius where the mean interior density equals a fixed multiple of the critical or mean density of the Universe, causing an artificial ‘pseudo-evolution’ of dark matter haloes by changing the radius of the halo (Diemer et al., 2013). Star formation, however, is governed by the amount of gas that enters these haloes and reaches their central regions. Wetzel & Nagai (2015) show that the growth of dark matter haloes is subject to this ‘pseudo-evolution’, whereas the accretion of gas is not. Because gas is able to cool radiatively, it decouples from the dark matter, tracking the accretion rate near a radius of $R_{200\mathrm{m}}$, the radius within which the mean density is 200 times the mean density of the universe, $\bar{\rho}$. As we try to connect the mass accretion rate of dark matter haloes to star formation, we define halo masses as the total mass within $R_{200\mathrm{m}}$,

[TABLE]

where $\bar{\rho}(t)={\rho}_{0}a(t)^{-3}$.

We assume that during gravitational collapse, the gas experiences strong shocks and thermalises its kinetic infall energy to the virial temperature of the halo,

[TABLE]

where we have assumed a uniform cloud of monatomic gas. $M_{h}$ is the mass of the halo, $\mu$ is the mean molecular weight of the gas in the halo, which we have assumed $\mu\approx 0.6$ for a fully ionized plasma of primordial composition, $k_{\mathrm{B}}$ is the Boltzmann constant, and $m_{p}$ the proton mass. Note from Eq. 9 that for a given halo mass, the radius of the halo, $R_{200\mathrm{m}}$, changes with time as the mean density of the Universe evolves.

2.3 The model

The analytical galaxy formation model developed here is comprised of three main components, which are summarised in the schematic diagram of Fig. 2:

• •

A cosmological model (blue block)

• •

An astrophysical model that sets the relation between galaxies and their haloes (green block)

• •

The galaxy formation outputs (orange block)

By using the Taylor expansion solution of the Friedmann equations introduced by Salcido et al. (2018), the formation and evolution of dark matter haloes can be described analytically. This component is shown as the blue block in Fig. 2. The growth rates of haloes depend on the cosmological parameters ${\rho}_{0}$, $\Lambda$, $H_{0}$ and the shape of the matter density fluctuation power-spectrum. We parametrise the variance of the spherically averaged smoothed density field, $S=\sigma^{2}$, as a power law $S\approx S_{0}(M_{h}/10^{12}\mathrm{M}_{\odot})^{-\gamma}$ with slope ${\gamma}$. Since we are interested in only a small range of halo mass, this is a sufficiently accurate description. The derivation of these equations are presented in Appendix A. For convenience, we define the cosmology dependent approximations for the equations that appear bellow using the function,

[TABLE]

where $t_{\Lambda}=\sqrt{{3}/{\Lambda c^{2}}}$, $\Lambda$ is the value of the cosmological constant, and the coefficients $A$ and $B$ are obtained by using the Taylor expansion solution of the Friedmann equations in Salcido et al. (2018).

Astrophysical processes (shown in green in Fig. 2), enter into the model through the effective star formation efficiency, which is fully described by the efficiency $\epsilon_{*}$ (Eq. 3). We consider two models, in which $\epsilon_{*}$ is a function of halo mass or virial temperature. This component of the model is described by four parameters, $M_{\mathrm{crit}}$ (or $T_{\mathrm{crit}}$), $\epsilon_{\mathrm{N}}$, $\alpha$ and $\beta$ following Eq. 4 or Eq. 8.

In order to simplify the numerical constants in the equations presented in this section, we have substituted the numerical values for the cosmological parameters for a standard $\Lambda$CDM universe as inferred by the Planck Collaboration et al. (2014), i.e. $\Omega_{\mathrm{m}}{=}0.307$, $\Omega_{\Lambda}{=}0.693$, $\Omega_{b}{=}0.04825$, $H_{0}{=}67.77\,\mathrm{km}\,\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$, $\rho_{0}=3.913\times 10^{10}\,{\hbox{M}_{\odot}}\,\hbox{Mpc}^{-3}$, $\gamma{=}0.3$, $t_{0}{=}13.8\,\hbox{Gyr}$ and $S_{0}{=}3.98$. The full cosmology dependence of the numerical constants is given in Appendix A, and are highlighted using a coloured superscript (^{{\color[rgb]{1,0,0}c^{*}}}).

• •

Halo mass history

The analytic form of the growth rate equations allows us to simply describe the growth of haloes as a function of their present-day mass, $M_{0}$:

[TABLE]

where $t_{\mathrm{m}}=\sqrt{{3}/{8\pi G{\rho}_{0}}}$, and ${\rho}_{0}$ is the mean matter density of the Universe at the present time. For the Planck Collaboration et al. (2014) cosmological parameters, $t_{m}{=}26.04\,\hbox{Gyr}$ and $t_{\Lambda}{=}17.33\,\hbox{Gyr}$. Figure 3 shows the individual mass histories for haloes of a given mass $M_{0}$ at the present cosmic time (represented by the colour coding).

• •

Halo mass function

Using the Press & Schechter formalism (Press & Schechter, 1974), the co-moving abundance of haloes of mass $M_{h}$ at time $t$ is given by,

[TABLE]

This equation consists of two parts, a low-mass power law dependence close to $M_{h}^{-1}$, and an exponential cut-off at high masses. For a given halo mass, the abundance initially increases as the exponential suppression is reduced, but at late times the halo abundance slowly decreases because of the power-law term.

• •

The galaxy stellar mass function and the origin of the Schechter function

The galaxy stellar mass function (GSMF) has been reasonably well measured over much of cosmic time, so that, for a know cosmology, the GSMF provides a good measurement of the efficiency by which haloes convert their baryons into stars. Typically, the GSMF (Eq. 2) is parameterised by a Schechter function (Schechter, 1976),

[TABLE]

where $\phi^{*}$ provides the normalisation, and $M^{*}$ is a characteristic galaxy stellar mass where the power-law form of the function cuts off. The form of this function was originally motivated by the halo mass dependence given in Eq. 13. Importantly, however, the shape GSMF is only indirectly related to the halo mass function (eg., Benson et al. (2003)), with observations showing that the power-law slope is much flatter than that expected for the halo mass function. Moreover, recent measurements of the GSMF at low redshift (e.g. Baldry et al., 2008; Li & White, 2009; Baldry et al., 2012; Moustakas et al., 2013), have proven that a single Schechter function is insufficient to describe the density of galaxies. Specifically, the low redshift GSMF shows a clear evidence for a low-mass upturn, or equivalently, a "pile-up" in the abundance of galaxies around $M^{*}$. Typically, a double Schechter function parametrisation has been used to better describe observational data,

[TABLE]

In the model presented here, the GSMF can be computed as a function of time, by combining the halo mass function (Eq. 13) and the efficiency of star formation through Eq. 2 and Eq. 7. These equations link the observed shape of the GSMF, to the underlying dark matter halo distribution, and hence to the cosmological background. They also link galaxies to their dark matter haloes, providing valuable information about the efficiency by which haloes convert their baryons into stars.

Further consideration shows that they provide a description of the non-trivial shape of the GSMF and the need for a double Schechter function to describe it. While the underlying distribution of dark matter haloes is theoretically predicted to be a single Schechter function (Press & Schechter, 1974), its transformation to the GSMF relies on Eq. 2. When the halo mass function is multiplied by the inverse of the logarithmic slope of the SHMR, the low-mass end is multiplied by a factor $1/(1+\alpha)$, while the high-mass end is multiplied by a factor of $1/(1-\beta)$. As both $\alpha$ and $\beta$ are positive, this creates an a kink in the gradient, shown with a red dot in the top-left panel of Fig. 4. These different changes in the normalisation cause a maximum in the logarithmic derivative of $(\phi_{h}/\varepsilon)$ shown top-right panel of the figure. Hence, there is an inflection point in the distribution, as the second derivative vanishes and changes sign at $\sim M_{\mathrm{crit}}$. At this point, the abundance of galaxies decreases slowly, or even rises, as a function of mass, creating a “bump” at the knee of the GSMF. The sharper the transition (ie., the larger $\alpha+\beta$), the more pronounced the bump at the knee of the GSM. These effects are illustrated in Fig. 4. Physically, this can interpreted as galaxies of similar masses “piling up” at the peak star formation efficiency, i.e. $M_{h}\approx M_{\mathrm{crit}}$, as they rapidly stop forming many more stars. In Section 3, we will systematically vary the four parameter in the efficiency model in Eq. 4 to investigate their effect on different galaxy formation outputs.

• •

The galaxy stellar mass growth

We now have all the necessary ingredients to calculate the stellar mass growth of individual galaxies through cosmic time. Substituting Eq. 34 into Eq. 3, the stellar mass is given by the integral of,

[TABLE]

where $\epsilon_{*}$ is given by either Eq. 4 or Eq. 8.

Assuming an instantaneous recycling approximation (Schmidt, 1963), the relation between the stellar mass growth due to star formation and the observed galaxy SFR is simply given by,

[TABLE]

where $R$ is the fraction of mass of gas that is instantaneously returned into the interstellar medium by an entire stellar generation. For a universal Chabrier (2003) initial mass function (IMF), $R=0.41$.

Furthermore, both in situ star formation and galaxy mergers contribute to the total stellar build up of galaxies. In low mass haloes, most of the stellar build up is expected to come from in situ star formation, while the most massive galaxies experience almost no internal star formation and grow mainly by mergers with smaller satellite galaxies. Hence, the fractional contribution of accreted stars to the total stellar mass build up of galaxies is a steep function of halo mass (e.g. Rodriguez-Gomez et al., 2016; Qu et al., 2017; Pillepich et al., 2018b). Assuming that all of the stellar growth of haloes of mass $M_{\mathrm{crit}}$ and bellow is due to internal star formation, we parametrise the fraction of stellar mass growth from in situ SFR by a broken power law as,

[TABLE]

where $M_{\mathrm{crit}}$ is the effective star formation peak efficiency defined in Eq. 4. For the virial temperature-dependent efficiency model (Model II), $M_{\mathrm{crit}}$ also varies with time, and can be calculated using the critical virial temperature in Eq. 10. We fix the value of $\eta$ by assuming that at redshift $z=0$, where $M_{\mathrm{crit}}\approx 10^{12}\mathrm{M}_{\odot}$, $f_{\mathrm{SFR}}(10^{13}\mathrm{M}_{\odot})$ is $\sim 50\%$ (Pillepich et al., 2018b), hence,

[TABLE]

Putting Eqs. 17 and 18 together, the fraction of stellar mass growth of central galaxies due to in situ formation is given by,

[TABLE]

• •

The cosmic star formation rate density

The total cosmic SFR density is given by the integral of all star formation in all haloes,

[TABLE]

Using the stellar mass growth rate from Eq. 16, the halo mass function from Eq. 13, together with the effective star formation efficiency from Eq. 3, the contribution to the cosmic SFR density from haloes of mass $M_{h}$ (the integrand of Eq. 21) is given by,

[TABLE]

where $\epsilon_{*}$ is modelled using Eq. 4 or Eq. 8. The differential form of Eq. 22 explicitly shows the contribution from haloes of different masses $M_{h}$, to the total cosmic SFR density.

Equations 2, 12, 13, 16, 20 and 22, together with a model of the effective star formation efficiency, Eq. 4 or Eq. 8, provide a full mathematical framework to explore the effects of cosmology and baryonic physics on galaxy formation. In the next section, we will explore the effect of the different efficiency parameters on the galaxy SFR, GSMF and the cosmic SFR density.

3 The impact of the effective star formation efficiency

We now use our model to explore the effect of the different efficiency parameters in the galaxy formation outputs in the orange block of Fig. 2. It is common to characterise galaxy properties over halo masses, and for simplicity, in this section we will only use a halo mass-dependent star formation efficiency model (i.e. $\epsilon_{*}$ is constant across cosmic time).

It has been estimated that the SHMR peaks at masses similar to Milky-Way sized haloes (${\sim}10^{12}\,{\hbox{M}_{\odot}}$). Typically, at low masses, the SFR is suppressed because of the efficiency of stellar feedback. On the other hand, at higher masses the cooling of the inflowing gas is suppressed by heating from supermassive BHs (e.g. White & Frenk, 1991; Bower et al., 2006; Benson, 2012). The “Fiducial” model captures this behaviour with both $\alpha$ and $\beta$ being positive and equal to 0.75.

We consider five alternative models varying the efficiency parameters systematically to explore the physics of galaxy formation. An extreme idealised case label as “Constant”, describes a model where a fixed fraction of the baryon budget is turned into stars, regardless of the halo mass. The “No AGN” model describes a scenario where the efficiency of feedback processes is weak for massive objects. Physically, this could be thought as a model where feedback from active galactic nuclei is inefficient. The “No SN” model describes a scenario where the efficiency of feedback processes is weak in small haloes. Physically, this could be thought as a model where feedback from supernovae is inefficient. While much more computationally expensive, a similar behaviour to these models is reproduced in full hydrodynamical simulations (see Appendix B for a couple of examples). In Fig. 5 we show the median stellar-halo mass ratio for three variations of the eagle simulations where the subgrid prescription for stellar and AGN feedback have been removed. Indeed, our models capture the overall behaviour attained in simulations. An additional model labelled “$M_{\mathrm{crit}}=10^{10}$” explores the effect of changing the critical, or transition, halo mass. A final model labelled “High efficiency” has the same slopes as the fiducial model, but with a different normalisation. We show in Fig. 6 the effective star formation efficiency for the six models, and their parameters are summarised in Table 1.

3.1 The build up of stellar mass

First, we explore the effect of the star formation efficiency in the build up of stellar mass in individual haloes for the six models, which can be calculated by integrating Eq. 16.

Fig. 7 shows an example of the evolution of the stellar mass in a halo of mass $M_{0}=10^{13}\mathrm{M}_{\odot}$ at the present time for the six efficiency models. For the constant efficiency model, the stellar mass grows steadily with time, and starts to slow down only at late times due to cosmic expansion, tracking the mass assembly of the dark matter halo. For the fiducial model, the build up of stellar mass is faster (steeper slope). Once the critical halo mass is reached ($M_{\mathrm{crit}}=10^{12}\mathrm{M}_{\odot}$, corresponding to $M_{*}\approx 10^{10}\mathrm{M}_{\odot}$ for this model), the stellar mass plateaus. The No AGN model has a similar behaviour at early times, but once the critical halo mass is reached, star formation does not slow down and the halo reaches a higher stellar mass at the present time. On the other hand, the No SN model produces much more stellar mass at early times, but once the critical halo mass is reached, star formation slows down, and the halo reaches a similar final stellar mass as the fiducial model. The $M_{\mathrm{crit}}=10^{10}\mathrm{M}_{\odot}$ model presents a similar behaviour to the fiducial model, i.e. once the halo reaches the critical mass it hardly produces any additional stellar mass. However, as the critical mass is lower for this model, the transition happens at earlier times. The high efficiency model has the same shape as the fiducial model, but as expected, a higher normalisation.

3.2 The stellar mass function

We use Eqs. 2, 5 and 13 to calculate the GSMF at redshift $z{=}0$ for the six efficiency models. We note that for a time evolving efficiency, Eq. 16 should be used to calculate the stellar mass of any halo as a function of time. These equations allow us to obtain the SHMR and convolve it with the halo mass function to calculate the GSMF.

Fig. 8 shows the GSMF at the present time for the six efficiency models. For the constant model, as it has been pointed out before (e.g Benson et al., 2003), if feedback process are inefficient both at the low mass and high mass end, i.e. a constant fraction of the baryon budget is turned into stars in every halo, the GSMF does not exhibit the characteristic knee obtained in observations and is identical to the halo mass function (dotted grey line) but shifted in mass by a constant value. Once feedback processes are implemented, the location of the knee of the GSMF is determined by the critical mass in the star formation efficiency (fiducial, $M_{\mathrm{crit}}=10^{10}$). Changing the critical mass also changes the normalisation of the distribution. All models with $\beta=0.75$, i.e. at high masses cooling is suppressed by AGN feedback, exhibit a sharp cut-off at the transition mass. Hence, AGN feedback is mainly responsible for the knee of the GSMF. The No AGN model has the same shallow slope at the faint end of the GSM function as the fiducial model, with a slight bend at high masses driven only by the exponential cut-off of the halo mass function. The No SN model presents the same knee as the fiducial model, but the slope of the faint end of the GSM function is much steeper. As discussed in Section 2, the low-mass and high-mass end slopes of the SHMR produce a “bump” at the knee of the GSMF. Finally, the high efficiency model, as $\alpha$ and $\beta$ are the same as for the fiducial model, the shape of the GSMF is the same. i.e. the relative abundance of galaxies to their haloes, and hence the shape, is independent of the normalisation (as $\varepsilon$ does not depend on $\epsilon_{\mathrm{N}}$ in Eq. 7). For a given halo mass, changing the normalisation maps that halo mass to a different galaxy mass. Hence, a change in the normalisation, $\epsilon_{\mathrm{N}}$, shifts the whole distribution only horizontally. In this case, the high efficiency model shifts the GSMF to the right compared to the fiducial model.

3.3 The cosmic SFR density

The cosmic history of star formation is perhaps one of the most fundamental observables of our Universe. It has been observed to peak approximately 3.5 Gyr after the Big Bang ($z\approx 2$), and decline exponentially thereafter (for a review see Madau & Dickinson, 2014). Different groups have tried to model the complex physics driving the cosmic SFR by using, for example, full hydrodynamical simulations (e.g. Schaye et al., 2015; Davé et al., 2016; Dubois et al., 2016; Pillepich et al., 2018a). Our analytic model disentangles the role of cosmology from the role of astrophysics, which in turn, allows us to examine the effect of the different efficiency parameters on the cosmic SFR density.

We begin by noting that the behaviour of Eq. 22 is governed by two main factors. First, a multiplier term that originates from both, the halo accretion rate, and the halo mass function, and is $\propto t^{-7/3}$. This, comes from the dynamical timescale of the universe getting larger. Second, an exponential term contribution due to the build up of haloes in the halo mass functions that is $\propto e^{-t^{-4/3}}$. For a given halo mass then, the exponential term dominates at early times, and the contribution to the cosmic SFR density is driven by the exponential build up of haloes. At late times, the exponential term asymptotically tends to a constant value, and the further evolution of the cosmic SFR is dominated by the multiplier term, i.e., it behaves as a power law. As discussed in Salcido et al. (2018), the contribution of dark energy is only relevant at late times, and at its observed value, it has a negligible impact on the star formation in the Universe.

Figure 9 shows the integrated cosmic SFR density for the six efficiency models computed using Eq. 22. For the fiducial model, while smaller haloes are more abundant than large objects, a smaller fraction of the inflowing material is converted into stars. As a result, the SFR density is dominated by the largest haloes in which star formation is able to proceed without generating efficient feedback. The smaller haloes only contribute significantly at very early times, when the abundance of larger objects is strongly suppressed by the exponential term in the mass function. We see therefore that the contribution of haloes of mass $\approx M_{\mathrm{crit}}=10^{12}\,{\hbox{M}_{\odot}}$, is representative of most of the SFR in the model.

If star formation is efficient at all halo masses (constant model), then the cosmic SFR behaves like a power law with time, which only deviates from this behaviour at late times due to the suppression due to the cosmological constant.

Examining the No SN model reveals the origin of the peak in the cosmic history of star formation is the efficient feedback in low mass galaxies. Without a mechanism to suppress star formation in small haloes, the history of the cosmic SFR density would not have its characteristic peak. Supernovae feedback is then mainly responsible for shaping the cosmic SFR density of the Universe. On the other hand, examining the No AGN model reveals that efficient feedback in high mass haloes only has a moderate effect on shaping the cosmic star formation. Without a mechanism to prevent star formation in massive galaxies, the cosmic SFR density would still exhibit a peak, only changing mildly its amplitude and localisation. However, the slope of the decline would be similar (orange vs green dashed lines).

Changing the transition mass $M_{\mathrm{crit}}$ has a great impact on the localisation of the SFR peak. As in the fiducial model, the contribution of haloes of mass $\approx M_{\mathrm{crit}}$, is representative of most of the SFR in the $M_{\mathrm{crit}}=10^{10}$. Hence, the peak happens at earlier times, but also has a higher normalisation, as $10^{10}\mathrm{M}_{\odot}$ haloes are more abundant than $10^{12}\mathrm{M}_{\odot}$ haloes.

Finally for the high efficiency model, the shape of the SFR of the Universe is identical to the fiducial model, but with a higher normalisation.

4 Fitting observations

In this section we compare the galaxy formation outputs from our analytic model with different observational datasets. We begin by calibrating our model to reproduce the GSMF at ${z}\sim 0$ using observations from the Galaxy And Mass Assembly (GAMA) survey (Baldry et al., 2012) and the Sloan Digital Sky Survey (SDSS) (Moustakas et al., 2013)333In this paper we used the standardised GSMF data from Behroozi et al. (2019), which assumes a Chabrier (2003) IMF, a Bruzual & Charlot (2003) stellar population synthesis model, dust corrections from Calzetti et al. (2000), and UV-stellar mass corrections.. We use the reduced chi-squared statistic to derive the best-fitting effective star formation efficiency $\epsilon_{*}({M_{h}})$. Because the model is fully analytic, this calibration process is fast and easy to perform. Figure 10 shows the best fit model in orange. Results from the eagle reference simulation are shown in green for reference. The figure shows that a constant halo mass-dependent efficiency model provides an excellent fit to the present-day GSMF (with reduced $\chi_{\nu}^{2}$ = 1.5). The best best fit efficiency parameters are shown in Table 2.

4.1 Contrasting halo mass and virial temperature efficiency models

Having established the best-fit efficiency parameters for the model, we can study the evolution of the model outputs. By construction, $\epsilon_{*}({M_{h}})$ is only a function of halo mass and is fixed in time. Hence, the evolution of the GSFM depends only on the evolution of the abundance of haloes of mass $M_{h}$ as a function of time, as described by the halo mass function. Figure 11 shows the evolution of the predicted GSMF for the halo mass-dependent star formation efficiency model in dashed lines. Different panels and colours represent different redshifts. Observational data from Baldry et al. (2012); Moustakas et al. (2013); Tomczak et al. (2014); Ilbert et al. (2013); Muzzin et al. (2013); Song et al. (2016)444In this paper we used the standardised GSMF data from Behroozi et al. (2019), which assumes a Chabrier (2003) IMF, a Bruzual & Charlot (2003) stellar population synthesis model, dust corrections from Calzetti et al. (2000), and UV-stellar mass corrections. are shown with coloured symbols.

Remarkably, a simple halo mass-dependent efficiency model reproduces very well the evolution of the GSMF up to redshift $z\approx 4$. While the observed data at higher redshifts is highly uncertain, the halo mass-dependent model significantly under predicts the abundance of distant galaxies. This may hint to the need of a time-evolving efficiency model. As discussed in • ‣ Section 2.1, it is reasonable to assume that the efficiency of star formation should be a function of the halo’s virial temperature, which naturally evolves with cosmic time. In. this section, we investigate an effective star formation efficiency model that depends on virial temperate.

4.1.1 The three stages of galaxy formation

An effective star formation efficiency model characterised by a time-independent critical virial temperature $T_{\mathrm{crit}}$ assumes that there exists a critical halo virial temperature at which there is a transition from where star formation driven outflows can escape, to where outflows stall inside the halo. Using the viral temperature of the halo to parameterise this tipping point, provides a natural evolution of the star formation efficiency. For a fixed halo mass, early collapsed haloes are more compact (denser), and one might expect a higher efficiency (for haloes with $T_{\mathrm{vir}}<T_{\mathrm{crit}}$).

In this simple picture, we can distinguish three stages of galaxy formation555While perhaps closely related, we distinguish these three stages of galaxy formation from the three phases in Clauwens et al. (2018), as the latter refer mainly to a morphological evolution, rather than the entropy state and buoyancy of the gas., characterised by the virial temperature of the halo:

• •

Stellar feedback regulated stage: star formation driven outflows effectively regulate the gas content of galaxies residing in haloes with virial temperature $T_{\mathrm{vir}}<T_{\mathrm{crit}}$. In this stage, efficient outflows prevent the density of central star forming gas building up, suppressing the growth of the central BH.

• •

Efficient star forming/rapid growing black hole stage: as haloes grow, the virial temperature increases to the point that the stellar outflows are no longer buoyant relative to their surroundings, and therefore stall (i.e $T_{\mathrm{vir}}\approx T_{\mathrm{crit}}$). The density of gas builds up within the halo triggering high star formation rates and rapid BH growth.

• •

Black hole feedback regulated stage: In haloes with $T_{\mathrm{vir}}>T_{\mathrm{crit}}$, the central BH is massive enough to produce efficient AGN feedback, in turn, regulating the gas content of the halo and preventing further star formation.

An additional advantage of using the virial temperature to characterise the star formation efficiency, is that we can add a proxy for the effect of cosmic reionisation. Ultraviolet radiation from the first stars formed reionised neutral hydrogen, raising its entropy to a temperature of $\approx 10^{4}\mathrm{K}$. This process prevented further cooling, hence preventing star formation in haloes with $T_{\mathrm{vir}}<10^{4}\mathrm{K}$ (Doroshkevich et al., 1967; Couchman & Rees, 1986; Rees, 1986; Efstathiou, 1992; Loeb & Barkana, 2001). As a result of this suppression of star formation, only a fraction of the haloes with present-day mass $\approx 10^{10}\mathrm{M}_{\odot}$ form a galaxy, and no galaxies form below a halo mass of $\approx 10^{7}\mathrm{M}_{\odot}$ (Sawala et al., 2013; Sawala et al., 2016; Fitts et al., 2017; Bose et al., 2018). We therefore include the effect of reionisation by setting $\epsilon_{*}(T_{\mathrm{vir}}<10^{4}\mathrm{K})=0$.

Of course, one could think of more complex ways in which the expected star formation efficiency might evolve with cosmic time. For instance, the evolution of cooling versus free-fall time of a cloud of gas, the evolution of metallicity and the UV background radiation, might all result in a more complex evolution. However, the aim here is to describe the main features of the universe as simply as possible, and so we leave exploration of more complex models for future work.

We calibrate the $\epsilon_{*}(T_{\mathrm{vir}})$ model to the GSMF at ${z}=0$ using the reduced chi-squared statistic to derive the best-fitting parameters. The best fit efficiency parameters are shown in Table 2.

It is important to highlight that the models were calibrated to reproduce only the observed GSMF at redshift $z\sim 0$. Figure 10 shows the best fit virial temperature model in blue. The figure shows that both the halo mass-dependent efficiency $\epsilon_{*}$, and the virial temperature-dependent efficiency $\epsilon_{*}(T_{\mathrm{vir}})$ models, provide a good fit to the present-day GSMF, both at the faint end and at the knee. Figure 11 shows that a star formation efficiency as a function of the virial temperature of the halo provides a good fit to the abundance of galaxies both at low and high redshift, but the evolution is too rapid at intermediate redshift ($z{=}1$ to $z{=}4$).

Figure 12 shows the build up of the stellar mass within haloes using both models (calculated integrating Eq. 16). Colour coding represents different present-day halo mass $M_{h}(t_{0})$. The transition of the star formation efficiency at $M_{\mathrm{crit}}$ can be clearly seen in very massive haloes, where there is a rapid rise of stellar mass, then, when the halo reaches the critical mass (or virial temperature), the build up of stellar mass slows down significantly. The change in slope is due to AGN feedback becoming efficient in those haloes, preventing any further star formation.

Figure 13 shows the predicted cosmic SFR for the two efficiency models. Using the analytic model, we can clearly see the contribution to the integrated SFR density from dark matter haloes of different masses (per dex) shown as coloured dashed lines (only shown for the halo mass-dependent efficiency model). The total SFR for the virial temperature efficiency model is shown in blue. The halo mass-dependent efficiency model is shown in orange. Results from the eagle simulation are shown in green for reference. Observational data compiled by Behroozi et al. (2013) are shown as grey symbols. The latest observational results from the GAMA survey from Driver et al. (2018) are shown as black symbols. The model using a halo mass efficiency reproduces the amplitude and shape of the observed SFR density remarkably well, while the virial temperature-dependent efficiency model, produces a higher SFR at high redshift.

Figure 14 shows the predicted SHMR from both efficiency models. Colour coding represents different redshifts. The virial temperature efficiency model is shown in solid lines. The halo mass-dependent efficiency model is shown with a dashed line (only shown for ${z}=0$ as the halo mass-dependent efficiency model is constant in time). The critical halo mass predicted in Bower et al. (2017) is shown in vertical dotted lines, which roughly coincide with the peak efficiency for a viral temperature efficiency model. Recently, McAlpine et al. (2018) showed that the critical halo mass predicted in Bower et al. (2017) agrees remarkably well with the triggering of a rapid black hole growth stage in the eagle simulations.

The model using a virial temperature efficiency predicts a SHMR relation that differs from observational contains using abundance and clustering properties of galaxy samples with predictions from a phenomenological halo models. For example, recently Cowley et al. (2018) calculated that the peak of the SHMR shifts to higher masses at earlier times. These methods however, depend heavily on the underlying modelling and assumptions. More sophisticated empirical models (e.g. Behroozi et al., 2019; Moster et al., 2018) find that the peak in the SHMR moves first to higher masses for low redshifts, and then to lower masses at high redshifts.

Finally, in Fig. 15 we show the sSFR of galaxies for different redshifts. The halo mass-dependent model is shown in dashed lines. The model using a virial temperature efficiency is shown in solid lines. Results for central galaxies from the eagle simulations are shown in dotted lines for reference. Observational data from Gilbank et al. (2010); Karim et al. (2011); Tomczak et al. (2016) are shown as symbols. While not calibrated to reproduce the sSFR of galaxies, the agreement of the halo mass-dependent model with the observational data is remarkable.

5 Discussion and Conclusions

In our current paradigm of galaxy formation, every galaxy forms within a dark matter halo. Due to the tight correlation observed between the properties of galaxies and their host haloes, it is natural to expect that individual galaxy assembly could be correlated with halo assembly (see Wechsler & Tinker 2018 for a review).

In this paper we developed a fully analytic model of galaxy formation that connects the growth of dark matter haloes in a cosmological background, with the build up of stellar mass within these haloes. The model restricts the role of baryonic astrophysics to setting the relation between galaxies and their haloes. We assume an effective star formation efficiency which captures all the physical processes involved in the conversions of gas into stars, i.e. cooling, star formation law, feedback mechanisms, etc.

We show that galaxy formation is revealed as a simple process where the effective star formation efficiency within haloes is only a function of their mass. We show that all the complex physics of galaxy formation, the interplay between cosmology and baryonic process can be understood as a simple set of equations. Despite its simplicity, the model reproduces self-consistently the shape and evolution of the cosmic star formation rate density, the specific star formation rate of galaxies, and the galaxy stellar mass function, both at the present time and at high redshift.

We use our model to investigate the origin of the characteristic shape of the GSMF and the need for a double Schechter function to describe it. Using the logarithmic slope of the SHMR, the model naturally explains an inflection point in the distribution causing the characteristic “bump” observed at the knee of the GSMF.

To demonstrate the flexibility and power of our mathematical framework, we introduced a physically motived model for the effective star formation efficiency, characterised by a time-independent critical virial temperature, $T_{\mathrm{crit}}$. The model assumes that there exists a critical halo virial temperature at which there is a transition from where star formation driven outflows can escape, to where outflows stall, triggering high star formation rates and rapid BH growth. We demonstrate that this model can reproduce the GSMF at high redshift ($z>4$) better than a simple halo mass-dependent model, but the evolution at intermediate redshifts is to rapid to reproduce observations.

While the aim of this paper is not to present a “perfect” model fitted to reproduce a large set of observational constraints, the two variations of an effective star formation efficiency presented here, already provide very valuable information about the average evolution of the galaxy population within a cosmological background. Furthermore, the model can be easily extended to include further modelling (such as time evolution of the model parameters, or a prescription for satellite galaxies, e.g. Grylls et al. 2019) or the use of advanced gradient-based minimization and Markov Chain Monte Carlo algorithms to fit to a larger number of datasets. Additionally, the model can be easily adapted to combine the equations developed here, with for example, halo merger trees from a dark matter simulation.

Our model is limited to the connection between haloes and central galaxies only. Sub-haloes and satellite galaxies are subject to complex processes, such as tidal and ram pressure stripping, which are not included.

Finally, one of the main advantages of the model is that by providing a set of analytic equations, the model can be easily “inverted” and allows for rapid experiments to be conducted, providing a great tool to explore the differential effects of baryonic physics, averaged over galaxy scales. We conclude therefore that there is a clear opportunity to use the analytic model developed in this paper to improve theoretical galaxy formation models.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 769130). This work was supported by the Science and Technology Facilities Council [ST/P000541/1]. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1 and ST/R002371/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Appendix A Derivation of the Model

A.1 Cosmological expansion

Here we provide a brief summary of the analytic solution of the Friedmann equation developed in Salcido et al. (2018). Equations with full cosmology dependence of the numerical constants are highlighted using a coloured superscript (^{{\color[rgb]{1,0,0}c^{*}}}), and can be used for arbitrary flat $\Lambda$CDM cosmologies.

Using a Taylor expansion, the expansion factor of the Universe can be written as,

[TABLE]

where the matter timescale is given by,

[TABLE]

and the dark energy timescale is given by,

[TABLE]

For the cosmological parameters given at the end of Section 1, $t_{\mathrm{m}}=26.04\,\hbox{Gyr}$ and $t_{\Lambda}=17.33\,\hbox{Gyr}$. At the present day, $t\equiv t_{0}=13.82\,\hbox{Gyr}$, so that $\frac{t_{0}}{t_{m}}=0.53$ and $\frac{t_{0}}{t_{\Lambda}}=0.8$. By convention, Eq. 23 is normalised so that $a(t_{0})=1$.

A.2 The growth of density perturbations and the halo accretion rates

Dark matter structures are assumed to have grown from small initial density perturbations. Expressing the density, $\rho$, in terms of the density perturbation contrast against a density background,

[TABLE]

the differential equation that governs the time dependence of the growth of linear perturbations in a pressureless fluid, such as e.g. dark matter, can be written as

[TABLE]

The growing mode of Eq. 27 can be written as,

[TABLE]

where $D(t)$ is the linear growth factor, which determines the normalisation of the linear matter power spectrum relative to the initial density perturbation power spectrum, and is computed by the integral

[TABLE]

Using the power-series approximation for $a(t)$ from Eq. 23, keeping the leading order terms and using the definition of $f_{\Lambda}$ in Eq. 11, we can obtain an analytic solution of Eq. 29,

[TABLE]

where $K_{D}$ is a normalisation constant with units of time*-2*. By convention, $K_{D}$ is chosen so that ${D}(t_{0})=1$. For the cosmological parameters inferred by the Planck Collaboration et al. (2014), $K_{D}=4.7\times 10^{-3}\,\mathrm{Gyr}^{-2}$. Collecting the numerical and cosmology dependent constants together gives,

[TABLE]

The growth rates of linear perturbations do not directly predict the growth rates of haloes; however, we can directly connect the two through the approach developed by Press & Schechter (1974). Correa et al. (2015) showed that the accretion rates of haloes can be written as (see also Neistein et al., 2006),

[TABLE]

where $M_{h}$ is the halo mass and $S$ is the variance of the density field on the length scale corresponding the halo mass. $\delta_{c}$ is a parameter that represents a threshold in the linearly extrapolated density field for halo collapse. We assume $\delta_{c}=1.68$ (Press & Schechter, 1974). The parameters, $q$ and $\gamma$, are related to the shape of the power-spectrum around the halo mass $M_{h}$. The scale dependence of the density field is approximated as a power-law around $10^{12}\,{\hbox{M}_{\odot}}$ haloes as $S=S_{0}(M_{h}/10^{12}\mathrm{M}_{\odot})^{-\gamma}$. Correa et al. (2015) find that this prescription works for different cosmologies because the halo mass histories are mainly driven by changes in $\sigma_{8}$ and $\Omega_{m}$. For the cosmological parameters inferred by the Planck Collaboration et al. (2014), $S_{0}\approx 3.98,\gamma\approx 0.3$ and $q\approx 3.16$. Collecting the numerical and cosmology dependent constants together yields,

[TABLE]

Using the series approximation Eq. 31, the specific growth rate of haloes can be written as,

[TABLE]

This differential equation can be solved by separation of variables to obtain the average mass history of dark matter haloes,

[TABLE]

where $M_{0}$ is the mass of a halo today. Integrating both sides and solving for $M(t)$ yields,

[TABLE]

Note that in the case $\gamma\rightarrow 0$ the LHS becomes the logarithm of the mass ratio $M_{h}/M_{0}$, and all haloes grow by the same factor in a given time interval. For realistic power spectra, however, the relative growth rate increases with mass because massive haloes arise from increasingly rare fluctuations in the initial density perturbation field. We can re-write Equation 36 as an explicit equation for the halo mass as a function of time. This form is useful for symbolic substitution into calculations that are driven by the halo mass.

[TABLE]

As $t\rightarrow 0$, the mass of the halo becomes small compared to the final mass so that we can write,

[TABLE]

where the first term in the RHS is much smaller than the second term. This shows that masses of early haloes depend very weakly on their average final mass, and that the halo mass initially grows roughly $\propto t^{4}$, (since $\gamma\approx 1/3$).

Finally, collecting the numerical and cosmology dependent constants together, we can write Eq. 37 as,

[TABLE]

A.3 The halo mass function

In the Press & Schechter analysis, the co-moving abundance of haloes of mass $M_{h}$ at time $t$ is given by (Press & Schechter, 1974),

[TABLE]

where we have assumed that the density power spectrum is a power law with exponent $\gamma$ and written the co-moving matter density of the Universe as ${\rho}_{0}$ following our convention. Using the evolution of the growth factor given by Eq. 31 and keeping the leading order terms we obtain,

[TABLE]

For the cosmological parameters adopted in this paper, $\rho_{0}=3.913\times 10^{10}\,{\hbox{M}_{\odot}}\,\hbox{Mpc}^{-3}$. Substituting for values of the constants and cosmological parameters, we can write Eq. 41 as,

[TABLE]

Appendix B Comparison with hydrodynamical simulations

In this section, we compare three variations of the eagle $(50\mathrm{cMpc})^{3}$ simulations to their equivalent analytic effective star formation efficiency model. The Ref-L050N0752 eagle model (Schaye et al., 2015; Crain et al., 2015; McAlpine et al., 2016), uses the same calibrated sub-grid parameters as the reference model $(100\mathrm{cMpc})^{3}$, ran with the same resolution, but in a smaller volume. The “No AGN” run uses the same calibrated sub-grid parameters as the reference model but removing feedback from BHs. For the “No SN” model (red, introduced in Bower et al., 2017), feedback from star formation has been removed. We note that, while the eagle “No SN” simulation removes the effect of star formation feedback, it still includes the effect of cosmic reionisation. Hence, there is a suppression of star formation is small haloes. In order to compare with the simulations, we have included the effect of cosmic reionisation, $\epsilon_{*}(T_{\mathrm{vir}}<10^{4}\mathrm{K})=0$ in both the “No SN” and “constant” star formation efficiency models. Finally, there is no eagle equivalent to the “constant” (or “no feedback”) model. In Fig. B16 we compare the GSMF at $z=0.1$, and in Fig. B17, we compare the SFR history of the Universe.

While much more computationally expensive, the behaviour of the full hydrodynamical simulations is well approximated by the analytic models introduced here.