The $I\kappa\epsilon\alpha$ model of feedback-regulated galaxy formation

TL;DR

The $Ipsilonlpha$ model explains galaxy formation through a feedback-regulated balance between stellar feedback energy and dark matter halo accretion, accurately matching simulation data and predicting galaxy properties over cosmic time.

Contribution

The paper introduces the $Ipsilonlpha$ model, a new framework that captures galaxy star formation regulation via feedback and accretion balance, independent of gas content or star formation law.

Findings

Accurately reproduces star formation rates in EAGLE simulations.

Predicts a star forming main sequence with weak mass dependence and redshift evolution.

Provides insights into mass-metallicity relation and galaxy mass function evolution.

Abstract

We present the $I\kappa\epsilon\alpha$ model of galaxy formation, in which a galaxy's star formation rate is set by the balance between energy injected by feedback from massive stars and energy lost by the deepening of the potential of its host dark matter halo due to cosmological accretion. Such a balance is secularly stable provided that the star formation rate increases with the pressure in the star forming gas. The $I\kappa\epsilon\alpha$ model has four parameters that together control the feedback from star formation and the cosmological accretion rate onto a halo. $I\kappa\epsilon\alpha$ reproduces accurately the star formation rate as a function of halo mass and redshift in the EAGLE hydrodynamical simulation, even when all four parameters are held constant. It predicts the emergence of a star forming main sequence along which the specific star formation rate depends weakly on stellar mass with an amplitude that increases rapidly with redshift. We briefly discuss the emerging mass-metallicity relation, the evolution of the galaxy stellar mass function, and an extension of the model that includes feedback from active galactic nuclei (AGN). These self-regulation results are independent of the star formation law and the galaxy's gas content. Instead, star forming galaxies are shaped by the balance between stellar feedback and cosmological accretion, with accurately accounting for energy losses associated with feedback a crucial ingredient.

Tables

Table 1. Table 1: Selected parameters of the eagle  simulations used here. From left-to-right the columns show: simulation name; co-moving box size; initial baryonic particle mass; maximum proper softening length, and comment.

Name	$L$	$m_g$	${ϵ}_{\mathrm{prop}}$	comment
	(Mpc)	(${10}^6{M}_{\odot }$)	(kpc)
REF	50	$1.81$	0.7	reference model
FBconst	50	$1.81$	0.7	$f_{\mathrm{th}}=1$
FBconstnoAGN	50	$1.81$	0.7	$f_{\mathrm{th}}=1$, no AGN
DMO	50	0	0.7	dark matter only

Equations

$$\dot{E}={\dot{E}}_{\mathrm{nucl}}-L.$$

$$E_{\mathrm{h}}$$

$${\Omega }_{\mathrm{h}}$$

$$U_{\mathrm{h}}$$

$$R_{\mathrm{h}}$$

$${v_{\mathrm{h}}}^2=\alpha (\mathrm{G}{M}_{\mathrm{h}}{)}^{2/3}(10H{)}^{2/3}.$$

$$\frac{d\ln |E_{\mathrm{h}}|}{dz}=\frac{5}{3}\frac{d\ln M_{\mathrm{h}}}{dz}+\frac{d\ln \alpha }{dz}+\frac{2}{3}\frac{d\ln H}{dz}.$$

$${\dot{E}}_{\mathrm{h}}=-\frac{\kappa }{2}{\dot{M}}_{\mathrm{h}}{v_{\mathrm{h}}}^2.$$

$$M_{\mathrm{h}}$$

$$m_{\mathrm{h}}(z)$$

$$\frac{d\ln M_{\mathrm{h}}}{dz}$$

$${\xi }_{\mathrm{h}}(z)$$

$$\bar{a}\approx 0.24,\bar{b}\approx 0.75.$$

$${v_{\mathrm{h}}}^2(z)=v_{\mathrm{h},0}^2{\left(m_{\mathrm{h}}(z)\mathcal{H}(z)\right)}^{2/3},$$

$$\mathcal{H}(z)\equiv \frac{H(z)}{H_0}.$$

$${\omega }_{\mathrm{b}}\equiv \frac{{\Omega }_{\mathrm{b}}}{{\Omega }_{\mathrm{dm}}+{\Omega }_{\mathrm{b}}}=\frac{{\Omega }_{\mathrm{b}}}{{\Omega }_{\mathrm{m}}}.$$

$${\dot{E}}_{\mathrm{g}}={\dot{E}}_{\star }-\frac{\kappa }{2}{\omega }_{\mathrm{b}}{\dot{M}}_{\mathrm{h}}{v_{\mathrm{h}}}^2.$$

$${\dot{E}}_{\star }=\frac{1}{2}{\dot{M}}_{\star }{v}_{\star }^2.$$

$$v_{\star }={\left(ϵ\eta \frac{2\times 10^{51}\mathrm{erg}}{100{\mathrm{M}}_{\odot }}\right)}^{1/2}\approx 400{\left(\frac{ϵ}{0.091}\times \frac{\eta }{1.74}\right)}^{1/2}{\mathrm{km}\mathrm{s}}^{-1}.$$

$$\frac{1}{2}{\dot{M}}_{\star }{v}_{\star }^2=\frac{\kappa }{2}{\omega }_{\mathrm{b}}{\dot{M}}_{\mathrm{h}}{v_{\mathrm{h}}}^2,$$

$${\dot{M}}_{\star }(z)$$

$${\dot{M}}_{\star ,0}$$

$${\Psi }_{\star }(z)$$

$$M_{\star }(t)=(1-\mathcal{R}){\int }_0^t{\dot{M}}_{\star }(t^{\prime })d{t}^{\prime },$$

$$M_{\star }(z)$$

$$M_{\star ,0}$$

$$m_{\star }(z)$$

$$m_{\star ,0}$$

$$\mathrm{sSFR}(z)$$

$$\frac{d\ln M_{\star }}{d\ln M_{\mathrm{h}}}{|}_{z=\mathrm{const}}=\frac{5}{3}.$$

$$\frac{d\ln M_{\star }}{d\ln M_{\mathrm{h}}}{|}_{M_{\mathrm{h},0}=\mathrm{const}}=-\frac{m_{\mathrm{h}}(z)}{m_{\star ,0}(a-b)m_{\star }(z)}(m_{\mathrm{h}}(z)\mathcal{H}(z){)}^{2/3}.$$

$$\frac{d{n}_{\mathrm{h}}}{d\log M_{\mathrm{h},0}}=n_0(\frac{M_{\mathrm{h},0}}{M_{\mathrm{ps}}}{)}^{-{\alpha }_{\mathrm{h}}}\exp (-\frac{M_{\mathrm{h},0}}{M_{\mathrm{ps}}}),$$

$$\frac{d{n}_{\mathrm{h}}}{d\log M_{\mathrm{h}}}=n_0(\frac{M_{\mathrm{h}}}{m_{\mathrm{h}}(z)M_{\mathrm{ps}}}{)}^{-{\alpha }_{\mathrm{h}}}\exp (-\frac{M_{\mathrm{h}}}{m_{\mathrm{h}}(z)M_{\mathrm{ps}}}),$$

$$\frac{d{n}_{\mathrm{g}}}{d\log M_{\star }}$$

$$M_{\star ,\mathrm{ps}}$$

$${\alpha }_{\star }$$

Full text

1 Introduction

The cold dark matter cosmogony links the small fluctuations detected in the cosmic microwave background (cmb) at redshift $z\sim 1000$ to the observed large-scale clustering of galaxies at all observable redshifts. The fluctuations in the cmb temperature correspond to density perturbations that grow in amplitude due to gravity, resulting in the formation of dark matter halos that host galaxies (see e.g. Springel et al. 2005 and reference therein for more background).

Whereas computer simulations can reliably predict virtually all properties of dark halos, the same can not be said for the properties of the galaxies that inhabit these halos. Even though our basic understanding of the underlying physics is probably correct - galaxies form as gas accretes onto a halo, cools, becomes self-gravitating and forms stars (White & Rees, 1978; White & Frenk, 1991) - numerous uncertainties remain. What sets the star formation rate of a galaxy in a given halo at a given redshift? How does the energetic feedback from stars and accreting black holes regulate star formation? What is the role of galaxy interactions such as mergers? Are there any other crucial processes, for example feedback from cosmic rays or reionisation, and what is the role of magnetic fields?

Models that are designed to reproduce a mock universe that looks and evolves like the one we observe may not care about the details of the relevant physical processes. Examples include halo occupation distribution models (hod, e.g. Peacock & Smith 2000) or subhalo abundance matching (sham, e.g. Vale & Ostriker (2004), see e.g. Wechsler & Tinker (2018) for recent reviews).

Semi-analytical models recognise that the physics of galaxy formation is complex, and use parametrizations to model poorly understood physical processes. Cosmological hydrodynamical simulations try to capture some of these physical processes as accurately as possible (cosmological accretion and cooling of gas onto halos for example), but also rely on more parametrised descriptions of physical processes to capture physics below the resolution scale (see Somerville & Davé 2015 and Naab & Ostriker 2017 for recent reviews).

Several of the semi-analytical models and recent hydrodynamical simulations yield mock universe that look impressively similar to the one observed. Even though these models typically all include the same ingredients, the details of how the processes are implemented may be quite different. It is therefore somewhat surprising that the resulting galaxy population is nevertheless very similar. At the very least this suggests some level of degeneracy in the modelling and that such calculations cannot be used to infer how the unresolved processes operate in detail. But it also suggests that many properties of galaxies do not actually depend on the details of many of these processes (see Hopkins et al. 2014 for a similar point of view).

Arguably one of the more striking features of the galaxy population as a whole is the emergence of a ‘star forming main sequence’ (or ‘blue cloud’), Noeske et al. (2007), on which galaxies form stars at a specific rate, $\dot{M}_{\star}/M_{\star}$, that depends weakly on stellar mass ($M_{\star}$), but increases rapidly with redshift. The scatter around the mean trend is small, of order 0.3 dex (see Schreiber et al. 2015 and references therein for more recent observational analysis and discussion).

The appearance of such a main sequence suggests that the rate at which a galaxy forms stars in a halo of given mass, is somehow self-regulating. Several papers argued just that (e.g. Bouché et al., 2010; Lilly et al., 2013; Dekel & Mandelker, 2014; Davé et al., 2012; Dayal et al., 2013). The aim of these models is not to be able to predict the properties of galaxies in great detail, but rather understand the origin of self-regulation. The current paper follows this philosophy, adopting simplifications to more clearly expose the feedback loop that operates on the star forming sequence.

This paper is organised as follows: section 2 exposes the basic physics behind self-regulation in our model and tests the central assumptions by comparing to galaxies from the eagle cosmological hydrodynamical simulation (Schaye et al., 2015). Section 3 explores consequences in terms of galaxy scaling relations (such as the galaxy stellar mass function and the mass-metallicity relation), compares these to simulations and data, and discusses successes and failures of the model. Section 4 puts our results into context by comparing to previous work, and discusses what we think are its main limitations. Section 5 summarises our findings and is followed by an appendix that contains a short overview of the eagle simulations, including a description of the reference model, ‘Ref-L100N1504’, in which the subgrid parameters are calibrated to reproduce redshift $z=0$ observations of the galaxy stellar mass function, the relation between galaxy size and mass, the relation between black hole mass and stellar mass, as described by Crain et al. (2015). The appendix also describes the eagle model ‘FbConstNoAGN’, in which the feedback parameters are kept constant and which does not include feedback from AGN, as well as another eagle variation, ‘FbConst’, in which the feedback parameters are kept constant and which does include AGN feedback.

2 Self-regulation of star formation in galaxies

The appearance of a star-forming sequence of galaxies is suggestive of the action of a feedback cycle. Such a feedback cycle is also important in understanding the main sequence of stars in a Hertzsprung-Russel diagram. Indeed: nuclear energy generation in main sequence stars is secularly stable - a prerequisite for their longevity. We begin this section by briefly describing the well-known reason behind this stability (see e.g. any text book on stellar structure, for example Prialnik 2009). We next investigate whether we can apply similar reasoning to star forming galaxies.

2.1 The secular evolution of main sequence stars

The total energy $E$ of a main sequence star of mass $M$ is the sum of its gravitational energy, $\Omega<0$, and its internal energy, $U=Mu$, where $u$ is its mean specific energy per unit mass. Stars are approximately in virial equilibrium, $E=\Omega/2=-Mu$, and as a consequence $dE/du<0$. Therefore, if a star loses energy for example through radiation so that $\dot{E}<0$, it will heat up, $\dot{u}>0$. The effective negative specific heat capacity of a star is a well-known but nevertheless intriguing feature of gravitationally bound systems, see e.g. Lynden-Bell & Lynden-Bell (1977), and is crucial for its longevity.

Indeed, consider a star losing energy through radiation (rate $L$), while gaining internal energy through nuclear fusion (rate $\dot{E}_{\rm nucl}$),

[TABLE]

In equilibrium, $\dot{E}=0$, however consider what happens for (small) deviations from equilibrium. Assuming $\dot{E}_{\rm nucl}<L$, say, $|E|$ increases since $E<0$, meaning $|u|$ increases and hence the temperature $T$ rises. The rate of energy generation through fusion is a rapidly increasing function of $T$, hence increasing $T$ increases $\dot{E}_{\rm nucl}$, so that $\dot{E}_{\rm nucl}<L$ results in an increase in $\dot{E}_{\rm nucl}$ towards equilibrium. Similarly, if $\dot{E}_{\rm nucl}>L$, the decrease in $T$ results in a decrease in the nuclear burning, until $\dot{E}_{\rm nucl}=L$. Clearly, the negative specific heat capacity of a star is not just an amusing feature of self-gravitating systems, but is key in understanding stability on the main sequence. As the star’s mean molecular weight changes due to fusion, $L$ and hence $\dot{E}_{\rm nucl}$ evolve secularly on a time scale which vastly exceeds $E/L$.

2.2 The evolution of a galactic halo

As a galactic halo111We will use the term ‘galactic halo’ to refer to a central galaxy (as opposed to a satellite galaxy) with gas and stars, together with its host dark matter halo. grows in mass due to cosmological accretion, its energy changes in time as well. At first sight there is little in common between the evolution of a galactic halo and that of a main sequence star. Indeed, the total energy of a star changes only secularly, $|\dot{E}|\ll L$, as self-regulation leads to a near balance between the energy generated by nuclear fusion and lost by radiation, but a galactic halo seems to have no equivalent channel for regulation. Does that mean that it is not secularly stable? The answer is partially yes: we show in the following that the dark matter halo is not secularly stable, in the sense that $\dot{E}_{\rm h}\neq 0$. However the same may not be true for the galaxy itself, because supernovae inject energy into the interstellar medium. Below we investigate whether that energy injection rate balances the loss of energy due to cosmological accretion, and if such a situation is a stable equilibrium - in analogy with the evolution of main sequence stars describe above. Before we do so we summarise some well known relations for the evolution of dark matter halos.

2.3 The growth of a dark matter halo

We begin by investigating the cosmological growth in mass and the associated change in energy of a dark matter halo. The concentration and mass of a dark matter halo may be affected by baryonic processes. Indeed, in the simulations presented by Duffy et al. 2010, strong cooling and inefficient feedback increases the central dark matter density of galaxy and group halos significantly, whereas strong feedback, for example from an AGN, decreases that density. Baryonic mass loss, associated with strong feedback, may lead to a decrease in the rate at which a dark matter halo increases its mass. These effects are relatively modest at the scale of galaxies in the eagle simulations, as shown by Schaller et al. 2015, and we will neglect them in this paper.

The total energy, $E_{\rm h}<0$, of a dark matter halo with mass $M_{\rm h}$ is the sum of its potential energy, $\Omega_{\rm h}<0$, and its internal energy, $U_{\rm h}$ (the total kinetic energy of all dark matter particles in the centre of mass rest frame, subscript $h$ for halo). Dark matter halos satisfy the virial theorem approximately, $E_{\rm h}\approx{\Omega_{\rm h}/2}\approx-U_{\rm h}$ (e.g. Neto et al., 2007), as we show in Fig. 1. There is clearly some evolution of the ratio $U_{\rm h}/\Omega_{\rm h}$ as the halo grows, but we will neglect this in what follows.

Assuming that the dark matter halo is in virial equilibrium, mass, radius and internal energy are related by,

[TABLE]

We used the standard way of assigning a ‘radius’, $R_{\rm h}$, to a halo, by requiring that the mean density within $R_{\rm h}$ is 200 times the critical density, $\rho_{c}=3H^{2}/(8\pi{\hbox{\rm G}})$, where $H(z)$ is the Hubble constant at redshift $z$. The value of the dimensionless parameter $\alpha$ depends on the halo’s density profile: $\alpha=3/5$ for constant density, \alpha={\hbox{R_{\rm h}}}/(6a) for the spherical profile with scale radius $a$ described by Hernquist (1990), and $\alpha$ is uniquely related to the concentration parameter, $c$, of a halo with an NFW (Navarro et al., 1997) profile. Equations (2.3) also define a characteristic ‘virial velocity ’ of the halo, $v_{\rm h}$, also given by

[TABLE]

If the accreting halo remains in virial equilibrium, then

[TABLE]

We will show below that the first term on the right hand side, |{5\over 3}{d\ln{\hbox{M_{\rm h}}}/dz}| is of order unity. How about the other terms? We tracked the evolution of the parameter $\alpha$ of halos in the eagle L0100N1504 dark matter only simulation along their merger tree. The result is plotted in Fig. 2, where different colours refer to halos in bins of their redshift $z=0$ mass, $M_{\rm h,0}$. As was the case of the virial ratio, there is clearly some evolution in $\alpha$ as a halo grows, but that evolution is relatively weak and we will neglect it. We also note that the term $(2/3)d\ln H/dz$ is always $<1/2$. Therefore the last two terms in Eq. (4) are small compared to the first term on the right hand side, therefore {d\ln|E_{\rm h}|/dz}\approx\kappa{d\ln{\hbox{M_{\rm h}}}/dz} with $\kappa\approx 5/3$. To test this approximation in more detail, we once more track halos along their merger tree to compute d\ln|{\hbox{E_{\rm h}}}|/d\ln{\hbox{M_{\rm h}}} directly, the result is plotted in Fig. 3; different colours refer to halos in bins of $M_{\rm h,0}$. As $M_{\rm h}$ increases, $|E_{\rm h}|$ increases, with d\ln|{\hbox{E_{\rm h}}}|/d\ln M_{\rm h}\approx 5/3. Combining this approximation with Eq. (3), motivates us to parametrize the rate of change of energy as a halo grows in mass by

[TABLE]

The variables $\alpha$ (Eq. 3) and $\kappa$ are two of the four parameters of the $I\kappa\epsilon\alpha$ model - and as we just showed, they are approximately independent of halo mass and redshift, and we will simply keep them constant at $\alpha=1$ and $\kappa=5/3$. We proceed by parametrizing the evolution of $M_{\rm h}$.

The increase with time of the halo mass in the extended Press-Schecher (EPS) or ‘excursion set’ formalism (Bond et al., 1991; Lacey & Cole, 1993) describes the growth measured in simulations very well. Here we will use the parametrisation described by Correa et al. (2015a, b), which we write in the form of the value of the halo mass at redshift $z=0$, $M_{{\rm h},0}$, times a dimensionless function $m_{\rm h}(z)$,

[TABLE]

The corresponding logarithmic growth rate is

[TABLE]

The dimensionless functions $m_{\rm h}$ and $\xi_{\rm h}(z)$ are both unity at $z=0$. Since halos grow in mass, {d\ln{\hbox{M_{\rm h}}}/dz}\leq 0, and in terms of the previous equation we have that the function $\xi_{\rm h}(z)>0$ but $a-b<0$. The parameters $a$ and $b$ depend on the mass of the halo at some reference redshift which we take to be $z=0$. Averaging over halo masses, Correa et al. (2015b) find

[TABLE]

We will use $\bar{a}$ and $\bar{b}$ and denote them by $a$ and $b$ in our derivations, but in our figures we will use the more elaborate but more accurate version discussed by Correa et al. (2015b) in which $a$ and $b$ are functions of $M_{\rm h,0}$ (except in Figure 5 and 15 in which we use the constant values from Eq. 8). Using $\bar{a}$ and $\bar{b}$, gives |(5/3)d\ln{\hbox{M_{{\rm h},0}}}/dz|=0.85 at $z=0$ and $1.25$ for $z\rightarrow\infty$, therefore |(5/3)d\ln{\hbox{M_{\rm h}}}/dz| is of order unity, as we used before.

The virial velocity’s evolution follows from Eq. (3),

[TABLE]

where the function $\mathcal{H}(z)$ is defined by

[TABLE]

2.4 The growth of a galaxy

2.4.1 Stability of feedback-regulated galaxy formation

A (central) galaxy too satisfies the equivalent of Eq. (5). We neglect any pre-processing of the accreted matter, so that the ratio of gas mass that accretes onto the galaxy to total mass accreted onto the galactic halo, is simply the cosmological ratio $\omega_{\rm b}$ of the baryon to the total matter density,

[TABLE]

Once more neglecting the effect of the growing galaxy on the dark matter halo leads us to deduce that cosmological accretion decreases the energy of a galaxy at a rate \dot{E}_{\rm g}=-(\kappa/2)\omega_{\rm b}{\hbox{\dot{M}{\rm h}}}{\hbox{v{\rm h}}}^{2} (subscript ‘$g$’ for galaxy).

However, unlike the case of the dark matter halo, the growing galaxy generates energy through feedback from stars (and AGN, discussed later), therefore

[TABLE]

In analogy with Eq. (1), we now speculate that \dot{E}_{\star}\approx{\kappa\over 2}\omega_{\rm b}\,{\hbox{\dot{M}{\rm h}}}{\hbox{v{\rm h}}}^{2}: feedback from star formation compensates the energy loss associated with cosmological accretion so that the galaxy grows at nearly constant energy. Figure 4 supports this Ansatz: it shows that, whereas the energy $E_{\rm h}$ of the dark matter halo (black curve) increases by almost 2 orders of magnitude from a look-back time of 10 Gyr to the present, the energy of the galaxy, $E_{\rm g}$, (blue curve) changes by less than $\sim 50$ per cent over the same time interval.

Most of the energy injected into the galaxy’s interstellar medium (ISM) is associated with star formation (i.e. supernovae and other processes associated with short-lived massive stars), therefore we write $\dot{E}_{\star}$ in terms of the star formation rate, $\dot{M}_{\star}$, and a characteristic velocity $v_{\star}$,

[TABLE]

We can obtain an order of magnitude estimate for $v_{\star}$ by assuming that most of the injected energy is from core collapse supernovae (SNe), which inject $10^{51}$ erg of energy each and occur once per $100/\eta$ solar masses worth of stars formed222$\eta=1.74$ for a Chabrier (2003) stellar initial mass function that consists of stars in the mass range of $[0.1,100]{\rm M}_{\odot}$, of which those with mass $6-100~{}{\rm M}_{\odot}$ explode as a core collapse SN., hence

[TABLE]

The factor $\epsilon$ accounts for radiative loses, with $\epsilon=1$ corresponding to no radiative losses and $\epsilon\ll 1$ when such losses are substantial. Numerical simulations of SNe going off in a range of gas densities (e.g. Thornton et al., 1998, and reference therein), and analytical models of the wind in M82 combined with simulations (e.g. Strickland & Heckman, 2009), suggest that a large fraction of the injected energy is radiated, $1-\epsilon\approx 90~{}\%$. The cooling rate of a radiating plasma also depends on its metallicity, therefore $\epsilon$ is is unlikely to be constant in all galaxies and at all times. In this paper we use $\epsilon$ as a fitting parameter when comparing to the simulations; we used a reference value of $\epsilon=0.091$ in Eq. (14) which is consistent with the expected radiative losses being substantial and yields a round number for $v_{\star}$.

Equation (12) that describes the rate of change of the energy of a galaxy is reminiscent of Eq. (1) that describes the rate of change of a main sequence star: whereas the star loses energy (becomes more bound) through radiative losses, the galaxy becomes more bound as the potential well of its host halo deepens due to cosmological accretion. While the star reacts by compensating the energy loss by nuclear fusion, the galaxy reacts by forming stars that inject energy in the galaxy’s star forming gas. For stars, this results in $\dot{E}=\dot{E}_{\rm nucl}-L=0$, and we propose here that the same is true for a galaxy, \dot{E}_{\rm g}=\dot{E}_{\star}-(\kappa/2)\omega_{\rm b}\dot{M}_{\rm h}{\hbox{v_{\rm h}}}^{2}\approx 0.

Why would the feedback from star formation be just so that $\dot{E}_{\rm g}\approx 0$? Is the equilibrium situation stable in the galaxy’s case, just as it was for the star? To examine this question, suppose that $\dot{E}_{\rm g}<0$ - i.e. $|E_{\rm g}|$ is increasing because the galaxy is currently undergoing too little star formation given the current cosmological accretion rate. With gas in the galaxy getting compressed by the deepening potential well, the internal energy $U_{\rm g}$ of the galaxy will increase. How does that affect the star formation rate?

In the eagle implementation of star formation, an increase in thermal energy per unit mass implies an increase in pressure, $P\propto u^{4}$ from Eq. (55), and hence an increase in star formation rate, $\dot{\rho}_{\star}\propto u^{4(n-1)/2}\approx u^{0.8}$ from Eq. (53) for $n=1.4$ from Eq. (54). Therefore an increase in the accretion rate results in an increase in the star formation rate (and conversely, a decrease in the accretion rate results in a decrease in the star formation rate), so that the equilibrium333If the dynamical time-scales are very short, then self-regulation may fail to keep the galaxy in equilibrium. This may happen for example at high redshift, (e.g. Duffy et al., 2010) situation, $\dot{E}_{\rm g}=0$, is secularly stable - just as in the case of nuclear fusion in a main sequence star, and for a similar reason444To take the analogy further, the galaxy in its galactic halo plays the role of the stellar core in the main sequence star.. We note in particular that the increase in star formation rate due to increased accretion, neither assumes nor requires that the gas mass - the fuel for star formation - increases. In our model, the gas reservoir is not regulating the star formation rate in a galaxy. We also note that stability requires that the star formation rate increases with the ISM’s pressure, but without requiring any detailed form of the dependence of $\dot{\rho}_{\star}$ on $P$: the details of exactly how star formation feedback operates are unimportant for the secular stability of the star formation rate in a star forming galaxy. Another consequence is that the star formation rate in a cosmological galaxy depends very little, if at all, on the star formation law that relates star formation rate to the gas mass555We note this is not true in simulations of an isolated galaxy, for which the simulator specifies the gas fraction..

The star formation rate in our model of feedback-regulated galaxy formation, depends on the stellar Initial mass function (through $\eta$ and the recycle fraction ${\cal R}$ discussed below) and the parameters $\mathbold{\kappa}$ (Eq. 5), $\mathbold{\epsilon}$ (Eq. 14) and $\mathbold{\alpha}$ (Eq. 2.3), which is why we called it $I\kappa\epsilon\alpha$. By computing the star formation rate and stellar mass as a function of halo mass, we next show that $I\kappa\epsilon\alpha$ galaxies lie on a star-forming main sequence.

2.4.2 The main sequence of star forming galaxies

Setting $\dot{E}_{\rm g}=0$ in Eq. (12) for a self-regulating galaxy results in a relation between a galaxy’s star formation rate and the cosmological accretion rate onto its host halo at a given redshift,

[TABLE]

which is the main result of this paper. The right-hand side is the cosmological energy accretion rate onto a halo of given mass. The left-hand side sets the corresponding star formation rate in the galaxy, in terms of the effective energy injection rate per stellar mass formed. Substituting the expressions for the accretion rate $\dot{M}_{\rm h}$ and virial velocity $v_{\rm h}$ from Eqs. (6) and (9) allows us to write the star formation rate as a product of its value at $z=0$, $\dot{M}_{\star,0}$, times a dimensionless function, {\hbox{\Psi_{\star}}}(z),

[TABLE]

The star formation rate scales $\propto M_{\rm h,0}^{5/3}\propto v_{\rm h,0}^{5}$; the function {\hbox{\Psi_{\star}}}(z=0)=1.

Since stars lose mass during stellar evolution, the time integral of the star formation rate does not equal the total stellar mass at some later time. In the ‘instantaneous recycling approximation’,

[TABLE]

where ${\cal R}$ is the fraction of mass originally in stars that is returned back to star forming gas through stellar mass loss; the stellar population models used in eagle have $1-{\cal R}\approx 0.55$ (Wiersma et al., 2009). The stellar mass is in this approximation

[TABLE]

with $m_{\star}(z=0)=1$. To evaluate $M_{\star,0}$ and $\dot{M}_{\star,0}$ we have used the cosmological parameters $\Omega_{\rm b}=0.0482519$, $\Omega_{\rm m}=0.307$, $\Omega_{\Lambda}=1-\Omega_{\rm m}$, $\omega_{\rm b}=0.157$ and $h=0.677$ from Planck Collaboration et al. (2014), and set $a=\bar{a}$ and $b=\bar{b}$ for the redshift dependence of the halo accretion rate from Eq. (8); numerical values in our figures correspond to the more general accretion histories from Correa et al. (2015b), for which $a$ and $b$ depend on $M_{\rm h,0}$.

The specific star formation rate, ${\rm sSFR}$, follows from combining Eqs. (16) and (18),

[TABLE]

This ratio would depend on halo mass and hence also on $M_{\star}$ if (i) we had taken into account that the halo accretion rate depends on halo mass rather than using the average accretion rate from Eq. (8) and (ii) if one or more of the $I\kappa\epsilon\alpha$ parameters were to depend on halo mass.

The expression for the sSFR at $z=0$ from Eq. (LABEL:eq:sSFR) looks suspiciously simple: what sets this numerical value? Tracing back the definitions of the dimensionless functions $\mathcal{H}(z)$ (Eq. 10), {\hbox{\Psi_{\star}}}(z) (Eq. 16) and $m_{\star}$ (Eq. 18), we see that these only depend on cosmology and the growth rate of dark matter halos. Changing the growth rate will change the value of the numerical constant $m_{\star,0}$ in Eq. (18). The only other $I\kappa\epsilon\alpha$ parameter that sets the sSFR is ${\cal R}$, the recycled mass fraction, which depends on the IMF. Therefore the value of the sSFR at $z=0$ depends on cosmology (through the accretion history of halos), and on the fraction of mass returned to the ISM during stellar evolution, ${\cal R}$ - and nothing else. This is of course a consequence of assuming that none of the $I\kappa\epsilon\alpha$ parameters evolve.

The dimensionless functions {\hbox{\Psi_{\star}}}(z), $m_{\star}(z)$, and {\hbox{\Psi_{\star}}}(z)/m_{\star}(z) provide the unique connection between the stellar properties of a galaxy and the properties of its host halo - they are plotted in Fig. 5. The star formation rate of an $I\kappa\epsilon\alpha$ galaxy varies over a factor of $\sim 6$ between $z=0$ and $z=6$, peaking at $z\sim 2$, with half the stellar mass forming below $z\sim 1$. The sSFR increases rapidly with redshift, and is higher than its $z=0$ value by factors of 4.6, 13.7 and 30 at redshifts 1, 2 and 3, respectively.

To summarise: $I\kappa\epsilon\alpha$ predicts a main sequence of star forming galaxies along which the specific star formation rate does not depend on $M_{\star}$, provided the $I\kappa\epsilon\alpha$ parameters themselves do not depend on halo mass. The value of this specific star formation rate increases rapidly with redshift.

2.4.3 Comparison to eagle

We test the ideas put forward in the previous section by comparing the star formation rate of galaxies as a function of halo properties and redshift to that of eagle galaxies. We emphasize that for a given assumed stellar IMF, the $\epsilon$ parameter of the the $I\kappa\epsilon\alpha$ model - a measure of the radiative loses in the ISM of the energy injected by SNe - is the central free parameter that sets the star formation rate in a cosmological halo. It does so by setting the characteristic velocity $v_{\star}$ through Eq. (14). The parameter $\epsilon$ likely depends on the properties of a galaxy’s ISM - presumably $\epsilon$ would be smaller (greater cooling losses) when the ISM is denser and more metal rich. Rather than proposing a more detailed model for this, at this stage we simply keep $\epsilon$ constant. However, the eagle reference simulation has a parameter $f_{\rm th}$ which explicitly changes the amount of energy injected into the ISM per solar mass of stars formed, depending on density and metallicity of the ISM (see Eq. (7) in Schaye et al. 2015). Therefore, to keep the comparison between $I\kappa\epsilon\alpha$ and eagle fair, we compare here to the ‘FBconstnoAGN’ eagle variation, in which $f_{\rm th}$ is kept constant (and which does not include AGN feedback either, see the Appendix for more details). We re-iterate though, that keeping $f_{\rm th}$ constant is not quite equivalent to keeping $\epsilon$ constant, because the cooling losses in eagle still depend on density and metallicity.

The star formation rate predicted by Eq. (16) is compared to the eagle FBconstnoAGN model in Fig. 6, taking $\alpha=1$, $\kappa=5/3$, $\eta=1.7$ and $\epsilon=0.2$ (so that $v_{\star}$ is constant, see Eq. 14); coloured lines are the $I\kappa\epsilon\alpha$ prediction at different redshifts, large dots are the median relation for eagle galaxies with the shaded region encompassing the $25^{\rm th}-75^{\rm th}$ percentile range. Even when keeping $v_{\star}$ constant, Eq. (16) captures accurately the increase in $\dot{M}_{\star}$ with the halo’s virial velocity $v_{\rm h}$ at fixed $z$, as well as the increase in $\dot{M}_{\star}$ with $z$ at fixed $v_{\rm h}$. With only one ‘free’ model parameter $\epsilon$ (which sets $v_{\star}$), we were astonished by the level of agreement between $I\kappa\epsilon\alpha$ and eagle.

In the case of Fig. 6, the increase in $\dot{M}_{\star}$ with $z$ at given $v_{\rm h}$ is due to the increase in the cosmological accretion rate onto a halo with given $v_{\rm h}$ at given $z$, as is apparent from Eq. (15). However, plotting $\dot{M}_{\star}$ as a function of $M_{\rm h}$ (Fig. 7), we see that the redshift dependence is stronger due to the $\mathcal{H}(z)^{2/3}$ dependence of Eq. (16). This is not surprising within the context of our self-regulation model: the star formation rate depends on virial velocity rather than halo mass.

2.5 The M_{\star}-{\hbox{M_{\rm h}}} relation

The stellar mass of a galaxy in Eq. (18) is the product of a dimensional number that depends on the galaxy’s halo mass at $z=0$, $M_{\rm h,0}$, times a dimensionless function $m_{\star}(z)$. This functional dependence allows us to answer the question of what is the M_{\star}-{\hbox{M_{\rm h}}} relation in $I\kappa\epsilon\alpha$ in two different ways, (i) ‘What is the M_{\star}-{\hbox{M_{\rm h}}} relation for a population of galaxies at a given redshift?’, and (ii) ‘How does the M_{\star}-{\hbox{M_{\rm h}}} ratio of a halo evolve?’ The answer to the first question follows from $M_{\star}\propto M_{\rm h,0}^{5/3}$ and {\hbox{M_{\rm h}}}\propto M_{\rm h,0}, therefore

[TABLE]

The value of the exponent can be traced back to the {\hbox{\dot{M}{\rm h}}}\,{\hbox{v{\rm h}}}^{2}\propto{\hbox{M_{\rm h}}}^{5/3} dependence of the star formation rate on halo mass, Eq. (16). We compare the predicted relation to that measured in eagle in Fig. 8: coloured lines are the theoretical predictions at different redshifts, large dots are the median relation for eagle galaxies with the shaded region encompassing the $25^{\rm th}-75^{\rm th}$ percentile range. Given that $I\kappa\epsilon\alpha$ predicts the dependence of $\dot{M}_{\star}$ on $M_{\rm h}$ as a function of $z$ in eagle variation FBconstnoAGN well, it is not very surprising that it also reproduces the relation between $M_{\star}$ and $M_{\rm h}$.

Although $I\kappa\epsilon\alpha$ galaxies lie along a line with M_{\star}/{\hbox{M_{\rm h}}}\propto{\hbox{M_{\rm h}}}^{2/3}, they do not evolve along this line. The M_{\star}-{\hbox{M_{\rm h}}} ratio for a given halo evolves as

[TABLE]

This logarithmic slope is $\approx 1.1$ at $z=0$ and increases with $z$ to become nearly constant at a value of 1.4 for $z\geq 4$. If this slope were $5/3$, then (star forming) galaxies would evolve along the $z=0$ $M_{\star}-M_{\rm h}$ relation so that the stellar mass in a halo of a given mass would be independent of redshift. Because the slope is less than 5/3, the $M_{\star}/M_{\rm h}$ versus $M_{\rm h}$ relation evolves with redshift, in the sense that the stellar mass increases with redshift at a constant halo mass, however, that evolution is not very strong. This is the redshift evolution seen in Fig. 8.

Summarising, we conclude that $I\kappa\epsilon\alpha$ reproduces the relation between halo mass, star formation rate, and stellar mass measured in the FBconstnoAGN eagle variation. The fact that $I\kappa\epsilon\alpha$ reproduces the dependence of $\dot{M}_{\star}$ on $M_{\rm h}$ is particularly encouraging, since it directly tests the very basis of the self-regulation argument of Eq. (15). Interestingly, the star formation rate predicted by Eq. (16) does not depend at all on the galaxy’s gas mass or indeed the assumed star formation law - as long as $\dot{\rho}_{\star}\propto u^{\zeta}$ for some sufficiently large and positive value of the exponent $\zeta$, so that the star formation rate increases with the pressure of the galaxy’s interstellar medium. Instead the star formation rate depends on the cosmological accretion rate, and on $v_{\star}$ - that is, on the efficiency of stellar feedback. We will return to this point in the discussion section.

2.6 The galaxy stellar mass function (GSMF)

We compute the GSMF by combining the M_{\star}-{\hbox{M_{\rm h}}} relation from $I\kappa\epsilon\alpha$ with a model for the evolution of the halo mass function. The Press-Schechter (PS, Press & Schechter 1974) approximation for the actual number density of halos per dex in halo mass (e.g. Reed et al., 2007), at $z=0$, is

[TABLE]

where $n_{0}\approx 1\times 10^{-4}{\rm Mpc}^{-3}$ is a normalisation constant, $M_{\rm ps}\approx 2\times 10^{14}{\rm M}_{\odot}$ a characteristic mass above which the number density of halos falls exponentially, and the exponent $\alpha_{\rm h}\approx 0.9$. In the approximation that all halos grow at the same logarithmic rate, {\hbox{M_{\rm h}}}(z)=M_{\rm h,0}\,m_{\rm h}(z), the halo mass function at redshift $z$ is

[TABLE]

where $n_{0}$ and $M_{PS}$ are redshift independent, and $n_{\rm h}$ is now the co-moving number density of halos per dex in halo mass. Provided the $I\kappa\epsilon\alpha$ parameters are constants, M_{\star}\propto{\hbox{M_{\rm h}}}^{5/3}, and the co-moving number density of galaxies per dex in stellar mass becomes

[TABLE]

In this approximation, the GSMF is just a scaled version of the halo mass function, with a power-law shape at low masses and an exponential cut-off at high masses. However it is well known that the ‘knee’ in the galaxy stellar mass function - above which the exponential sets in - does not correspond to the knee in the halo mass function, but rather is a consequence of AGN feedback (Bower et al., 2006; Croton et al., 2006). We discuss how this can be incorporated in the model in §3.4 below.

It is interesting to note that we can make the same argument that lead to Eq. (24) to the star formation rate of a galaxy and compute the ‘galaxy star formation rate function’, GSRF, the number density of galaxies per dex of star formation rate. Since \dot{M}_{\star}\propto{\hbox{M_{\rm h}}}^{5/3}, the GSRF has the same shape as the GSMF,

[TABLE]

The constants $\dot{M}_{\star,0}(M_{\rm h,0}=10^{12}{\rm M}_{\odot},z=0)$ and $M_{\star,0}(M_{\rm h,0}=10^{12}{\rm M}_{\odot},z=0)$ are the $z=0$ star formation rate and stellar mass of a galaxy in a halo of mass $10^{12}{\rm M_{\odot}}$; the numerical values for these are taken from Eqs. (16) and (18), respectively.

At sufficiently low halo mass, these functions are power laws with slope $3\alpha_{\rm h}/5\approx 0.54$. The co-moving number density of galaxies with a given stellar mass increases monotonically with decreasing redshift $\propto m_{\star}(z)^{\alpha_{\star}}$. The corresponding evolution of the co-moving number density of galaxies with a given star formation rate is $\propto\Psi_{\star}(z)^{\alpha_{\star}}$: this function is not monotonic but peaks around $z\sim 2$. It falls to approximately $0.5$ and $0.54$ times its $z=2$ value at redshifts $z=5.5$ and $z=0$, respectively

The $I\kappa\epsilon\alpha$ star formation in a halo with low $v_{\rm h}$ is much less than the rate at which that halo accretes gas. Indeed, according to Eq.(15), only a fraction \kappa{\hbox{v_{\rm h}}}^{2}/v_{\star}^{2} of the accreted gas goes into stars. What happens to the remaining gas? Also, self-regulation due to feedback from star formation must eventually fail for sufficiently high values of {\hbox{v_{\rm h}}}\approx v_{\star}/\kappa^{1/2}\approx 310~{}{\rm km}~{}{\rm s}^{-1}, since then the star formation rate required to self-regulate would exceed the gas accretion rate. To investigate the consequence of these considerations in more detail, we next examine the gas properties in $I\kappa\epsilon\alpha$.

3 Galactic winds and the failure of self-regulating stellar feedback

3.1 Galaxy sizes and gas fractions

The star formation rate in the $I\kappa\epsilon\alpha$ model does not depend on the gas mass. Instead, we compute the mass of star forming gas by assuming a star formation law. Taking the Kennicutt-Schmidt (Kennicutt, 1998) star formation law and assuming that star forming gas is in an exponential disk with scale-length $R_{\rm gas}$, the (total) star formation rate of a galaxy is related to its gas mass by

[TABLE]

where $A$ and $n$ are the parameters of the Kennicutt-Schmidt law (Eq. 54). We follow Mo et al. (1998) (see also Kravtsov 2013) by assuming666Mo et al. (1998) apply this reasoning to the stellar disk; Navarro et al. (2017) show that $R_{\rm gas}$ scales better with the scale radius of the halo, but since we neglect variations in halo concentration by taking $\alpha=1$, these are equivalent. that disk size, $R_{\rm gas}$, is a constant fraction, $\lambda$, of the halo’s virial radius. Using Eq. (2.3), this yields

[TABLE]

where $\lambda=0.01$ yields a reasonable reference scale-length of $R_{\rm gas,0}=2$ kpc for the galaxy inhabiting a $10^{12}{\rm M}_{\odot}$ halo at $z=0$. Using the $\dot{M}_{\star}-M_{\rm h,0}$ relation from Eq. (16) and the $M_{\star}-M_{\rm h,0}$ relation from Eq. (18) allows us to relate galaxy size to star formation rate and stellar mass,

[TABLE]

Sizes of a galaxies with a given $M_{\star}$ depend on redshift $\propto m_{\rm h}^{1/3}/(\mathcal{H}^{2/3}m_{\star}^{1/5})$. The ratio $m_{\rm h}^{1/3}/m_{\star}^{1/5}$ varies by less than a factor 0.75 below $z=6$, meaning that the size scales approximately as $1/\mathcal{H}(z)^{2/3}=1/(1+z)$ for $z\gg 1$, and slower than that at lower $z$. This agrees rather well with the observed scaling: Allen et al. (2017) quotes a scaling $\propto(1+z)^{-0.97}$ for redshifts 5-7 and van der Wel et al. (2014) quotes a scaling $(1+z)^{-0.75}$ for redshifts 0-3. At a given value of $\dot{M}_{\star}$, $R_{\rm gas}\propto(1+z)^{-1.7}$ for $z\gg 1$, which is steeper than the $(1+z)^{-1.1}$ quoted by Shibuya et al. (2015). The weak dependence of size on mass, $R_{\rm gas}\propto M_{\star}^{1/5}$, is consistent the scaling $R_{\star}\propto M_{\star}^{0.22}$ for the stellar size-$M_{\star}$ relation found by van der Wel et al. (2014).

Substituting Eq. (27) into Eq. (26) yields

[TABLE]

Using the $\dot{M}_{\star}-M_{\rm h,0}$ relation from Eq. (16) allows us to compute the $M_{\rm gas}-\dot{M}_{\star}$ relation, and the result is compared to the eagle simulation in Fig. 9, where we used the values of $A$ and $n$ from Kennicutt (1998). The $I\kappa\epsilon\alpha$ prediction reproduces very well the slope of the relation and the normalisation at $z=0$. The simulated evolution is somewhat weaker than predicted. Although pleasing, the excellent agreement between the theoretical prediction and the simulation is not surprising: galaxies in eagle follow the Kennicutt-Schmidt relation of Eq. (54), and that relation results in galaxies following Eq. (26) at least approximately.

The evolution of the gas mass, as governed by the dimensionless function $m_{\rm gas}(z)$ from Eq. (29), is plotted in Fig. 5. The ratio of $m_{\rm gas}(z)$ over its value at $z=0$, is 0.5, 0.8, 1.2 and 1.4 at $z=4$, 3, 2 and 1, respectively, meaning that the gas mass of a forming galaxy changes by slightly more than a factor of 2 since $z=4$. Therefore assuming that galaxies form stars at nearly constant gas mass is a relatively good approximation below $z\sim 4$; it forms the basis of the equilibrium model of Davé et al. (2012), see also Bouché et al. (2010); Krumholz & Dekel (2012).

3.2 Galactic winds

Galactic winds are a natural outcome of a model in which cosmological accretion sets the star formation rate but a star formation law sets the gas mass. Indeed, conservation of baryon mass requires that

[TABLE]

where $\dot{M}_{\rm w}$ is the rate at which the galaxy loses mass through a galactic wind, and the ratio $\beta\equiv\dot{M}_{\rm w}/\dot{M}_{\star}$ is usually called the ‘mass loading factor’ of the wind. Solving for $\beta$ gives

[TABLE]

where we used the main $I\kappa\epsilon\alpha$ relation of Eq. (15) to relate $\dot{M}_{\star}$ and \dot{\hbox{M_{\rm h}}}. We can compute $\beta$ as a function of redshift and $v_{\rm h}$ or halo mass by integrating this equation using the relation between the gas mass and $\dot{M}_{\star}$ (Eq. 26), the result is shown in Figs. 10 and 11. The velocity of the outflow can be estimated by assuming that the wind conserves energy once launched777This assumption may not be unreasonable because $\epsilon$ already accounts significant radiative losses before the wind is launched.,

[TABLE]

${\cal M}=v_{\rm w}/c_{s}$ is the wind’s Mach number, $\gamma=5/3$ is the adiabatic index, $v_{\phi}^{2}/2$ is the change in potential of a Navarro-Frenk-White (NFW) halo (Navarro et al., 1997) between the launch cite, $R_{L}$, and the location $R$ where it is measured (e.g. Łokas & Mamon, 2001), $c$ is the halo’s concentration parameter which depends on $M_{\rm h}$ and $z$ (e.g. Ludlow et al., 2014) and we assume the launch radius, $R_{\rm L}=R_{\rm gas}$. $v_{\phi}\approx 0$ if the wind speed is measured very close to the launch site, and $v_{\phi}$ equals the escape speed from the halo if the wind speed is measured at infinity. This expression also neglects any ram-pressure the outflow may suffer. If the outflow is cold, ${\cal M}\rightarrow\infty$, and

[TABLE]

Clearly this treatment of the wind is quite approximate and in particular it is not obvious how one should compare our value of $v_{\rm w}$ to observations, in which the wind speed is often expressed in terms of the full width half maximum of an emission line. Fortunately, the behaviour of the mass loading $\beta$ is independent of these considerations, although here it is not so clear whether $\beta$ refers to gas leaving the galaxy or gas leaving the halo.

Given these limitations, we plot $\beta$, and the wind speed, $v_{\rm w}$, at a distance of 5 times the gas scale radius, $R_{\rm gas}$, as a function of halo mass, virial velocity, and redshift in Fig.10. The \beta-{\hbox{M_{\rm h}}} relation evolves with redshift, as is clear from the left panels of the figure, basically because the relation between $\dot{M}_{\star}$ and $\dot{M}_{\rm h}$ depends on virial velocity according to Eq. (15). Most of that redshift dependence is removed if we plot $\beta$ as a function of $v_{\rm h}$, as is seen from the right panels in the figure. As $v_{\rm h}$ increases, $\beta$ decreases and $v_{\rm w}$ increases. Also notice that as $v_{\rm h}$ tends to a critical value of around $v_{\rm h,c}\approx 180$ km s*-1*, $\beta$ drops precipitously whereas the wind speed increases rapidly.

Winds in low $v_{\rm h}$ galaxies are slow and strongly mass-loaded, $\beta\gg 1$, as can be seen from Fig. 10. When $\beta\gg 1$ and making the further approximation that $|\dot{M}_{\rm gas}|\ll\dot{M}_{\star}$, Eqs (30) and (33) combine to

[TABLE]

Therefore, the wind speed tracks the halo’s virial velocity (in low $v_{\rm h}$ galaxies at $z<4$), as is apparent from Fig. 10.

The relation between gas mass and star formation rate that results from the Kennicutt-Schmidt star formation law, Eq. (29), and the equation for the mass-loading of winds, Eq. (31), have interesting consequences, namely (i) the emergence of a mass-metallicity relation, and (ii) the existence of a characteristic value of $v_{\rm h}$ above which self-regulation due to feedback from stars fails. We investigate these next.

3.3 The mass-metallicity relation

The metal mass of the star forming gas, $M_{Z}\equiv Z\,M_{\rm gas}$, changes due to metals synthesised and released by stars, metals accreted, metals lost in a galactic wind, and metals locked-up in long-lived stars. Its rate of change is therefore

[TABLE]

where $y$ is the stellar yield, $Z_{0}$ is the metallicity of accreted gas and $Z_{\rm w}$ is the metallicity of the wind which may differ from that of the gas, for example because enriched gas is more like to be ejected by feedback (see e.g. Creasey et al., 2015). Combining this relation with Eq. (30), which expresses baryon mass conservation, and the main $I\kappa\epsilon\alpha$ relation between $\dot{M}_{\rm h}$ and $\dot{M}_{\star}$ from Eq. (15), we find that provided $Z_{\rm w}=Z$ and $Z_{0}=0$,

[TABLE]

The recycled fraction ${\cal R}$ does not affect $\dot{Z}$ in the instantaneous recycling approximation, and the wind’s mass loading $\beta$ does not affect $\dot{Z}$ provided $Z_{\rm w}=Z$. Integrating this equation in time, we compare the relation between $Z$ and $v_{\rm h}$ as a function of redshift to the results from eagle (simulation FbConstnoAGN) in Fig. 12; the agreement is quite good, with $I\kappa\epsilon\alpha$ showing a somewhat steeper dependence of $Z$ on $M_{\rm h}$ and a lower normalisation at Milky Way-like values of {\hbox{v_{\rm h}}}\sim 140 km s*-1*.

Interestingly though, both $I\kappa\epsilon\alpha$ and eagle show very little evolution of the Z-{\hbox{v_{\rm h}}} relation in Fig. 12. Indeed, Eq. (36) shows that the metallicity of a galaxy tends to a value Z\approx\kappa y{\hbox{v_{\rm h}}}^{2}/v_{\star}^{2} that in fact only depends on a halo’s virial velocity and not explicitly on redshift. In this approximation, the metallicity of a galaxy changes only secularly, tracking the evolution of {\hbox{v_{\rm h}}}^{2}. Such a behaviour is an attractor of Eq. (36): given that $\dot{M}_{\star}/M_{\rm gas}>0$, $\dot{Z}$ is positive (negative) when $Z<y\kappa v_{\rm h}^{2}/v^{2}_{\star}$ (when $Z>y\kappa v^{2}_{\rm h}/v^{2}_{\star}$). Therefore $Z$ approaches the secular value,

[TABLE]

on the gas consumption time-scale, $M_{\rm gas}/\dot{M}_{\star}$. This secular value reproduces the evolution from Eq. (36) very well, as shown by the dashed lines in Fig. 13.

Using this secular expression for Z({\hbox{v_{\rm h}}}), taking $v_{\star}=400$ km s*-1*, and $y=0.04$, $Z_{\odot}=0.0127$ for the total metal yield and solar metallicity as done in the the stellar evolution models collected from the literature by Wiersma et al. (2009), we obtain

[TABLE]

The reference values of $\dot{M}_{\star}$ and $M_{\star}$ for the star formation rate and stellar mass, are taken from Eqs. (16) and (18). The normalisation of this relation, $0.68Z_{\odot}$ for $M_{\rm h,0}=10^{12}{\rm M}_{\odot}$, depends on $I\kappa\epsilon\alpha$ parameters $\propto\left({\kappa\alpha/(\epsilon\eta)}\right)^{3/5}$. The observed normalisation is uncertain but at face value higher than what we find by a factor of two (e.g. Tremonti et al., 2004).

The dependence of $Z$ on $v_{\rm h}$ implies that $I\kappa\epsilon\alpha$ galaxies fall on a mass-metallicity relation, as well as on a star-formation rate-metallicity relation. Similarly to the stellar mass-halo mass relation, we can compute how $Z$ depends on $M_{\star}$ at a given redshift,

[TABLE]

independent of redshift, with the value of the exponent resulting from the {\hbox{v_{\rm h}}}^{2}\propto M_{\star}^{2/5} dependence of Eq. (18). As a galaxy grows in mass, its metallicity increases as

[TABLE]

The evolution of $Z$ at a given stellar mass or star formation rate, is $\propto(m_{\rm h}\mathcal{H})^{2/3}$ according to Eq. (38). With increasing $z$, $m_{\rm h}(z)$ decreases whereas $\mathcal{H}(z)$ increases, resulting in little evolution in the $Z-M_{\star}$ relation. At a given value of $M_{\star}$, $Z$ decreases with increasing $z$ by factors 0.9 and 0.76 compared to its $z=0$ value for $z=2$ and $3$, respectively. The observed evolution is somewhat stronger and better reproduced by the eagle reference model in which $\epsilon$ varies with the local gas properties (De Rossi et al., 2017).

Why does $Z$ depend on $M_{\star}$ in $I\kappa\epsilon\alpha$? The $I\kappa\epsilon\alpha$ metallicity of a galaxy is $Z\approx y\dot{M}_{\star}/(\omega_{\rm b}\dot{M}_{\rm h})$, the ratio of the rate at which stars metal enrich the ISM over the rate at which these metals are being diluted by accreting (primordial) gas. The reason this ratio depends on $M_{\star}$ is that the star formation efficiency depends on $v_{\rm h}$: $\dot{M}_{\star}/(\omega_{\rm b}\dot{M}_{\rm h})\propto v_{\rm h}^{2}\propto M_{\star}^{2/5}$, given that $M_{\star}\propto v_{\rm h}^{5}$. In $I\kappa\epsilon\alpha$ , the origin of the mass-metallicity relation is the dependence of the star formation efficiency on the halo’s virial velocity. The $M_{\star}-Z$ relation evolves because the $M_{\star}-v_{\rm h}$ relation evolves. The first part of this claim agrees with Davé et al. (2012), but the second part does not: in their model, evolution is caused by the increase in metallicity of accreting gas. Note that, as long as the galaxy self-regulates888Clearly this would not true in case of a recent merger which might increase $\dot{M}_{\star}$ and dilute $Z$ by gas accretion. its gas metallicity is set by the instantaneous star formation rate rather than a consequence of the build-up of metals that fail to escape from the potential well of its host halo. In other words, the reason that $Z$ depends on $v_{\rm h}$ is because $\dot{M}_{\star}/\dot{M}_{\rm h}$ depends on $v_{\rm h}$, rather than that it is ‘easier for metals to escape from halos with low $v_{\rm h}$’, as is often claimed. Indeed we have assumed that $Z_{\rm w}=Z$ so that an outflow by itself does not affect $Z$ at all. Instead, low $v_{\rm h}$ halos have galaxies with low $Z$ because they are inefficient at forming stars.

3.4 When stellar feedback fails

The basic $I\kappa\epsilon\alpha$ relation of Eq. (15) between the halo accretion rate and the star formation rate results in $\dot{M}_{\star}\propto v_{\rm h,0}^{5}$, where $v_{\rm h,0}$ is the virial velocity of the halo at redshift $z=0$, so that halos with a large virial velocity form stars at a greater rate. For low values of $v_{\rm h,0}$, only a very small fraction of the accreted baryons are converted into stars with the majority of the accreted gas expelled in an outflow, as discussed in §3.1. The rate of gas accretion increases $\propto v^{3}_{\rm h,0}$ but the star formation rate increases $\propto v_{\rm h,0}^{5}$. Since obviously the star formation rate cannot be higher than when all accreted gas is converted to stars, \dot{M}_{\star}\leq\omega_{\rm b}\dot{\hbox{M_{\rm h}}}, it eventually becomes impossible to satisfy Eq. (15) when $v_{\rm h}$ is larger than the critical value that results from inserting \dot{M}_{\star}=\omega_{\rm b}\dot{\hbox{M_{\rm h}}} in Eq. (15),

[TABLE]

Equation 31 gives a slightly weaker limit when requiring that the mass-loading factor $\beta\geq 0$ so that any outflow decreases the baryon fraction of the halo rather than spuriously increasing it.

As a halo grows and $v_{\rm h}$ increases, $\beta$ starts to drop rapidly to values below $1$, as seen in Fig. 10. At $z\leq 2$, $\beta\rightarrow 0$ for {\hbox{v_{\rm h}}}\rightarrow v_{\rm h,{\rm max}} but $\beta$ already plunges to values $\beta\leq 1$ as $v_{\rm h}$ approaches a somewhat smaller critical velocity. At higher redshifts, this critical halo virial velocity decreases, basically because it is no longer a good approximation to neglect $\dot{M}_{\rm gas}$.

What is the consequence of this failure of self-regulation for halos with too high $v_{\rm h}$? The $I\kappa\epsilon\alpha$ GSMF discussed in §2.6 is a power-law that tracks the power-law shape of the halo mass function. In contrast, the observed GSMF has an exponential cut-off at stellar masses above a characteristic stellar mass. It is thought that feedback from accreting black holes (AGN) suppresses star formation in such massive galaxies and this is the cause of the observed break in the GSMF (e.g. Bower et al., 2006; Croton et al., 2006).

This motivates us to associate the critical velocity above which stars cannot self-regulate galaxy formation with those halos in which AGN regulate galaxy formation instead. Using the subscript ‘agn’ as a mnemonic, we see from Fig. 10 that the onset of AGN activity takes place at a nearly redshift-independent value of $v_{\rm h}$ of order

[TABLE]

for which the corresponding virial temperature is999$m_{p}$ is the proton mass, $k_{\rm B}$ is Boltzmann’s constant and $\mu$ the mean molecular weight.

[TABLE]

In the model described by Bower et al. (2017), seed black holes start to grow exponentially in mass when the outflow that is powered by feedback from star formation ceases to be buoyant in the hot corona that fills the dark matter halo. This causes a build-up of gas that fuels the growth of the black hole. The episode of exponential growth ends when the black hole is sufficiently massive that its feedback regulates the forming galaxy. In practise this results in a significant decrease in $\dot{M}_{\star}/M_{\star}$. This model describes well the behaviour of galaxies in the eagle simulation, with the transition between star formation and AGN feedback regulated galaxies occurring in halos with a virial temperature nearly identical to that of Eq. (43) (McAlpine et al., 2018).

At first sight it seems that the reasoning that led to Eq. (43), ‘stellar feedback fails because $v_{\star}^{2}$, a measure of the thermal energy of feedback-heated gas, is too low compared to \kappa{\hbox{v_{\rm h}}}^{2}’  is very different from that of Bower et al. (2017), ‘stellar feedback fails because outflows are no longer buoyant in the hot corona’. However, the build-up of the hot halo is itself depending on the efficiency of stellar feedback (Correa et al., 2018). Put in terms of $I\kappa\epsilon\alpha$: the higher $\epsilon$, the higher the value of $v_{\rm h}$ above which a hot corona develops (see in particular Fig. 14 in Correa et al. 2018). Within the current interpretation, the failure of stellar feedback is not due to the formation of a hot corona, but rather the formation of a hot halo is facilitated by failing stellar feedback.

The results from previous sections allow us to compute other properties of the halo and the galaxy when {\hbox{v_{\rm h}}}=v_{\rm h,{\rm agn}}, the onset of AGN activity. The halo mass, stellar mass, and star formation rate in a halo with {\hbox{v_{\rm h}}}=v_{\rm h,{\rm agn}} at $z=0$, are

[TABLE]

and the corresponding values at redshift $z$ are

[TABLE]

We do not expect the $I\kappa\epsilon\alpha$ GSMF to be correct for halos with {\hbox{v_{\rm h}}}\geq v_{\rm h,{\rm agn}}. We therefore plot the GSMF discussed in §2.6 up to halos of mass $M_{\star,{\rm agn}}(z)$, and compare to the eagle GSMF (simulation FbConst, in which the stellar feedback efficiency is a constant and which does include feedback from AGN) in Fig. 15. The $I\kappa\epsilon\alpha$ model reproduces the power-law shape of the eagle mass function up to $M_{\rm h,{\rm agn}}(z)$ well, getting the evolution of the normalisation approximately correct as well. The value of $M_{\star,{\rm agn}}(z)$ is close to where eagle predicts a rapid decrease in the number density of galaxies, which is due to the action of AGN feedback in the simulation. The number density of galaxies at the knee decreases with increasing $z$. The previous equations elucidate the reason for this in $I\kappa\epsilon\alpha$. Consider two redshifts $z_{1}$ and $z_{2}$, with $z_{1}<z_{2}$, say. Haloes with {\hbox{v_{\rm h}}}=v_{\rm h,{\rm agn}} at a redshift $z_{2}$ will be more massive at $z=0$ than those that have {\hbox{v_{\rm h}}}=v_{\rm h,{\rm agn}} at a redshift $z_{1}$, by the factor $\mathcal{H}(z_{2})/\mathcal{H}(z_{1})$, which is $\approx\left((1+z_{2})/(1+z_{1})\right)^{3/2}$ for $z_{1}\geq 1$. The corresponding ratio of number densities then follows from the slope of the PS halo-mass function, $\left(\mathcal{H}(z_{2})/\mathcal{H}(z_{1})\right)^{\alpha_{\rm h}}$. For example the co-moving number density at $z=6$ is lower than at $z=1$ by a factor 4.8.

3.5 Reality check

Up to now we have compared $I\kappa\epsilon\alpha$ to an eagle simulation in which the feedback parameters are kept constant (simulation FbConst). That simulation does not reproduce the observed properties as well as the eagle reference simulation. So, how well does $I\kappa\epsilon\alpha$ describe the observations?

For a fiducial value of $v_{\star}=400~{}{\rm km}~{}{\rm s}^{-1}$, $I\kappa\epsilon\alpha$ predicts that a $z=0$ galaxy with stellar mass $M_{\star}=5\times 10^{10}{\rm M}_{\odot}$ has a star formation rate of $\dot{M}_{\star}=3.5~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ and is hosted in a dark matter halo of mass $M_{\rm h}=1.9\times 10^{12}{\rm M}_{\odot}$. For the Milky Way, the inferred values are $M_{\star}=(5\pm 1)\times 10^{10}{\rm M}_{\odot}$, $\dot{M}_{\star}=(1.65\pm 0.19){\rm M}_{\odot}~{}{\rm yr}^{-1}$ and {\hbox{M_{\rm h}}}=(1.1\pm 0.3)\times 10^{12}{\rm M}_{\odot} (Bland-Hawthorn & Gerhard, 2016), respectively. However, the scatter in $M_{\star}$ and $\dot{M}_{\star}$ for a halo with given $M_{\rm h}$ is substantial, and the $I\kappa\epsilon\alpha$ value for $M_{\star}$ is consistent with the abundance matching analysis by Guo et al. (2010) and the star formation rate of $\dot{M}_{\star}=3.5~{}{\rm M}_{\odot}~{}{\rm yr}^{-1}$ falls well within the blue cloud for a galaxy with that $M_{\star}$ in the MPA-JHU DR7101010https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/ catalogue. This reasonable level of agreement is of course not surprising: we chose $I\kappa\epsilon\alpha$’s feedback efficiency parameter $\epsilon$ which sets $v_{\star}$ by comparing to these data sets.

The $I\kappa\epsilon\alpha$ specific star formation rate is $\dot{M}_{\star}/M_{\star}\approx 0.07~{}{\rm Gyr}^{-1}$ at $z=0$, independent of $\epsilon$, as compared to an observed value of 0.1 Gyr*-1* at $M_{\star}=10^{10}{\rm M}_{\odot}$ (see the discussion of the data compilation by Behroozi et al. 2018). The observed sSFR increases to a value of $1~{}{\rm Gyr}^{-1}$ ($2~{}{\rm Gyr}^{-1}$) by redshift $z=1$ ($z=2$, Behroozi et al. 2018), as compared to the $I\kappa\epsilon\alpha$ values of 0.3 (1). The $I\kappa\epsilon\alpha$ values are actually very close to those in eagle (simulation FbConstNoAGN). The faster observed evolution might signal that $\epsilon$ does evolve.

The $M_{\star}\propto M_{\rm h}^{5/3}$ dependence of stellar mass on halo mass according to Eq. (18) results in a redshift-independent low-mass slope of the galaxy stellar mass function of $dn/d\log(M_{\star})\propto M_{\star}^{\approx-0.54}$. The faint-end slope of the Schechter luminosity function (Schechter, 1976),

[TABLE]

is $\alpha_{\rm g}\approx 0.48$ at redshift $z=0$ in the gama ‘z’-band (Loveday et al., 2012), a long enough wavelength so that stellar mass is approximately proportional to z-band luminosity. The level of agreement between the two slopes, 0.54 versus 0.48, is encouraging, but not surprising given that $I\kappa\epsilon\alpha$ reproduces the eagle GSMF at the low mass-end well (Fig. 15). Observationally there is no convincing evolution of this slope out to $z\sim 3$ in the $K$-band (Mortlock et al., 2017), also consistent with the $I\kappa\epsilon\alpha$ prediction of no evolution.

The observed evolution in the location of the knee of the Schechter luminosity function is claimed to be consistent with little or no evolution in the value of the stellar mass at which the transition occurs (e.g. Song et al. 2016) but an alternative interpretation is that the transition occurs at a nearly constant star formation rate. Indeed, according to Parsa et al. (2016), the absolute 1500Å magnitude of galaxies at the knee of the Schechter luminosity function occurs at $M_{1500,c}=-19.6$, -20.3, -20.6 and -20.68 for redshifts $z=1$, 2, 3, 4, respectively. If we make the reasonable assumption that UV-luminosity is proportional to star formation rate, then the star formation rate $\dot{M}_{\star}$ of those galaxies increases compared to the value at $z=1$ by factors $\dot{M}_{\star}(z)/\dot{M}_{\star}(z=1)=1.9$, 2.5, 2.8 and 2.9 at $z=2$, 3, 4 and 5. The prediction from $I\kappa\epsilon\alpha$ follows from Eq. (LABEL:eq:Magn), $\dot{M}_{\star,{\rm agn}}(z)/\dot{M}_{\star,{\rm agn}}(z=1)=1.6$, 2.2, 2.8 and 3.4, respectively, impressively close to the observations.

We conclude from this brief comparison to data that $I\kappa\epsilon\alpha$ reproduces observations of the observed galaxy population and its evolution rather well, although there are some differences too.

3.6 Incorporating AGN feedback

An important limitation of the model as described so far is the absence of AGN feedback. Following the arguments that led us to stellar-feedback self-regulation, an obvious way to include AGN in the model is by modifying Eq. (15) to

[TABLE]

with the understanding that AGN feedback sets in when111111$v_{\rm agn}$, which characterises the energy input by the AGN per unit of mass accreted onto the BH, is not to be confused with $v_{\rm h,{\rm agn}}$ - the virial velocity of the halo above which stellar feedback fails. {\hbox{v_{\rm h}}}\gtrapprox v_{\rm h,{\rm agn}}. Here, $\epsilon_{r}\approx 0.1$ is the radiative efficiency of the AGN and $\epsilon_{f}\approx 0.15$ the fraction of radiated energy that couples to the gas (see the discussion in Schaye et al. 2015, their section 4.6). As in the case of feedback from star formation discussed in §2.4.1, AGN feedback will be self-regulating provided that the black hole accretion rate increases with the pressure in the ISM.

We can examine what to expect for the black hole mass of a halo with {\hbox{v_{\rm h}}}\gg v_{\rm h,{\rm agn}} by integrating Eq. (47) for $v_{\star}=0$,

[TABLE]

where the second line is at redshift $z=0$; this scaling is plotted as a blue line in Fig.16. How does this compare to observations? Ferrarese (2002) claim that black hole mass scales with the circular speed as $M_{\rm BH}\propto v_{c}^{5.5}$. This scaling is shown as a black dashed line and is close to Eq. (LABEL:eq:BH). Kormendy & Ho (2013) argue that, because the scatter in the $M_{\rm BH}\hbox{--}v_{c}$ relation is large at low $v_{c}$, the Magorrian relation (Magorrian et al., 1998) between black hole mass and bulge mass, is more fundamental. We would argue instead that low-mass black holes are not in a self-regulating regime.

The observed relation between black hole mass and (3D)121212We have assumed that $\sigma_{\star}^{2}=3\sigma^{2}$, where $\sigma$ is the line-of-sight stellar velocity dispersion. stellar velocity dispersion $\sigma_{\star}$, is $M_{\rm BH}=1.1\times 10^{7}\,(\sigma_{\star}/200~{}{\rm km}~{}{\rm s}^{-1})^{5.12}~{}{\rm M}_{\odot}$ (McConnell et al., 2011). Provided {\hbox{v_{\rm h}}}\sim\sigma_{\star}, the observed dependence on velocity is close to our prediction while the normalisation requires reasonable values for $\epsilon_{r}$ and $\epsilon_{f}$. The scaling of the $M_{\rm BH}-\sigma_{\star}$ relation in the model by Silk & Rees (1998) is identical to ours basically because both are based on energy arguments, however, our normalisation is significantly more realistic, as shown in Fig. 16. The model by King (2003) is based on momentum arguments; their scaling, $M_{\rm BH}\propto\sigma_{\star}^{4}$, is shallower than observed. Booth & Schaye (2010) obtain a $M_{\rm BH}\propto M_{\rm h}^{1.55}$ scaling by arguing that the net total energy injected by an AGN is of order of the binding energy of a halo. This is somewhat similar to our reasoning, except that we argue that it is the rate of energy injection by the AGN that tracks the rate of energy accretion by the halo due to self-regulation. The secular growth rate of a black hole - and hence the time-averaged luminosity of the AGN - therefore depends on the cosmological accretion rate onto its host halo and therefore on redshift, and not just on halo properties.

4 Discussion

4.1 Comparison to previous work

The paper by Bouché et al. (2010) sparked interest in trying to understand the basic physics underlying self-regulation of galaxies. That paper, and several that followed, contain equations that resemble those of Section 3 - but the underlying assumptions are sometimes strikingly different, as we discuss here. The starting point of Bouché et al. (2010) is their realisation that the dependence of $\dot{M}_{\star}$ on stellar mass and redshift, resembles that of the cosmological accretion rate, suggesting that the gas accretion rate $\dot{M}_{\rm gas,acc}\propto\dot{M}_{\rm h}$. The proportionality constant is argued to be less than $\omega_{b}$, the cosmological gas to total matter density, because only cold accreted gas is assumed to be eligible for star formation. The resulting star formation rate, is then determined by the efficiency with which gas is converted into stars - that is - by the star formation law.

This reasoning results in $\dot{M}_{\star}\propto\omega_{b}\dot{M}_{h}$, as in our Eq. (16), with the important distinction that the efficiency of galaxy formation,

[TABLE]

is set by the efficiency of star formation,

[TABLE]

where $\tau_{d}$ is a characteristic time that still needs to be determined. The onus of getting the observed M_{\star}/{\hbox{M_{\rm h}}} relation is now wholly on the star formation law, Eq. (50). The solution advocated by Bouché et al. (2010) , is to assume that halos do not form any stars as long as their halo mass is below some minimum value, $M_{h,min}\approx 10^{10}-10^{11}{\rm M}_{\odot}$, which conspires to result in $\epsilon_{g}$ increasing with $M_{\rm h}$. They stress repeatedly that their results are completely independent of the efficiency of feedback.

Lilly et al. (2013) build on this work, and in their ‘gas regulator’ frame work, $\dot{M}_{\star}$ is regulated by the gas reservoir of the galaxy, $M_{\rm gas}$ in our notation. Rather than assuming a minimum halo mass $M_{h,min}$ below which no stars form, the model introduces two main fitting parameters, which in our notation are the product $\epsilon_{\star}\,\tau_{d}$ (their variable $\epsilon$) and $\beta$ (their variable $\lambda$). In the follow-up paper by Birrer et al. (2014), they show how the evolution of galaxies over cosmic time can be modelled well once $\epsilon$ and $\lambda$ are parameterised as functions of $M_{\star}$. Note that these cannot be independent of $M_{\star}$, since otherwise the ratio M_{\star}/{\hbox{M_{\rm h}}} is constant as well, since a constant fraction of the accreted gas is converted into stars.

The ‘minimum bathtub’ model described by Dekel & Mandelker (2014) has very similar ingredients, in that $\dot{M}_{\star}$ is also regulated by $M_{\rm gas}$ through the star formation law. These authors stress that many properties of galaxies follow from this model if it is assumed that the system is in a quasi-steady state, $\dot{M}_{\rm gas}=0$.

These models ‘self-regulate’ in the sense that the star formation rate is determined by the gas mass by mass conservation, in our notation \dot{M}_{\rm gas}=\omega_{b}\dot{\hbox{M_{\rm h}}}-(1-{\cal R}+\beta)\dot{M}_{\star} (Eq. 30), so that too much star formation depletes the gas reservoir which ultimately decreases $\dot{M}_{\star}$. Conversely too little star formation leads to a build-up of $M_{\rm gas}$, and through the star formation law, this increases $\dot{M}_{\star}$. A very nice feature of these models, in addition to prediction correctly the rapid increase in $\dot{M}_{\star}/M_{\star}$ with redshift because the gas accretion rate \propto{\hbox{\dot{M}_{\rm h}}}, is that they correctly predict secondary parameter dependencies, for example the fact that galaxies that lie above the main sequence are more gas rich and more metal poor, see also Dayal et al. (2013).

What all these models have in common is that the star formation rate is set by the gas reservoir through the star formation law. The origin of that law is not discussed in detail, but presumably it results from a balance between cooling and heating from star formation, as originally envisioned by White & Frenk (1991). In these models, feedback from star formation is only important in setting the star formation law, basically parameterised by $\epsilon_{\star}$. Combined with a model for the build-up of dark matter halos, or using dark matter-only simulations that follow the growth of halos, these ‘self-regulation’ models are very successful in building realistic looking mock universes (see e.g. Moster et al., 2018; Tacchella et al., 2018).

In our opinion, there are two major weaknesses to this basic model: (i) to be predictive the model needs to be able to predict how the efficiency of star formation, $\epsilon_{\star}$, and the mass loading factor, $\beta$, depend on halo (or stellar mass), a formidable task. More worryingly, (ii) there is evidence that one of the main assumptions - that the star formation rate depends on the gas mass through the star formation law - is not quite correct.

At first sight it seems impossible that the rate of star formation in a galaxy is not dependent on the star formation law - and in fact it would be if the galaxy were isolated. However a galaxy in a cosmological setting can gain mass through accretion and lose it through winds - and therefore the amount of gas in the reservoir is not some constant, rather $M_{\rm gas}$ too is set by the physics of galaxy formation. Demanding that $\dot{M}_{\star}$ depends on $M_{\rm gas}$ through a star formation law, and vice versa, results in a ‘chicken and egg’ problem.

Numerical simulations can be very helpful in distinguishing cause from effect. The owls simulations described by Schaye et al. (2010) are cosmological hydrodynamical simulations performed with gadget (Springel, 2005), but the parameters of sub-grid models are varied over a very wide range and not calibrated to observations as in eagle. In particular, the owls simulation suite includes parameter variations in which the efficiency of feedback from stars (i.e. the value of $v_{\rm\star}$ in our notation) and the star formation law (the values of $A$ and $n$ in Eq. 54), are varied separately. By plotting $M_{\rm gas}$ and $\dot{M}_{\star}$ versus a variable that does not depend on either $v_{\star}$ or the star formation law, such as halo mass, $M_{\rm h}$, it becomes possible to test the very core assumption of the gas-regulator or bath-tub models.

Haas et al. (2013a) compares models with the same star formation law (same value of $A$ and $n$) but different values of the feedback efficiency. Compare in particular their models REF and WML4: these have identical numerical parameters, except that the value of $v_{\star}^{2}$ in simulation WML4 is twice that of REF. Maybe not surprisingly, \dot{M}_{\star}/{\hbox{M_{\rm h}}} in the simulation with the stronger feedback is about half as large as in REF (their Fig. 4). Because the star formation law in these simulations is the same, this also implies that M_{\rm gas}/{\hbox{M_{\rm h}}} is also approximately half in WLM4 compared to REF, as is also born out by the same figure.

However now compare models REF and SFAMPLx3 in Haas et al. (2013b): these have identical feedback parameters, but the value of $A$ (from Eq. 54) in simulation SFAMPLx3 is three times that in simulation REF. Figure 5 in Haas et al. (2013b) shows that nevertheless the ratio \dot{M}_{\star}/{\hbox{M_{\rm h}}} is nearly identical in the two simulations: the star formation rate in a halo of given mass is not, or only very weakly dependent, on $A$: a direct violation of the main assumption in the ‘gas-regulator’ models. Given that the star formation rates are the same in these models, but the star formation law differs, this must imply that the gas reservoir in SFAMPLx3 is less than that in REF at a given value of $M_{\rm h}$: the same figure 5 shows that indeed M_{\rm gas}/{\hbox{M_{\rm h}}} is about a factor of three lower in model SFAMPLx3 compared to REF. As stressed by Haas et al. (2013b) and confirming what was found by Schaye et al. (2010): stellar feedback regulates the star formation rate by determining the amount of (star forming) gas. In this interpretation, $\dot{M}_{\star}$* regulates $M_{\rm gas}$ through stellar feedback, rather than $M_{\rm gas}$ setting $\dot{M}_{\star}$ through a star formation law*.

The model presented by Davé et al. (2012) incorporate self-regulation through feedback, as envisioned here. Because they limit their analysis to equilibrium states defined by $\dot{M}_{\rm gas}=0$, their results are actually very similar to the various incarnations of the bath-tub models.

In our interpretation, self-regulation follows from energy conservation, Eq. (15), and in particular the fact that $\dot{E}_{g}=0$ is a secularly stable equilibrium (provided that $\dot{\rho}_{\star}$ increases with pressure of the star forming gas). Therefore accretion sets the star formation rate, once the net energy input generated by forming stars is known. This sets the ‘efficiency of galaxy formation’ (the ratio of the star formation rate over the cosmological baryon accretion rate onto a halo) to be

[TABLE]

which does not depend on the star formation law but on the properties of the halo (through $v_{\rm h}$) and the efficiency of feedback (through $v_{\star}$) This is in contrast to Eq. (50). The star formation law then determines the gas reservoir in the $I\kappa\epsilon\alpha$ model, with any excess accreted gas expelled in a wind.

Combining the main $I\kappa\epsilon\alpha$ relation between the star formation rate and the accretion rate on a halo of Eq. (15), with the relation between halo virial velocity and halo mass (Eq. 3) and the equation for the growth of a halo (Eq. 6), allows us to write the star formation rate in terms of the halo accretion rate in the form

The relation between stellar mass and halo mass (Eq. 18) can be cast in the form

[TABLE]

which has the form of Eq. (1) in the paper by Salcido et al. (2019), with their $\epsilon(M_{h},z)=5/3$. These author show that a halo mass - stellar mass of this form can be integrated to give analytical relations for the galaxy stellar mass function and the evolution of the cosmic star formation rate density.

4.2 Limitations of the model

A forming galaxy can fail to be able to attain its equilibrium star formation rate given by Eq. (15) for several reasons. Consider for example what happens if $\dot{M}_{\rm h}$ suddenly decreases - for example because the galaxy becomes a satellite. Star formation will nevertheless continue in accordance with the star formation law, depleting the gas reservoir. In such galaxies, the star formation rate is set by the gas consumption time scale, rather than regulated by feedback. A less extreme version of the same phenomenon occurs when $\dot{M}_{\rm h}$ for a particular halo is unusually small compared to the ensemble average. The $I\kappa\epsilon\alpha$ model does not correctly describe this situation and in particular is not applicable to satellite galaxies.

We have neglected the finite lifetimes of massive stars. We think this is unlikely to be a major limitation at lower redshifts when the dynamical time of any galaxy is much larger than the lifetimes of massive stars. However, the limitation may affect the onset of star formation in small galaxies at high redshift. When $v_{\rm h}$ is very low, gas cannot cool and our self-regulation argument will not correctly predict $\dot{M}_{\star}$. When the halo grows in mass it may pass the threshold where gas can cool on a short time scale, and star formation may be unable to self-regulate because of the finite lifetimes of massive stars. This may lead to a star burst which $I\kappa\epsilon\alpha$ does not model correctly.

Not unrelated is what happens at high values of $v_{\rm h}$ at low redshift. The $I\kappa\epsilon\alpha$ model predicts that feedback becomes inefficient for {\hbox{v_{\rm h}}}\approx 180 km s*-1* following similar reasoning to Bower et al. (2017). We argued, as did Bower et al. (2017), that the resulting increase in gas mass triggers the AGN, which, once the black hole mass has increased sufficiently, will regulate the galaxy. However, by construction this occurs in the same halos that develop a hot halo of gas, so that it becomes unlikely that the right hand side of Eq. (15) describes correctly the rate at which gas enters the galaxy: it may simply add to the hot halo instead (see the discussion in Bouché et al. (2010) on hot versus cold accretion). We think therefore that it is unlikely that $I\kappa\epsilon\alpha$ models such galaxies accurately. Moreover, galaxy-galaxy mergers contribute significantly to the mass growth of such galaxies, and we have not attempted to include these in the model either.

We also neglected that stars may form from gas lost by previous generations of stars - such recycling may affect the star formation rate of galaxies at late times when their stellar masses are high but the cosmological accretion rate low (e.g. Oppenheimer et al., 2010; van de Voort, 2017). Gas lost from galaxies by winds may re-accrete later - again we have neglected this effect. More in general, we have neglected the possibility that the accretion rate differs from \omega_{b}{\hbox{\dot{M}_{\rm h}}}.

If Eq. (15) is indeed applicable, then it might be possible to estimate the scatter around the main sequence of star forming galaxies from the scatter of $\dot{M}_{\rm h}$ around the ensemble average. This would provide a good test of the basic assumption in our model.

5 Summary and conclusions

We have presented a model for star formation in galaxies that is motivated by the origin of the stability of nuclear fusion in main sequence (MS) stars. The energy generated by nuclear fusion in a MS star equals the rate at which energy is lost through radiation. This equilibrium is secularly stable because if the star loses energy, it heats up, which increases the rate at which fusion occurs. The analogy with a star forming galaxy is that the rate of energy injection by supernovae (and winds from their massive progenitor stars) equals the rate at which energy is lost due to cosmological accretion. This equilibrium is stable provided the star formation rate increases with the pressure of the star forming gas.

Equation (15), (1/2)\dot{M}_{\star}v_{\star}^{2}=(\kappa/2)\omega_{b}{\hbox{\dot{M}{\rm h}}}\,{\hbox{v{\rm h}}}^{2}, encapsulates this energy balance. Here, $v_{\star}^{2}$ is a measure of the effective energy injected per unit mass of star formed by feedback, so that the left hand side is the rate at which feedback increases the galaxy’s energy. The right hand side of the equation is the energy loss term due to cosmological accretion ($\omega_{b}$ is the cosmological baryon to total mass fraction), with {\hbox{v_{\rm h}}}^{2} a measure of the depth of the dark halo’s potential. In our ‘$I\kappa\epsilon\alpha$’ model, the star formation rate is set by the cosmological accretion rate by energy balance. The predicted dependence of $\dot{M}_{\star}$ on redshift and virial velocity, $v_{\rm h}$, or halo mass, $M_{\rm h}$, agrees very well with that measured in the eagle cosmological hydrodynamical simulation (Schaye et al., 2015), as shown in Figs. 6 and 7, respectively.

The $I\kappa\epsilon\alpha$ model has four parameters (I, $\kappa$, $\epsilon$, and $\alpha$; hence the name), which together shape the star forming sequence of galaxies. The parameter ‘I’ stands for the (stellar) Initial mass function (IMF), which sets how much energy is available for feedback from star formation, in particular from the supernovae (SN) associated with star formation, as well as the recycled fraction ${\cal R}$ that relates the time integral of star formation to the stellar mass formed. We have kept the IMF constant in this paper. The dimensionless parameters $\kappa$ and $\alpha$ quantify the rate of cosmological accretion onto a halo ($\kappa$), and the concentration of such halos ($\alpha$, see Eq. 2.3). We find that $\kappa\approx 5/3$ and $\alpha\approx 1$, and have kept these constant as well.

We think that the main numerical parameter that affects our results is $\epsilon$, which is a measure of the fraction of the energy that is injected by SNe that effectively increases the energy of the star forming gas, rather than being radiated away. It relates $v_{\star}^{2}$ to the energy produced by SNe per unit mass (or more generally to the energy injected in the ISM as a result of recent star formation), see Eq. (14). If feedback is efficient, $\epsilon$ is large, and $\dot{M}_{\star}$ is small. The eagle simulation has a parameter, $f_{\rm th}$, that controls what fraction of the available supernova energy is injected into the star forming gas. This means that $f_{\rm th}\approx\epsilon$, provided radiative loses are small. Because feedback is efficient131313Gas heated by SNe has its temperature increased by $\Delta T=10^{7}$ K where its cooling rate is minimal and mostly independent of metallicity. in eagle, radiative loses in SN-heated gas are mostly small, which explains why the $I\kappa\epsilon\alpha$ model reproduces the eagle model with $f_{\rm th}$ held constant relatively well. However, in the eagle reference model, $f_{\rm th}$ is allowed to vary as a function of density and metallicity in a way that is calibrated so that the simulation reproduces (some) observations. Therefore to improve the agreement of $I\kappa\epsilon\alpha$ with data, we would need to understand how radiative loses depend on the interstellar medium of a star forming galaxy. It seems unlikely that there is a simple way to do so.

A striking feature of the model is that $\dot{M}_{\star}$ does not depend on the gas mass, $M_{\rm gas}$, unlike what is assumed in many self-regulating models (e.g. Bouché et al., 2010; Lilly et al., 2013). We use a star formation law (in our case the Kennicutt-Schmidt law, Kennicutt 1998) to infer $M_{\rm gas}$ from $\dot{M}_{\star}$ - rather than the other way around. Doing so allows us to reproduce the $M_{\rm gas}-\dot{M}_{\star}$ relation in eagle  (Fig. 9) as well as the mass-metallicity relation (Fig. 12).

We tried to incorporate feedback from accreting black holes (AGN) by (i) identifying when feedback from star formation fails so that a black hole can grow, and (ii) include AGN in the self-regulation process. Stellar feedback fails in galaxies with deep enough potential wells, so that energy injected by stars cannot compensate for energy lost through accretion even if all accreted gas is converted into stars. We showed that this occurs in halos with virial velocity above a nearly redshift independent critical value of $\sim 180~{}{\rm km}~{}{\rm s}^{-1}$. Demanding that the AGN regulates galaxy formation results in a relation between the black hole mass and the virial velocity of the halo of the form M_{\rm BH}\propto{\hbox{v_{\rm h}}}^{5}, which closely follows the observed relation.

In the Introduction we discussed how gas cooling is thought to play an important role in determining the rate at which a galaxy forms stars, to the extent that it may even be the main property that determines the location of the peak in the redshift evolution of the star formation rate density of the Universe (Hernquist & Springel, 2003). Numerical simulations at first sight support this claim directly: a simulation where the contribution from metals is not included when calculating the cooling rate - and hence where the cooling rate is lower - yields lower values of M_{\star}/{\hbox{M_{\rm h}}} than when metals are included (compare models NOZCOOL and model REF in Fig.3 of Haas et al. 2013a). However, a lower metallicity of star forming gas reduces cooling loses of injected feedback energy, increasing $\epsilon$ and hence reducing $\dot{M}_{\star}$: that sequence of events is also consistent with the findings from Haas et al. 2013a). The main impact of metallicity on the cooling rate of the gas may be on the efficiency of feedback, rather than on the accretion rate. Of course this argument breaks down in halos where the virial temperature is so high that most of the gas is and remains hot.

We think that $I\kappa\epsilon\alpha$ provides a simple way of calculating the properties of a galaxy in terms of those of its host halo - and the results so obtained agree reasonably well with those from much more sophisticated models and importantly also with data. We suggest that a better description of how cooling losses depend on the properties of a galaxy through its history would improve the quality of the theoretical prediction.

Acknowledgements

We thank our colleagues (J. Schaye, M. Schaller, R. Crain and R. Bower) for sharing with us the data from the eagle simulation; and we are grateful to L. Heck and J. Helly for providing the computing support. This study was funded by the Science and Technology Facilities Council [grant number ST/F001166/1]. The study made use of the DiRAC Data Centric system at Durham University which is run by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk); the equipment was funded by BIS National E-Infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC; as a part of the National E-Infrastructure. This research was also supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. The International Centre for Radio Astronomy Research (ICRAR) is a Joint Venture of Curtin University and The University of Western Australia, funded by the Western Australian State government. M.S. is supported by an ASTRO 3D fellowship at the ICRAR, Curtin University.

Appendix A Simulation details

We compare the results of the model described in this paper to galaxies from the eagle simulation, which we briefly describe here. ‘Evolution and Assembly of Galaxies and their Environments’  (eagle) is a suite of cosmological, hydrodynamical simulations, performed using an evolution of the gadget smoothed particle hydrodynamics code described by Springel (2005). eagle uses a set of sub-grid modules to encode unresolved physics, described in detail by Schaye et al. (2015), which we briefly summarise here.

The sub-grid modules contain a set of numerical parameters, whose values are calibrated in a reference run to reproduce a small number of $z\approx 0$ observables, namely the galaxy stellar mass function, the relation between galaxy stellar mass, $M_{\star}$, and size, and between $M_{\star}$ and black hole mass, as detailed by Crain et al. (2015). Given these calibrated values, the simulation also reproduces several observable relations that were not part of the calibration, in particular yielding a ‘main sequence’ of blue star forming galaxies in which $\dot{M}_{\star}$ depends on $M_{\star}$ and redshift as observed (Furlong et al., 2015), as well as a ‘red sequence’ of quenched galaxies (Trayford et al., 2015, 2017). The $z=0$ galaxy colours correlate with galaxy morphology as observed (Correa et al., 2017; Trayford et al., 2018).

Most relevant for the comparisons in this paper are the implementation of star formation, of stellar feedback, and of feedback from accretion blackholes (AGN) in eagle:

• •

Star formation: sufficiently dense gas in eagle is converted into star particles at a rate per unit volume, $\dot{\rho}_{\star}$, that depends on the gas pressure, $P$, as

[TABLE]

The normalisation of this relation and the exponent $n$ are set by the Kennicutt-Schmidt law (Kennicutt, 1998) that relates the surface density of star formation, $\dot{\Sigma}_{\star}$ and of gas, $\Sigma_{\rm g}$,

[TABLE]

The underlying assumption connecting these relations is that volume and surface densities are related by the local Jeans length, as motivated by Schaye & Dalla Vecchia (2008). The simulation does not resolve the multi-phase nature of the interstellar medium and star forming gas is assumed to have a minimum pressure (Schaye et al., 2015),

[TABLE]

where $u$ is the thermal energy per unit mass.

• •

Stellar feedback is implemented as described by Dalla Vecchia & Schaye (2012): a newly formed star particle increases the temperature of surrounding gas by an amount $\Delta T$. The quantity of gas heated depends on the effective energy injected by star formation, $f_{\rm th}\Delta E$, where $\Delta E$ is the total energy released by the winds from massive stars and core-collapse supernovae which in turn depends on the assumed stellar initial mass function (IMF). The value of $\Delta T$ is chosen such that gas is heated to a temperature where its cooling rate is small: this makes the feedback efficient. The value of $1-f_{\rm th}$ quantifies the fraction of injected energy that is lost from the star forming region, for example through radiative cooling; $f_{\rm th}$ is one of the main calibration parameters in eagle.

• •

Black holes and AGN: The seeding, merging, accretion, and feedback from black holes (BHs) in eagle is described by Rosas-Guevara et al. (2015). Seed BHs are inserted in each dark matter halo once it becomes sufficiently well resolved. When a BH accretes mass and becomes an AGN, it injects thermal energy in the surrounding gas.

The origin of red galaxies in eagle is investigated by Trayford et al. (2016). Ram-pressure stripping and ‘strangulation’ dramatically decreases the star formation rate of satellite galaxies, causing them to leave the blue cloud of star forming galaxies and settle onto the red sequence. The simple self-regulating model described in this paper does not attempt to describe these effects, and we will therefore only compare to central, i.e. non-satellite, eagle galaxies. Similarly, AGN feedback suppresses star formation in massive galaxies, causing them to become passive. Since that mechanism is also not included in the model, most of the comparison in this paper i to eagle variation FbConstNoAGN, in which $f_{\rm th}$ is a constant, and which does not include AGN feedback. We also use variation FbConst, in which $f_{\rm th}$ is kept constant and which does include AGN.

Table 1 contains a list of parameters of the eagle runs that we used. Simulation ‘REF’ is the default eagle model from Table 2 of Schaye et al. (2015). The simulation FBconst with $f_{\rm th}=1$ appears in Table 1 of Crain et al. (2015), simulation. Simulation DMO is a dark matter-only version of the same volume. All simulations are initialised from the same Gaussian initial conditions, so that halo masses are nearly identical in all runs.

Galaxies of the eagle reference model look like observed galaxies in many of their properties. Keeping $f_{\rm th}$ constant, the simulated galaxies have similar stellar masses and star formation rates, but are typically smaller than in the reference model. Therefore this model is not as good a representation of the real galaxy population, but we believe its physics is still reasonable - and it is much easier to compare to our simple model. Many of the properties of the population of eagle galaxies can be extracted directly from the public database141414http://icc.dur.ac.uk/Eagle/database.php. (McAlpine et al., 2016), which we used extensively in preparing the figures.