How sensitive are predicted galaxy luminosities to the choice of stellar population synthesis model?

Abstract

1 INTRODUCTION

The luminosity function is a basic observable of the galaxy population and is of fundamental importance since it encodes information about the physics of galaxy formation and evolution. If galaxies populate their host dark matter haloes in a simple way such that the mass-to-light ratio of all the galaxies within a halo is constant, we would expect to observe more galaxies at both bright and faint luminosities than are seen (Baugh 2006). The formation of galaxies at extreme luminosities is prevented by different mechanisms related to the processes affecting the cooling of gas (White & Rees 1978; Kauffmann, White & Guiderdoni 1993; Benson et al. 2003; De Lucia et al. 2010; Lacey et al. 2011). At high redshift, photoionization of the intergalactic medium (IGM) prevents the formation of faint galaxies, through the suppression of both the accretion of baryons and gas cooling in low circular velocity haloes (Benson et al. 2002; Okamoto, Gao & Theuns 2008). However, this is not enough by itself to explain the observed faint end of the galaxy luminosity function at the present day. Stellar feedback, the reheating of the gas by supernovae (SNe) in the discs of galaxies forming stars, also has an impact on the faint end of the luminosity function (Cole et al. 1994; Benson et al. 2003). Models of galaxy formation and evolution suggest that the bright end of the luminosity function is mostly affected by the feedback from active galactic nuclei (AGN) and the duration of starbursts (Bower et al. 2006; Croton et al. 2006, Lacey et al., in preparation).

Galaxy formation models are useful tools to understand which physical processes shape the evolution of observed galaxies properties, such as their luminosities. Stellar population synthesis (SPS) models are used to calculate the luminosity of galaxies from the predicted star formation and metal-enrichment histories directly predicted by the galaxy formation models. Thus, the comparison between model predictions and observations is subject to the accuracy of the SPS models (Conroy, Gunn & White 2009; Chen et al. 2010). The performance of SPS models is mainly limited by the uncertainties in poorly constrained stellar evolutionary phases, in particular, the asymptotic giant branch and later phases (e.g. Charlot, Worthey & Bressan 1996). Mancone & Gonzalez (2012) studied the colours modelled using different SPS models, and found that they agree best in the optical for old ages and solar metallicities. In their study, they pointed out that the dependence of SPS uncertainties on both age and wavelength will translate into uncertainties which change with redshift in a non-trivial way.

Over the last decade, SPS models have incorporated improved descriptions of stars in a TP-AGB phase (Bruzual & Charlot 2003; Maraston 2005; Conroy & Gunn 2010), and also of horizontal branch (HB) stars which are particularly relevant for the ultraviolet (UV) emission from old stellar populations (Conroy & Gunn 2010; Smith, Lucey & Carter 2012). The contribution from TP-AGB stars is notably controversial, since some observations, in particular those of the characteristics of globular clusters, favour models with stronger emission from TP-AGB stars (Goudfrooij et al. 2006; Eminian et al. 2008; MacArthur et al. 2010; Riffel et al. 2011; Lyubenova et al. 2012; Andreon 2013), while other observations, such as those focusing on post-starburst galaxies, favour models with a smaller contribution from TP-AGB stars (Kriek et al. 2010; Melbourne et al. 2012; Zibetti et al. 2013). Stars in the TP-AGB phase can dominate the near-infrared (NIR) luminosity (Maraston 2005). Rest-frame NIR luminosities are less affected by dust than shorter wavelengths and thus, have been used to derive stellar masses. However, the uncertainties regarding the TP-AGB phase could hinder these determinations and better constraints are needed (Kelson & Holden 2010; McGaugh & Schombert 2013).

Semi-analytical galaxy formation models combined with an SPS model with a strong TP-AGB phase are able to reproduce the observed numbers of extremely red galaxies in place at z > 1 (Tonini et al. 2009, 2010; Fontanot & Monaco 2010; Henriques et al. 2011, 2012). However, the need for a strong TP-AGB phase to match this galaxy population is not found in all models (Gonzalez-Perez et al. 2009). Here, we assess the robustness of the predicted luminosities at different wavelengths by combining a new version of the galform semi-analytical model with seven different SPS models. Among these, there are two which include a strong TP-AGB phase (Maraston 2005; Bruzual 2007) and one in which the intensity of the TP-AGB phase can be modified (Conroy et al. 2009). All of the SPS models studied cover a large range in wavelength, from the UV to the NIR and beyond in some cases, and most are publicly available.

The model used in this study is a new development of galform, the semi-analytical model for the formation and evolution of galaxies developed mainly at Durham University (Cole et al. 2000). Much of the previous work with galform has used the Millennium Simulation. The cosmology assumed for the Millennium Simulation, which is close to that derived from the Wilkinson Microwave Anisotropy Probe (WMAP1) observations (Spergel et al. 2003), is inconsistent with current observations. Here, we have updated the Lagos et al. (2012) model to the WMAP7 cosmology (Komatsu et al. 2011), which is also close to Planck cosmology [Ade et al. (Planck Collaboration) 2013] on the scales relevant to galaxy formation. In order to do this update, we have implemented the galform semi-analytical model using the dark matter halo trees derived from the MS-W7 N-body simulation (Lacey et al. in preparation), which is based on a WMAP7 cosmology. A similar MS-W7 simulation has previously been used in combination with a semi-analytical model of galaxy formation to explore how to rescale the model predictions from those for an underlying cosmology similar to WMAP1 to the WMAP7 one (Guo et al. 2013). In order to update the galform semi-analytical model to the WMAP7 cosmology, we have found the best-fitting model parameters to reproduce a set of observations as discussed in Section 2.

Lagos et al. (2011b) introduced a new implementation of the star formation law which follows the atomic and molecular hydrogen in the interstellar medium (ISM) based on the empirical relation between the star formation rate and the surface density of molecular hydrogen inferred by Blitz & Rosolowsky (2006). This new implementation is more realistic than the simplified one used before in the majority of semi-analytical models, in which the star formation was assumed to be proportional to the total mass of cold gas in a galaxy (e.g. Cole et al. 2000).

The new model presented here uses a single initial mass function (IMF) in all modes of star formation and can be considered as an alternative to the Lacey et al. (in preparation) model, which is also implemented in a WMAP7 cosmology but which does not assume a universal IMF. The universality of the IMF has been challenged by different observations (e.g. Wilkins et al. 2008; Cappellari et al. 2012; Geha et al. 2013; La Barbera et al. 2013). In galform, a top-heavy IMF has been found to be needed in order to reproduce the observed numbers and redshift distribution of submillimeter galaxies, i.e. luminous dusty galaxies at redshifts around 1 < z < 3 (see Fig. A8 in the Appendix for further details). However, it is very uncertain how and in which environments the IMF changes. The comparison of the model predictions using different SPS models is simpler for a model with a single IMF, as is the case in this paper.

In Section 2, we describe the new model of galaxy formation and evolution used in this study. Section 3 summarizes the differences between the SPS models considered. In Sections 4 and 5, we compare the luminosity functions and number counts predicted using different SPS models with observations. Our conclusions are presented in Section 6. The Appendix contains further comparisons between observations and predictions of the new model presented in this work. This model is available in the Millennium Archive Data Base.¹ There, it is also possible to get the post-processing calculation introduced by Lagos et al. (2012) for estimating the emission of the most widely used tracer of molecular hydrogen: the carbon monoxide, ¹²CO.

2 GALAXY FORMATION MODEL

We present a new development of galform (Cole et al. 2000), a semi-analytical galaxy formation model set in a Λ cold dark matter universe. Semi-analytical models use simple, physically motivated equations to follow the fate of baryons in a universe in which structure grows hierarchically through gravitational instability (see Baugh 2006 and Benson 2010, for an overview of hierarchical galaxy formation models).

galform models the main processes which shape the formation and evolution of galaxies. These include: (i) the collapse and merging of dark matter haloes; (ii) the shock-heating and radiative cooling of gas inside dark matter haloes, leading to the formation of galactic discs; (iii) quiescent star formation in galactic discs; (iv) feedback from SNe, from AGN and from photoionization of the IGM; (v) chemical enrichment of the stars and gas; (vi) galaxy mergers driven by dynamical friction within common dark matter haloes, leading to the formation of stellar spheroids, which also may trigger bursts of star formation. The model also computes the scale size of the disc and bulge component of the galaxy. Galaxy luminosities are obtained by combining an SPS model with the star formation and metal enrichment histories predicted for each model galaxy. The attenuation of starlight by dust is modelled in a physically self-consistent way, based on the results of a radiative transfer calculation (Ferrara et al. 1999) for a realistic geometry in which stars are distributed in a disc plus bulge system (Cole et al. 2000; Lacey et al. 2011; Gonzalez-Perez et al. 2013). The end product of the calculation is a prediction of the number and properties of galaxies that reside within dark matter haloes of different masses. The free parameters in the semi-analytical model are chosen in order to reproduce the set of observations described in Section 2.2.

The model presented here is an updated version of the Lagos et al. (2012) model. The Lagos et al. model adopts the cosmological parameters of the Millennium Simulation (Springel et al. 2005), which correspond approximately to the results of the first year of the WMAP satellite (WMAP1 Spergel et al. 2003). For the work presented here, we have updated this model to the best-fitting cosmological parameters of the WMAP 7 year data set (WMAP7 Komatsu et al. 2011). In Section 2.2, we give details of how this updated model differs from the original one.

Lagos et al. introduced a new calculation of the star formation rate in galaxies compared with that used in the Bower et al. (2006) model. Lagos et al. also increased the starburst duration, without further changes to other parameters. The Lagos et al. model was not actually returned to obtain a good match to observational data, but nevertheless, gave reasonable agreement with the observed b_j and K-band luminosity functions at z = 0, the H i and H₂ mass functions at z = 0, the observed cold gas content of galaxies up to z = 2, the bimodality in the star formation rate versus stellar mass plane for local galaxies and the cosmic density evolution of H i (Lagos et al. 2011a,b, 2012).

In Bower et al., the star formation rate is assumed to be simply proportional to the mass of cold gas present in the galaxy and inversely proportional to the dynamical time. Recent observations have motivated a more sophisticated calculation in which the quiescent surface density of the star formation rate is assumed to be proportional to the surface density of molecular hydrogen mass in the ISM (Blitz & Rosolowsky 2006; Bigiel et al. 2008; Leroy et al. 2008). In this calculation, the ratio between the molecular and total gas is determined by a pressure law (Blitz & Rosolowsky 2006). Lagos et al. (2011b) implemented a self-consistent calculation into galform in which H i and H₂ are tracked explicitly and the star formation in discs is assumed to depend on the amount of molecular gas, H₂, rather than on the total mass of cold gas. The introduction of this observationally motivated model of the star formation rate leads to a reduction in the size of the available parameter space, as the star formation parameters are determined empirically from the observations. Furthermore, the constraints on the model are extended through the ability to make new predictions such as the H i and H₂ mass functions.

2.1 Dark matter halo merger trees

The model that we present here is implemented in a new Millennium Simulation run with the WMAP7 cosmology (MS-W7; Lacey et al. in preparation). This cosmology is very similar to the best fit for the recently released Planck data [Ade et al. (Planck Collaboration) 2013]; the power spectra corresponding to the best-fitting models to WMAP7 and Planck data are very close to one another on the scales relevant to galaxy formation. The MS-W7 N-body simulation uses 2160³ particles, each with a mass of 9.35 × 10⁸ h⁻¹ M_⊙, in a box of side 500 h⁻¹ Mpc and with a starting redshift of z = 127 (see Springel et al. 2005, for details of the original Millennium Simulation). The MS-W7 simulation uses the WMAP7 cosmological parameters (Komatsu et al. 2011), that are also summarized in Table 1: matter density, Ω_m0 = 0.272, cosmological constant, Ω_Λ0 = 0.728, baryon density, Ω_b0 = 0.045, a normalization of density fluctuations given by σ₈ = 0.807 and a Hubble constant today of H₀ = 100 h km s⁻¹ Mpc⁻¹, with h = 0.704. The 61 outputs of the simulation are approximately spread evenly in the logarithm of the expansion factor and so do not correspond to round numbers in redshift.

galform calculates the evolution of galaxies in halo merger trees which describe the assembly and merger histories of cold dark matter haloes. For the model presented here, the halo merger trees have been extracted from the MS-W7 N-body simulation.

Briefly, the construction of the merger trees starts with two consecutive steps: (1) groups of dark matter particles are identified in each simulation snapshot using the friends-of-friends algorithm (Davis et al. 1985); (2) self-bound, locally overdense sub-groups are identified using the algorithm subfind (Springel et al. 2001). In some cases, this procedure identifies structures as groups that might better to be considered as distinct haloes for the purpose of implementing the semi-analytical galaxy formation model. The merger tree algorithm used in this work deals with such cases and ensures that the resulting trees are strictly hierarchical, i.e. once two haloes are considered to have merged they remain merged at all later times. A full description of the construction method of the halo merger trees is presented in Jiang et al. (2013). This new method is similar to that described in Merson et al. (2013) but differs in two points.

• When determining the main progenitor for each subfind sub-group, the new scheme tries to identify the most bound ‘core’ of dark matter particles, rather than the most massive one.

• When looking for descendants, the new scheme looks for the same most bound ‘core’ of particles. As with (i), before this was done by looking for the most massive progenitor.

These changes were mainly driven by issues that arise in the construction of trees in very high resolution N-body simulations and do not have a noticeable impact here.

2.2 galform in the WMAP7 cosmology

The new model presented here is essentially an update of the Lagos et al. model to the WMAP7 cosmology. Table 1 summarizes the cosmological parameters used in the Lagos et al. model, which are close to those obtained from WMAP1 (Spergel et al. 2003), and the WMAP7 parameters used here (Komatsu et al. 2011). The growth of structure in the WMAP1 and WMAP7 cosmologies can be compared by looking at the parameter combination |$\Omega _{\text{m}0}^{0.55} \sigma _{8}$| (Eke, Cole & Frenk 1996; Linder 2005). The differences between the best-fitting values of |$\Omega _{\text{m}0}$| and σ₈ in these two cosmologies almost compensate one another in this parameter combination, so that |$\Omega _{\text{m}0}^{0.55} \sigma _{8}$| differs by only 6 per cent. Although this change result in only a little difference in the mass function of dark matter haloes in these cosmologies, it does have a visible impact on the predicted luminosity functions at z = 0. Fig. 1 shows the luminosity functions at z = 0 in b_J and K-bands predicted using the original set of galaxy formation parameters from the Lagos et al. model but implemented in the WMAP7 cosmology. This figure shows the effect that changing the cosmological parameters has on the model of galaxy evolution which, when implemented in the WMAP1 cosmology, satisfactorily reproduced the observations.

Figure 1.

The predicted luminosity functions at z = 0 (solid lines), in the b_J-band (λ_eff = 4500 Å, top) and in the K-band (λ_eff = 2.2 μm, bottom), compared with observations from Cole et al. (2001) (triangles), Kochanek et al. (2001) (crosses), Norberg et al. (2002) (squares) and Driver et al. (2012) (circles). The blue lines show the predictions from the Lagos et al. model run in the WMAP7 cosmology, while the red lines show the predictions from the new model presented here. The dashed lines show the predicted luminosity function without dust attenuation.

Open in new tabDownload slide

In this section, we describe the parameters that control the physical processes in the model which were changed in order to reproduce the b_J and K-band luminosity functions at z = 0 and to give reasonable evolution of the predicted rest-frame UV and K-band luminosity functions.

Before starting to tune the model parameters, for this work we decided not to consider the metal yield,²y, and recycled fraction,³R, as free parameters. These two parameters are set to those values calculated using the Padova stellar evolution models for the Kennicutt (1983) IMF: yield = 0.021 and recycled fraction = 0.44 (for a discussion of the uncertainties on these values we refer the reader to Cole et al. 2000). Note that these values are different from those used in Bower et al. (2006) and Lagos et al. (2012).

In galform, it is assumed that the available reservoir of cold gas is consumed during a starburst event with a finite duration. In our model, as well as in Bower et al. and Lagos et al., starbursts can be triggered by both mergers and discs becoming dynamically unstable. The time-scale of a starburst event is considered to be proportional to the bulge dynamical time, τ_bulge, through a parameter f_dyn, except when this time-scale is too small, then a minimal duration time is considered, τ_min (for further details see Lacey et al. 2011). Thus, due to the scaling with the dynamical time of the bulge, the duration of starbursts will change with redshift. Starting from the values advocated by Lagos et al., we find a compromise by decreasing τ_min from 0.1 Gyr to 0.05 Gyr, while increasing f_dyn from 2 to 10. A model including these further changes slightly overpredicts the bright end of the K-band luminosity function at z = 0. Therefore, we increased very slightly the strength of the AGN feedback.

The onset of the AGN suppression of the cooling flow in the model is governed by a comparison of the cooling time of the gas, t_cool with the free-fall time for the gas to reach the centre of the halo, t_ff. Galaxies with t_cool > t_ff/α_cool, where α_cool is a model parameter, are considered to be hosted by haloes undergoing quasi-hydrostatic cooling. In such haloes, cooling is suppressed if the luminosity released by gas accreted on to a central supermassive black hole balances or exceeds the cooling luminosity (see Bower et al. 2006; Fanidakis et al. 2011, for further details). Increasing α_cool slightly from 0.58 to 0.6 was found to be enough to compensate the changes in the burst duration which affected the K-band luminosity function. This change implies that slightly more haloes will be quasi-static, leading to the possible quenching of gas cooling and hence star formation.

Table 1 lists the parameters changed from the values used in the Lagos et al. (2012) model in order to update it to the WMAP7 cosmology, once the value of the yield and recycled fractions were fixed to be consistent with the predictions from stellar evolution models for the Kennicutt IMF used here. The predicted b_J- and K-band luminosity function at z = 0 can be seen in Fig. 1 compared with observations. Further predictions of the model introduced here are presented in the Appendix.

3 STELLAR POPULATION SYNTHESIS MODELS

SPS models provide a link between the predictions of galaxy formation and evolution models and observations. Models like galform predict the intrinsic properties of galaxies, such as their stellar masses, star formation and metal enrichment histories, gas content, etc. Given this information, the SPS models provide a way to estimate the spectral energy distribution (SED) of the model galaxies. In order to make this link, the SED of simple stellar populations (SSPs), i.e. of a group of coeval stars with the same metallicity, must be known.

The SPS models compute the characteristics of SSPs by first adopting a stellar evolutionary model. For a given initial stellar mass and metallicity, these models calculate the stellar isochrones, i.e. the luminosity of a star at different ages, from the main sequence zero-point to its death. Then stellar spectral libraries are used to assign a full spectrum to each set of values for the initial stellar mass, metallicity and age. The SED of a SSP is then obtained by integrating the light from coeval stars weighted by the chosen stellar IMF. The IMF describes the initial number of stars formed in a given stellar mass bin. The IMF is usually approximated by one or more power laws (though, other parametrizations have also been used in the literature, in particular for masses below 1 M_⊙, see Chabrier 2003; van Dokkum 2008). This approximation is such that, in a given stellar mass interval, the IMF can be expressed as: dN(m)/d ln m ∝ m^−x.

3.1 Comparing SPS models

As described above, the main ingredients of an SPS model are: (i) the stellar evolutionary tracks, (ii) stellar spectral libraries and (iii) the IMF. By default, galform assumes a Kennicutt (1983) IMF for which dN(m)/d ln m ∝ m^−0.4 if m ≤ 1 M_⊙ and dN(m)/d ln m ∝ m^−1.5 otherwise. We will be using the different SPS models in combination with the Kennicutt (1983) IMF with a mass range set to |$0.1 \le m(\text{M}_{\odot }) \le 100$|⁠, except when this IMF is not provided in the public release of the SPS model (only for Figs 2, 3 and 7). In such cases, we will use a Salpeter IMF for which dN(m)/d ln m ∝ m^−1.35, at all masses.

Figure 2.

An illustration of the differences between the single stellar population luminosity obtained for a Salpeter IMF and solar metallicity, normalized to 1 M_⊙, as a function of wavelength for fixed ages using the BC99, BC03, MN05, CB07, CW09, PD01 and P2 SPS models as indicated in the legend. From left to right and top to bottom, the panels show the spectra of the SSPs at ages: 0.5 Myr, 10 Myr, 500 Myr and 1.5 Gyr. Note that the y-axis limits change from 10 to 0.01, going from the first to second row.

Open in new tabDownload slide

Figure 3.

An illustration of the evolution with age of different single stellar population luminosity obtained for a Salpeter IMF with solar metallicity, normalized to 1 M_⊙, using the BC99, BC03, MN05, CB07, CW09, PD01 and P2 SPS models, as indicated in the legend. The luminosities are shown as a function of age convolved with a top-hat filter covering a fixed range in wavelength: [1450,1550] Å (left), [5350,5450] Å (centre) and [21900,22 000] Å (right). Note that the scale of the y-axis changes from the UV panel to the other two. The lower panels show the single stellar population luminosity for each SPS model, divided by that from the BC99 model on a linear scale.

Open in new tabDownload slide

Table 2 summarizes the main ingredients of the following SPS models: PD01 (Silva et al. 1998), P2⁴ (Fioc & Rocca-Volmerange 1997, 1999), BC99 (an updated version of Bruzual & Charlot 1993), BC03 (Bruzual & Charlot 2003), MN05 (Maraston 2005), CB07 (Bruzual 2007) and CW09 (FSPS v2.3, Conroy et al. 2009, 2010; Conroy & Gunn 2010). BC99 is the default SPS model used in the present model and also in other galform developments such as: Cole et al. (2000), Bower et al. (2006) and Lagos et al. (2012). PD01 is the SPS model used in the version of galform described by Baugh et al. (2005) and Lacey et al. (2011).

As can be seen in Table 2, all of these SPS models cover similar ranges in stellar metallicities and ages. In general, these SPS models can be used in combination with a variety of low- and high-resolution stellar spectral libraries. In order to cover the largest possible wavelength and metallicity range, here we use the SPS models combined with a low-resolution stellar spectra library. The PD01 and BC99 models are run in combination with an extended version of the spectra from Kurucz (1993) (see also Bressan et al. 1998). The other models make use of BaSeL, a semi-empirical library of stellar spectra (Lejeune et al. 1997, 1998; Westera et al. 2002), with extensions and modifications as indicated in Table 2.

All the SPS models used here, except the MN05 one, use isochrones to compute the SSPs for stars in different evolutionary stages. The PD01, P2, BC99 and BC03 models are based on the stellar tracks from the Padova group (mostly described in Bertelli et al. 1994), including modifications that differ from model to model, as indicated in Table 2. Both the CB07 and CW09 models are based on a later version of the stellar tracks from the Padova group which include a prescription for the TP-AGB evolution of low- and intermediate-mass stars (Marigo & Girardi 2007; Marigo et al. 2008). The CW09 SPS model also includes a modification for stars with very low masses with respect to the Pd08 stellar tracks.

The MN05 SPS model is the only one among those used in this study that computes the characteristics of the SSPs of post-main-sequence stars by integrating the amount of hydrogen and/or helium that is consumed during a given post-main-sequence phase. The fuel at a given age depends on the mass of stars completing their hydrogen burning phase. Thus, this integration neglects the dispersion of stellar masses in post-main-sequence phases, an aspect supported by observations. This approach has the advantage of minimizing the use of uncertain theoretical models for describing post-main-sequence phases of star evolution. For main-sequence stars, this model uses the tracks and isochrones from Cassisi, Castellani & Castellani (1997a), Cassisi, degl'Innocenti & Salaris (1997b), Cassisi et al. (2000) complemented with those from the Geneva group (Schaller et al. 1992; Meynet et al. 1994) for very young evolutionary stages and from the Padova group for systems with high metallicities. The MN05 model allows for blue and red populations on the HB, recommending the use of red populations for metallicities Z ≥ 0.01.

3.1.1 Comparing simple stellar populations

In order to gain some insight into the differences between the SPS models summarized in Table 2, in this section we compare the luminosities of SSPs. These are groups of coeval stars with a fixed metallicity born in a single instantaneous or short burst, with their luminosities weighted by the IMF.

To produce SSP spectra, the SPS models can be run in most cases with a small number of different IMFs. However, the IMF is not a free parameter for most of them. We have run the SPS models with a Salpeter IMF since it is one of the few IMFs available in all cases.

Fig. 2 shows the SSPs computed with the SPS models listed in Table 2 for a Salpeter IMF and solar metallicity. Each panel shows a fixed age. The SSPs plotted in Fig. 2 show similar global trends. At very early ages, the SSP luminosity is dominated by the strong UV emission from hot, massive stars. The emission at λ < 1000 Å declines with time, becoming almost non-existent after about 10 Myr. At this point in time, the emission in the infrared range starts to increase. After a few hundred Myr there is again feeble emission at λ ≤ 1000 Å coming from UV-bright old stars, such as hot stars on the HB or the post asymptotic giant branch stars (e.g. Renzini & Fusi Pecci 1988). After 1 Gyr, the NIR light is dominated by red giant branch stars and then the luminosity in the optical and infrared ranges declines with time.

The trends seen in Fig. 2 are found to be similar for both sub-solar and supersolar metallicities.

Though the general trends described above are followed by SSPs from all the SPS models, it is worth noting the cases where they differ. Fig. 2 shows that at ages below 1 Myr, the PD01 model produces more luminous SSPs at λ > 10⁴ Å than the other SPS models. The PD01 model also predicts a larger fraction of ionizing photons shortwards of 912 Å for young SSPs. These two differences decrease rapidly with age. The difference seen for NIR wavelengths is related to thermal continuum emission from H ii regions, whose contribution is included in the PD01 model (Bressan, Silva & Granato 2002) but not in the other models.

After ∼100 Myr, the SSP spectra become very different at wavelengths λ < 10³ Å. At these ages, this wavelength region is mostly sensitive to the light coming from stars in a hot post red giant branch phase, which is far less well understood than earlier stellar phases, in particular at ages above 5 Gyr (e.g. Brown et al. 2008; Conroy et al. 2009; Smith et al. 2012). This region of ages and wavelengths is the only one where the BC99 and BC03 models appear to differ noticeably from one another and, thus, the predicted evolution of the UV (∼1500 Å), optical and IR luminosity is practically indistinguishable for these two SPS models.

For the MN05 SPS model, Fig. 2 shows the SSPs of red HB stars, since we are considering the case of solar metallicity. Above 1.5 Gyr, the luminosity of the SSPs for blue HB stars from MN05 model extends bluewards, compared with the red HB ones shown in Fig. 2, though this is still far from the spectra obtained with the other SPS models. Nevertheless, we will compute luminosity function at wavelengths longer than 10³ Å, and thus these strong differences will not affect our later discussion. It is worth noting that for younger ages than ∼1.5 Gyr, the SSPs derived from the MN05 SPS model with blue HB stars are much fainter, by about an order of magnitude in the optical range, than the SSPs for the rest of the models. The use of blue HB stars within the MN05 SPS model is recommended for very metal poor cases.

The left-hand panel in Fig. 3 compares the SSP (Salpeter IMF, solar metallicity) luminosities obtained as a function of age and convolved with a top-hat filter defined over the wavelength range 1450 ≤ λ(Å) ≤ 1550. In this figure, it can be seen clearly that the UV luminosities for the different SSPs are in good agreement up to 1 Gyr, beyond which noticeable differences arise which are likely to be related to the poorly constrained advanced phases of stellar evolution, as mentioned before.

Both Fig. 2 and the middle panel of Fig. 3 show that the different SSPs are in remarkably good agreement in the optical range, with differences of less than a factor of 2 in luminosity at λ ≈ 5400 Å. The optical luminosity is dominated by main-sequence, sub-giant and red giant branch stars (e.g. Bruzual & Charlot 1993), stellar evolutionary phases which are well understood.

The right-hand panel in Fig. 3 compares the SSP luminosities as a function of age, convolved with a top-hat filter defined between 21 900 ≤ λ(Å) ≤ 22 000. Around 10 Myr both the MN05 and PD01 models produce SSPs with stronger emission in the infrared than in the other models. This difference is probably related to the modelling of supergiant stars. The difference lasts less than 50 Myr and thus will have a negligible effect on the predicted global luminosity functions.

At ages between 0.1 and 1 Gyr, the P2, CB07 and MN05 models produce SSPs with a noticeably higher NIR luminosity than the rest of the models, due to stars that are in a TP-AGB phase. Stars with initial masses from 0.8 to 8 M_⊙ evolve on to the asymptotic giant branch and can enter the TP-AGB phase, though only those with M > 2 M_⊙ will contribute significantly to integrated NIR luminosities (e.g. Frogel, Mould & Blanco 1990). Stars in the TP-AGB phase can be very luminous at wavelengths above 6000 Å, as can be seen in the right-hand panels of Fig. 2. However, observed TP-AGB stars have unknown ages, metallicities and mass-losses rates, which makes it very difficult to model them and to constrain their characteristics (e.g. Maraston 2005; Conroy et al. 2009). As can also be seen in Fig. 3, the time and duration at which the TP-AGB stars becomes dominant depends on the SPS model.

Between 1 and 2 Gyr, all SPS present an increase in NIR emission related to the evolution of low-mass stars, M_* < 2 M_⊙, through the helium flash (e.g. Bruzual & Charlot 2003).

4 THE EVOLUTION OF THE LUMINOSITY FUNCTION

In this section, we present the predicted galaxy luminosity function computed using different SPS models within galform, at three rest-frame wavelengths: UV (∼1500 Å), optical (∼5400 Å) and NIR (∼22 000 Å). As shown in the previous section, over the ranges in which we are interested, both the BC99 and BC03 models give very similar SSP spectra and thus, we will show results only for the BC99 model. The CB07 model differs from the previous two mainly in the infrared and thus, we will limit the comparison for this model to that wavelength range.

4.1 The rest-frame UV wavelength range

Fig. 4 shows the evolution of the predicted luminosity function in the rest-frame UV (1500 Å) for all galaxies at z = 0, 3.1 and 6.2. The predicted luminosity functions obtained starting from different SPS models are very similar (we note that this is also the case when making the comparison for a Salpeter IMF, including the BC03 and CB07 SPS models). At z = 0, the model overpredicts the number of UV bright galaxies. The same trend is found when increasing the resolution of the dark matter merger trees by using a Monte Carlo approach (Lagos et al. 2013a).

Figure 4.

The predicted rest-frame far-UV (∼1500 Å) luminosity function at z = 0 (top panel), z ∼ 3 (central panel) and z ∼ 6 (bottom panel), obtained with the SPS models from BC99 (black lines), CW09 (red lines), MN05 (blue lines), PD01 (yellow lines) and P2 (purple lines). The dotted lines show the luminosity function for the default SSP model (BC99) without attenuation by dust. The observational data is as follows: z = 0 – Wyder et al. (2006) (filled circles, 1500 Å) and Driver et al. (2012) (empty circles, 1500 Å); z = 3.1 – Sawicki & Thompson (2006) (filled circles, 1700 Å) and Reddy & Steidel (2009) (empty circles, 1700 Å); z = 6.2 – Bouwens et al. (2007) (filled circles, 1350 Å) and McLure et al. (2009) (empty circles, 1500 Å).

Open in new tabDownload slide

At z = 3.1 and z = 6.2, the predicted luminosity function turns down at luminosities fainter than those plotted in Fig. 4, due to the finite halo mass resolution of the MS-W7 simulation. For brighter magnitudes, those shown in Fig. 4, the match between predictions and observations is very good at z = 3.1, while at z = 6.2 the model underpredicts the number of UV faint galaxies. At these high redshifts there is a significant contribution to the bright end of the rest-frame UV luminosity function from galaxies experiencing a starburst. The mismatch seen at z = 6.2 could be alleviated if longer duration starbursts are allowed. However, increasing the burst duration increases the tension between the model and observations for the bright end of the K-band luminosity function at high redshift.

Similar trends to those shown in Fig. 4 are expected for colour selected Lyman Break Galaxies, since, as shown in Gonzalez-Perez et al. (2013), the predicted far-UV luminosity function in this case is almost identical to the total UV luminosity function up to the corresponding observational limits.

The rest-frame UV luminosity is sensitive mainly to both the recent star formation history and dust attenuation in a galaxy (see Gonzalez-Perez et al. 2013, for a discussion of the impact that the dust attenuation modelling has on the rest-frame UV luminosity). However, these two aspects depend on our modelling of the formation and evolution of galaxies and are treated independently of the particular SPS model adopted. The main uncertainties in the SPS models affecting the UV luminosity function are the modelling of the HB stars and the blue stragglers, which contribute to the galaxy light after about 5 Gyr (Conroy et al. 2009). While the characteristics and fractions of these two types of stars can have a strong impact on the colours of optically red galaxies (Brown et al. 2008; Smith et al. 2012), Conroy et al. reported that changing the fraction of stars in the blue HB has a negligible effect on the UV luminosity function.

4.2 The rest-frame optical wavelength range

Fig. 5 shows the evolution of the luminosity function in the rest-frame V band. This is remarkably similar for all of the SPS models considered, including the BC03 and CB07 SPS models. In fact, the predicted luminosity function obtained with the PD01 and P2 SPS models are practically indistinguishable from that with the BC99 SPS model. This was expected since the SPS models are mainly constrained by observations in the optical and, thus, mainly by stars on the main sequence (e.g. Conroy & Gunn 2010).

Figure 5.

From top to bottom, the predicted rest-frame V-band (λ_eff = 5404 Å) luminosity function at z = 0.5, 1.3, 1.8 and z = 3.1, obtained using the default SPS model (BC99, black lines) and the SPS models from CW09 (red lines), MN05 (blue lines), PD01 (yellow lines) and P2 (purple lines). The dotted lines show the BC99 luminosity function without attenuation by dust. The filled circles show the observations from Marchesini et al. (2012).

Open in new tabDownload slide

The predicted faint-end normalization increases by less than a factor of 2 from z = 3 to z = 1.3, and then remains practically constant. The model predicts stronger evolution at the bright end, where the break steepens at lower redshifts since at later times there are fewer bright stars on the main sequence. Both trends are in good agreement with the observations from Marchesini et al. (2012).

4.3 The rest-frame NIR wavelength range

The predicted rest-frame NIR luminosity function shows little evolution from z = 2.5 to z = 0, for all the implemented SPS models except the MN05 and CB07 ones. At higher redshifts, z = 3.1, the rest-frame NIR luminosity function normalization declines. These trends can be seen in Fig. 6 for the J and K bands and similar results are found for the H band. As for the rest-frame optical luminosity functions, the predicted rest-frame NIR luminosity function obtained with the PD01 and P2 SPS models are practically indistinguishable from the default SPS one, BC99.

Figure 6.

The predicted rest-frame J-band (λ_eff = 1.2 μm, left-hand panels) and K-band (λ_eff = 2.2 μm, right-hand panels) luminosity function at z = 0.5, 1, 2.4 and 3.1 (from top to bottom), compared to observational data. The lines show the predictions using a Kennicutt (1983) IMF and the SPS models from BC99 (black lines), CW09 (red lines), MN05 (blue lines), PD01 (yellow lines) and P2 (purple lines). The dotted lines show the BC99 luminosity function without attenuation by dust. The symbols correspond to observational data from Drory et al. (2003) (crosses), Pozzetti et al. (2003) (diamonds), Saracco et al. (2006) (triangles), Caputi et al. (2006) (circles), Cirasuolo et al. (2010) (squares) and Stefanon & Marchesini (2013) (stars).

Open in new tabDownload slide

At z = 1 and for all the SPS models studied, there are too many NIR faint galaxies (i.e. around ∼2 mag fainter than L_*). Lagos, Lacey & Baugh (2013b) argued that this could be related to the need for a more realistic modelling of the mass loading in SNe outflows. Alternatively, Henriques et al. (2013) proposed solving this problem by changing the reincorporation time-scale of gas heated by SNe and ejected from haloes.

The predicted rest-frame NIR luminosity functions without taking into account the attenuation by dust are also shown in Fig. 6 for the BC99 SPS model. Comparing these with the attenuated ones, we can see that there is an appreciable attenuation and that it increases with redshift, as noted by Mitchell et al. (2013). For a model galaxy with a given extinction optical depth and inclination, the attenuation is calculated in a self-consistent way taking into account the results of a radiative transfer model (Ferrara et al. 1999). In our model, the extinction optical depth is assumed to be proportional to the mass of cold gas and inversely proportional to the square of radial scale-length of the disc (Gonzalez-Perez et al. 2013). Thus, the increase of attenuation with redshift is probably related to galaxies having both larger amounts of cold gas and smaller disc sizes at high redshifts.

Fig. 6 shows that at z > 1, the rest-frame NIR luminosity function using the MN05 SPS model agrees with the observed bright end of the luminosity function, while overpredicting the number of faint galaxies. Besides, at z ≥ 1, the rest-frame NIR luminosity function predicted using this SPS model exhibits a global increase in normalization with respect to that at z = 0. This behaviour is also seen when the CB07 model is used. In order to be able to compare the predictions obtained using the CB07 SPS model with the other, we use the Salpeter IMF which is available for all seven considered SPS models. Note that changing the IMF affects the number of stars of a given mass present in a galaxy and thus, differences in the predicted luminosity functions are expected for the same SPS model. Fig. 7 shows the predicted rest-frame luminosity function in the K band (note that similar results are found in the J and H bands). This figure illustrates that, due to the intrinsic uncertainties in the SPS models, the predicted luminosity functions cannot match observations better than the spread due to the use of different SPS models.

Figure 7.

The predicted rest-frame K-band luminosity function at z = 1 (top panel) and z = 3.1 (bottom panel), using a Salpeter IMF combined with the default SPS model (BC99, black lines) and the SPS models from CW09 (red lines), MN05 (blue lines), PD01 (yellow lines), P2 (purple lines) and CB07 (green lines).

Open in new tabDownload slide

Recall that Fig. 3 showed how, compared with the other considered SPS models, both MN05 and CB07 SPS models produce more luminous stars in the NIR for stellar populations with ages over 0.5 Gyr. The rest-frame K-band luminosity of stars with ages 0.3 <age(Gyr)<2 is dominated by TP-AGB stars (Maraston 2005). For both the MN05 and CB07 models, this translates into both a shallower slope in the bright end and a change in the global normalization of the luminosity function at z > 1 with respect to that at z = 0. This clear change in normalization is not seen in the other SPS models, even if they explicitly consider the contribution of TP-AGB stars. We find similar trends when using the Lacey et al. (in preparation) model, whose free parameters are tuned using the MN05 SPS model.

We note that the shift in luminosity due to changing the BC99 SPS model by the MN05 one is larger for fainter model galaxies. The shallower slope of the luminosity function in the faint end makes this larger change far less visible than in the bright end.

In order to explore this point further, we have run the CW09 SPS model varying the properties of the TP-AGB stars. Conroy et al. (2009) modelled the uncertainties associated with the TP-AGB phase by allowing a shift in both their luminosity, Δ_L = Δ log(L_bol), and temperature, Δ_T = Δ log(T_eff), with respect to the default tracks (see Table 2). Their comparison to star cluster colours and the fitting of broad-band photometry of blue ∼L_* galaxies at z ∼ 0 suggested that the allowed ranges are −0.2 ≤ Δ_L ≤ 0.2 and −0.1 ≤ Δ_T ≤ 0.1. We have used these four extremes to produce new SSPs to be fed into the semi-analytical model and find that only the bright end of the predicted K-band luminosity function at z > 1 is affected by these changes. Fig. 8 shows the K-band luminosity function at both z = 0 and z=1.5 for the default CW09 SPS model and that modified by increasing the luminosity of the TP-AGB stars by Δ_L = 0.2. It is clear from this plot that the CW09 SPS model modified in this way reproduces better the observed evolution of the K-band luminosity function. However, Conroy et al. (2009) found that a shift of Δ_L = −0.4 was preferred in order to reproduce the colours of star clusters with their SPS model. These different preferred values depending on the data set used for comparison illustrates the uncertainties derived from the difficulty modelling the TP-AGB phase.

Figure 8.

The predicted rest-frame K-band luminosity function at z = 0 (red lines) and z = 1.5 (blue lines) compared to observations. The solid lines show the predictions using the MN05 SPS model, the dotted lines correspond to the predictions from the default CW09 SPS model and the dashed lines show the predictions when the luminosity of the TP-AGB stars is increased in the CW09 model. The observational data are as follows: z = 0 – Driver et al. (2012) (filled circles), Cole et al. (2001) (triangles), Kochanek et al. (2001) (crosses); z = 1.5 – Pozzetti et al. (2003) (filled squares), Caputi et al. (2006) (open circles), Cirasuolo et al. (2010) (diamonds).

Open in new tabDownload slide

Fig. 8 also shows the evolution of the K-band luminosity function predicted using the MN05 SPS model, which is similar to that predicted with the CB07 SPS model, as seen in Fig. 7. There is a clear increase in the luminosity function normalization from z = 0 to z = 1.5 with these two SPS models. The other models, including the modified CW09 SPS one, predict a much smaller change in normalization; for these models most of the evolution occurs at the bright end of the luminosity function. The observations appear to agree with the latter type of evolution.

As discussed in Section 3.1, the TP-AGB stars dominate the NIR light of stars with ages between 0.2 and 4 Gyr, depending on the SPS model. We have explored the predicted mean ages of galaxies when the SPS model is varied. Only small variations are found when we compare the distribution of both stellar mass-weighted ages for galaxies with M_* > 10⁸ h⁻¹ M_⊙ and V-luminosity weighted ages for galaxies with M_AB(V) − 5 log h ≤ −18. At z = 0, when the luminosity functions from both BC99 and MN05 SPS models agree very well, galaxies with M_AB(K) − 5 log h ≤ −18 have median K-weighted ages about 25 per cent younger when using the MN05 SPS model, compared to the BC99 model. This difference increases to 40 per cent for galaxies with M_AB(K) − 5 log h ≤ −22. However, the offset is reduced at higher redshifts. At z = 1.5, the median K-weighted age obtained with the MN05 model is less than 20 per cent younger than for the BC99 model. These differences are likely to be related to the different contribution of intermediate age stars to the integrated NIR luminosity in the different SPS models. The difference in ages cannot completely explain the variation in the NIR luminosity function found at z > 1.

Constraining the TP-AGB phase using the predicted rest-frame K-band luminosity function evolution is not possible due to the uncertainties in other physical processes which can also change this prediction. However, it is still interesting to note the difference between the observed evolution and that predicted with either the MN05 or the CB07 SPS models.

Due to the uncertainties in the modelling of the TP-AGB phase, our results advice against using the rest-frame K-band luminosity function at z > 1 to tune the free parameters in models of galaxy formation and evolution. We have found that, at z > 1, the rest-frame optical luminosity functions are robust against uncertainties in the SSPs models and thus, favoured for tuning the model free parameters.

5 NUMBER COUNTS

The number counts of galaxies depend on how the galaxy luminosity function varies with redshift. Thus, although they are a simple quantity to measure, the number counts contain information about the processes that affect the evolution of galaxies (e.g. Guiderdoni & Rocca-Volmerange 1990; Barro et al. 2009).

Fig. 9 shows the predicted number counts in R, K and [3.6] bands compared with observations The optical and NIR number counts predicted using different SPS models are practically indistinguishable and give a very good match to the observations in the R and K bands (similar results are also found for the i band, λ_eff = 7481 Å), though they are underpredicting by less than a factor of 2 the bright [3.6] number counts.

Figure 9.

Main panels: differential number counts in R band (λ_eff = 6490 Å, left), K band (λ_eff = 2.2 μm, centre) and [3.6] band (λ_eff = 3.6 μm, right). Bottom panels: differential number counts normalized by the expected counts in a Euclidean universe, i.e. log (dn/dm/deg²) − 0.6(m − 12.), with m being the corresponding magnitude. The symbols in the left-hand and central panel correspond to the compilation from Nigel Metcalfe (http://astro.dur.ac.uk/ nm/pubhtml/counts/counts.html). The symbols in the right-hand panel correspond to observations from the three fields in Fazio et al. (2004), the error bars show the Poisson uncertainties. K-band observations: Mobasher, Ellis & Sharples (1986) (dots), Gardner, Cowie & Wainscoat (1993) (diamonds), Soifer et al. (1994) (crosses), Djorgovski et al. (1995) (plus signs), McLeod et al. (1995) (squares), Gardner et al. (1996) (rotated triangles), Moustakas et al. (1997) (filled rotated triangles), Bershady, Lowenthal & Koo (1998) (filled squares), Szokoly et al. (1998) (dilled diamonds), McCracken et al. (2000) (asterisks), Väisänen et al. (2000) (stars), Huang et al. (2001) (inverted triangles), Kümmel & Wagner (2001) (triangles), Maihara et al. (2001) (filled triangles), Saracco et al. (2001) (filled stars), Iovino et al. (2005) (hexagons), Metcalfe et al. (2006) (filled inverted triangles), Imai et al. (2007) (pentagons), Keenan et al. (2009) (filled pentagons) and Cristóbal-Hornillos et al. (2009) (filled hexagons). R-band observations: Koo (1986) (empty diamonds), Infante, Pritchet & Quintana (1986) (rotated empty triangle), Yee & Green (1987) (filled triangles), Tyson (1988) (empty squares), Metcalfe et al. (1991) (asteriks), Picard (1991) (empty pentagons), Jones et al. (1991) (empty triangles), Steidel & Hamilton (1993) (rotated filled triangle), Couch, Jurcevic & Boyle (1993) (upsided-down filled triangles), Driver et al. (1994) (filled pentagons), Metcalfe et al. (1995) (dots), Metcalfe et al. (2001) (filled squares), Huang et al. (2001) (empty stars), Kümmel & Wagner (2001) (crosses), Yasuda et al. (2001) (filled circles), McCracken et al. (2003) (empty circles), Capak et al. (2004) (filled stars).

Open in new tabDownload slide

The faint end of the number counts is expected to be dominated by high-redshift galaxies while the bright end is determined by nearby ones. Thus, the differences seen in the rest-frame K-band luminosity function at z > 1 will start to be appreciable at longer wavelengths in the observer frame.

The faint [3.6]-band number counts are expected to be dominated by the contribution of galaxies at high redshifts, z ∼ 1, and thus, by their rest-frame NIR light. The right-hand panel in Fig. 9 shows that around [3.6]_Vega = 18, there is a difference of a factor of ∼1.6 between the number counts using the MN05 model and the BC99 one. This is in agreement with the different evolution seen in the rest-frame K-band luminosity function at z ≥ 0.8, i.e. where the TP-AGB stars start to dominate the global luminosity.

5.1 The number counts of EROs

Since the predicted number counts of K-selected galaxies are in excellent agreement with observations, we further explore the nature of galaxies that are bright in the NIR and that have been claimed to be particularly sensitive to the modelling of the TP-AGB phase (Fontanot & Monaco 2010; Henriques et al. 2011): the Extremely Red Objects (EROs).

EROs are K-selected galaxies with very red optical to NIR colours, for example (R − K)_Vega > 5. Bright EROs with m_Vega(K) < 20 have inferred stellar masses above 10^10.5  h⁻¹ M_⊙ (Conselice et al. 2008), they appear to be in place at z ∼ 1 (Cimatti et al. 2002) and about half of them are dominated by an old stellar component (Väisänen & Johansson 2004). Such characteristics posed a problem to hierarchical models of galaxy formation and evolution for about a decade (Smith et al. 2002). In hierarchical models, more massive haloes tend to form more recently. However, the formation of galaxy mass does not follow that of the host dark matter haloes, since processes affecting the star formation, gas cooling and galaxy merging can also have an impact on the mass assembly of galaxies.

In Gonzalez-Perez et al. (2009) and Gonzalez-Perez et al. (2011), we showed that the predicted numbers and characteristics of the EROs depend strongly on the feedback processes which suppress the formation of massive galaxies, in particular the AGN feedback, and on the star formation time-scale. We presented a model, Bower et al. (2006) with the BC99 SPS, that matched the observed ERO numbers. Both Fontanot & Monaco (2010) and Henriques et al. (2011) found that a galaxy formation model coupled with the MN05 SPS model could also reproduce the observed number counts of EROs.

Fig. 10 shows the predicted cumulative K-band number counts for EROs selected by their (R − K) colour, compared with different observations. The ERO number counts predicted using the BC99, PD01, P2 or CW09 SPS models are very similar (for clarity we have included in the figure only the default case, the BC99 SPS model). The semi-analytical model coupled with all of these SPS models underpredicts the observed number counts by a factor of ∼2.5 for EROs with (R − K)_Vega > 5 and by more than an order of magnitude for EROs with (R − K)_Vega > 7, though these differences can be alleviated by modifying the colour cut bluewards by 0.2 and 0.5 mag, respectively. Note that a large fraction of the dispersion found among observations is likely related to the sensitivity of the colour cut to the choice of filters and to aperture effects when estimating colours (Simpson et al. 2006). For the intermediate colour cuts, these models match some of the observational data sets, though not those covering the largest observed areas (Lawrence et al. 2007; Palamara et al. 2013).

Figure 10.

The predicted cumulative number counts of EROs selected by their (R − K) colour, as indicated in each panel. The solid black lines correspond to the BC99 SSP model and the blue solid lines to the MN05 one. The dotted lines show the BC99 SPS model counts without attenuation by dust. The dashed black lines show the predicted number counts when the BC99 SPS model is used in combination with the old prescription for the star formation law. The symbols correspond to observations from Daddi, Cimatti & Renzini (2000) (squares), Smith et al. (2002) (filled triangles), Smail et al. (2002) (empty circles), Roche et al. (2002) (empty triangles), Väisänen & Johansson (2004) (diamonds), Kong et al. (2006) (crosses), the UKIDSS survey (Simpson et al. 2006; Lawrence et al. 2007, filled circles) and Palamara et al. (2013) (stars). The errors shown are Poisson.

Open in new tabDownload slide

The mismatch between the default predictions and observations is larger than expected, considering that the AGN feedback in the model presented in this work is very similar to that in the Bower et al. (2006) model. The time-scale of starbursts has been increased from the Bower et al. model to the one presented here. Nevertheless, this change has a minimal impact on the EROs number counts.

One of the main changes introduced in the model by Lagos et al. (2011b) is an improved treatment of star formation in discs. In our work, we make use of this calculation of the star formation rate in discs, in which the star formation rate surface density is assumed to be proportional to the molecular hydrogen surface density in galaxies. When we run the model using the old implementation for calculating the star formation rate in discs, we obtain number counts that adequately reproduce the observational data and even overpredict the number of EROs with (R − K)_Vega > 7. These results are closer to those reported in Gonzalez-Perez et al. The K-band luminosity function predicted at z = 0 with the old star formation law in discs although consistent with observation, it has a break at fainter luminosities than observed.

As in Gonzalez-Perez et al. (2009), we find that using the old star formation law most of the predicted EROs are passively evolving and had their last burst of star formation more than 1 Gyr ago for both star formation implementations (the exact split varies with magnitude). However, with the new star formation implementation the fraction of EROs that have experienced a recent burst of star formation increases.

The comparison between the number counts calculated with the two implementations of the star formation rate law shows that, besides the AGN feedback, the ERO counts depend on the way star formation is modelled, since this has a strong impact on the numbers of galaxies that experience significant bursts of star formation at different cosmic times. A basic reason for the difference might be the average slower depletion of gas by quiescent star formation in the new star formation prescription compared with the old prescription. It is important to note that there are compelling observations from z = 0 to z ∼ 3 supporting the new star formation implementation (Bigiel et al. 2011; Tacconi et al. 2013). Moreover, some observations indicate that massive and passively evolving galaxies (as most of the predicted EROs are) could have gas depletion time-scales even longer than the 2 Gyr assumed in our model (Saintonge et al. 2011; Martig et al. 2013) and, thus, even further away from the old one.

When using the MN05 SPS model, the match between the predicted and observed ERO number counts is somewhat improved, with respect to the other SPS models, except for the most extreme case of EROs with (R − K)_Vega > 7. This is in agreement with previous results from Fontanot & Monaco (2010) and Henriques et al. (2011). We have also explored the prediction obtained when using a modified version of the CW09 SPS model similar to that described in Section 4.3, for which we have increased the luminosity of the TP-AGB stars. This change has a visible impact increasing the predicted cumulative number counts of EROs. However, this increase is still too small to match the predictions using the MN05 SPS model, in agreement with the different evolution of the K-band luminosity function seen in Fig. 8 using the MN05 and modified CW09 SPS models.

The predicted number counts without including the attenuation by interstellar dust agree remarkably well with observations, except for EROs with (R − K)_Vega > 7. Thus, uncertainties in the treatment of dust attenuation might also be partially responsible for the mismatch between predictions and observations.

6 CONCLUSIONS

We have presented a new hierarchical model of galaxy formation and evolution based on the MS-W7 simulation (Lacey et al. in preparation). This run has the same specifications as the Millennium Simulation but uses a WMAP7 cosmology (Komatsu et al. 2011), which is close to Planck cosmology for the scales relevant to galaxy evolution [Ade et al. (Planck Collaboration) 2013]. This new model is a development of the galform semi-analytical model, which includes physical treatments of the hierarchical assembly of dark matter haloes, shock-heating and cooling of gas, star formation, feedback from SNe, AGN and photoionization of the IGM, galaxy mergers and chemical enrichment. The luminosities of galaxies are calculated from an SPS model and dust attenuation is then included using a self-consistent theoretical model based on the results of a radiative transfer calculation.

The model presented here is an extension of the Lagos et al. (2011b) model which introduced a new treatment of star formation which follows the atomic and molecular hydrogen in the ISM. The free parameters in the model have been chosen in order to reproduce the rest-frame luminosity functions in b_J and K bands at z = 0 with the WMAP7 cosmology, and to predict a reasonable evolution of the luminosity function in both the rest-frame UV and K bands.

The variation in the model predictions arising from the choice of SPS model gives an indication of how accurately different properties can be predicted. Thus, with this new model we have explored how sensitive the predicted luminosity function evolution from the rest-frame UV to the NIR wavelength range is to the particular choice of SPS model. In order to do this, we have run our galaxy formation and evolution model coupled with seven different SPS models: PD01 (Silva et al. 1998), P2 (Fioc & Rocca-Volmerange 1997, 1999), BC99 (an updated version of Bruzual & Charlot 1993), BC03 (Bruzual & Charlot 2003), MN05 (Maraston 2005), CB07 (Bruzual 2007) and CW09 (Conroy et al. 2009, 2010; Conroy & Gunn 2010). Our default model uses the BC99 SPS model and a Kennicutt IMF. All of the studied SPS models cover a large range in wavelength, at least from the UV to the NIR, and most are publicly available. In the wavelength range of interest here, and for stellar ages above 10 Myr, the main differences between the spectra of SSPs obtained with different SPS models arise from stars in post-main-sequence phases, in particular, those in a thermally pulsating asymptotic giant branch (TP-AGB) phase. This is a phase which is particularly difficult to model due to its variability, including uncertain mass-losses by winds. There are observations both supporting and disfavouring the modelling of a strong TP-AGB phase (e.g. Lyubenova et al. 2012; Zibetti et al. 2013). From the SPS models explored, the MN05 and CB07 ones include a strong TP-AGB phase and the CW09 one allows changes in the intensity of this phase.

We found that the predicted rest-frame UV luminosity function is insensitive to the choice of SPS model. Compared with observations, we have found that our model overpredicts the number of bright rest-frame UV galaxies at z = 0. At higher redshifts, z ≤ 6, the model reproduces the observed rest-frame UV luminosity function reasonably well.

The predicted rest-frame optical luminosity function is also found to be insensitive to the choice of SPS model. In this case, the model agrees remarkably well with observations up to z ∼ 3.

We have found that the evolution of the rest-frame NIR luminosity function strongly depends on the particular treatment of the TP-AGB phase done in the SPS models. When we couple the galaxy formation and evolution model with the two SPS models which include the strongest contribution from TP-AGB stars, the MN05 and CB07 models, the predicted luminosity function strongly evolves from z = 0 to z = 1.5, with a change in the magnitude corresponding to L_* galaxies. For the other SPS models, even when we artificially increase the rest-frame luminosity of TP-AGB stars in the CW09 model, we only find a change at the bright end of the rest-frame NIR luminosity function. Similar trends are found when using the Lacey et al. (in preparation) model, whose free parameters are tuned using the MN05 SPS model. The observations of the luminosity function evolution are in better agreement with an evolution of the rest-frame NIR luminosity function mainly affecting its bright end. When using the MN05 model, we predict galaxies at z = 1.5 to have slightly smaller K-weighted ages. The difference in K-weighted ages increases towards lower redshifts, being 20 per cent at z = 0 for galaxies brighter than M_K − 5 log h < −18 (the percentage increases for brighter cuts). At z = 0, despite this difference in K-weighted ages, the NIR luminosity functions predicted with the different SPS models agree. Thus, the difference in ages appears to be related with the specific contribution of intermediate age stars to the global NIR luminosity when using either the BC99 or the MN05 SPS models, rather than the origin of the difference of the predicted rest-frame NIR luminosity function. Our results suggest that, once other set of observations are also used to constrain the model, the predicted evolution of the rest-frame NIR luminosity function might help in constraining the strength of the TP-AGB phase.

The number counts of galaxies depend on the evolution of the observed luminosity function. Thus, as a further exploration we have compared the predicted number counts from the optical to the NIR wavelength range. We find very good agreement between observations and predictions. Only at wavelengths as long as 3.6 μm, we start to see the impact of TP-AGB stars on the global galactic number counts.

Galaxies selected by their red optical to NIR colours are particularly sensitive to the implementation of the TP-AGB in the SPS models. However, we find that the number counts of galaxies with (R − K) > 5 are more sensitive to the treatment of the star formation in discs than to the strength of the TP-AGB phase.

In conclusion, we find that we can use the predicted rest-frame UV to optical luminosity functions without being hindered by uncertainties in the particular choice of SPS model done within the galaxy formation and evolution model. The rest-frame NIR luminosity function is more problematic, since it depends on the particular treatment done in the SPS model to account for TP-AGB stars.

We would like to thank Andrew Benson, Richard Bower, Shaun Cole and Carlos Frenk for their help developing galform; Stephane Charlot, Charlie Conroy, Michel Fioc, Claudia Maraston and Brigitte Rocca-Volmerange for their helpful comments on the SPS models; Olivier Ilbert and David Palamara for providing us with tabulated observational data; Peder Norberg for making the model runs for the Millennium Archive Data Base. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. We acknowledge support from the Durham STFC rolling grant in theoretical cosmology. VGP acknowledges financial support from the Agence Nationale de la Recherche OMEGA grant ANR-11-JS56-003-01, past support from the UK Space Agency and logistic support from St Andrews University. DJRC acknowledges support from the Royal Astronomical Society Grant for a summer studentship. VGP and DJRC thank Alastair Edge for his help getting this studentship.

REFERENCES