We present predictions for the abundance of submillimetre galaxies (SMGs) and Lyman-break galaxies (LBGs) in the Λ cold dark matter cosmology. A key feature of our model is the self-consistent calculation of the absorption and emission of radiation by dust. The new model successfully matches the LBG luminosity function, as well as reproducing the properties of the local galaxy population in the optical and infrared. The model can also explain the observed galaxy number counts at 850 μm, but only if we assume a top-heavy initial mass function for the stars formed in bursts. The predicted redshift distribution of SMGs depends relatively little on their flux over the range 1–10 mJy, with a median value of z≈ 2.0 at a flux of 5 mJy, in good agreement with the recent measurement by Chapman et al. The counts of SMGs are predicted to be dominated by ongoing starbursts. However, in the model these bursts are responsible for making only a few per cent of the stellar mass locked up in massive ellipticals at the present day.
The detection of populations of star-forming galaxies at high redshift (z≳ 2) in the late 1990s opened a window on the process of galaxy formation when the Universe was less than 20 per cent of its current age. Two main techniques have been used to discover such objects. (1) Lyman-break galaxies (LBGs) are detected in optical bands from their stellar emission in the rest-frame ultraviolet, using the spectral break around 912 Å produced by absorption by intervening neutral hydrogen (Steidel et al. 1996). The largest samples of LBGs have been obtained for z∼ 3–4 (Steidel et al. 1999). The dust extinction corrections required to derive star formation rates of LBGs from their rest-frame ultraviolet (UV) luminosities are probably quite large (a factor of ∼3–10) (e.g. Meurer, Heckman & Calzetti 1999; Adelberger & Steidel 2000), but remain uncertain, this being a major obstacle preventing an accurate determination of the contribution of LBGs to the global star formation rate. (2) Submillimetre galaxies (SMGs) are detected by emission from warm dust in the rest-frame far-infrared (IR)/submillimetre bands (Smail, Ivison & Blain 1997; Hughes et al. 1998). From their discovery using the Submillimetre Common User Bolometer Array (SCUBA) instrument on the James Clerk Maxwell Telescope (JCMT), the faint mJy-flux SMGs were suspected to lie at high redshifts (e.g. Yun & Carilli 2002) and this has now been firmly established by measurements of spectroscopic redshifts, which give a median redshift z≈ 2.2–2.4 for sources with 850-μm fluxes around 5 mJy (Chapman et al. 2003, 2004). The SMGs are generally interpreted as being dust-enshrouded starbursts, in which the dust is being heated by UV radiation from young stars, but estimates of their star formation rates, based on the measured 850-μm fluxes, are uncertain because of a lack of knowledge of the shapes of their spectral energy distributions (SEDs) at the shorter wavelengths which presumably dominate the bolometric luminosity of the dust. There also remains the possibility that a significant fraction of the submillimetre emission may be powered by heating of dust by active galactic nuclei (AGN) rather than young stars. Recent X-ray observations suggest that the AGN contribution to the submillimetre emission is small in most SMGs (Alexander et al. 2003), but this issue will only be finally settled by further observational data.
A third technique that has been used to search for high-redshift star-forming galaxies is to look for the Lyα emission line (e.g. Hu, Cowie & McMahon 1998). The usefulness of Lyα emission as a quantitative star formation indicator is, however, limited by the fact that many star-forming galaxies do not show Lyα in emission (e.g. Steidel et al. 2000), and that resonant scattering of Lyα photons by neutral hydrogen makes corrections for dust extinction even more uncertain than for the UV continuum. Thus, more reliable estimates of the star formation rates in these systems typically use measurements of the rest-frame UV luminosity, as in LBGs. We present predictions from our semi-analytical model for Lyα-emitting galaxies in a separate paper (Le Delliou et al. 2004), but do not consider them further here.
Within the uncertainties, galaxies detected as LBGs and as SMGs appear to contribute comparably to the global star formation rate at high redshifts (e.g. Smail et al. 2002). In combination, they should, in principle, provide a census of all high-mass star formation at redshifts z≳ 2, as one should either see the UV radiation from these stars directly (as in LBGs) or after reprocessing by dust (as in SMGs). It is therefore of fundamental importance to understand both types of galaxy, in order to understand the process of galaxy formation as a whole. In this paper, we describe the main results of our attempt to do this within the framework of the cold dark matter (CDM) model of structure formation; subsequent papers will describe different aspects of our model in more detail.
We build on our previous work using semi-analytical models of galaxy formation to calculate the evolution of the galaxy population ab initio. These models predict the mass accretion, merging and star formation history of each galaxy, along with the evolution of its gas mass and metallicity, and the mass and size of the disc and bulge. In Baugh et al. (1998) we interpreted the first observational results on LBGs in terms of a hierarchical model of galaxy formation, but without including any treatment of dust. In Granato et al. (2000, hereafter G00), we presented a new, multiwavelength version of our model which included a comprehensive treatment of the reprocessing of stellar radiation by dust (using the grasil code developed by Silva et al. 1998, hereafter S98), and showed that it reproduced the luminosity functions and other properties of galaxies at z= 0 all the way from the far-UV to the far-IR. However, as we show in this paper, when the G00 model in its original form is applied to the high-z universe, it dramatically underpredicts both the number density of bright LBGs and the number counts of SMGs at 850 μm. Of these, the latter problem seems the most serious, as the observed submillimetre fluxes of SMGs are usually interpreted as requiring objects with star formation rates in excess of 100 M⊙ yr−1 (for a standard initial mass function, IMF) to be commonplace; such vigorous star formation rates at high redshifts are difficult to accommodate in the CDM model. To see if it is possible to achieve such high submillimetre luminosities in high-z galaxies, we have therefore made significant modifications to our original model, in regard to both star formation and the behaviour in bursts, in order to arrive at a new self-consistent model, which simultaneously reproduces the properties of the LBGs and SMGs at high redshift and the main properties of the galaxy population at low redshift. Here we present the first results from this new model.
The challenge of modelling galaxies in the high-redshift universe has received much attention over the past decade. Simple theoretical and empirical models have been presented which use the number counts of galaxies observed at different wavelengths to try to constrain the global star formation history of the Universe (e.g. Pearson & Rowan-Robinson 1996; Blain et al. 1999a,b; Chary & Elbaz 2001; Balland, Devriendt & Silk 2003). More sophisticated models that attempt to follow the physics of galaxy formation in greater detail have stressed the importance of bursts of star formation at high redshift in explaining the abundance of SMGs (Guiderdoni et al. 1998; Devriendt & Guiderdoni 2000) and LBGs (Somerville, Primack & Faber 2001). Most recently, Granato et al. (2004) have developed a coupled model for the evolution of galactic spheroids and quasi-stellar objects (QSOs), which successfully reproduces the counts of SMGs and the evolution of the K-band luminosity function with redshift (e.g. Pozzetti et al. 2003; Drory et al. 2004).
Our methodology represents a number of significant advances over earlier attempts to model the SMGs and LBGs within the framework of hierarchical galaxy formation. Previous work has tried to model either SMGs (e.g. Guiderdoni et al. 1998; Blain et al. 1999b) or LBGs (e.g. Baugh et al. 1998; Somerville et al. 2001), but not both populations together. This allowed them to tune their model parameters to fit either one population or the other. In this paper, we try to explain both populations together. Our model has the following important features that improve on earlier work.
It directly predicts the merger histories of dark matter haloes and of the galaxies contained within them. This allows us to predict the redshifts and gas masses of bursts of star formation triggered by galaxy mergers.
The absorption of starlight by dust is calculated self-consistently by radiative transfer, using the dust mass obtained from the gas mass and metallicity, and the disc and bulge scalelengths predicted by the model. This is particularly important when computing the amount of extinction at short wavelengths. In turn, this determines the amount of energy absorbed by the dust and re-radiated at long wavelengths.
The spectrum of dust emission is calculated using a distribution of dust temperatures within each galaxy, which is obtained by solving for the radiative equilibrium temperatures of individual dust grains.
In contrast to our detailed treatment of dust emission, previous works either treated the dust temperature as a free parameter in fitting the SMG number counts (e.g. Kaviani, Haehnelt & Kauffmann 2003), or used empirical IR/submillimetre SEDs estimated from present-day galaxies (e.g. Pearson & Rowan-Robinson 1996; Blain et al. 1999a; Devriendt & Guiderdoni 2000). The task of fitting the observed number counts of SMGs is much simpler if the dust temperature is treated as a free parameter: for a dust emissivity εν∝νβ, bolometric dust luminosity Ld and temperature Td, the luminosity per unit frequency at long wavelengths scales as Lν∝ν2+βLd/T3+βd, which varies roughly as T−5d at fixed Ld for β≈ 2. This might suggest that it is easy to reconcile any theoretical model with the observed submillimetre counts simply by making modest adjustments to the assumed dust temperature. However, treating the dust temperature as a free parameter is physically inconsistent, as it is actually determined by the condition that the dust grains should be in thermal equilibrium between radiative heating and cooling. For dust mass Md, with a single dust temperature, thermal equilibrium implies Ld∝MdTβ+4d, so that Lν∝ν2+βM(3+β)/(4+β)dL1/(4+β)d if Md and Ld are treated as the physical parameters, rather than Td and Ld. This argument shows that once one imposes the condition that dust temperatures should be physically consistent, the predicted submillimetre luminosities in a model can only be changed by large factors by making large changes to the dust masses.
Additional shortcomings of previous works are that they either did not include halo and galaxy mergers (e.g. Guiderdoni et al. 1998), or treated the amount of dust extinction in high-z galaxies as a free parameter (e.g. Somerville et al. 2001).
We now give a brief overview of our new galaxy formation model. In Section 2.1, for completeness, we list some general differences from the model presented by G00, which were motivated by changes in the preferred values of the cosmological parameters and by improvements in the modelling of feedback processes (Benson et al. 2003; hereafter B03). We demonstrate in Section 3 that neither the G00 model nor the model presented by B03 can reproduce the number counts of submillimetre galaxies. The additional changes required to the model of B03 to rectify this are set out in detail in Section 2.2. The impact of these changes on the predicted submillimetre counts is demonstrated in Fig. 5 in Section 3. Full details of the semi-analytic galaxy formation model can be found in Cole et al. (2000) and B03. The grasil model for the emission from stars and dust is described by S98, and its implementation in the semi-analytic model is described in G00. A full description of our new model will be given in Lacey et al. (Paper 1, in preparation).
The basic model
The galaxy formation model from which we start is the one described by B03. This contains a number of changes from the model of G00. We assume a flat, ΛCDM cosmology, with density parameter Ω0= 0.3, a Hubble constant of h=H0/(100 km s−1 Mpc) = 0.7 and with a perturbation amplitude set by the linear rms fluctuation in spheres of radius 8 h−1 Mpc of σ8= 0.9.
One of the major differences is motivated by a revision in the value of the cosmological baryon density, Ωb, to be consistent with recent constraints from big bang nucleosynthesis models (O'Meara et al. 2001) and from measurements of fluctuations in the cosmic microwave background radiation (e.g. Spergel et al. 2003). We adopt a value of Ωb= 0.04, which is twice that used by G00. To prevent the formation of too many luminous galaxies, B03 found it necessary to invoke a ‘superwind’ mode of feedback, in which cold gas is ejected from a galaxy in proportion to the star formation rate. The superwind operates in addition to the feedback used by G00 (see Cole et al. 2000). There is observational support for the existence of superwinds from the detection of high-velocity outflows in massive galaxies at both low and high redshift (e.g. Martin 1999; Heckman et al. 2000; Pettini et al. 2002; Adelberger et al. 2003; Smail et al. 2003).
The B03 model also includes a simple treatment of the feedback from photoionization of the intergalactic medium, which suppresses the collapse of gas into low-mass dark haloes after reionization. We model this by assuming that no gas cools in haloes with circular velocities Vc < 60 km s−1 at z < 6. This has a similar effect on the predicted form of the galaxy luminosity function to the more realistic and sophisticated treatment developed by Benson et al. (2002b).
The grasil code computes the emission from both the stars and dust in a galaxy, based on the star formation and metal enrichment histories predicted by the semi-analytical model. It includes radiative transfer through a two-phase dust medium, with a diffuse component and giant molecular clouds, and a distribution of dust grain sizes. Stars are assumed to form inside the clouds and then gradually leak out. The output from grasil is the galaxy SED from the far-UV to the submillimetre. The motivation for the choice of dust parameters in the grasil model and the impact of varying these choices are discussed in detail in G00 (see also Prouton et al. 2004).
Key features of the new model
In Section 3, we will show that the predictions of the G00 and B03 models do not match the observed number counts of submillimetre sources. The positive effect on the number counts of increasing the baryon fraction by a factor of 2 between G00 and B03 is offset by the inclusion of superwinds in the latter. It is therefore necessary to make a number of important changes to the B03 model in order to reproduce the SMG counts. Our freedom to change the B03 model is, however, restricted for two reasons. First, we retain the cosmological parameters adopted by B03. This sets the rate at which structure grows in the dark matter, which in turn determines the halo merger trees in which galaxies form and evolve. Secondly, we follow the philosophy set out by Cole et al. (2000), and demand that any new model should reproduce basic properties of the local galaxy population, such as the luminosity functions in the optical, near-IR and far-IR, the distribution of disc scale sizes, the gas and metal contents of discs and spheroids, and the mix of morphological types. The task of finding a successful model is made even more demanding by our use of a self-consistent calculation of the absorption and re-emission of radiation by dust. We emphasize again that the temperature distribution of the dust grains is not a free parameter in our model, but is determined by the stellar luminosity and the mass, distribution and composition of the dust. We largely retain the same parameters in the grasil code as used by G00, with modifications to the escape time of young stars from dust clouds and to the long-wavelength emissivity of the dust grains, which we discuss in the next sections.
The model set out below was reached after an extensive exploration of parameter space, within the boundaries set by the stringent constraints outlined above, with the goal of improving the model predictions for the high-redshift galaxy population. The consequences of systematically varying key model assumptions will be illustrated in Section 3.
Star formation time-scale in discs. The star formation rate is given by ψ=Mgas/τ★, where Mgas is the cold gas mass. G00 and B03 assumed that the star formation time-scale in discs, τ★, incorporated a scaling with the dynamical time of the disc. In this case, the star formation time-scale becomes shorter with increasing redshift. Consequently, mergers at high redshift tend to be between gas-poor discs, with little fuel available to power a starburst. In order to move more star formation from the quiescent to the burst mode at high redshifts, we adopt a prescription for the star formation time-scale in discs without the dynamical time scaling:Cole et al. 1994; Baugh, Cole & Frenk 1996; Baugh et al. 1998). There are no compelling theoretical arguments for preferring one of these prescriptions over the other. (We test the star formation recipe in our models against observational data in Bell et al. 2003.) In both cases, appropriate parameter values can be chosen to reproduce the observed gas fraction–luminosity relation at z= 0 (see fig. 9 of Cole et al. 2000); here, we choose τ★0= 8 Gyr and α★=−3. In the new model, galactic discs at low redshift have roughly the same star formation time-scales as in G00, but discs at high redshifts have much longer star formation time-scales than before. Consequently, in the new model discs at high redshift are much more gas-rich than in G00, resulting in more gas being available for star formation in starbursts triggered by mergers between these discs. This type of picture is hinted at observationally by the very large gas fractions implied by CO measurements of SMGs (e.g. Neri et al. 2003). Similar schemes resulting in gas-rich mergers at high redshift have previously been invoked to explain the luminosity function of Lyman-break galaxies (Somerville et al. 2001) and the evolution of the quasar population (Kauffmann & Haehnelt 2000).
Star formation bursts triggered by galaxy mergers. Bursts of star formation are assumed to be triggered only by galaxy mergers. We define a galaxy merger as major or minor according to whether the ratio of masses of the two galaxies exceeds fellip or not. In major mergers, the stellar discs are assumed to be transformed into a new stellar spheroid, while in minor mergers, the stellar disc of the larger galaxy is preserved. Bursts are assumed to be triggered in all major mergers (as in G00), but also (unlike in G00) in minor mergers that satisfy both of the following conditions: (i) the galaxy mass ratio exceeds fburst (where fburst < fellip); (ii) the gas fraction in the primary galaxy exceeds fgas,crit. The latter condition is motivated by the suggestion of Hernquist & Mihos (1995) that gas-rich discs should be more susceptible to bursts triggered by the accretion of small satellites. We adopt values of fellip= 0.3, fburst= 0.05 and fgas,crit= 0.75. Another modification is that the time-scale of star formation in a burst is taken to be τ⋆burst= max [fdynτdyn, τ⋆burst,min], where τdyn is the dynamical time of the newly formed spheroid (G00 assumed τ⋆burst=fdynτdyn). We adopt fdyn= 50 and τ⋆burst,min= 0.2 Gyr; these values were found to produce the best match to the present-day 60-μm luminosity function and the abundance of SMGs.
The IMF in bursts. Perhaps the most dramatic change from our previous work is the adoption of a top-heavy IMF for stars formed in bursts. This is the single most important change to the model, which has the biggest impact on the predicted counts of submillimetre sources. For stars formed quiescently in discs, we use a solar neighbourhood IMF, with the form proposed by Kennicutt (1983): dN/d ln m∝m−x, with x= 0.4 for m < 1 M⊙ and x= 1.5 for m > 1 M⊙. In bursts, we instead use a single power law with x= 0. In both cases, the IMF covers a mass range 0.15 < m < 125 M⊙. (In G00, we used a Kennicutt IMF for all star formation.) There are two factors motivating the adoption of a top-heavy IMF in starbursts. First, the total energy radiated in the UV per unit mass of stars formed is increased due to the larger fraction of high-mass stars, which increases the amount of radiation heating the dust. For example, the energy per unit mass emitted at 1500 Å is around four times larger for the top-heavy x= 0 IMF compared with the Kennicutt IMF. Secondly, a top-heavy IMF also results in a higher yield of metals from Type II supernovae, so that more dust is produced to absorb this energy. For stars formed with the top-heavy IMF, more than twice the mass is recycled to the interstellar medium and over six times as many metals are produced, compared with the case of stars formed with a Kennicutt IMF. The increased production of dust for a top-heavy IMF is essential for boosting the luminosity of galaxies in the submillimetre range, as discussed in the Introduction. If the stellar UV luminosity were to be increased without changing the mass of dust, the dust would be heated to higher temperatures, shifting the peak of the dust emission spectrum to shorter wavelengths, and resulting in relatively little increase in the luminosity at wavelengths longwards of the peak (including the submillimetre). Including the increased dust production, using a top-heavy IMF results in the most efficient production of luminosity at submillimetre wavelengths for a given rate of star formation.
We first make some general comparisons between the predictions of the G00 model and our new model, before focusing on the properties of galaxies at high redshift.
Cosmic star formation history. We compare the global star formation histories in G00 and the new model in Fig. 1. The overall star formation rate in the new model (Fig. 1a) is similar to that in G00 (Fig. 1b). There is a broad peak in the total star formation density at z≈ 2–3, corresponding to a look-back time of ∼10 Gyr. However, in contrast to G00, where star formation in bursts is around 5 per cent of the total at all redshifts, in the new model the fraction of star formation occurring in bursts increases sharply with redshift, and exceeds that in quiescent discs for z≳ 3. This increase in the importance of bursts at high z relative to G00 results both from the change in the star formation time-scale in quiescent discs, which causes them to be more gas-rich at high z, and from allowing bursts to be triggered by minor mergers. However, the fraction of star formation in bursts summed over the history of the Universe remains modest. Only 30 per cent of star formation takes place in the burst mode in the new model. The contribution of bursts to the mass locked up in stars is even less impressive. Taking into account the large fraction of mass recycled by dying stars with the top-heavy IMF, less than 7 per cent of the mass locked up in stars today is predicted to have been made in bursts.
Present-day galaxy luminosity function. The present-day galaxy luminosity function plays an essential role in constraining the parameters of the galaxy formation model. Our philosophy is to consider only those models that produce reasonable matches to the optical galaxy luminosity function, tracing the stellar population, and the far-IR luminosity function, tracing the emission from dust. Fig. 2 compares the predictions of the new model with the observed luminosity functions at z= 0 in the optical (B band) and in the far-IR (60 μm). In Fig. 2 (and Fig. 3), the errorbars on the model curves represent the statistical uncertainties resulting from the finite sample of galaxies computed; they are calculated as described in G00. The new model reproduces the local optical and far-IR luminosity function data almost as well as the model of G00. As in G00, ongoing bursts dominate the bright end of the 60-μm luminosity function but are less important in the optical range. The new model also reproduces other properties of the local galaxy population, such as the ratio of gas mass to luminosity as a function of luminosity and the stellar metallicity versus luminosity relation (see figs 9 and 10 from Cole et al. 2000, which are very similar to the results from our new model). Further comparisons of the new model with observed properties of present-day galaxies will be given in Lacey et al. (Paper 1, in preparation).
Lyman-break galaxies. A very important prediction of the model is the galaxy luminosity function in the rest-frame UV at high redshift. This depends both on the distribution of galactic star formation rates and on the amount of extinction by dust. We plot the rest-frame UV luminosity function at z= 3 in Fig. 3, comparing the model predictions with observational data for LBGs. The top panel shows the predictions of our new model, while the lower panel shows the fiducial model of G00. The figure shows that the G00 model would match the observed UV luminosity function if there were no dust extinction, but once the effects of dust have been included in a consistent way, the G00 model is typically 2 mag too faint at a given galaxy abundance. In the new model, the dust-extinct UV luminosity function matches the observations well. The luminosity function would be brighter by 1.5–2.5 mag or a factor 4–10 in luminosity in the absence of dust. The level of extinction in our model galaxies is similar to that inferred by Steidel et al. (1999) using a simple extinction law and the observed G−R colours of their Lyman-break sample. However, a detailed comparison requires reproducing the observational LBG selection in our model, and we defer this to a future paper (Lacey et al., Paper 2, in preparation).
As, on average, only a small fraction of the UV light escapes from galaxies with large intrinsic UV luminosities, the dust-extinct UV luminosity function is sensitive to the parameters of the dust model. In the new model, we have adopted a value for the escape time of stars from molecular clouds of tesc= 1 Myr for both bursts and quiescent star formation in order to match the observed luminosity function of LBGs. For this choice, ongoing bursts dominate the extinct UV luminosity function at high z. However, our current choice of dust parameters may not be the only one that allows a match to the observational data. The submillimetre emission from these galaxies is much less sensitive than the UV emission to the details of the dust model, such as cloud masses, radii and escape times, provided that most of the UV light is absorbed. The impact of varying key assumptions of the new model on the abundance of Lyman-break galaxies will be explored in a forthcoming paper, in which we will also explore in detail the connection between the submillimetre fluxes and rest-frame UV luminosities of high-z galaxies, and between SMGs and LBGs (Lacey et al., Paper 2, in preparation).
Submillimetre galaxies. The other crucial constraint on high-redshift galaxies comes from the number counts of faint SMGs. Fig. 4 shows the cumulative number counts as a function of flux for galaxies selected by their emission at 850 μm. The model predictions are compared with a compilation of data obtained by different surveys using SCUBA. The model of G00, shown in Fig. 4(b), is seen to underpredict the observed counts, by a factor ∼20 at 3 mJy. The new model presented in this paper is much more successful (Fig. 4a), agreeing spectacularly well with the observations. To obtain a better fit to the number counts of SMGs, we have modified the emissivity of dust grains in bursts at wavelengths λ > 100 μm, from the canonical εν∝ν2 (which we retain for quiescent star formation, and which G00 used for both modes of star formation), to εν∝ν3/2. There is independent evidence that the dust emissivity in some starbursts may be flatter than εν∝ν2: S98 found that modelling the SED of Arp 220 favoured εν∝ν1.6. As reviewed by Draine (2003), simple theoretical models of dust grains predict that the dust emissivity should vary as εν∝ν2 at long enough wavelengths, independently of the grain size distribution. However, laboratory studies of some possible dust materials find a much shallower dependence than this even at submillimetre wavelengths. Assuming an emissivity such as εν∝ν3/2 thus remains controversial. We illustrate the impact on our predictions of reverting to the canonical choice of εν∝ν2 later on in this section.
In the new model, the counts are dominated by ongoing bursts at fluxes 0.1 ≲Sν(850 μm) ≲ 30 mJy. Bursts account for a factor of 103 times more sources at 3 mJy in the new model than they did in G00; the counts from sources that are not experiencing an ongoing burst have increased by a more modest factor, ≈1.5–6, depending on the flux.
We now investigate the impact on the counts of SMGs of systematically switching off, one at a time, the changes made in the new model. In Fig. 5, the total counts predicted by the new model are shown by the solid line, which is reproduced in each panel. The total counts in the variant models are shown by the dashed line in each panel. Note that in this experiment, we have not altered any additional parameters beyond those stated; however, in all cases, the variant models predict optical luminosity functions at z= 0 that are in reasonably good agreement with the observational estimates. In Fig. 5(a), we illustrate the effect of using a standard Kennicutt IMF in both modes of star formation. This has the biggest impact on the model predictions, reducing the counts by a factor of ≈60 at 3 mJy. The resulting total counts are similar to those predicted by the G00 model. Changing the emissivity of the dust in bursts to be the same as in the quiescent mode, εν∝ν2, has a much less dramatic effect, reducing the counts by less than a factor of 2 at 3 mJy (Fig. 5b). Given the observational scatter, this variant could be a viable alternative model, especially if dust-obscured active galactic nuclei made a comparable contribution to the 850-μm emission to that of dust-obscured starbursts (although Alexander et al. 2003 argue that this is unlikely to be the case).
The importance of bursts of star formation triggered by minor mergers is illustrated in Fig. 5(c), in which we show how the predicted counts change when minor merger bursts are turned off by setting fburst=fellip. The counts at 3 mJy are down by a factor of ≈7 without the contribution of minor merger bursts. Finally, in Fig. 5(d) we show how the predicted counts are affected when we relax the condition for triggering bursts in minor mergers that the gas fraction in the primary galaxy should exceed fgas,crit, by setting fgas,crit= 0. This variant model slightly overpredicts the 850-μm source counts. However, this model greatly overpredicts the number of bright 60-μm galaxies at the present day. As a direct consequence of the parametrization of the star formation time-scale (discussed in Section 2.2 above), galactic discs at high redshift tend to be gas-rich in the new model, whereas they become progressively more gas-poor as z= 0 is approached, through the depletion of gas by quiescent star formation. The gas fraction threshold therefore has the effect of suppressing most minor merger bursts at low redshift, but has relatively little effect at high redshift.
In Fig. 6 we show what the new model predicts for the redshift distributions of sources selected at 850 μm. The median redshift is seen to vary little over the flux range 1 < Sν(850 μm) < 8 mJy, from z50= 2.3 at the faint end to z50= 1.7 at the bright end. At a flux of 5 mJy, the median redshift is predicted to be 2.0, in good agreement with the recent estimate by Chapman et al. (2003, 2004) from a sample of SMGs with spectroscopic redshifts. Fig. 6 also shows that the redshift distributions for the quiescent galaxies peak at lower redshifts than the bursts. The quiescent galaxy peak is at very low redshifts for the brighter part of the flux range shown. Low-z quiescent galaxies are predicted to dominate the total counts at fluxes Sν(850 μm) ≳ 30 mJy.
We will discuss the predictions of the new models for the properties of the SMGs and LBGs at length in Lacey et al. (Paper 2, in preparation). In summary, at z= 2, SMGs with S850 μm≥ 5 mJy have median stellar masses of ∼1010h−1 M⊙ and reside in haloes with a median mass of ∼1012h−1 M⊙. The SMGs are thus predicted to reside in the more massive haloes in place at z= 2, and therefore will be more strongly clustered than the dark matter at this epoch (Baugh et al., in preparation), consistent with tentative observational constraints (Blain et al. 2004).
Discussion and Conclusions
We have combined two powerful theoretical techniques, semi-analytical modelling of galaxy formation and radiative transfer modelling of the reprocessing of stellar radiation by dust, to address the issue of how faint submillimetre galaxies fit into current theories of galaxy formation. Our approach to this problem represents a significant advance over previous work in this area for two reasons. First, our overriding philosophy is to require that any acceptable model must reproduce a basic set of observational data for present-day galaxies, including the optical, near-IR and far-IR luminosity functions, gas fractions and metallicities, and the distribution of disc scale sizes. Secondly, the dust extinction and emission are computed self-consistently, rather than being put in by hand. Our adopted methodology severely restricts the range of available parameter space, making a successful reproduction of the high-redshift galaxy population more challenging, but, at the same time resulting in a model with more predictive power.
The new model presented in this paper extends the successes of the G00 model to high redshift. Our model matches the measured rest-frame UV luminosity function of Lyman-break galaxies at z≈ 3 (Steidel et al. 1999). At the same time, it reproduces the number counts of 850-μm sources, and predicts redshifts for these sources in very good agreement with recent measurements from Chapman et al. (2003, 2004). This is the first time that any realistic CDM-based galaxy formation model has successfully matched the properties of both major populations of star-forming galaxies at high redshift.
In order to achieve these successes while still remaining within the framework of the hierarchical assembly of galaxies, we found it necessary to enhance the importance of bursts at high redshift, by: (i) modifying the quiescent star formation time-scale in galactic discs, in order to make mergers at high redshift more gas-rich; (ii) allowing triggering of bursts by minor and major mergers; and (iii) adopting a top-heavy IMF in bursts.
Our new model shares some features with those of Guiderdoni et al. (1998) and Devriendt & Guiderdoni (2000). They also proposed that the fraction of star formation occurring in bursts increased strongly with redshift, and suggested that the IMF in bursts was top-heavy, in order to explain the cosmic far-IR background and the submillimetre number counts. However, in other respects, their model is very different from ours. In particular, their model does not include galaxy mergers, so the rate of starbursts as a function of redshift is put in by hand. Also, they use empirical SEDs for the IR/submillimetre emission, rather than calculating dust temperatures from radiative equilibrium.
The most controversial of our changes is the adoption of a top-heavy IMF in bursts of star formation. For simplicity, we have chosen to use a flat IMF in starbursts. We have explored a range of tilted and truncated IMFs. A flat IMF slope is not essential to match the number counts of SMGs; we found that reasonable matches to the counts could be obtained with slopes given by x < 0.35. The critical point is that the relative proportion of high- to low-mass stars produced in a starburst should be greater than in quiescent star formation.
We are by no means the first to suggest the abandonment of a universal IMF. Some theoretical calculations predict a form for the IMF that is biased towards more massive stars, compared with the Kennicutt (1983) IMF, under the conditions that exist in starbursts (e.g. Padoan, Nordlund & Jones 1997; Larson 1998). However, the bulk of the support for a variation in the IMF is indirect, resulting from the challenge of matching observations of early-type galaxies and explaining the metal content of the hot gas in clusters (e.g. Renzini et al. 1993; Gibson & Matteucci 1997; Chiosi et al. 1998; Thomas, Greggio & Bender 1999). The metal content of the intracluster medium and stellar absorption-line strengths in ellipticals can be more readily explained if significant star formation occurred with an IMF which had a higher fraction of massive stars compared with that in the solar neighbourhood. This boost in the relative proportion of high-mass stars can be achieved either by tilting the IMF or by truncating it below a fixed stellar mass. Nagashima et al. (2004) have taken the new model described in this paper and incorporated a detailed treatment of chemical enrichment by Type Ia and Type II supernovae. The new model reproduces the observed abundances of O, Mg, Si and Fe in the intracluster medium of X-ray clusters. Nagashima et al. find that reverting to a standard IMF in bursts leads to the model underpredicting the abundances of these elements by a factor of 2–3.
Granato et al. (2004) have recently presented an alternative model which can explain the submillimetre number counts without assuming either a top-heavy IMF or a modified dust emissivity, based on spheroid formation by monolithic collapse at high redshift. They posit that the formation of very massive stellar spheroids occurs much earlier than in the model presented here, through the joint effects of very strong feedback from QSOs suppressing the cooling of gas in small haloes at high redshift, and a rate of gas cooling in high-mass haloes that is greatly enhanced by clumping effects. This model also reproduces the inferred evolution of the K-band luminosity function (e.g. Pozzetti et al. 2003; Drory et al. 2004). However, unlike the model presented here, the Granato et al. model does not include the effects of either halo or galaxy mergers, and does not currently treat the formation of galactic discs. The differences between these two models highlight the need to develop a better physical understanding of the processes of gas cooling and feedback, in order to understand whether the assembly history of galactic spheroids was closer to the ‘monolithic’ or ‘hierarchical’ pictures.
The relative contributions to the 850-μm fluxes of SMGs from dust heating by AGNs and by starbursts have been the subject of much recent debate. The first Chandra X-ray images of areas which had been mapped in the submillimetre range with SCUBA revealed little overlap between the X-ray and submillimetre sources, suggesting that SMGs were not powered by obscured QSOs (e.g. Almaini et al. 2003). Deeper X-ray images have revealed a much closer correlation between X-ray and submillimetre sources, with 30–50 per cent of bright SMGs (S850 μm>5 mJy) having X-ray-detected AGN counterparts, and some others showing X-ray emission consistent with star formation (Alexander et al. 2003). However, despite the apparent commonness of AGNs in SMGs, Alexander et al. argue, based on an analysis of SED shapes, that even in the submillimetre sources which host AGNs, the AGN contributes on average <2 per cent of the 850-μm flux. The SMGs thus appear to be powered by star formation, as assumed in our model.
In our new model, star formation in bursts is very important at high redshift; however, bursts are still responsible for only 30 per cent of the total star formation when integrated over all redshifts, compared with 5 per cent in G00. Once the recycling of mass from dying stars is taken into account, less than 7 per cent of the mass locked up in stars today is predicted to have been produced in bursts. The new model and the G00 model produce similar stellar mass fractions in discs and spheroids; in the new model 57 per cent of the stellar mass is in discs and 43 per cent is in spheroids. This is in very good agreement with the observational estimate by Benson, Frenk & Sharples (2002a). The most massive galaxies in the model tend to be bulge-dominated, in agreement with observations (e.g. Kochanek et al. 2001; Nakamura et al. 2003). However, typically only a few per cent of the present-day stellar mass in these galaxies was formed in bursts at z∼ 2, around the peak in the predicted redshift distribution of 850-μm sources. In our model, the overwhelming bulk of the stellar mass of ellipticals is built up in discs and then rearranged into spheroids during mergers (Baugh et al. 1996; Kauffmann 1996).
Our model is able to address many other issues relating to submillimetre galaxies and to galaxy evolution generally, which will we pursue in subsequent papers in this series. These include: a more detailed comparison with galaxy properties in the local Universe (Lacey et al., Paper 1, in preparation); the environments and clustering of SMGs and LBGs (Baugh et al., in preparation); a more detailed analysis of the properties of LBGs and SMGs and of the relationship between them (Lacey et al., Paper 2, in preparation); the descendants of SMGs (Frenk et al., in preparation); and a comprehensive set of predictions for number counts at different wavelengths in the IR and submillimetre ranges, and for the extragalactic background light (Lacey et al., Paper 3, in preparation).
CMB and AJB acknowledge receipt of Royal Society University Research Fellowships. This work was supported in part by a PPARC rolling grant at Durham. We thank Ian Smail for his helpful comments and encouragement, and the referee for providing a detailed and helpful report.