Bayesian inferences of galaxy formation from the K-band luminosity and H i mass functions of galaxies: constraining star formation and feedback

Abstract

We infer mechanisms of galaxy formation for a broad family of semi-analytic models (SAMs) constrained by the K-band luminosity function and H i mass function of local galaxies using tools of Bayesian analysis. Even with a broad search in parameter space the whole model family fails to match to constraining data. In the best-fitting models, the star formation and feedback parameters in low-mass haloes are tightly constrained by the two data sets, and the analysis reveals several generic failures of models that similarly apply to other existing SAMs. First, based on the assumption that baryon accretion follows the dark matter accretion, large mass-loading factors are required for haloes with circular velocities lower than 200 km s⁻¹, and most of the wind mass must be expelled from the haloes. Second, assuming that the feedback is powered by Type II supernovae with a Chabrier initial mass function, the outflow requires more than 25 per cent of the available supernova kinetic energy. Finally, the posterior predictive distributions for the star formation history are dramatically inconsistent with observations for masses similar to or smaller than the Milky Way mass. The inferences suggest that the current model family is still missing some key physical processes that regulate the gas accretion and star formation in galaxies with masses below that of the Milky Way.

methods: numerical, methods: statistical, galaxies: evolution, galaxies: formation, galaxies: luminosity function, mass function

1 INTRODUCTION

The formation and evolution of galaxies is regulated by mass accretion, star formation, and associated physical processes, such as feedback from supernova (SN) explosions. Unravelling this physical complexity is one of the great challenges in modern astrophysics (see Mo, van den Bosch & White 2010, for an overview). An important hurdle for any theoretical model of galaxy formation is a natural explanation for the differing shapes of the dark matter halo mass function and the galaxy luminosity and stellar mass function. The mass function of dark matter haloes, n(M), scales with halo mass as n(M) ∝ M⁻² at the low-mass end, in contrast to the observed luminosity function of galaxies at low z, Φ(L), which has a shallow faint-end slope: Φ(L) ∝ L⁻¹. This suggests that star formation in low-mass haloes must be inefficient, but the physical cause remains unexplained. Similarly, the dependence of the baryonic mass fraction on halo mass requires explanation. One might expect that each halo has a baryonic mass fraction that is close to the universal value of 17 per cent. In reality, the total inventory of baryons in Milky Way-sized haloes is only 5–10 per cent of the halo mass, and the baryon mass fraction decreases rapidly with decreasing halo mass at the low-mass end (e.g. Yang, Mo & van den Bosch 2003; Papastergis et al. 2012). Therefore, some processes must keep a large fraction of baryons outside of galaxies.

The process most often considered to suppress star formation in low-mass haloes is SN feedback; the total amount of energy released by SNe can be larger than the binding energy of the gas in low-mass haloes (Dekel & Silk 1986; White & Zaritsky 1992). Therefore, it is energetically possible to expel large amounts of baryonic matter from low-mass haloes and, thereby, reduce the efficiency of star formation. Feedback may do this in multiple ways. For example, the feedback energy may heat a fraction of the cold interstellar medium (ISM) in the disc, either causing it to mix with hot halo gas or ejecting it directly without mixing (Somerville & Primack 1999; Benson et al. 2003). The feedback energy may also interact with the halo gas, heat it, and then eject it from the halo (Croton et al. 2006). Alternatively, galactic winds may be driven by radiation pressure associated with massive stars and SN explosions (e.g. Murray, Quataert & Thompson 2005). The total amount of energy in radiation is two orders of magnitude larger than the available kinetic energy in SN explosions, and so the energy budget can be increased significantly. Some simulations of isolated galaxies demonstrate that radiation pressure-driven (also called momentum-driven) winds are effective at ejecting cold gas from galaxies (e.g. Hopkins, Quataert & Murray 2011; Hopkins et al. 2013), but other simulations question their efficiency (e.g. Krumholz & Thompson 2012). Furthermore, cosmological simulations that employ simplified versions of the momentum-driven model have not successfully reproduced the observed galaxy mass function (e.g. Oppenheimer & Dave 2006; Puchwein et al. 2012), although in the latest simulations the discrepancies are predominantly at high masses where SN winds may not be the dominant feedback mechanism (Davé et al. 2013). In summary (1) none of the feedback processes is understood in detail; and (2) the effect of feedback on the star formation efficiency is poorly determined; the efficiency parameter adopted in the models ranges from a few per cent to about 100 per cent (Bower, Benson & Crain 2012; Mutch, Poole & Croton 2013), with a large range of uncertainty.

Observationally, significant outflows are observed only for starburst galaxies, and there is no clear evidence for massive winds from normal spiral galaxies (e.g. Heckman 2003). The mass loading and velocity of galactic winds are poorly constrained (Martin 2005; Weiner et al. 2009; Chen et al. 2010). Martin (2006) estimated that mass outflow rates (OFRs) in the cold component of the wind is on the order of 10 per cent of the star formation rate (SFR) in the starburst, while other observations suggest that the rate may be as high as or even higher than the SFR of the galaxy (Martin, Kobulnicky & Heckman 2002).

The feedback processes described above assume that baryons accrete into a halo at the cosmic baryon fraction and that a large fraction of the cooled baryons are subsequently ejected by feedback. Most semi-analytic models (SAMs) have focused on gas ejection through winds and we will follow suit here. Owing to the uncertainty in the physical details, we use the SAM developed by Lu et al. (2011) that both incorporates sufficiently general parametrizations of galaxy formation mechanisms so that its results will apply to a large class of SAMs and uses advanced Markov chain Monte Carlo (MCMC) techniques to effectively explore the high-dimensional parameter space. Similar approaches have been adopted by Kampakoglou, Trotta & Silk (2008), Henriques et al. (2009, 2013), and Mutch et al. (2013). In addition, Gaussian process model emulators have also been adopted to explore the parameter space of SAMs by Bower et al. (2010) and Gómez et al. (2012). More generally, since the Bayesian inference approach places the model on a probabilistic footing, comparisons between models, probabilities within models, and the consistency of the models given the data (i.e. goodness of fit) can all be assessed quantitatively.

In Lu et al. (2012), we applied the Bayesian SAM to make model inferences conditional on the observed K-band luminosity function of local galaxies. We found that star formation and SN feedback are degenerate when one only uses the luminosity function to constrain the model. The model implied two modes for low-mass haloes: one where the star formation efficiency is increasingly lower for lower halo masses and one where the SN feedback is increasingly effective in reheating the cold gas. The dominating low star formation mode hides a large fraction of baryonic mass in the cold gas phase, thereby vastly overpredicting the H i mass function of galaxies. In this paper, we use both the K-band luminosity function and the H i mass function of galaxies to constrain a model family based on Lu et al. (2012), with additional feedback processes for ejecting baryons out of haloes, to see if there are interesting regions in the extended parameter space that can explain both data sets simultaneously. We choose the K-band luminosity function because it is less affected by dust extinction than are optical bands. The cold gas component, although a relatively small fraction of the total baryon gas (e.g. Bell et al. 2003b; Keres, Yun & Young 2003; Giovanelli et al. 2005; Zwaan et al. 2005), depends crucially on gas cooling and accretion, star formation, and feedback (Lu et al. 2012; Wang, Weinmann & Neistein 2012) and thereby complements the information from the stellar component. We first use our SAM with 17 free parameters to obtain the posterior distribution from the constraining data. Our philosophy in this paper is to take a broad prior range for the parameters, a range that is sometimes broader than a conventional interpretation of observations and theory would allow, but one that is broad enough to contain values assumed in published SAMs. We find that using these very broad ranges for the priors of the model parameters, our posterior distribution recovers many parameter combinations of previously published SAMs, but the model is still unable to match both of the constraining data sets simultaneously. If we restrict our prior ranges to those suggested by other observations (and theory), then the mismatch with the constraining observation becomes much worse, as we discuss later. We then analyse the posterior distribution of the key parameters that govern star formation and feedback in an attempt to gain insight into the underlying physical processes and to investigate the shortcomings of the model. Finally, we use the constrained posterior to make predictions and compare them with existing data to further test the models.

The paper is organized as follows. Section 2 describes the extended model based on Lu et al. (2011, 2012). We present the data, its error model, and the definition of the likelihood function for the Bayesian inference in Section 3. Section 4 presents the results of our model inference using the constrained posterior distribution of model parameters (Section 4.1) and model predictions (Section 4.2). Finally, we discuss our results in Section 5. Throughout the paper, we assume a Λ cold dark matter (ΛCDM) cosmology with Ω_{M, 0} = 0.26, Ω_{B, 0} = 0.044, h = 0.71, n = 0.96, and σ₈ = 0.79, which are consistent with the 5-year Wilkinson Microwave Anisotropy Probe (WMAP5) data (Dunkley et al. 2009; Komatsu et al. 2009).

2 THE SEMI-ANALYTIC MODEL

We adopt the SAM developed by Lu et al. (2011; see also Lu et al. 2012), which includes star formation, SN feedback, galaxy mergers, and active galactic nucleus (AGN) feedback. The processes are parametrized with standard analytic functions so that the model encompasses a large range of previously used prescriptions. We only describe the parametrizations relevant to SN feedback and star formation in this section; see Lu et al. (2011) for a comprehensive description of the other details.

Our model assumes that quiescent star formation occurs in a high-density, exponentially distributed cold gas disc. The characteristic radius for the exponential disc, r_disc, is related to the virial radius and the spin parameter of the host halo as follows: |$r_{\rm disc}= \lambda r_{\rm vir}/\sqrt{2}$|⁠, where λ is the spin parameter of the halo (Mo, Mao & White 1998). We choose a single value λ = 0.035 for all haloes; this value is roughly the average value measured in cosmological N-body simulations (e.g. Macciò et al. 2007). The fiducial size of the disc is then |$r_{\rm disc,0}=0.035r_{\rm vir}/\sqrt{2}$|⁠. Observationally, gas forms stars when its surface density exceeds a threshold, Σ_SF; the observed threshold surface density is ∼10 M_⊙ pc⁻² (e.g. Kennicutt 1998; Kennicutt et al. 2007; Bigiel et al. 2008). For an exponential disc, the cold gas mass with a surface density higher than the threshold surface density Σ_SF is

\begin{equation} m_{\rm sf}=m_{\rm cold}\left[1-\left(1+\ln {m_{\rm cold} \over 2\pi f_{\rm SF} r_{\rm disc,0}^2}\right) {2\pi f_{\rm SF} r_{\rm disc,0}^2 \over m_{\rm cold}}\right], \end{equation}

(1)

where m_cold is the total cold gas mass in the galaxy, and f_SF is a model parameter defined as

\begin{equation} f_{\rm SF}=\left({\lambda \over 0.035}\right)^2 \Sigma _{\rm SF}. \end{equation}

(2)

The model parameter f_SF describes both the size of the disc and the threshold surface density of the cold gas. Thus, the parameters λ and Σ_SF are intrinsically degenerate through the combined parameter f_SF. If f_SF is significantly lower than 10 M_⊙ pc⁻², then either the cold gas disc is more compact than the fiducial size, or the threshold surface density is required to be lower than 10 M_⊙ pc⁻². We assign a generous prior for the parameter f_SF to accommodate the broad uncertainties in the implementation of star formation in the model (see Table 1), and we test the sensitivity of our results on the choice of the prior, which will be discussed in Section 4.1.

We assume that the SFR, ϕ, is proportional to the mass of star-forming gas and inversely proportional to the dynamical time-scale of the disc, τ_disc = r_disc/V_vir:

\begin{equation} \phi =\epsilon _{\rm sf}{m_{\rm sf} \over \tau _{\rm disc}}, \end{equation}

(3)

where ϵ_sf is an overall star formation efficiency factor. To further generalize the dependence of the SFR on halo properties, the model further assumes that ϵ_sf has a broken power-law dependence on the circular velocity of the host halo:

\begin{equation} \epsilon _{\rm sf} = \left\lbrace \begin{array}{ll}\alpha _{\rm SF} &\quad {V_{\rm vir}\ge V_{\rm SF}}, \\ \alpha _{\rm SF}\left({V_{\rm vir}\over V_{\rm SF}}\right)^{\beta _{\rm SF}} &\quad {V_{\rm vir}< V_{\rm SF}}, \end{array} \right. \end{equation}

(4)

where α_SF, β_SF, and V_SF characterize the star formation efficiency. For haloes with a circular velocity higher than the critical velocity scale, V_SF, the star formation efficiency is a constant, α_SF; for haloes with a circular velocity lower than V_SF, the star formation efficiency follows a power-law function of halo circular velocity defined by α_SF and β_SF. The pivot point V_SF in the parametrization is needed to avoid unrealistically high star formation efficiencies in high V_vir haloes when β_SF is positive. However, we note that it introduces an intrinsic degeneracy into the parametrization. For example, when β_SF = 0, V_SF does not have any impact on the model. This affects our study, as we discuss in Section 4.1. We note that allowing β_SF to have non-zero values conflicts with some observational results, which find that the SFR is largely determined by the available supply of cold dense gas and not on the mass (or circular velocity) of the host halo (Schmidt 1959; Kennicutt 1989).

Star formation in our SAM can both consume cold gas by forming stars and affect the cold gas supply for subsequent star formation through feedback. In general, star formation feedback can influence the cold gas in three different ways: first, feedback can heat up the cold gas in the disc so that it joins the hot halo; second, feedback can directly eject cold gas out of the gravitational potential of the host halo; and third, feedback can expel some of the hot halo gas out the host halo. Our model attempts to capture all these processes. In our previous studies (Lu et al. 2011), we did not explore the dependence on cold gas ejection owing to its strong degeneracy with other mechanisms when the luminosity function or stellar mass function alone is used as a constraint. The addition of the H i mass function breaks this degeneracy, and we will see that cold gas ejection is required by the H i data. Our previous papers enforced conservation of the mechanical energy produced by SN explosions for the chosen stellar initial mass function (IMF) when modelling the mass loading in galaxy winds. In this paper, we allow the possibility of other energy sources, e.g. radiatively driven winds.

We parametrize the mass-loss rate of the cold gas to be proportional to the SFR, ϕ, as

\begin{equation} \dot{m}_{\rm out}=\alpha _{\rm ld} \phi , \end{equation}

(5)

where the coefficient α_ld is the loading factor of the outflow. Because the quantity of energy or momentum produced during stellar evolution is finite, the mass loading depends on the initial wind velocity and the physics of the coupling with the multiphase ISM. Two types of mass-loading models have been extensively explored in the literature: (1) an energy-driven wind that couples a fixed fraction of the feedback energy with the wind, which motivates |$\alpha _{\rm ld}\propto v_{\rm w}^{-2}$| (e.g. Okamoto et al. 2010); (2) a momentum-driven wind that couples a fixed fraction of momentum (e.g. in a radiation field) with the wind, which motivates |$\alpha _{\rm ld}\propto v_{\rm w}^{-1}$| (e.g. Oppenheimer & Dave 2006; Oppenheimer et al. 2010; Hopkins et al. 2013). Attempts to observe v_w and infer its dependence on galaxy properties (e.g. Martin 1999) have been inconclusive. Some theoretical models assume a constant wind velocity (e.g. Croton et al. 2006), while others assume a halo-mass-dependent velocity (e.g. Puchwein et al. 2012). To parametrize the various possibilities, we assume the loading factor to be a power-law function of halo circular velocity, V_vir,

\begin{equation} \alpha _{\rm ld}=\alpha _{\rm LD} \left({V_{\rm vir}\over V_0} \right)^{-\beta _{\rm LD}}, \end{equation}

(6)

where the power index β_LD and the normalization α_LD are model parameters, and V₀ is an arbitrary velocity scale that we fix at 220 km s⁻¹.

In this model, a given value of β_LD does not determine a particular wind launching mechanism (energy driven or momentum driven). For example, β_LD = 0 implies independence of V_vir, but this does not affect the α_ld–v_w relation predicted by a particular wind driving mechanism, because the v_w–V_vir relation is undetermined. Only in the special case when v_w ∝ V_vir is the value of β_LD related to the wind driving mechanism: β_LD = 2 for an energy-driven wind and β_LD = 1 for a momentum-driven wind. Because of all the uncertainties, we assign a large prior range for β_LD, allowing β_LD to take any value between 0 and 8 in our inference.

After the gas is removed from the galaxy, it can be trapped within the potential well of the dark matter halo or ejected from the halo into the intergalactic medium (IGM). We assume that the fraction of removed gas that is ejected from the halo is

\begin{equation} f_{\rm ej}=\left[1+\left({ V_{\rm vir}\over V_{\rm EJ}} \right)^{\beta _{\rm EJ}}\right]^{-1}, \end{equation}

(7)

where V_EJ and β_EJ are free parameters. Thus, for a halo with a circular velocity lower than V_EJ, most of the outflow gas is ejected from the halo, while for haloes with circular velocities much larger than V_EJ, the ejected fraction follows a power-law function of the halo circular velocity. Our model also includes an outflow of hot halo gas. The strength of the effect is controlled by the parameter α_SN, which describes the fraction of the total kinetic energy released by Type II SN that is needed to power all the feedback processes in the model. The excess energy over the amount used for reheating and ejecting the cold gas is used to power an outflow of hot halo gas. The amount of halo gas mass expelled is written as

\begin{equation} \Delta m_{\rm wind} = \epsilon _{\rm W} \left\lbrace \alpha _{\rm SN}\frac{V_{\rm SN}^2}{V_{\rm esc}^2} - \alpha _{\rm ld} \left[\left(\frac{V_{\rm vir}}{V_{\rm esc}}\right)^2+f_{\rm ej}\right]\right\rbrace \phi \Delta t, \end{equation}

(8)

where V_esc is the escape velocity of the halo. For a Navarro–Frenk–White (NFW) halo with a concentration c (Navarro, Frenk & White 1996), |$V_{\rm esc}=V_{\rm vir}\sqrt{ 2c \over \ln (1+c) -{c \over 1+c}}$| (Puchwein & Springel 2013).

The model assumes that the ejected gas can re-collapse into the halo at later times as hot halo gas. As in Springel et al. (2001) and De Lucia & Blaizot (2007), the rate of re-infall of ejected gas is given by

\begin{equation} \dot{m}_{\rm re\hbox{-}infall}=f_{\rm RI} \left( { M_{\rm ej} \over \tau _{\rm dyn}} \right). \end{equation}

(9)

Here, f_RI is a free parameter, M_ej is the total mass of gas in the ‘ejected’ reservoir, including both ejected gas from the cold gas disc and expelled gas from the hot halo, and τ_dyn = r_vir/V_vir is the dynamical time of the halo.

To convert the stellar mass to stellar light, we apply the stellar population synthesis (SPS) model of Bruzual (2007, BC07), which includes an improved treatment of thermally pulsing asymptotic giant branch (AGB) stars. The SPS model is used to predict the K-band luminosity of a galaxy from its star formation history. To implement the SPS model, we use a Chabrier IMF (Chabrier 2003), and create a look-up table for the luminosities on a grid with 220 values of age evenly spaced from 0 to 15 Gyr, and with Z = 0.0001, 0.0004, 0.004, 0.008, 0.02, and 0.05.

Finally, our SAM uses Monte-Carlo-generated halo merger trees tuned to match the conditional mass functions found in N-body simulations following Parkinson, Cole & Helly (2008) with final masses in the range from 10⁹ to 10¹⁵ h⁻¹ M_⊙ (see Lu et al. 2011, for additional details). Since haloes and their merger trees are randomly sampled from the halo mass function, a model prediction based on a finite merger tree sample suffers from sampling variance. We, therefore, choose the number of halo merger trees in each mass bin such that the standard deviation in the predictions, namely the K-band luminosity function and the H i mass function, are at least two times smaller than the error in the observational data in all the bins. Specifically, we use 1000 merger trees for haloes with present masses in the range 10¹¹–10^12.5 h⁻¹ M_⊙, 1500 in 10^12.5–10^13.5 h⁻¹ M_⊙, 400 in 10⁹–10¹⁰ h⁻¹ M_⊙, 400 in 10¹⁰–10¹¹ h⁻¹ M_⊙, and 100 in 10^13.5–10¹⁵ h⁻¹ M_⊙. The logarithmic halo masses are evenly distributed in each range. Since massive haloes are rare in the assumed cosmology, their contribution to scatter in the stellar mass function is negligible. We choose a mass resolution for the merger trees that varies with the final halo mass. For haloes with final masses smaller than 10¹⁰ h⁻¹ M_⊙, the mass resolution is 10^8.5 h⁻¹ M_⊙; for haloes with final mass between 10¹⁰ and 10¹² h⁻¹ M_⊙, the mass resolution is 10^9.2 h⁻¹ M_⊙; for haloes with final mass between 10¹² and 10¹⁴ h⁻¹ M_⊙, the mass resolution is 10^10.2 h⁻¹ M_⊙; and for haloes with final masses larger than 10¹⁴ h⁻¹ M_⊙, it is 10¹¹ h⁻¹ M_⊙. All merger trees are sampled at 100 redshifts equally spaced in log (1 + z) from z = 7 to 0. To make predictions for the total galaxy population we weight each predicted galaxy by the halo mass function of Sheth, Mo & Tormen (2001). The merger tree set we adopt allows us to run each model rapidly with sufficient accuracy so that we can explore the parameter space within a reasonable amount of time. After we obtained the posterior, we checked our model predictions using a larger number of trees and higher time and mass resolutions. The deviations of the results from those obtained with the smaller tree sample and lower resolutions described above were found to be well within the 1σ range of the posterior predictive distributions, demonstrating that the number of trees and mass resolution that we adopt are adequate for our purposes.

Table 1.

Open in new tab

Model parameters.

#	Parameter	Meaning	Prior
1	log M_CC ( M_⊙)	Cooling cut-off halo mass	[1.5, 4.5]
2	log α_SF	Star formation efficiency power-law amplitude	[−3, 0]
3	β_SF	Star formation efficiency power-law index	[−1, 12]
4	log V_SF (km s⁻¹)	Star formation law turn-over halo circular velocity	[1.5, 3.0]
5	log f_SF	Star formation threshold gas surface density	[0, 1.2]
6	log α_SN	SN feedback energy fraction	[−3, 1]
7	log α_LD	SN feedback reheating power-law amplitude	[−3, 2]
8	β_LD	SN feedback reheating power-law index	[0, 14]
9	log V_EJ	SN feedback ejection pivot halo circular velocity	[1.6, 3.0]
10	β_EJ	SN feedback ejection power-law index	[0, 10]
11	ϵ_EJ	SN feedback ejection escaping velocity factor	[1, 5]
12	log ϵ_W	Fraction of surplus SN feedback energy used for powering wind	[−3, 0]
13	log f_RI	Fraction of re-infall ejected hot gas	[−2, 0]
14	log f_DF	Merging time-scale in dynamical friction time-scale	[0, 2]
15	log α_SB	Merger triggered starburst efficiency power-law amplitude	[−2, 0]
16	β_SB	Merger triggered starburst efficiency power-law index	[0, 2]
17	\|$\arctan (\alpha _{\rm IN})$\|	K-band luminosity function faint-end incompleteness	[0, 0.177]

#	Parameter	Meaning	Prior
1	log M_CC ( M_⊙)	Cooling cut-off halo mass	[1.5, 4.5]
2	log α_SF	Star formation efficiency power-law amplitude	[−3, 0]
3	β_SF	Star formation efficiency power-law index	[−1, 12]
4	log V_SF (km s⁻¹)	Star formation law turn-over halo circular velocity	[1.5, 3.0]
5	log f_SF	Star formation threshold gas surface density	[0, 1.2]
6	log α_SN	SN feedback energy fraction	[−3, 1]
7	log α_LD	SN feedback reheating power-law amplitude	[−3, 2]
8	β_LD	SN feedback reheating power-law index	[0, 14]
9	log V_EJ	SN feedback ejection pivot halo circular velocity	[1.6, 3.0]
10	β_EJ	SN feedback ejection power-law index	[0, 10]
11	ϵ_EJ	SN feedback ejection escaping velocity factor	[1, 5]
12	log ϵ_W	Fraction of surplus SN feedback energy used for powering wind	[−3, 0]
13	log f_RI	Fraction of re-infall ejected hot gas	[−2, 0]
14	log f_DF	Merging time-scale in dynamical friction time-scale	[0, 2]
15	log α_SB	Merger triggered starburst efficiency power-law amplitude	[−2, 0]
16	β_SB	Merger triggered starburst efficiency power-law index	[0, 2]
17	\|$\arctan (\alpha _{\rm IN})$\|	K-band luminosity function faint-end incompleteness	[0, 0.177]

Table 1.

Open in new tab

Model parameters.

#	Parameter	Meaning	Prior
1	log M_CC ( M_⊙)	Cooling cut-off halo mass	[1.5, 4.5]
2	log α_SF	Star formation efficiency power-law amplitude	[−3, 0]
3	β_SF	Star formation efficiency power-law index	[−1, 12]
4	log V_SF (km s⁻¹)	Star formation law turn-over halo circular velocity	[1.5, 3.0]
5	log f_SF	Star formation threshold gas surface density	[0, 1.2]
6	log α_SN	SN feedback energy fraction	[−3, 1]
7	log α_LD	SN feedback reheating power-law amplitude	[−3, 2]
8	β_LD	SN feedback reheating power-law index	[0, 14]
9	log V_EJ	SN feedback ejection pivot halo circular velocity	[1.6, 3.0]
10	β_EJ	SN feedback ejection power-law index	[0, 10]
11	ϵ_EJ	SN feedback ejection escaping velocity factor	[1, 5]
12	log ϵ_W	Fraction of surplus SN feedback energy used for powering wind	[−3, 0]
13	log f_RI	Fraction of re-infall ejected hot gas	[−2, 0]
14	log f_DF	Merging time-scale in dynamical friction time-scale	[0, 2]
15	log α_SB	Merger triggered starburst efficiency power-law amplitude	[−2, 0]
16	β_SB	Merger triggered starburst efficiency power-law index	[0, 2]
17	\|$\arctan (\alpha _{\rm IN})$\|	K-band luminosity function faint-end incompleteness	[0, 0.177]

#	Parameter	Meaning	Prior
1	log M_CC ( M_⊙)	Cooling cut-off halo mass	[1.5, 4.5]
2	log α_SF	Star formation efficiency power-law amplitude	[−3, 0]
3	β_SF	Star formation efficiency power-law index	[−1, 12]
4	log V_SF (km s⁻¹)	Star formation law turn-over halo circular velocity	[1.5, 3.0]
5	log f_SF	Star formation threshold gas surface density	[0, 1.2]
6	log α_SN	SN feedback energy fraction	[−3, 1]
7	log α_LD	SN feedback reheating power-law amplitude	[−3, 2]
8	β_LD	SN feedback reheating power-law index	[0, 14]
9	log V_EJ	SN feedback ejection pivot halo circular velocity	[1.6, 3.0]
10	β_EJ	SN feedback ejection power-law index	[0, 10]
11	ϵ_EJ	SN feedback ejection escaping velocity factor	[1, 5]
12	log ϵ_W	Fraction of surplus SN feedback energy used for powering wind	[−3, 0]
13	log f_RI	Fraction of re-infall ejected hot gas	[−2, 0]
14	log f_DF	Merging time-scale in dynamical friction time-scale	[0, 2]
15	log α_SB	Merger triggered starburst efficiency power-law amplitude	[−2, 0]
16	β_SB	Merger triggered starburst efficiency power-law index	[0, 2]
17	\|$\arctan (\alpha _{\rm IN})$\|	K-band luminosity function faint-end incompleteness	[0, 0.177]

3 DATA AND LIKELIHOOD

Our SAM inference is conditional on both the K-band luminosity function and H i mass function of galaxies to constrain models. We use the K-band luminosity function from Bell et al. (2003a); the details of the data and how it is used to constrain the model can be found in Lu et al. (2012). Different from the luminosity function, whose errors can be safely approximated by Poisson errors, the H i mass function contains errors that are correlated across different bins. Here, we give a detailed description for our use of the H i mass function from Zwaan et al. (2005), derived from the HI Parkes All Sky Survey (HIPASS) Bright Galaxy Catalogue (hereafter HIPASS BGC; see Koribalski et al. 2004), and show how we derive a full covariance matrix for the likelihood function utilized in our Bayesian analysis.

3.1 Uncertainties and covariance in the observed H i mass function

A galaxy's H i mass is estimated from its 21-cm flux density and distance as

\begin{equation} M_{\rm H\,{\small i}}=2.36 \times 10^5 d^2F_{\rm H\,{\small i}} \,{\rm M_{{\odot }}}, \end{equation}

(10)

where d is the distance of the source in Mpc, and |$F_{\rm H\,{\small i}}$| is the integrated 21-cm flux density in units of Jansky km⁻¹ s⁻¹. The value of |$F_{\rm H\,{\small i}}$| is provided for each galaxy in the catalogue, along with its uncertainty, |$\sigma (F_{\rm H\,{\small i}})$|⁠. We use the Hubble distance,

\begin{equation} d=v_{\rm LG}/H_0, \end{equation}

(11)

where H₀ is the Hubble constant in units of km s⁻¹ Mpc⁻¹, v_LG is the recession velocity of the source, v_sys, corrected to the Local Group rest frame:

\begin{equation} v_{\rm LG}=v_{\rm sys} + 300 \sin l \cos b, \end{equation}

(12)

where l and b are the Galactic longitude and latitude of the source, respectively. The observational error in v_sys is provided in the source catalogue and is used in our analysis.

The detectability of a H i galaxy is affected by both its peak flux density, S_p, and its linewidth, w₂₀, defined to be the wavelength difference between the two points where the flux density is 20 per cent of S_p. Following Koribalski et al. (2004), we estimate the uncertainty in the peak flux density as

\begin{equation} \sigma (S_{\rm p})^2= {\rm rms}^2 + (0.05 S_{\rm p})^2, \end{equation}

(13)

where the rms = 13 mJy. The uncertainty in w₂₀ is estimated as σ(w₂₀) ≈ 3σ(v_sys).

Using these values of v_LG, |$F_{\rm H\,{\small i}}$|⁠, S_p, and w₂₀ and their standard deviations, we randomly generate 1000 replicas containing 1000 galaxies, each assuming independent Gaussian distributions. For each of the replicas, we adopt the same procedure used by Zwaan et al. (2005) to obtain the H i mass function. Specifically, we apply the two-dimensional, stepwise maximum likelihood method proposed by Zwaan et al. (2003) to find the maximum likelihood estimator of the mass function. These ensembles yield the average mass function and the corresponding covariance matrix, |${\boldsymbol {\Sigma}}_{\rm obs}$|⁠. This covariance matrix includes both H i-measurement variance and the sampling variance owing to the finite survey volume. We will assume that both types of variance are independent and use this covariance matrix to specify the likelihood function in Section 3.3.

3.2 Corrections for molecular fraction and incompleteness

Our SAM predicts a galaxy's total cold mass, M_cold, which includes atomic H i gas, molecular H₂ gas, as well as heavier elements. We make the following assumptions to predict the H i mass function. To start, we write

\begin{equation} M_{\rm cold}={ M_{\rm H\,{\small i}} + M_{\rm H_2} \over \beta }, \end{equation}

(14)

where β ≈ 0.74 is used to correct for the mass in helium (He) and in the small fraction of heavier elements in a cosmic gas. Given the molecular to atomic hydrogen mass ratio, |$\eta =M_{\rm H_2}/M_{\rm H\,{\small i}}$|⁠, we have

\begin{equation} M_{\rm H\,{\small i}}={\beta M_{\rm cold} \over 1+\eta }. \end{equation}

(15)

The value of η depends on galaxy properties. Following Obreschkow et al. (2009), we assume that η depends only on the total mass of cold gas:

\begin{equation} \log \eta = q + k \log \left({M_{\rm cold} \over 10^9 \,h^{-2} \,{\rm M_{{\odot }}}} \right) + \sigma , \end{equation}

(16)

where |$q=-0.51^{+0.03}_{-0.04}$|⁠, |$k=-0.24^{+0.05}_{-0.05}$|⁠, and σ = 0.30. The H i mass derived from the observed 21-cm flux is only a fraction of the total H i mass owing to self-absorption. Although the details of self-absorption depend on inclination, the average correction expected for a sample of randomly oriented galaxies is about 15 per cent (Zwaan et al. 2003). This gives a relation between the predicted, observed H i mass, |$M^{\rm p}_{\rm H\,{\small i}}$|⁠, and the real H i mass |$M_{\rm H\,{\small i}}$|⁠:

\begin{equation} M^{\rm p}_{\rm H\,{\small i}}=(1-f_{\rm H\,{\small i}})M_{\rm H\,{\small i}}=\beta {1-f_{\rm H\,{\small i}} \over 1+\eta } M_{\rm cold}. \end{equation}

(17)

We take the self-absorption correction factor, |$f_{\rm H\,{\small i}}$|⁠, to be a random number between 0 and 0.3.

The detection of H i galaxies is also limited by a minimal velocity width w_m, below which galaxy signals cannot be reliably distinguished from radio-frequency interference. Thus, galaxies with low inclinations and low H i mass, which thereby have low velocity widths, are undersampled. In the HIPASS BGC, this effect is specified in terms of the velocity width at 50 per cent of the peak flux:

\begin{equation} w_{50}=(w_0^2 \sin ^2 i + w_{\rm t}^2)^{1/2} + w_{\rm i}, \end{equation}

(18)

where i is the inclination, w₀ is the intrinsic velocity spread owing to rotation, w_t is the velocity width owing to turbulence, and w_i is the broadening owing to the instrument. Following Zwaan et al. (2003), we set w_t = 6 km s⁻¹ and w_i = 2.3 km s⁻¹. Galaxies included in the HIPASS BGC have w₅₀ > w_m = 26.4 km s⁻¹, so the minimal reliably detected inclination is

\begin{equation} i_{\rm m} = \arcsin \left[{(w_{\rm m}-w_{\rm i})^2 - w_{\rm t}^2 \over w_0^2} \right]^{1/2}. \end{equation}

(19)

The fraction of galaxies potentially missed at a given w₀ is then

\begin{eqnarray} \zeta =\int _0^{i_{\rm m}} \sin i\, {\rm d}i = 1-\cos i_{\rm m} = 1-\left[ 1- {(w_{\rm m}-w_{\rm i})^2 - w_{\rm t}^2 \over w_0^2} \right]^{1/2}. \nonumber\\ \end{eqnarray}

(20)

Zwaan et al. (2003) used the H i Tully–Fisher relation from Lang et al. (2003) to infer |$w_0=0.35 M_{\rm H\,{\small i}}^{0.3}$|⁠, ignoring the scatter. This implies that ζ ≈ 10 per cent for |$M_{\rm H\,{\small i}}=2\times 10^7 \,{\rm M_{{\odot }}}$| and ζ ≈ 5 per cent for |$M_{\rm H\,{\small i}}=5\times 10^7 \,{\rm M_{{\odot }}}$|⁠. We include this effect in our model predictions for the H i mass function as

\begin{equation} \Phi ^{\rm p}\left(M^{\rm p}_{\rm H\,{\small i}}\right) = \Phi \left(M^{\rm p}_{\rm H\,{\small i}}\right) (1-\zeta ), \end{equation}

(21)

with ζ given by equation (20).

3.3 The likelihood function

For the K-band luminosity function, we use the same likelihood function as described in Lu et al. (2012). For the H i mass function, we construct a likelihood function that independently represents the H i-measurement error and the sample variance. To start, we neglect the covariance between measurement error and sampling variance. Such covariance may be important for small sample sizes but is unimportant for the current sample where the variance is measurement dominated (Papastergis, private communication). The uncertainties in the process leading to the H i mass estimate for each galaxy (Sections 3.1 and 3.2) induce correlations between the H i mass bins. Ideally, the error model describing this covariance would be used to predict the observed data. However, our galaxy formation model is not capable of predicting all the necessary observables (e.g. w₂₀ and S_p). Instead, we proceed by making two additional assumptions: (1) the likelihood function is a multidimensional normal distribution; and (2) the sampling is a pure Poisson process. Then, the errors in the simulated data presented in Section 3.1 can be modelled as follows. We define the equivalent volume, |${\boldsymbol V}_{\rm eq}$|⁠, implied by the H i measurement process including all the selection effects described in Section 3.2. This equivalent volume can be estimated from the data as |${\boldsymbol V}_{\rm eq}=N_{\rm obs}/{\boldsymbol {\Phi} }_{\rm obs}$|⁠. The Poisson variance in each modelled bin is then |${\boldsymbol {\Phi} }_{\rm mod}/{\boldsymbol V}_{\rm eq}$| and only affects the diagonal terms of the covariance matrix. We, therefore, write the diagonal terms as

\begin{equation} \sigma _{\rm mod}^{2} = {\boldsymbol\sigma} _{\rm mst}^{2} + \left(\frac{{\boldsymbol {\Phi}} _{\rm mod}}{{\boldsymbol V}_{\rm eq}} \right), \end{equation}

(22)

where |${\boldsymbol \sigma} _{\rm mst}^{2}$| is some unknown measurement variance. Under the same assumptions, the diagonal elements of the simulated covariance matrix from Section 3.1 are

\begin{equation} {\boldsymbol \sigma} _{\rm obs}^{2} = {\boldsymbol \sigma} _{\rm mst}^{2} + \left(\frac{{\boldsymbol {\Phi}} _{\rm obs}}{{\boldsymbol V}_{\rm eq}} \right). \end{equation}

(23)

We now use equations (22) and (23) to eliminate |${\boldsymbol \sigma} _{\rm mst}^{2}$|⁠:

\begin{equation} \sigma _{\rm mod}^2 = {\boldsymbol \sigma} _{\rm obs}^2 + {{\boldsymbol {\Phi}} _{\rm mod} - {\boldsymbol {\Phi} }_{\rm obs} \over {\boldsymbol V}_{\rm eq}}. \end{equation}

(24)

This replaces the variance in the diagonal terms of the simulated covariance matrix to obtain the predicted covariance |$\boldsymbol {\Sigma }_{\rm pred}$|⁠. Thus, the likelihood function is written as

\begin{eqnarray} &&\!\!\! L({\boldsymbol {\Phi} }_{\rm obs}| \boldsymbol {\Theta }) = \frac{L_0}{(2 \pi )^{I/2} |{\rm det}(\boldsymbol {\Sigma }_{\rm pred})|^{1/2}} \nonumber \\ &&\quad \times \exp \left[-{1 \over 2} ({\boldsymbol {\Phi} }_{\rm obs} - {\boldsymbol {\Phi} }_{\rm mod})^{\rm T} \cdot {\boldsymbol {\Sigma} }^{-1}_{\rm pred} \cdot ({\boldsymbol {\Phi} }_{\rm obs} - {\boldsymbol {\Phi} }_{\rm mod})\right], \end{eqnarray}

(25)

where |$\boldsymbol {\Theta }$| denotes the model parameter vector, |${\boldsymbol {\Phi} }_{\rm obs}$| and |${\boldsymbol {\Phi} }_{\rm mod}$| are, respectively, the vectors of the observed and predicted H i mass functions over the H i mass bins, and I is the rank of the covariance matrix. We assume that the K-band luminosity function and the H i mass function are statistically independent and, therefore, the total likelihood is simply the product of the two:

\begin{equation} L(D|\boldsymbol {\Theta }) =L(D_{K}|\boldsymbol {\Theta })\times L(D_{\rm H\,{\small i}}|\boldsymbol {\Theta }), \end{equation}

(26)

where D_K and |$D_{\rm H\,{\small i}}$| denote the observed data of the K-band luminosity function and the H i mass function, respectively.

4 MODEL INFERENCES

4.1 Inference from the posterior distribution

We used a tempered version of the differential evolution algorithm (ter Braak 2006) implemented in the Bayesian Inference Engine (BIE; see Weinberg 2013) to sample the posterior density distribution with 256 chains. All simulations run for more than 7000 iterations and the convergence is monitored with the Gelman–Rubin |${\hat{R}}$| test (Gelman & Rubin 1992). We declare convergence when |${\hat{R}}\lesssim 1.2$|⁠. All simulations typically converge after 3000 iterations, and the converged states are used to sample the posterior distribution. After removing 12 outlier chains, we obtain ∼10⁶ chain states. As the autocorrelation length of the chains is typically ∼10, there are about 10⁵ independent chain states to sample the posterior distribution.

Fig. 1 shows the posterior predictions for the K-band luminosity function and H i mass function of local galaxies. We account for potential incompleteness in the faint end of the K-band luminosity function by marginalizing over the incompleteness parameter α_IN (see Lu et al. 2012, for details). The orange and yellow bands in each panel encompass 66 and 95 per cent of the posterior probability, respectively, and the red line is the median. To quantify the goodness of fit, we use posterior predictive check (PPC; Rubin 1984; Gelman, Meng & Stern 1996; Gelman et al. 2003). PPC compares the value of one or more discrepancy statistics for the observed data to a reference distribution based on the posterior (which is conditioned on the observed data). The reference distribution is obtained by generating predicted data from the posterior distribution. If a model fits the data well, the observed data should be likely under the posterior predictive distribution, and the discrepancy statistics of the observation should be likely under the distribution of the discrepancy statistics of the posterior models. On the other hand, large discrepancies between the observed data and the posterior predictive distribution indicate that the model performs poorly. Following our previous paper (Lu et al. 2012), we adopt the Bayesian p-value to quantify the goodness of fit, with p defined as

\begin{equation} p = {1 \over N} \sum _{i=1}^N I_{\chi _{i, {\rm mod}}^2 \ge \chi _{\rm obs}^2}, \end{equation}

(27)

where |$\chi _{i, {\rm mod}}^2$| and |$\chi _{\rm obs}^2$| are the discrepancy statistics for the model and the data. I_q is the indication function for the condition q, with I_q = 1 if q is true and I_q = 0 otherwise. In other words, the fraction of posterior model samples with |$\chi _{i, {\rm mod}}^2 \ge \chi _{\rm obs}^2$| is an estimate of p. A small value of the p reflects the implausibility of the data under the model (and, hence, the lack of fit of the model to the data) and, therefore, suggests problems with the model. For the data and posterior distribution shown in Fig. 1, the Bayesian p-value from the PPC for the H i mass function is |$p_{\rm H\,{\small i}}=0.25$|⁠, which indicates an adequate fit, but for the K-band luminosity function it is only p_K = 0.011. This is much lower than the value of p_K = 0.66 obtained in Lu et al. (2012), which used only the K-band luminosity function as the data constraint. Hence, the model constrained using both the H i mass function and the K-band luminosity function does not fit the observational data well. In the remainder of this paper we will explore the causes for this discrepancy and their implication for models of galaxy formation generally.

Figure 1.

The Bayesian posterior predictions for the K-band luminosity function and the H i mass function. The black solid lines with error bars are the observational data. The orange and yellow bands encompass the 66 and 95 per cent credible ranges of the predictions, respectively, while the red lines are the median values.

Open in new tab Download slide

Fig. 2 shows marginalized posterior distributions for 10 parameters characterizing star formation and feedback and we now describe those with marginal distributions that differ significantly from those in Lu et al. (2012). In the upper left-hand panel of Fig. 2, the vertical axis shows the model parameter V_EJ: the halo circular velocity below which all outflowing mass reheated from the disc must be ejected out of the host halo. The horizontal axis shows the parameter M_CC: the critical halo mass above which radiative cooling of the halo gas is shut off to mimic radio-mode AGN feedback (Bower et al. 2006; Croton et al. 2006). The mass M_CC is well constrained, with a value about 10^12.5 M_⊙. This result also agrees with the value obtained by Lu et al. (2012) using the K-band luminosity function alone as a data constraint, and with other existing SAMs that implement ‘halo quenching’ (e.g. Cattaneo et al. 2006; Somerville et al. 2008). Including the H i mass function as an additional constraint does not affect the inferred value of M_CC significantly because M_CC controls gas cooling and star formation in massive haloes while the H i mass function mainly constrains low-mass haloes, which contain most of the H i gas.

$The marginalized posterior distributions of 10 model parameters. Panel 1: AGN feedback shuts off gas cooling for haloes with Mvir > MCC and winds expel gas from haloes with circular velocities Vvir > VEJ. Panel 2: αSF parametrizes the normalization of star formation efficiency, which is a constant for haloes with a circular velocity higher than parameter VSF, and a power-law function of Vvir with βSF being the power-law index of the relation. Panel 3: fSF is the cold gas surface density threshold for star formation and αSF is as above. Panel 4: αLD parametrizes the loading factor of the SN feedback and βLD is the power index of the power-law dependence of the loading factor on the halo circular velocity. The blue lines denote the kinetic energy limit based on Type II SN for outflows in haloes with a circular velocity of 220 km s−1 for two IMFs: Salpeter (solid) and Chabrier (dashed). The red lines show the extreme parameter combination of the two mass-loading parameters for haloes with a circular velocity of 100 km s−1. Models in regions to the right of the limiting lines require more than 100 per cent of the total available energy for the cold gas ejection. The magenta dotted line denotes models with a loading factor equal to 4 for haloes with a circular velocity of 70 km s−1. Panel 5: the joint marginalized posterior distribution in the βLD–βSF plane. Panel 6: the joint marginalized posterior distribution in the βLD–VSF plane. Panel 7: αSN characterizes the fraction of SN energy that is required to power the ejection of hot halo gas. Panel 8: the joint marginalized posterior distribution in the βLD–log fDF plane, where fDF is a coefficient characterizing the satellite dynamical friction time-scale. The grey-scale (colour in on-line version) denotes the confidence levels as shown by the colour bar at the top of the figure. Symbols with different colours in the various panels denote the corresponding parameters adopted in some earlier SAMs, as indicated. They do not appear in every panel for all models because some parameters are absent in the parametrization of those models.$

Figure 2.

The marginalized posterior distributions of 10 model parameters. Panel 1: AGN feedback shuts off gas cooling for haloes with M_vir > M_CC and winds expel gas from haloes with circular velocities V_vir > V_EJ. Panel 2: α_SF parametrizes the normalization of star formation efficiency, which is a constant for haloes with a circular velocity higher than parameter V_SF, and a power-law function of V_vir with β_SF being the power-law index of the relation. Panel 3: f_SF is the cold gas surface density threshold for star formation and α_SF is as above. Panel 4: α_LD parametrizes the loading factor of the SN feedback and β_LD is the power index of the power-law dependence of the loading factor on the halo circular velocity. The blue lines denote the kinetic energy limit based on Type II SN for outflows in haloes with a circular velocity of 220 km s⁻¹ for two IMFs: Salpeter (solid) and Chabrier (dashed). The red lines show the extreme parameter combination of the two mass-loading parameters for haloes with a circular velocity of 100 km s⁻¹. Models in regions to the right of the limiting lines require more than 100 per cent of the total available energy for the cold gas ejection. The magenta dotted line denotes models with a loading factor equal to 4 for haloes with a circular velocity of 70 km s⁻¹. Panel 5: the joint marginalized posterior distribution in the β_LD–β_SF plane. Panel 6: the joint marginalized posterior distribution in the β_LD–V_SF plane. Panel 7: α_SN characterizes the fraction of SN energy that is required to power the ejection of hot halo gas. Panel 8: the joint marginalized posterior distribution in the β_LD–log f_DF plane, where f_DF is a coefficient characterizing the satellite dynamical friction time-scale. The grey-scale (colour in on-line version) denotes the confidence levels as shown by the colour bar at the top of the figure. Symbols with different colours in the various panels denote the corresponding parameters adopted in some earlier SAMs, as indicated. They do not appear in every panel for all models because some parameters are absent in the parametrization of those models.

Open in new tab Download slide

The value of V_EJ is also well constrained, with log V_EJ ≈ 2.3 (V_EJ ≈ 200 km s⁻¹), suggesting that all haloes with circular velocities lower than this are required to have efficient mass ejection to fit the data, and the same ejection process can be inefficient in haloes with higher circular velocities. There is also a non-negligible tail to even larger V_EJ; models in this tail require efficient mass ejection at even higher halo circular velocities. This result agrees with both the model of Somerville et al. (2008) and a recent analysis of the Munich model (Henriques et al. 2013): efficient gas ejection from haloes is required to match multiple data sets simultaneously. Without the H i constraint, the observed K-band luminosity function of galaxies can be reproduced by making the star formation efficiency a strong function of circular velocity, resulting in a massive cold gas component. However, this makes the amplitude of the corresponding predicted H i mass function an order of magnitude higher than that observed (Lu et al. 2012; Wang et al. 2012). In our model, the strong degeneracy between outflow and star formation efficiency is broken by the constraint provided by the H i mass function. Consequently, one requires efficient gas ejection from haloes to prevent a large reservoir of cold gas and only requires the star formation efficiency to be a weak function of halo circular velocity (see below). If we fix V_EJ = 0, we find that the model can no longer match the observational data; the best-fitting model, i.e. in this case the model with the maximum a posteriori (MAP) probability, has a likelihood that is almost 9000 000 times lower than our fiducial MAP. This suggests that efficient gas ejection is unavoidable in our model family.

4.1.1 Star formation

With the H i mass function as an additional constraint, the posterior probability of the star formation parameters is now constrained to a fairly narrow range (see Panel 2 of Fig. 2). The posterior of log (α_SF) peaks sharply at −1.8, implying that a galaxy converts about ≈1.6 per cent of its cold gas into stars in every dynamical time. This is consistent with the observed efficiency in galaxies (Elmegreen 1997; Silk 1997). The posterior for β_SF is now peaked near zero, much lower than the value of β_SF ≈ 10 obtained in Lu et al. (2012) constrained using the K-band luminosity function alone. This low value for β_SF means that the SFR is largely determined by the available supply of cold dense gas and does not depend on the mass (or circular velocity) of the host halo. This is consistent with observed star formation laws (Schmidt 1959; Kennicutt 1989), in which the SFR depends only on gas properties. However, there still remains a weaker mode – its posterior odds relative to the main mode is 2.5 times smaller – where α_SF ≈ 2 per cent and β_SF ≈ 2.5. We will show that the parameter β_SF is covariant with feedback parameters in another panel of the figure.

In Lu et al. (2012), the star formation efficiency parameter α_SF was found to be strongly covariant with the threshold of cold gas surface density for star formation, f_SF: a high star formation efficiency with a high threshold value was equivalent to a low star formation efficiency with a low threshold value. As expected, the addition of the H i mass function breaks this degeneracy. Panel 3 of Fig. 2 shows the marginalized posterior probability distribution in the log f_SF–log α_SF plane. Remember that in our parametrization, f_SF is equal to the surface gas density threshold for star formation, Σ_SF, under the assumption that r_disc = r_{disc, 0}, i.e. that the disc size is given by the model of Mo et al. (1998). The combined data sets require that the threshold surface density for star formation is small; we find Σ_SF ≈ 1 M_⊙ pc⁻². This value is about one order of magnitude lower than the observed threshold surface density (Kennicutt 1998; Bigiel et al. 2008; Leroy et al. 2008) and is also right against the prior bound, so even lower values may be preferred by the model, although such low values would be in even greater conflict with observations. In the model, a high gas surface density threshold would result in too much cold gas in low-mass haloes at large radii and too little star formation to power enough outflows. This problem was highlighted by Mo et al. (2005), who found that the standard threshold surface density for star formation would lead to an H i mass function that was too high, even if each low-mass halo only contained a cold gas disc at the threshold surface density.

Our SAM requires both powerful, efficient feedback to remove the gas from the halo and a threshold gas surface density that is much lower than the observed star formation threshold (and still fails to match the data). To explore explicitly the sensitivity to a higher surface density threshold, we constrained the prior of f_SF to range between 3 and 30 M_⊙ pc⁻². We find that this greatly reduces the goodness of fit and that the best fits are always obtained with f_SF at the lower boundary of the prior. The MAP using this restricted prior has a likelihood that is smaller than the MAP in the fiducial run by a factor of over 3000 000. If we further restrict f_SF to match the observed value of 10 M_⊙ pc⁻², the likelihood decreases by a further factor of over 60 000 000, and the p-value for the H i mass function becomes very small, |$p_{\rm H\,{\small i}}<0.005$|⁠, indicating an unacceptable fit. Thus, the current model family cannot explain the observed H i mass function using the observed gas density threshold for star formation.

4.1.2 Feedback

Panel 4 of Fig. 2 shows the marginalized posterior probability distribution of the two SN gas ejection parameters: α_LD and β_LD. For comparison, we also plot the combinations of these parameters used in a number of existing SAMs. These models include Cole et al. (1994), Kauffmann & Charlot (1998), Somerville & Primack (1999), Somerville et al. (2008), De Lucia, Kauffmann & White (2004), Kang et al. (2005), Menci et al. (2005), Cattaneo et al. (2006), and Guo et al. (2011). The parameter choices are converted into our parametrization, where the amplitude of the mass-loading factor is normalized at a halo circular velocity scale of 220 km s⁻¹. As expected, because of the similarity in our phenomenological parametrizations, most of the other SAMs are located in the high-probability region of our inferred posterior. These two parameters are covariant and populate three main modes: β_LD ≈ 0 that is equivalent to the Croton model (Croton et al. 2006), β_LD ≈ 2 that is similar to the Somerville model (Somerville et al. 2008) and a number of other models, and a large value of β_LD ≈ 6 (e.g. Cole et al. 1994).

Because we parametrize the outflow loading factor as a power-law function of halo circular velocity, |$\alpha _{\rm ld}=\alpha _{\rm LD}(V_{\rm vir}/220\,{\rm km\,s^{-1}})^{-\beta _{\rm LD}}$|⁠, a covariance represented by a straight line in the log α_LD–β_LD plane corresponds to a group of models that have a fixed loading factor for a specific circular velocity. The magenta dotted line in Panel 4 of Fig. 2 roughly captures the covariance of the two parameters and corresponds to a loading factor of 4 for haloes with a circular velocity of 70 km s⁻¹ or a mass of about 1.5 × 10¹¹ M_⊙ at the present time. The covariance suggests that these haloes may be critical for fitting the K-band luminosity function and the H i mass function of galaxies simultaneously.

Our fiducial model does not limit the total energy available for feedback but, of course, nature does impose limits. The total energy available for outflows powered by Type II SNe is

\begin{equation} E_{\rm SN, max}= \phi \eta _{\rm SN} \epsilon _{\rm SN}, \end{equation}

(28)

where ϵ_SN ∼ 10⁵¹ erg = 5 × 10⁷ M_⊙ km² s⁻² is the total energy output of a SN, and η_SN is the number of SNe expected from the formation of 1 M_⊙ of stars. Assuming that Type II SN progenitors are stars with masses 8 < m < 20 M_⊙, then η_SN = 5.88 × 10⁻³ M_⊙⁻¹ for a Salpeter IMF (Salpeter 1955), and η_SN = 9.76 × 10⁻³ M_⊙⁻¹ for a Chabrier IMF (Chabrier 2003). The total energy consumed by the outflow depends on the detailed processes and the final state of the outflow material. If the outflow is thermalized and settles into an equilibrium state in the gravitational potential of the host halo, the energy required is |${5 \over 4} m_{\rm rh} V_{\rm vir}^2$|⁠, where m_rh is the mass of the reheated gas (e.g. Kauffmann & Charlot 1998). Since only a fraction of the energy is expected to heat the gas:

\begin{equation} \eta _{\rm SN} \epsilon _{\rm SN} \ge {5 \over 4} \alpha _{\rm LD} \left({220\,{\rm km\,s}^{-1} \over V_{\rm vir}}\right)^{\beta _{\rm LD}} V_{\rm vir}^2. \end{equation}

(29)

Thus, for haloes with V_vir = 220 km s⁻¹, the maximum normalization of the loading factor, |$\alpha _{\rm LD, max}(V_{\rm vir}=220) =4.13 \times {\eta _{\rm SN} \over 5\times 10^{-3}}$|⁠, which is 4.8 for a Salpeter IMF and 8.0 for a Chabrier IMF. These values are plotted in the lower right-hand panel of Fig. 2 as the blue solid line and the blue dashed line, respectively.

If the outflow is ejected directly from the centre of the host halo, the energy required to escape is |${1 \over 2} m_{\rm ej} V_{\rm esc}^2$|⁠, where V_esc is the escape velocity. For a NFW halo (Navarro et al. 1996) with a concentration c = 10, the escape velocity from the centre is |$V_{\rm esc}^2\approx 13\,V_{\rm vir}^2$|⁠. Thus, to eject the gas from such a halo, the lower limit on the required energy is |$E={13 \over 2} m_{\rm LD} V_{\rm vir}^2$|⁠, where m_LD is the total gas mass loaded in the wind. The mass-loading factor in such haloes is, therefore, limited by the total energy provided by SN explosions, and the solid red line and the red dashed line in Panel 4 of Fig. 2 show these limits for a Salpeter and a Chabrier IMF, respectively. Models lying above such a limiting line need more energy than is available from haloes with circular velocities equal to or below 100 km s⁻¹. The posterior contours for all the modes with β < 4 are ∼0.3–0.6 dex to the left of but not far from the energy limits for haloes with V_vir ≤ 100 km s⁻¹, indicating that a large fraction, i.e. 25–50 per cent, of all the SN energy has to be put into the outflow in such galaxies.

The required high loading factor results from the need to simultaneously match the shallow low-end slope of both the K-band luminosity function and the H i mass function. The total baryonic mass that can be accreted into a halo is f_bm_vir, where f_b is the universal baryon fraction. Suppose a fraction f_c of the gas can cool and collapse on to the central galaxy. For low-mass haloes with V_vir ≲ 150 km s⁻¹, f_c is expected to be close to unity because of the high cooling efficiency (Thoul & Weinberg 1995; Lu & Mo 2007). If the expelled gas can fall back on to the galaxy and go through multiple cycles, the effective f_c can be larger than 1. However, to explain the shallow low-mass end of the H i mass function and the faint end slope of the K-band luminosity function, the cold gas and stellar mass have to be very low relative to f_bm_vir. Indeed, the stellar mass to halo mass ratio for haloes with masses ∼10¹¹ M_⊙ is only about 3 × 10⁻³ (Papastergis et al. 2012; Yang et al. 2012; Behroozi, Wechsler & Conroy 2013), and the cold gas mass to stellar mass ratio is no more than 10 (Kannappan 2004; Papastergis et al. 2012). Thus, the gas associated with low-mass haloes either is never accreted into the galaxy or has been expelled by galactic winds. In the latter case, if we conservatively assume that the cold gas mass is 10 times the stellar mass, then the total cold baryonic mass in a 10¹¹ M_⊙ halo is ∼0.033 of the halo mass. If we make the further conservative assumption that for every unit of stellar mass formed at early times only half of it remains in stars today owing to stellar mass loss, the total mass that was ever involved in star formation is ∼0.006 of the final halo mass. Based on these conservative assumptions, one finds that the mass-loading factor is (f_b − 0.033)/0.006 = 23, which is close to the upper limit obtained from the total kinetic energy expected from SN.

The posterior marginals (Fig. 2) imply that 25–50 per cent of the available SN energy, depending on the assumed IMF, is required to power the wind. This required high efficiency of energy loading poses a severe challenge for the standard SN feedback scenario (e.g. Dekel & Silk 1986; White & Zaritsky 1992). Indeed, detailed hydrodynamical simulations by Mac Low & Ferrara (1999) and Strickland & Stevens (2000) demonstrated that SN feedback is very inefficient in expelling mass because the onset of Rayleigh–Taylor instabilities severely limits the mass-loading efficiency of galactic winds. There are a couple of possible solutions to this dilemma. First, a top-heavy or bottom-light IMF would yield a larger η_SN. Our data constraints are from the local galaxy population, which by themselves do not distinguish different IMFs. Second, the ultraviolet (UV) photons from massive stars could effectively transfer radiation momenta to the gas via dust opacity, thereby producing momentum-driven winds (Murray et al. 2005). The total radiation energy available is about 100 times as high as the kinetic energy of SN explosions, vitiating the mechanical energy limit. The viability of this mechanism is uncertain: while some simulations show that it is effective (e.g. Hopkins et al. 2011), others find that the development of Rayleigh–Taylor-type instabilities can suppress the formation of outflows (Krumholz & Thompson 2012). Thus, although the required high efficiency may rule out mechanisms based on SN energy-driven winds alone, the viability including momentum-driven winds is still unclear.

The star formation feedback models implemented in SAMs are crude. In our model, the outflow reduces the total cold gas mass of the galaxy and the remaining cold gas is redistributed in an exponential profile at each time step. This is similar to some cosmological hydrodynamical simulations (e.g. Davé et al. 2013) where the effective viscosity of the disc gas, which could include actual gas viscosity and/or gravitational torques, causes the cold gas from the outer disc to move rapidly inwards to replenish the star-forming gas, and so star formation can continue near the disc centre even if the total amount of gas in the disc is reduced. However, since the gas viscosity in such simulations is artificial and the resolution is still limited, it is unclear whether the effect is physical or numerical. For example, if the effective viscosity of the disc gas were unimportant, much of the cold gas would stay in the outer region of the disc without forming stars. Some simulations of individual discs (Forbes, Krumholz & Burkert 2012) do find that star formation feedback is localized to regions where the star formation is active, producing a disc with a reduced cold gas density only in the central part of the disc, leaving the cold gas in the outskirts intact. Others confirm the cosmological simulation results (e.g. Hirschmann et al. 2013). To explore this in our SAM, we have considered a model in which star formation feedback is assumed to only affect the gas within the radius where star formation can happen, i.e. the cold density is above the threshold surface density. Without a lower limit on f_SF, the posterior distribution is nearly unchanged owing to the low surface density threshold for star formation and the high OFR required by the shallow low-mass end of the galaxy luminosity function. However, if the value of f_SF is forced to match the observed star formation threshold surface density of 10 M_⊙ pc⁻², then the model fails to fit the data (⁠|$p_{\rm H\,{\small i}}<0.005$|⁠), consistent with the results obtained by Mo et al. (2005).

As we remarked before, we find that there are three major modes in the log α_LD–β_LD plane located at β_LD ≈ 0, 2, and 6, as one can see in Panel 4 of Fig. 2. The relative posterior odds in each mode are 4:1:2, respectively. Hence, the β = 0 mode is favoured by the combined data sets in terms of marginalized posterior probability. Since our prior distribution is intentionally very broad, as we describe in Section 2, and this will quantitatively affect these odds. That said, the maximum likelihood is located in the β_LD ≈ 6 mode, and the highest likelihood values in the β_LD ≈ 2 and 0 modes are 30 and 1400 times lower, indicating that those two modes produce a worse fit to the data than the β_LD ≈ 6 mode. The PPC p-values calculated for the K-band luminosity function predicted by each of the modes also indicate that the β_LD ≈ 6 mode produces a better, but still inadequate, fit to the data, p_K = 0.016 for the β_LD ≈ 6 mode, and p_K < 0.005 for the other two modes. In Fig. 3, we show the K-band luminosity function produced by the three modes. The shaded blue region shows the posterior predictions of models with β_LD < 1.25, the magenta region shows the posterior predictions of models with 1.25 < β_LD < 4, and the green region shows the posterior predictions of models with β_LD > 4. Note that the three modes produce indistinguishable predictions for the H i mass function, as the predictive distribution of the H i mass function produced by the entire posterior shown in Fig. 1 is narrow. The predicted faint-end slope of the K-band luminosity function is steeper for a smaller value of β_LD, and the current data apparently favour the mode with the highest β_LD. This trend confirms recent hydro-dynamical simulation results of Puchwein & Springel (2013) and Davé et al. (2013). Moreover, this demonstrates that discrimination between these modes could be possible given improved data for the faint end of the luminosity function.

Figure 3.

The predicted K-band luminosity function for the three different modes present in the lower right-hand panel of Fig. 2 represented by different values of β_LD. The red line is the median prediction for the entire posterior.

Open in new tab Download slide

Even though the β_LD ≈ 0 mode has relatively lower likelihood, it has the highest marginal posterior probability. This implies that the modes with high values of β_LD are supported by proportionately smaller volumes in parameter space relative to the β_LD ≈ 0 mode. This could be investigated in detail using the distribution of marginal likelihood, but the required resolution to compute the volume in the high-dimensional parameter space is computationally prohibitive. Nonetheless, relying on the posterior marginals shown in Fig. 2, we believe that the complex posterior results from the fine-tuning necessary to simultaneously match the K-band luminosity and H i mass function. Specifically, Panel 5 indicates that β_SF correlates with β_LD. As β_SF increases, the increasingly smaller star formation efficiency in low circular velocity haloes leaves more cold gas in the disc. To simultaneously match the H i mass function, the model thus needs to increase the feedback efficiency in lower circular velocity haloes by increasing β_LD.

In addition, sometimes a particular parametrization can also introduce specific features into the posterior. Panel 6 of Fig. 2 shows an example; the parameter V_SF is not constrained when β_LD ≈ 0, but is strongly constrained when β_LD ≈ 6. When β_LD ≈ 0, β_SF is also about 0, meaning a constant star formation efficiency. In this case, as described in Section 2, V_SF no longer has any effect on the model predictions, because the parametrization for the star formation efficiency no longer has a transition scale that V_SF describes. Thus, the likelihood is the same no matter what value V_SF takes when β_LD ≈ 0.

Moreover, the high β_LD mode requires two unusual conditions. First, as shown in Panel 7 of Fig. 2, for β_LD ≈ 6, the parameter α_SN, which characterizes the fraction of SN kinetic energy used to eject the hot halo gas, is always required to be higher than one, implying that the required energy is higher than the kinetic energy provided by Type II SN. The reason is that β_LD ≈ 6 leads to a rapid decrease of cold gas ejection with increasing halo circular velocity, and a strong outflow of hot halo gas in relatively massive haloes is needed by the mode to fit the data. This conclusion agrees with the result of Mutch et al. (2013) that one needs more energy than that available from Type II SN to explain the evolution of the galaxy stellar mass function from z ≈ 0.8 to 0. Second, the β_LD ≈ 6 mode also requires a long time-scale, about 10 times as long as the fiducial dynamical friction time-scale, for satellite galaxies to orbit in the host halo before merging into central galaxies (see Panel 8 of Fig. 2). This can be understood as follows. When β_LD is large, the effect of star formation feedback is increasingly weaker for higher mass haloes. In this case, the merging time-scale for satellite galaxies is required to be long to prevent them from merging into the halo centre and hence significantly increasing the mass of the central galaxy. Since a longer merger time also leads to a higher satellite fraction, observational data of satellite fraction can tighten the constraint. We will come back to this topic in a future paper.

4.2 Posterior predictions

The model family generally requires efficient outflows from low-mass galaxies for its best fit, but even that is not enough to match both of the constraining observational data sets. We now use the posterior distribution to predict the OFR and other observables to investigate the failure of the model family. Since the marginalized posterior in the α_LD–β_LD plane is multimodal, we distinguish the modes according to their value of β_LD when making model predictions. We marginalize the posterior probability to generate predictions as described in Lu et al. (2012).

Fig. 4 shows the posterior prediction for the cold baryonic mass fraction (cold gas plus stars) of central galaxies normalized by the cosmic baryon fraction, f_b, as a function of halo mass. The red line shows the median of the whole posterior prediction, while the shaded regions with different colours encompass the 95 per cent credible ranges from the modes with β_LD ∼ 0, 2, and 6, as indicated in the figure. The black bars show the results of abundance matching obtained by Papastergis et al. (2012) using the Arecibo Legacy Fast ALFA (ALFALFA) and Sloan Digital Sky Survey (SDSS) data.

$The posterior predictive distribution of the cold baryon (stars plus cold gas) mass fraction as a function of halo mass. The shaded bands with different colours encompass the 95 per cent credible range of the prediction for the three modes represented by different values of βLD, while the red line is the predicted median of the entire posterior predictive distribution. The black points with error bars are from Papastergis et al. (2012). The vertical yellow dashed line marks the lowest halo mass constrained by the K-band luminosity function.$

Figure 4.

The posterior predictive distribution of the cold baryon (stars plus cold gas) mass fraction as a function of halo mass. The shaded bands with different colours encompass the 95 per cent credible range of the prediction for the three modes represented by different values of β_LD, while the red line is the predicted median of the entire posterior predictive distribution. The black points with error bars are from Papastergis et al. (2012). The vertical yellow dashed line marks the lowest halo mass constrained by the K-band luminosity function.

Open in new tab Download slide

The observed K-band luminosity function only determines absolute magnitudes up to K = Mag_K − 5log h = −20, which only constrains the halo virial mass above ∼10¹¹ M_⊙ (indicated as the vertical line in Fig. 4). The 95 per cent range of the posterior prediction contains the abundance matching result for haloes more massive than a few times 10¹¹ M_⊙, but all three modes overpredict the cold baryon fraction for lower mass haloes. Among the three modes, models with larger β_LD produce lower cold baryon mass fractions for low-mass haloes. The mode with β_LD ≈ 0 vastly overpredicts the cold baryon mass fraction in haloes with masses ≤10^11.5 M_⊙.

To explore the physical implications of feedback, we predicted the average SFR and OFR in each mass bin for all central galaxies at z = 0. Fig. 5 shows the ratio OFR/SFR, which is commonly referred to as the mass-loading factor of the wind. Again the red lines show the medians of the whole posterior prediction, while the shaded regions with different colours encompass the 95 per cent credible range from the three modes defined by their values of β_LD. The predicted OFR/SFR ratio generally decreases with stellar mass; it is approximately 10 for galaxies with stellar masses ≲ 10^8.5 M_⊙ and close to 1 for Milky Way mass galaxies. Observationally, there is no strong evidence for the existence of such a strong stellar mass dependence for outflows in normal disc galaxies at the present day (see Rupke, Veilleux & Sanders 2002; Heckman 2003; Veilleux, Cecil & Bland-Hawthorn 2005; Martin 2006). In starburst galaxies, the mass OFRs are roughly comparable to the SFR, so that OFR/SFR ∼1, and lower mass galaxies do not show stronger outflows. Thus, the high OFR/SFR ratio predicted for the general dwarf galaxy population does not have any observational support, but then again there is no strong observational evidence against it. Comparing the predictions for the different posterior modes, we see that different choices of β_LD produce very different trends of the loading factor with galaxy stellar mass. The β_LD ≈ 0 mode produces a constant mass-loading factor, OFR/SFR ≈ 3 over a large stellar mass range, while the other two modes predict a rapidly decreasing OFR/SFR with increasing stellar mass. The roughly constant OFR/SFR predicted by the β_LD ∼ 0 mode is consistent with the observational results of Martin (2006), but this mode significantly overpredicts the cold baryon mass fraction for low-mass haloes. Clearly, as one needs a large mass-loading factor to achieve better fits to the constraining data, observations of gas outflows in low-mass galaxies, and how they scale with galaxy mass can provide important constraints on the feedback model family considered here.

Figure 5.

The predicted ratio of mass OFR to SFR as a function of stellar mass for central galaxies at z = 0. The colour bands and lines are as described in Fig. 4.

Open in new tab Download slide

Fig. 6 shows the predicted cosmic SFR density as a function of redshift using the posterior samples selected from the three modes of β_LD. The grey shaded region in the figure denotes a compilation of observational estimates of the cosmic SFR density as a function of redshift from Behroozi et al. (2013). Compared to the data, our model family predicts a different shape. The cosmic SFR density in the model reaches its maximum at higher redshifts (z ≈ 3.5), while data suggest it more sharply peaks around z ≈ 2, although the observed drop at high redshifts could partly owe to missing low-mass galaxies (Lu et al. 2014). The model also underpredicts the SFR density around z = 2, as the median of the posterior predictive distribution is only marginally enclosed by the lower bound of the observational data. The predictions from the β_LD ≈ 0 mode are systematically higher and might be more consistent with the observations. The high β_LD mode predicts a SFR density that is too low to match the observational data at z > 1. This stems from the powerful SN-driven outflows in low-mass haloes required to fit the local galaxy luminosity function and H i mass function. Because the parametrization for feedback does not have an explicit redshift dependence, such strong outflows work throughout all cosmic time in the model. This feedback suppresses star formation at high redshift where such low-mass haloes dominate.

Figure 6.

The posterior prediction of the comoving SFR density versus redshift. Points with error bars are observational data. The grey band denotes observational estimates compiled in Behroozi et al. (2013). The other colour bands and lines are as described in Fig. 4.

Open in new tab Download slide

Fig. 7 shows the SFR as a function of redshift for haloes with masses of M_vir = 10¹² M_⊙ today (z = 0) for the three modes in β_LD. The model predictions are compared with the recent results obtained by Lu et al. (2014) and Behroozi et al. (2013) using empirical models tuned to fit the observed stellar mass functions at different redshifts. The empirical models predict a peak near z ≈ 1 while our SAM family predicts a much weaker peak at higher z. This is consistent with the underprediction of the SFR density at z ≈ 1 in Fig. 6. The β_LD = 0 mode also overpredicts the SFR at z > 2, although this mode better matches the observed cold gas fraction in haloes with masses of 10¹² M_⊙ (see Fig. 4). Overall, the predicted star formation histories for Milky-Way-sized galaxies by our SAM family are inconsistent with the results obtained from the empirical models.

The posterior prediction of the star formation rate history for haloes with present-day mass 1012 M⊙. Points with error bars show results obtained with two recent empirical models. The colour bands and lines are as described in Fig. 4.

Figure 7.

The posterior prediction of the star formation rate history for haloes with present-day mass 10¹² M_⊙. Points with error bars show results obtained with two recent empirical models. The colour bands and lines are as described in Fig. 4.

Open in new tab Download slide

Weinmann et al. (2012) found that SAMs and numerical simulations cannot generally reproduce the observed number density evolution of low-mass galaxies. For comparison, we predict the number densities of galaxies with stellar masses in the range 9.27 ≤ log M_⋆/ M_⊙ ≤ 9.77 from z = 0 to 6, and compare them with the observational data compiled by Weinmann et al. (2012). Fig. 8 shows that the model family can accommodate the observational data at z = 0, but vastly overpredicts the number density of low-mass galaxies at higher redshift. This confirms the results of Weinmann et al. (2012), and suggests that the current model family cannot match the data even if a large parameter space is explored.

Figure 8.

The posterior prediction for the evolution of the number density of galaxies with stellar masses in the range log (M_⋆/ M_⊙) = 9.27–9.77 as a function of redshift. Points with error bars show data compiled by Weinmann et al. (2012). The colour bands and lines are as described in Fig. 4.

Open in new tab Download slide

Fig. 9 shows the prediction of the comoving mass density of cold gas as a function of redshift. The observational data at high redshifts, estimated using damped Lyman α systems (Péroux et al. 2003; Prochaska & Wolfe 2009), suggest that the cold gas density at high redshifts is higher than in the local universe. Our predicted increase of the cold gas mass density with redshift in the range from z = 0 to 2 is consistent with the observations. The model also predicts a decrease of the cold gas density with increasing redshift at z > 3. Unfortunately the current data are still too uncertain to provide meaningful constraints on the model.

$The posterior prediction for the comoving cold gas mass density normalized by the critical density of the universe at the present time. The blue dashed line shows the total mass in haloes with mvir ≥ 5 × 109 h−1 M⊙ times the universal baryon fraction fb. Points with error bars are observational data from Zwaan et al. (1997, 2005), Martin et al. (2010), Rao & Briggs (1993), Rao, Turnshek & Nestor (2006), Lah et al. (2007), Péroux et al. (2003), and Prochaska & Wolfe (2009). The colour bands and lines are as described in Fig. 4.$

Figure 9.

The posterior prediction for the comoving cold gas mass density normalized by the critical density of the universe at the present time. The blue dashed line shows the total mass in haloes with m_vir ≥ 5 × 10⁹ h⁻¹ M_⊙ times the universal baryon fraction f_b. Points with error bars are observational data from Zwaan et al. (1997, 2005), Martin et al. (2010), Rao & Briggs (1993), Rao, Turnshek & Nestor (2006), Lah et al. (2007), Péroux et al. (2003), and Prochaska & Wolfe (2009). The colour bands and lines are as described in Fig. 4.

Open in new tab Download slide

In summary, the model family considered here, which has physical prescriptions similar to most published SAMs, in addition to not fitting one of the constraining data sets, the observed K-band luminosity function, also fails to predict additional observations as follows.

It requires very efficient feedback to eject gas from low-mass haloes. The energy required is probably too high if the outflows are driven by the kinetic energy of SN explosions alone, although this problem may be alleviated by including radiation pressure.
It requires a density threshold for star formation that is much lower than current observations suggest, owing to the energy needed to eject cold gas from the disc.
The high mass-loading factor required by the present-day distribution of low-mass galaxies predicts mass OFRs that might be too high to match observation.
The predicted SFR history of Milky-Way-sized galaxies is inconsistent with that obtained from empirical models.
It overpredicts the number density of low-mass galaxies at high redshift.

These results together suggest that some important physics may still be missing in the current SAMs.

5 DISCUSSION AND CONCLUSIONS

We use the observed low-redshift K-band luminosity function and H i gas mass function of galaxies as data constraints to explore star formation and feedback in a SAM of galaxy formation, which uses MCMC-based Bayesian methodology. The combination of these two data sets breaks the degeneracy in the parametrized phenomenology based on using the K-band luminosity function data constraint alone (Lu et al. 2012). Our model assumes that baryons accrete into the assembling halo where they cool, form stars, and the resulting feedback ejects some fraction of the gas. Our model has many features in common with most published SAMs and our results focus on conclusions generic to a large class of published models. We marginalize over all the uncertainties in the model family that are not constrained by the data sets.

Our broad parameter space covers and extends the ranges of the important parameters for star formation and feedback adopted in many existing SAMs. Our exhaustive parameter-space search results in the robust identification of several compact modes that best match the constraining data. These modes are likelihood dominated, i.e. the existence of these modes does not depend sensitively on the prior distribution as long as the modes are supported. The analyses presented in this paper focus on the behaviour of the model constrained to those modes, namely on the physical implications of the constrained parameters and the predictions of these modes. These conclusions are not expected to depend sensitively on the prior. However, the choice of the prior distribution is generally important in complex high-dimensional problems, especially when the likelihood is not sufficiently informative. In such cases, strong prior belief easily sways weak observational constraints. Conversely, a thorough investigation of the poorly constrained model motivates the acquisition of data that robustly discriminates between competing hypotheses of interest.

Overall, our model fails to reproduce the joint data sets and, therefore, more work will be required to understand the nature of the failure and to propose a remedy before attempting a quantitative Bayesian analysis. To assess the fit, we use the PPC method with a residual sum of squares for each bin, for both the K-band luminosity function and the H i gas mass function, as a discrepancy statistic. In this method, the discrepancy with the data is compared with the distribution predicted from the posterior distribution. We find that the data are in the tail of the distribution defined by the posterior, suggesting the implausibility of the data given the model.

The primary cause for the failure is that to simultaneously produce the observed stellar component and the observed cold gas fraction at the low-mass end requires the SN feedback to be extremely efficient in expelling gas. This requires a high mass-loading factor that consumes much of the available energy budget. The low p-value we found by exhaustively exploring the parameter space has its roots in the over simplicity of the phenomenological model and suggests that key physical processes are missing or misrepresented.

In addition, we find that the mass-loading factor must depend strongly on halo circular velocity to obtain the shallow faint-end slope of the luminosity function. The model requires that winds from haloes with circular velocities lower than ∼200 km s⁻¹ are ejected from the halo, completely. However, in many mass-loading models the mass-loading factor is a weak function of the halo circular velocity. Recent simulations concur: Davé et al. (2013) show that a steeper relation between mass loading and halo circular velocity tends to produce a shallower stellar and H i mass function at the low-mass end. This also suggests that accurate observational data for the faint-end slope of the galaxy luminosity function are crucial to constrain feedback models.

The inferred requirement for efficient feedback is consistent with the recent results of Mutch et al. (2013) and Henriques et al. (2013) who also found that a high efficiency of SN feedback is needed to fit the stellar mass and/or luminosity functions of galaxies at multiple redshifts. However, the high efficiency implied by the posterior distribution obtained here is not supported by hydrodynamical simulations (e.g. Mac Low & Ferrara 1999), which showed that the feedback efficiency of SN kinetic energy is usually quite low. Thus, either some other important energy sources are missing or the way feedback works is not correctly modelled in the current model family. For example, radiation pressure associated with star formation and SN explosions may provide a viable solution to the energy problem found here (e.g. Stinson et al. 2013).

Finally, the model predicts the star formation history over all time, although it is only constrained by data at z = 0. The observed star formation histories of Milky-Way-sized haloes and the observed redshift evolution in the number density of low-mass galaxies at high z are inconsistent with our posterior predictions. This suggests that star formation and feedback may have a more complex redshift dependence than assumed here.

In summary, our analysis shows that we require a high wind mass-loading factor and mass ejection rate of galactic winds to match the observed luminosity and cold gas mass functions. This result applies to many if not all ejective feedback scenarios that depend on winds, largely independent of the mechanistic details. This is a direct consequence of a short radiative cooling time-scale in low-mass haloes combined with the small fraction of baryons in stars and cold gas relative to the cosmic baryon fraction. Moreover, the predicted high mass-loading factor and mass-loss rates are not observed, suggesting that a feedback-generated wind may not be the agent suppressing star formation and cold gas accretion in low-mass haloes, at least at low z.

Alternatively, this suggests exploring mechanisms that prevent the gas accretion in the first place. For example, the IGM may be pre-heated by some processes so that low-mass haloes can accrete baryonic matter only at a reduced rate, either by properly including the effects of already considered processes like SN winds or by additional processes like gravitational pancaking or blazar heating (e.g. Mo & Mao 2002; Mo et al. 2005; Lu & Mo 2007; Anderson & Bregman 2010; Zhu, Feng & Fang 2011; Pfrommer, Chang & Broderick 2012). This may also introduce a characteristic time before and after which the star formation and feedback proceed differently, as seems to be required to match the star formation history and number density evolution of low-mass galaxies (see Lu et al. 2014). We will investigate the effects of pre-heating on galaxy formation in a future paper.

This work is partly supported by grants NSF AST-1109354 (HJM and MDW), NSF AST-0908334 (HJM and NK), and NASA ATP NNX10AJ95G (NK). Part of this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. YL thanks Tom Abel, Peter Behroozi, Eric Bell, Andrew Benson, Darren Croton, Romeel Davé, Andrey Kravtsov, Emmanouil Papastergis, Joel Primack, Rachel Somerville, Frank van den Bosch, and Risa Wechsler for useful discussion.

REFERENCES

Anderson

M. E.

Bregman

J. N.

. ,

ApJ

2010

, vol.

714

pg.

320

Month:	Total Views:
January 2017	2
February 2017	2
April 2017	2
June 2017	1
July 2017	5
August 2017	1
September 2017	2
October 2017	3
November 2017	2
December 2017	9
January 2018	9
February 2018	11
March 2018	14
April 2018	7
May 2018	9
June 2018	2
July 2018	13
August 2018	10
September 2018	3
October 2018	5
November 2018	5
December 2018	9
January 2019	4
February 2019	9
March 2019	8
April 2019	12
May 2019	10
June 2019	5
July 2019	9
August 2019	4
September 2019	5
October 2019	7
November 2019	6
December 2019	2
January 2020	1
February 2020	2
March 2020	5
April 2020	2
May 2020	1
June 2020	9

Article Contents

Abstract

1 INTRODUCTION

2 THE SEMI-ANALYTIC MODEL

3 DATA AND LIKELIHOOD

3.1 Uncertainties and covariance in the observed H i mass function

3.2 Corrections for molecular fraction and incompleteness

3.3 The likelihood function

4 MODEL INFERENCES

4.1 Inference from the posterior distribution

4.1.1 Star formation

4.1.2 Feedback

4.2 Posterior predictions

5 DISCUSSION AND CONCLUSIONS

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only