Mon. Not. R. Astron. Soc. 000, 000–000 (0000) Printed 8 November 2021 (MN LATEX style file v2.2)
Galaxy evolution in the infra-red: comparison of a hierarchical
galaxy formation model with SPITZER data
C. G. Lacey ⋆,1 C. M. Baugh,1 C.S. Frenk,1 L. Silva,2 G.L. Granato,3 and A. Bressan,3
1Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham, DH1 3LE, UK
2INAF, Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy
3INAF, Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 2, I-35122 Padova, Italy.
8 November 2021
ABSTRACT
We present predictions for the evolution of the galaxy luminosity function, number counts and
redshift distributions in the IR based on the ΛCDM cosmological model. We use the combined
GALFORM semi-analytical galaxy formation model and GRASIL spectrophotometric code to
compute galaxy SEDs including the reprocessing of radiation by dust. The model, which is
the same as that in Baugh et al. (2005), assumes two different IMFs: a normal solar neigh-
bourhood IMF for quiescent star formation in disks, and a very top-heavy IMF in starbursts
triggered by galaxy mergers. We have shown previously that the top-heavy IMF seems to be
necessary to explain the number counts of faint sub-mm galaxies. We compare the model
with observational data from the Spitzer Space Telescope, with the model parameters fixed
at values chosen before Spitzer data became available. We find that the model matches the
observed evolution in the IR remarkably well over the whole range of wavelengths probed by
Spitzer. In particular, the Spitzer data show that there is strong evolution in the mid-IR galaxy
luminosity function over the redshift range z ∼ 0 − 2, and this is reproduced by our model
without requiring any adjustment of parameters. On the other hand, a model with a normal
IMF in starbursts predicts far too little evolution in the mid-IR luminosity function, and is
therefore excluded.
Key words: galaxies: evolution – galaxies: formation – galaxies: high-redshift – infrared:
galaxies – ISM: dust, extinction
1 INTRODUCTION
In recent years, the evolution of galaxies at mid- and far-infrared
wavelengths has been opened up for direct observational study by
infrared telescopes in space. Already in the 1980s, the IRAS satel-
lite surveyed the local universe in the IR, showing that much of
present-day star formation is optically obscured, revealing a pop-
ulation of luminous and ultra-luminous infrared galaxies (LIRGs
with total IR luminosities LIR ∼ 10
− 1012L⊙ and ULIRGs
with LIR & 10
12L⊙), and providing the first hints of strong evo-
lution in the number density of ULIRGs at recent cosmic epochs
(e.g. Wright et al. 1984; Soifer et al. 1987a; Sanders & Mirabel
1996). The next major advance came with the discovery by COBE
of the cosmic far-IR background which has an energy density
comparable to that in the optical/near-IR background (Puget et al.
1996; Hauser et al. 1998). This implies that, over the history of the
universe, as much energy has been emitted by dust in galaxies as
reaches us directly in starlight, after dust extinction is taken into
account. This discovery made apparent the need to understand the
IR as much as the optical emission from galaxies in order to have
⋆ E-mail: Cedric.Lacey@durham.ac.uk (CGL)
a complete picture of galaxy evolution. In particular, it is essen-
tial to understand IR emission from dust in order to understand the
cosmic history of star formation, since most of the radiation from
young stars must have been absorbed by dust over the history of
the universe, in order to account for the far-IR background (e.g.
Hauser et al. 1998).
Following these early discoveries, the ISO satellite enabled the
first deep surveys of galaxies in the mid- and far-IR. The deep-
est of these surveys were in the mid-IR at 15µm, and probed the
evolution of LIRGs and ULIRGs out to z ∼ 1, showing strong
evolution in these populations, and directly resolving most of the
cosmic infrared background at that wavelength (Elbaz et al. 1999,
2002; Gruppioni et al. 2002). Deep ISO surveys in the far-IR at
170µm (Dole et al. 2001; Patris et al. 2003) probed lower red-
shifts, z ∼ 0.5. Around the same time, sub-mm observations using
the SCUBA instrument on the JCMT revealed a huge population of
high-z ULIRGs (Smail, Ivison & Blain 1997; Hughes et al. 1998)
which were subsequently found to have a redshift distribution peak-
ing at z ∼ 2 (Chapman et al. 2005), confirming the dramatic evo-
lution in number density for this population seen at shorter wave-
lengths and lower redshifts. The sub-mm galaxies have been stud-
c© 0000 RAS
http://arxiv.org/abs/0704.1562v2
2 Lacey et al.
ied in more detail in subsequent SCUBA surveys (e.g. SHADES,
Mortier et al. 2005).
Now observations using the Spitzer satellite (Werner et al.
2004), with its hugely increased sensitivity and mapping speed
are revolutionizing our knowledge of galaxy evolution at IR wave-
lengths from 3.6 to 160 µm. Spitzer surveys have allowed direct
determinations of the evolution of the galaxy luminosity func-
tion out to z ∼ 1 in the rest-frame near-IR and to z ∼ 2 in
the mid-IR (Le Floc’h et al. 2005; Perez-Gonzalez et al. 2005;
Babbedge et al. 2006; Franceschini et al. 2006). Individual galax-
ies have been detected by Spitzer out to z ∼ 6 (Eyles et al. 2005).
In the near future, the Herschel satellite (Pilbratt 2003) should
make it possible to measure the far-IR luminosity function out to
z ∼ 2, and thus directly measure the total IR luminosities of galax-
ies over most of the history of the universe.
Accompanying these observational advances, various types of
theoretical models have been developed to interpret or explain the
observational data on galaxy evolution in the IR. We can distinguish
three main classes of model:
(a) Purely phenomenological models: In these models,
the galaxy luminosity function and its evolution are described
by a purely empirical expression, and this is combined with
observationally-based templates for the IR spectral energy dis-
tribution (SED). The free parameters in the expression for the
luminosity function are then chosen to obtain the best match
to some set of observational data, such as number counts
and redshift distributions in different IR bands. These pa-
rameters are purely descriptive and provide little insight into
the physical processes which control galaxy evolution. Exam-
ples of these models are Pearson & Rowan-Robinson (1996);
Xu et al. (1998); Blain et al. (1999); Franceschini et al. (2001);
Chary & Elbaz (2001); Rowan-Robinson (2001); Lagache et al.
(2003); Gruppioni et al. (2005).
(b) Hierarchical galaxy formation models with phenomeno-
logical SEDs: In these models, the evolution of the luminosity func-
tions of the stellar and total dust emission are calculated from a
detailed model of galaxy formation based on the cold dark mat-
ter (CDM) model of structure formation, including physical mod-
elling of processes such as gas cooling and galaxy mergers. The
stellar luminosity of a model galaxy is computed from its star
formation history, and the stellar luminosity absorbed by dust,
which equals the total IR luminosity emitted by dust, is calculated
from this based on some treatment of dust extinction. However,
the SED shapes required to calculate the distribution of the dust
emission over wavelength from the total IR dust emission are ei-
ther observationally-based templates (e.g. Guiderdoni et al. 1998;
Devriendt & Guiderdoni 2000) or are purely phenomenological,
e.g. a modified Planck function with an empirically chosen dust
temperature (e.g. Kaviani et al. 2003). In this approach, the shape
of the IR SED assumed for a model galaxy may be incompatible
with its other predicted properties, such as its dust mass and radius.
(c) Hierarchical galaxy formation models with theoretical
SEDs: These models are similar to those of type (b), in that the
evolution of the galaxy population is calculated from a detailed
physical model of galaxy formation based on CDM, but instead
of using phenomenological SEDs for the dust emission, the com-
plete SED of each model galaxy, from the far-UV to the radio, is
calculated by combining a theoretical stellar population synthesis
model for the stellar emission with a theoretical radiative transfer
and dust heating model to predict both the extinction of starlight
by dust and the IR/sub-mm SED of the dust emission. The ad-
vantages of this type of model are that it is completely ab initio,
with the maximum possible theoretical self-consistency, and all of
the model parameters relate directly to physical processes. For ex-
ample, the typical dust temperature and the shape of the SED of
dust emission depend on the stellar luminosity and the dust mass,
and evolution in all of these quantities is computed self-consistently
in this type of model. Following this modelling approach thus al-
lows more rigorous testing of the predictions of physical models
for galaxy formation against observational data at IR wavelengths,
as well as shrinking the parameter space of the predictions. Ex-
amples of such models are Granato et al. (2000) and Baugh et al.
(2005). (An alternative modelling approach also based on theoret-
ical IR SEDs but with a simplified treatment of the assembly of
galaxies and halos has been presented by Granato et al. (2004) and
Silva et al. (2005).)
In this paper, we follow the third approach, with physical mod-
elling both of galaxy formation and of the galaxy SEDs, includ-
ing the effects of dust. This paper is the third in a series, where
we combine the GALFORM semi-analytical model of galaxy forma-
tion (Cole et al. 2000) with the GRASILmodel for stellar and dust
emission from galaxies (Silva et al. 1998). The GALFORM model
computes the evolution of galaxies in the framework of the ΛCDM
model for structure formation, based on physical treatments of the
assembly of dark matter halos by merging, gas cooling in halos,
star formation and supernova feedback, galaxy mergers, and chem-
ical enrichment. The GRASILmodel computes the SED of a model
galaxy from the far-UV to the radio, based on theoretical models of
stellar evolution and stellar atmospheres, radiative transfer through
a two-phase dust medium to calculate both the dust extinction and
dust heating, and a distribution of dust temperatures in each galaxy
calculated from a detailed dust grain model. In the first paper in
the series (Granato et al. 2000), we modelled the IR properties of
galaxies in the local universe. While this model was very success-
ful in explaining observations of the local universe, we found sub-
sequently that it failed when confronted with observations of star-
forming galaxies at high redshifts, predicting far too few sub-mm
galaxies (SMGs) at z ∼ 2 and Lyman-break galaxies (LBGs) at
z ∼ 3. Therefore, in the second paper (Baugh et al. 2005), we
proposed a new version of the model which assumes a top-heavy
IMF in starbursts (with slope x = 0, compared to Salpeter slope
x = 1.35), but a normal solar neighbourhood IMF for quiescent
star formation. In this new model, the star formation parameters
were also changed to force more star formation to happen in bursts.
This revised model agreed well with both the number counts and
redshift distributions of SMGs detected at 850µm, and with the
rest-frame far-UV luminosity function of LBGs at z ∼ 3, while
still maintaining consistency with galaxy properties in the local uni-
verse such as the optical, near-IR and far-IR luminosity functions,
and gas fractions, metallicities, morphologies and sizes.
This same model of Baugh et al. (2005) was found by
Le Delliou et al. (2005a, 2006) to provide a good match to the ob-
served evolution of the population of Lyα-emitting galaxies over
the redshift range z ∼ 3−6. Support for the controversial assump-
tion of a top-heavy IMF in bursts came from the studies of chem-
ical enrichment in this model by Nagashima et al. (2005a,b), who
found that the metallicities of both the intergalactic gas in galaxy
clusters and the stars in elliptical galaxies were predicted to be sig-
nificantly lower than observed values if a normal IMF was assumed
for all star formation, but agreed much better if a top-heavy IMF in
bursts was assumed, as in Baugh et al. . In this third paper in the
series, we extend the Baugh et al. (2005) model to make predic-
tions for galaxy evolution in the IR, and compare these predictions
with observational data from Spitzer. We emphasize that all of the
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 3
model parameters for the predictions presented in this paper were
fixed by Baugh et al. prior to the publication of any results from
Spitzer, and we have not tried to obtain a better fit to any of the
Spitzer data by adjusting these parameters1.
Our goals in this paper are to test our model of galaxy evo-
lution with a top-heavy IMF in starbursts against observations of
dust-obscured star-forming galaxies over the redshift range z ∼
0−2, and also to test our predictions for the evolution of the stellar
populations of galaxies against observational data in the rest-frame
near- and mid-IR. The plan of the paper is as follows: In Section 2,
we give an overview of the GALFORM and GRASIL models, fo-
cusing on how the predictions we present later on are calculated.
In Section 3, we compare the galaxy number counts predicted by
our model with observational data in all 7 Spitzer bands, from 3.6
to 160 µm. In Section 4, we investigate galaxy evolution in the
IR in more detail, by comparing model predictions directly with
galaxy luminosity functions constructed from Spitzer data. In Sec-
tion 5, we present the predictions of our model for the evolution
of the galaxy stellar mass function and star formation rate distribu-
tion, and investigate the insight our model offers on how well stel-
lar masses and star formation rates can be estimated from Spitzer
data. We present our conclusions in Section 6. In the Appendix, we
present model predictions for galaxy redshift distributions in the
different Spitzer bands, to assist in interpreting data from different
surveys.
2 MODEL
In this paper use the GALFORM semi-analytical model to predict
the physical properties of the galaxy population at different red-
shifts, and combine it with the GRASIL spectrophotometric model
to predict the detailed SEDs of model galaxies. Both GALFORM
and GRASIL have been described in detail in various previous
papers, but since the descriptions of the different model compo-
nents, as well as of our particular choice of parameters, are spread
among different papers, we give an overview of both of these here.
GALFORM is described in §2.1, and GRASIL in §2.2. Particularly
important features of our model are the triggering of starbursts by
mergers (discussed in §2.1.4) and the assumption of a top-heavy
IMF in starbursts (discussed in §2.1.7). We further discuss the
choice of model parameters in §2.3. Readers who are already fa-
miliar with the Baugh et al. (2005) model can skip straight to the
results, starting in §3.
2.1 GALFORM galaxy formation model
We compute the formation and evolution of galaxies within the
framework of the ΛCDM model of structure formation using the
semi-analytical galaxy formation model GALFORM. The general
methodology and approximations behind the GALFORM model are
set out in detail in Cole et al. (2000) (see also the review by Baugh
(2006)). In summary, the GALFORM model follows the main pro-
cesses which shape the formation and evolution of galaxies. These
include: (i) the collapse and merging of dark matter halos; (ii) the
1 A closely related model of galaxy formation obtained by applying
GALFORM principles to the Millennium simulation of Springel et al. (2005)
has recently been published by Bower et al. (2006). This model differs
from the current one primarily in that it includes feedback from AGN activ-
ity, but does not have a top-heavy IMF in bursts. We plan to investigate the
IR predictions of this alternative model in a subsequent paper.
shock-heating and radiative cooling of gas inside dark halos, lead-
ing to the formation of galaxy disks; (iii) quiescent star formation
in galaxy disks; (iv) feedback both from supernova explosions and
from photo-ionization of the IGM; (v) chemical enrichment of the
stars and gas; (vi) galaxy mergers driven by dynamical friction
within common dark matter halos, leading to the formation of stel-
lar spheroids, and also triggering bursts of star formation. The end
product of the calculations is a prediction of the numbers and prop-
erties of galaxies that reside within dark matter haloes of different
masses. The model predicts the stellar and cold gas masses of the
galaxies, along with their star formation and merger histories, their
sizes and metallicities.
The prescriptions and parameters for the different processes
which we use in this paper are identical to those adopted by
Baugh et al. (2005), but differ in several important respects from
Cole et al. (2000). All of these parameters were chosen by com-
parison with pre-Spitzer observational data. The background cos-
mology is a spatially flat CDM universe with a cosmological con-
stant, with “concordance” parameters Ωm = 0.3, ΩΛ = 0.7,
Ωb = 0.04, and h ≡ H0/(100km s
−1Mpc−1) = 0.7. The am-
plitude of the initial spectrum of density fluctuations is set by the
r.m.s. linear fluctuation in a sphere of radius 8h−1Mpc, σ8 = 0.93.
For completeness, we now summarize the prescriptions and param-
eters used, but give details mainly where they differ from those in
Cole et al. (2000), or where they are particularly relevant to pre-
dicting IR emission from dust.
2.1.1 Halo assembly histories
As in Cole et al. (2000), we describe the assembly histories of dark
matter halos through halo merger trees which are calculated using
a Monte Carlo method based on the extended Press-Schechter ap-
proach (e.g. Lacey & Cole 1993). The process of galaxy forma-
tion is then calculated separately for each halo merger tree, follow-
ing the baryonic physics in all of the separate branches of the tree.
As has been shown by Helly et al. (2003), the statistical properties
of galaxies calculated in semi-analytical models using these Monte
Carlo merger trees are very similar to those computed using merger
trees extracted directly from N-body simulations.
2.1.2 Gas cooling in halos
The cooling of gas in halos is calculated using the same simple
spherical model as in Cole et al. (2000). The diffuse gas in halos
(consisting of all of the gas which has not previously condensed
into galaxies) is assumed to be shock-heated to the halo virial tem-
perature when the halo is assembled, and then to cool radiatively by
atomic processes. The cooling time depends on radius through the
gas density profile, which is assumed to have a core radius which
grows as gas is removed from the diffuse phase by condensing into
galaxies. The gas at some radius r in the halo then cools and col-
lapses to the halo centre on a timescale which is the larger of the
cooling time tcool and the free-fall time tff at that radius. Thus, for
tcool(r) > tff(r), we have hot accretion, and for tcool(r) < tff(r),
we have cold accretion 2. In our model, gas only accretes onto the
central galaxy in a halo, not onto any satellite galaxies which share
2 Note that contrary to claims by Birnboim & Dekel (2003), the process
of “cold accretion”, if not the name, has always been part of semi-analytical
models (see Croton et al. (2006) for a detailed discussion)
c© 0000 RAS, MNRAS 000, 000–000
4 Lacey et al.
that halo. We denote all of the diffuse gas in halos as “hot”, and all
of the gas which has condensed into galaxies as “cold”.
2.1.3 Star formation timescale in disks
The global rate of star formation ψ in galaxy disks is assumed to be
related to the cold gas mass, Mgas, by ψ = Mgas/τ∗,disk, where
the star formation timescale is taken to be
τ∗,disk = τ∗0
Vc/200 km s
, (1)
where Vc is the circular velocity of the galaxy disk (at its half-
mass radius) and τ∗0 is a constant. We adopt values τ∗0 = 8Gyr
and α∗ = −3, chosen to reproduce the observed relation between
gas mass and B-band luminosity for present-day disk galaxies. As
discussed in Baugh et al. (2005), this assumption means that the
disk star formation timescale is independent of redshift (at a given
Vc), resulting in disks at high redshift that are much more gas-rich
than at low redshift, and have more gas available for star formation
in bursts triggered by galaxy mergers at high redshift.
2.1.4 Galaxy mergers and triggering of starbursts
In the model, all galaxies originate as central galaxies in some halo,
but can then become satellite galaxies if their host halo merges into
another halo. Mergers can then occur between satellite and central
galaxies within the same halo, after dynamical friction has caused
the satellite galaxy to sink to the centre of the halo. Galaxy mergers
can produce changes in galaxy morphology and trigger bursts. We
classify galaxy mergers according to the ratio of masses (including
stars and gas) M2/M1 6 1 of the secondary to primary galaxy
involved. We define mergers to be major or minor according to
whether M2/M1 > fellip or M2/M1 < fellip (Kauffmann et al.
1993). In major mergers, any stellar disks in either the primary or
secondary are assumed to be disrupted, and the stars rearranged into
a spheroid. In minor mergers, the stellar disk in the primary galaxy
is assumed to remain intact, while all of the stars in the secondary
are assumed to be added to the spheroid of the primary. We adopt a
threshold fellip = 0.3 for major mergers, consistent with the results
of numerical simulations (e.g. Barnes 1998), which reproduces the
observed present-day fraction of spheroidal galaxies. We assume
that major mergers always trigger a starburst if any gas is present.
We also assume that minor mergers can trigger bursts, if they sat-
isfy both M2/M1 > fburst and the gas fraction in the disk of the
primary galaxy exceeds fgas,crit. Following Baugh et al. (2005),
we adopt fburst = 0.05 and fgas,crit = 0.75. The parameters for
bursts in minor mergers were motivated by trying to explain the
number of sub-mm galaxies. An important consequence of assum-
ing eqn.(1) for the star formation timescale in disks, combined with
the triggering of starbursts in minor mergers, is that the global star
formation rate at high redshifts is dominated by bursts, while that at
low redshifts it is dominated by quiescent disks (see Baugh et al.
for a detailed discussion of these points).
In either kind of starburst, we assume that the burst consumes
all of the cold gas in the two galaxies involved in the merger, and
that the stars produced are added to the spheroid of the merger rem-
nant. During the burst, we assume that star formation proceeds ac-
cording to the relation ψ =Mgas/τ∗,burst. For the burst timescale,
we assume
τ∗,burst = max [fdynτdyn,sph; τ∗,burst,min] , (2)
where τdyn,sph is the dynamical time in the newly-formed spheroid.
We adopt fdyn = 50 and τ∗,burst,min = 0.2Gyr (these parame-
ters were chosen by Baugh et al. (2005) to allow a simultaneous
match to the sub-mm number counts and to the local 60µm lu-
minosity function). The star formation rate in a burst thus decays
exponentially with time after the galaxy merger. It is assumed to be
truncated after 3 e-folding times (where the e-folding time takes ac-
count of stellar recycling and feedback - see Granato et al. (2000)
for details), with the remaining gas being ejected into the galaxy
halo at that time.
2.1.5 Feedback from photo-ionization
After the intergalactic medium (IGM) has been reionized at redshift
zreion, the formation of low-mass galaxies is inhibited, both by the
effect of the IGM pressure inhibiting collapse of gas into halos, and
by the reduction of gas cooling in halos due to the photo-ionizing
background. We model this in a simple way, by assuming that for
z < zreion, cooling of gas is completely suppressed in halos with
circular velocities Vc < Vcrit. We adopt Vcrit = 60 kms
−1, based
on the detailed modelling by Benson et al. (2002). We assume in
this paper that reionization occurs at zreion = 6, for consistency
with Baugh et al. (2005), but increasing this to zreion ∼ 10 in
line with the WMAP 3-year estimate of the polarization of the mi-
crowave background (Spergel et al. 2006) has no significant effect
on the model results presented in this paper.
2.1.6 Feedback from supernovae
Photo-ionization feedback only affects very low mass galaxies.
More important for most galaxies is feedback from supernova ex-
plosions. We assume that energy input from supernovae causes gas
to be ejected from galaxies at a rate
Ṁej = β(Vc)ψ = [βreh(Vc) + βsw(Vc)] ψ (3)
The supernova feedback is assumed to operate for both quiescent
star formation in disks and for starbursts triggered by galaxy merg-
ers, with the only difference being that we take Vc to be the circular
velocity at the half mass radius of the disk in the former case, and
at the half-mass radius of the spheroid in the latter case. For sim-
plicity, we keep the same feedback parameters for starbursts as for
quiescent star formation.
The supernova feedback has two components: the reheating
term βrehψ describes gas which is reheated and ejected into the
galaxy halo, from where it is allowed to cool again after the halo
mass has doubled through hierachical mass build-up. For this, we
use the parametrization of Cole et al. (2000):
βreh = (Vc/Vhot)
−αhot , (4)
where we adopt parameter values Vhot = 300 kms
−1 and αhot =
2. The reheating term has the largest effect on low-mass galaxies,
for which ejection of gas from galaxies flattens the faint-end slope
of the galaxy luminosity function.
The second term βswψ in eqn.(3) is the superwind term, which
describes ejection of gas out of the halo rather than just the galaxy.
Once ejected, this gas is assumed never to re-accrete onto any halo.
We model the superwind ejection efficiency as
βsw = fsw min
1, (Vc/Vsw)
based on Benson et al. (2003). We adopt parameter values fsw = 2
and Vsw = 200 km s
−1, as in Baugh et al. (2005). The superwind
term mainly affects higher mass galaxies, where the ejection of gas
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 5
from halos causes an increase in the cooling time of gas in halos
by reducing the gas densities. This brings the predicted break at
the bright end of the local galaxy luminosity function into agree-
ment with observations, as discussed in Benson et al. (2003). The
various parameters for supernova feedback are thus chosen in or-
der to match the observed present-day optical and near-IR galaxy
luminosity functions, as well as the galaxy metallicity-luminosity
relation.
We note that the galaxy formation model in this paper, un-
like some other recent semi-analytical models, does not include
AGN feedback. Instead, the role of AGN feedback in reducing the
amount of gas cooling to form massive galaxies is taken by su-
perwinds driven by supernova explosions. The first semi-analytical
model to include AGN feedback was that of Granato et al. (2004),
who introduced a detailed model of feedback from QSO winds
during the formation phase of supermassive black holes (SMBHs),
with the aim of explaining the co-evolution of the spheroidal com-
ponents of galaxies and their SMBHs. The predictions of the
Granato et al. model for number counts and redshift distributions
in the IR have been computed by Silva et al. (2005) using the
GRASIL spectrophotometric model, and compared to ISO and
Spitzer data. However, the Granato et al. (2004) model has the
limitations that it does not include the merging of galaxies or of
dark halos, and does not treat the formation and evolution of galac-
tic disks. More recently, several semi-analytical models have been
published which propose that heating of halo gas by relativistic
jets from an AGN in an optically inconspicuous or “radio” mode
can balance radiative cooling of gas in high-mass halos, thus sup-
pressing hot accretion of gas onto galaxies (Bower et al. 2006;
Croton et al. 2006; Cattaneo et al. 2006; Monaco et al. 2007).
However, these AGN feedback models differ in detail, and all are
fairly schematic. None of these models has been shown to repro-
duce the observed number counts and redshifts of the faint sub-mm
galaxies.
The effects of our superwind feedback are qualitatively quite
similar to those of the radio-mode AGN feedback. Both superwind
and AGN feedback models contain free parameters, which are ad-
justed in order to make the model fit the bright end of the ob-
served present-day galaxy luminosity function at optical and near-
IR wavelengths. However, since the physical mechanisms are dif-
ferent, they make different predictions for how the galaxy lumi-
nosity function should evolve with redshift. Current models for the
radio-mode AGN feedback are very schematic, but they have the
advantage over the superwind model that the energetic constraints
are greatly relaxed, since accretion onto black holes can convert
mass into energy with a much higher efficiency than can supernova
explosions. We will investigate the predictions of models with AGN
feedback for the IR and sub-mm evolution of galaxies in a future
paper.
2.1.7 The Stellar Initial Mass Function and Chemical Evolution
Stars in our model are assumed to form with different Initial Mass
Functions (IMFs), depending on whether they form in disks or in
bursts. Both IMFs are taken to be piecewise power laws, with slopes
x defined by dN/d lnm ∝ m−x, with N the number of stars and
m the stellar mass (so the Salpeter slope is x = 1.35), and covering
a stellar mass range 0.15 < m < 120M⊙. Quiescent star forma-
tion in galaxy disks is assumed to have a solar neighbourhood IMF,
for which we use the Kennicutt (1983) paramerization, with slope
x = 0.4 for m < M⊙ and x = 1.5 for m > M⊙. (The Kennicutt
(1983) IMF is similar to other popular parametrizations of the solar
neighbourhood IMF, such as that of Kroupa (2001).) Bursts of star
formation triggered by galaxy mergers are assumed to form stars
with a top-heavy IMF with slope x = 0. As discussed in detail
in Baugh et al. (2005), the top-heavy IMF in bursts was found to
be required in order to reproduce the observed number counts and
redshift distributions of the faint sub-mm galaxies. Furthermore,
as shown by Nagashima et al. (2005a,b), the predicted chemical
abundances of the X-ray emitting gas in galaxy clusters and of the
stars in elliptical galaxies also agree better with observational data
in a model with the top-heavy IMF in bursts, rather than a universal
solar neighbourhood IMF.
A variety of other observational evidence has accumulated
which suggests that the IMF in some environments may be top-
heavy compared to the solar neighbourhood IMF. Rieke et al.
(1993) argued for a top-heavy IMF in the nearby starburst M82,
based on modelling its integrated properties, while Parra et al.
(2007) found possible evidence for a top-heavy IMF in the ultra-
luminous starburst Arp220 from the relative numbers of super-
novae of different types observed at radio wavelengths. Evidence
has been found for a top-heavy IMF in some star clusters in in-
tensely star-forming regions, both in M82 (e.g. McCrady et al.
2003), and in our own Galaxy (e.g. Figer et al. 1999; Stolte et al.
2005; Harayama et al. 2007). Observations of both the old and
young stellar populations in the central 1 pc of our Galaxy
also favour a top-heavy IMF (Paumard et al. 2006; Maness et al.
2007). Fardal et al. (2006) found that reconciling measurements of
the optical and IR extragalactic background with measurements of
the cosmic star formation history also seemed to require an average
IMF that was somewhat top-heavy. Finally, van Dokkum (2007)
found that reconciling the colour and luminosity evolution of early-
type galaxies in clusters also favoured a top-heavy IMF. Larson
(1998) summarized other evidence for a top-heavy IMF during the
earlier phases of galaxy evolution, and argued that this could be
a natural consequence of the temperature-dependence of the Jeans
mass for gravitational instability in gas clouds. Larson (2005) ex-
tended this to argue that a top-heavy IMF might also be expected
in starburst regions, where there is strong heating of the dust by the
young stars.
In our model, the fraction of star formation occuring in the
burst mode increases with redshift (see Baugh et al. (2005)), so the
average IMF with which stars are being formed shifts from being
close to a solar neighbourhood IMF at the present-day to being very
top-heavy at high redshift. In this model, 30% of star formation
occured in the burst mode when integrated over the past history of
the universe, but only 7% of the current stellar mass was formed
in bursts, because of the much larger fraction of mass recycled by
dying stars for the top-heavy IMF. We note that our predictions for
the IR and sub-mm luminosities of starbursts are not sensitive to
the precise form of the top-heavy IMF, but simply require a larger
fraction of m ∼ 5− 20M⊙ stars relative to a solar neighbourhood
In this paper, we calculate chemical evolution using the instan-
taneous recycling approximation, which depends on the total frac-
tion of mass recycled from dying stars (R), and the total yield of
heavy elements (p). Both of these parameters depend on the IMF.
We use the results of stellar evolution computations to calculate
values of R and p consistent with each IMF (see Nagashima et al.
(2005a) for details of the stellar evolution data used). Thus, we use
R = 0.41 and p = 0.023 for the quiescent IMF, and R = 0.91 and
p = 0.15 for the burst IMF. Our chemical evolution model then
predicts the masses and total metallicities of the gas and stars in
each galaxy as a function of time.
c© 0000 RAS, MNRAS 000, 000–000
6 Lacey et al.
2.1.8 Galaxy sizes and dust masses
For calculating the extinction and emission by dust, it is essential to
have an accurate calculation of the dust optical depths in the model
galaxies, which in turn depends on the mass of dust and the size
of the galaxy. The dust mass is calculated from the gas mass and
metallicity predicted by the chemical enrichment model, assuming
that the dust-to-gas ratio is proportional to metallicity, normalized
to match the local ISM value at solar metallicity. The sizes of galax-
ies are computed exactly as in Cole et al. (2000): gas which cools
in a halo is assumed to conserve its angular momentum as it col-
lapses, forming a rotationally-supported galaxy disk; the radius of
this disk is then calculated from its angular momentum, includ-
ing the gravity of the disk, spheroid (if any) and dark halo. Galaxy
spheroids are built up both from pre-existing stars in galaxy merg-
ers, and from the stars formed in bursts triggered by these mergers;
the radii of spheroids formed in mergers are computed using an
energy conservation argument. In calculating the sizes of disks and
spheroids, we include the adiabatic contraction of the dark halo due
to the gravity of the baryonic components. This model was tested
for disks by Cole et al. (2000) and for spheroids by Almeida et al.
(2007) (see also Coenda et al. in preparation, and Gonzalez et al.
in preparation). During a burst, we assume that the gas and stars
involved in the burst have a distribution with the same half-mass
radius as the spheroid (i.e. η = 1 in the notation of Granato et al.
(2000), who used a value η = 0.1).
2.2 GRASIL model for stellar and dust emission
For each galaxy in our model, we compute the spectral energy dis-
tribution using the spectrophotometric model GRASIL (Silva et al.
1998; Granato et al. 2000). GRASIL computes the emission from
the stellar population, the absorption and emission of radiation by
dust, and also radio emission (thermal and synchrotron) powered
by massive stars (Bressan et al. 2002).
2.2.1 SED model
The main features of the GRASIL model are as follows:
(i) The stars are assumed to have an axisymmetric distribution in a
disk and a bulge. Given the distribution of stars in age and metal-
licity (obtained from the star formation and chemical enrichment
history), the SED of the stellar population is calculated using a
population synthesis model based on the Padova stellar evolution
tracks and Kurucz model atmospheres (Bressan et al. 1998). This
is done separately for the disk and bulge.
(ii) The cold gas and dust in a galaxy are assumed to be in a 2-phase
medium, consisting of dense gas in giant molecular clouds embed-
ded in a lower-density diffuse component. In a quiescent galaxy,
the dust and gas are assumed to be confined to the disk, while for
a galaxy undergoing a burst, the dust and gas are confined to the
spheroidal burst component.
(iii) Stars are assumed to be born inside molecular clouds, and then
to leak out into the diffuse medium on a timescale tesc. As a result,
the youngest and most massive stars are concentrated in the dustiest
regions, so they experience larger dust extinctions than older, typ-
ically lower-mass stars, and dust in the clouds is also much more
strongly heated than dust in the diffuse medium.
(iv) The extinction of the starlight by dust is computed using a ra-
diative transfer code; this is used also to compute the intensity of
the stellar radiation field heating the dust at each point in a galaxy.
(v) The dust is modelled as a mixture of graphite and silicate grains
with a continuous distribution of grain sizes (varying between 8Å
and 0.25 µm), and also Polycyclic Aromatic Hydrocarbon (PAH)
molecules with a distribution of sizes. The equilibrium temperature
in the local interstellar radiation field is calculated for each type
and size of grain, at each point in the galaxy, and this information
is then used to calculate the emission from each grain. In the case of
very small grains and PAH molecules, temperature fluctuations are
important, and the probability distribution of the temperature is cal-
culated. The detailed spectrum of the PAH emission is obtained us-
ing the PAH cross-sections from Li & Draine (2001), as described
in Vega et al. (2005). The grain size distribution is chosen to match
the mean dust extinction curve and emissivity in the local ISM, and
is not varied, except that the PAH abundance in molecular clouds
is assumed to be 10−3 of that in the diffuse medium (Vega et al.
2005).
(vi) Radio emission from ionized gas in HII regions and from syn-
chrotron radiation from relativistic electrons accelerated in super-
nova remnant shocks are calculated as described in Bressan et al.
(2002).
The output from GRASIL is then the complete SED of a
galaxy from the far-UV to the radio (wavelengths 100Å . λ .
1m). The SED of the dust emission is computed as a sum over the
different types of grains, having different temperatures depending
on their size and their position in the galaxy. The dust SED is thus
intrinsically multi-temperature. GRASIL has been shown to give an
excellent match to the measured SEDs of both quiescent (e.g. M51)
and starburst (e.g. M82) galaxies (Silva et al. 1998; Bressan et al.
2002).
The assumption of axisymmetry in GRASIL is a limitation
when considering starbursts triggered by galaxy mergers. However,
observations of local ULIRGs imply that most of the star formation
happens in a single burst component after the galaxy merger is sub-
stantially complete, so the assumption of axisymmetry for the burst
component may not be so bad.
2.2.2 GRASIL parameters
The main parameters in the GRASIL dust model are the fraction
fmc of the cold gas which is in molecular clouds, the timescale tesc
for newly-formed stars to escape from their parent molecular cloud,
and the cloud masses Mc and radii rc in the combination Mc/r
which determines the dust optical depth of the clouds. We assume
fmc = 0.25, Mc = 10
6M⊙ and rc = 16pc as in Granato et al.
(2000), and also adopt the same geometrical parameters as in that
paper. We make the following two changes in GRASIL parame-
ters relative to Granato et al. , as discussed in Baugh et al. (2005):
(a) We assume tesc = 1Myr in both disks and bursts (instead
of the Granato et al. values tesc = 2 and 10Myr respectively).
This value was chosen in order to obtain a better match of the pre-
dicted rest-frame far-UV luminosity function of galaxies at z ∼ 3
to that measured for Lyman-break galaxies. (b) The dust emissivity
law in bursts at long wavelengths is modified from ǫν ∝ ν
−2 to
ǫν ∝ ν
−1.5 for λ > 100µm. This was done in order to improve
slightly the fit of the model to the observed sub-mm number counts.
In applying GRASIL to model the SEDs of a sample of nearby
galaxies, Silva et al. (1998) found that a similar modification (to
ǫν ∝ ν
−1.6) seemed to be required in the case of Arp220 (the only
ultra-luminous starburst in their sample), in order to reproduce the
observed sub-mm data for that galaxy. This modification in fact has
little effect on the IR predictions presented in the present paper, but
we retain it for consistency with Baugh et al. (2005).
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 7
2.2.3 Interface with GALFORM
For calculating the statistical properties of the galaxy population
from the combined GALFORM+GRASIL model, we follow the
same strategy as described in Granato et al. (2000). We first run
the GALFORM code to generate a large catalogue of model galaxies
at any redshift, and then run the GRASIL code on subsamples of
these. For the quiescent galaxies, we select a subsample which has
equal numbers of galaxies in equal logarithmic bins of stellar mass,
while for the bursting galaxies, we select a subsample with equal
numbers of galaxies in equal logarithmic bins of burst mass. For the
burst sample, we compute SEDs at several different representative
stages in the burst evolution, while for the quiescent sample, we
only compute SEDs at a single epoch. Using this sampling strat-
egy, we obtain a good coverage of all the different masses, types
and evolutionary stages of galaxies, while minimizing the compu-
tational cost of running the GRASIL code. The statistical properties
of the galaxy population are then obtained by assigning the model
galaxies appropriate weights depending on their predicted number
density in a representative cosmological volume.
The outputs from the GALFORM galaxy formation model re-
quired by GRASIL to calculate the galaxy SEDs are: the combined
star formation history and metallicity distribution for the disk and
bulge, the radii of both components, and the total mass of dust. The
dust mass is calculated from the mass and metallicity of the cold
gas in the galaxy, assuming that the dust-to-gas ratio is proportional
to the metallicity. Since the gas mass and metallicity both evolve,
so does the dust mass, and this evolution is fully taken into account
in GRASIL. For simplicity, we assume that the size distribution of
the dust grains and PAH molecules does not evolve, apart from the
normalization.
Once we have calculated the SEDs for the model galaxies, we
compute luminosities in different observed bands (e.g. the optical
B-band or the Spitzer 24µm band) by convolving the SED with the
filter+detector response function for that band. For computing the
predicted fluxes from galaxies in a fixed observer-frame band, we
redshift the SED before doing the convolution.
The GRASIL code is quite CPU-intensive, requiring several
minutes of CPU time per galaxy. Consequently, we are limited to
running samples of a few thousand galaxies at each redshift. As a
result, quantities such as luminosity functions and redshift distri-
butions still show some small amount of noise, rather than being
completely smooth curves, as can be seen in many of the figures in
this paper.
2.3 Choice of parameters in the GALFORM+GRASIL model
The combined GALFORM+GRASIL model has a significant num-
ber of parameters, but this is inevitable given the very wide range
of physical processes which are included. The parameters are con-
strained by requiring the model predictions to reproduce a limited
set of observational data - once this is done, there is rather little
freedom in the choice of parameters. We have described above how
the main parameters are fixed, and more details can be found in
Cole et al. (2000) and Baugh et al. (2005). For both of these pa-
pers, large grids of GALFORM models were run with different pa-
rameters, in order to decide which set of parameters gave the best
overall fit to the set of calibrating observational data. These papers
also show the effects of varying some of the main model parameters
around their best-fit values. The parameters in the standard model
for which we present results in this paper were chosen to reproduce
the following properties for present-day galaxies: the luminosity
functions in the B- and K-bands and at 60µm, the relations between
gas mass and luminosity and metallicity and luminosity, the size-
luminosity relation for galaxy disks, and the fraction of spheroidal
galaxies. In addition, the model was required to reproduce the ob-
served rest-frame far-UV (1500Å) luminosity function at z = 3,
and the observed sub-mm number counts and redshift distribution
at 850µm (Baugh et al. 2005). The sub-mm number counts are the
main factor driving the need to include a top-heavy IMF in bursts.
The parameters for our standard model are exactly the same as
in Baugh et al. (2005), which were chosen before Spitzer data be-
came available. Since these parameters were not adjusted to match
any data obtained with Spitzer, the predictions of our model in the
Spitzer bands are genuine predictions. We could obviously have
fine-tuned our parameters in order to match better the observational
data we considered in this paper, but this would have conflicted
with our main goal, which is to present predictions for a wide set of
observable properties based on a single physical model in a series
of papers.
Since our assumption of a top-heavy IMF in bursts is a con-
troversial one, we will also show some predictions from a variant
model, which is identical to the standard model, except that we
assume the same solar neighbourhood (Kennicutt) IMF in bursts
and in disks. Comparing the predictions for the standard and vari-
ant models then shows directly the effects of changing the IMF
in bursts. We note that the variant model matches the present-day
optical and near-IR luminosity functions almost as well as the stan-
dard model, though it is a poorer fit to the local 60µm luminosity
function for the brightest galaxies (see Fig. 9). The variant model
underpredicts the 850µm counts by a factor of 10–30.
3 NUMBER COUNTS
We begin our comparison of the predictions of our galaxy forma-
tion model against Spitzer data with the galaxy number counts.
Fig. 1 shows number counts in the four IRAC bands (3.6, 4.5, 5.8
and 8.0 µm), and Fig. 2 does the same for the three MIPS bands
(24, 70 and 160 µm). Each panel is split in two: the upper sub-
panel plots the counts per logarithmic flux interval, dN/d lnSν ,
while the lower sub-panel instead plots SνdN/d lnSν . The latter
is designed to take out much of the trend with flux, in order to
show more clearly the differences between the model and the on-
servational data. In each case we plot three curves for our standard
model: the solid blue line shows the total number counts includ-
ing both extinction and emission by dust, the solid red line shows
the contribution to this from galaxies currently forming stars in a
burst, and the solid green line shows the contribution from all other
galaxies (star-forming or not), which we denote as “quiescent”. In
Fig. 1, we also plot a dashed blue line which shows the predicted
total counts if we ignore absorption and emission from interstellar
dust (emission from dust in the envelopes of AGB stars is still in-
cluded in the stellar contribution, however). In the MIPS bands, the
predicted counts are negligible in the absence of interstellar dust,
so we do not plot them in Fig. 2. In the lower sub-panels, we also
show by a dashed magenta line the prediction from a variant model
which assumes a normal (Kennicutt) IMF for all star formation, but
is otherwise identical to our standard model (which has a top-heavy
IMF in bursts). This variant model fits the local B- and K-band and
60 µm luminosity functions about as well as our standard model,
but dramatically underpredicts the 850 µm number counts. The ob-
served number counts are shown by black symbols with error bars.
Overall, the agreement between the predictions of our stan-
c© 0000 RAS, MNRAS 000, 000–000
8 Lacey et al.
Figure 1. Galaxy differential number counts in the four IRAC bands. The curves show model predictions, while the symbols with error bars show observational
data from Fazio et al. (2004) (with different symbols for data from different survey fields). Each panel is split in two: the upper sub-panel plots the counts
as dN/d lnSν vs Sν , while the lower sub-panel plots SνdN/d lnSν (in units mJy deg
−2) on the same horizontal scale. The upper sub-panels show four
different curves for our standard model - solid blue: total counts including dust extinction and emission; dashed blue: total counts excluding interstellar dust;
solid red: ongoing bursts (including dust); solid green: quiescent galaxies (including dust). The lower sub-panels compare the total counts including dust for
the standard model (solid blue line) with those for a variant model with a normal IMF for all stars (dashed magenta line). The vertical dashed line shows the
estimated confusion limit for the model. (a) 3.6 µm. (b) 4.5 µm. (c) 5.8 µm. (d) 8.0 µm.
dard model and the observed counts is remarkably good, when one
takes account of the fact that no parameters of the model were ad-
justed to improve the fit to any data from Spitzer. Consider first the
results for the IRAC bands, shown in Fig. 1. Here, the agreement
of the model with observations seems best at 3.6 and 8.0 µm, and
somewhat poorer at 5.8 µm. The model predicts somewhat too few
objects at fainter fluxes in all of the IRAC bands. Comparing the red
and green curves, we see that quiescent galaxies rather than bursts
dominate the counts at all observed fluxes in all of the IRAC bands,
but especially at the shorter wavelengths, consistent with the expec-
tation that at 3.6 and 4.5 µm, we are seeing mostly light from old
stellar populations. Comparing the solid and dashed blue lines, we
see that the effects of dust are small at 3.6 and 4.5 µm, with a small
amount of extinction at faint fluxes (and thus higher average red-
shifts), but negligible extinction for brighter fluxes (and thus lower
redshifts). On the other hand, dust has large effects at 8.0 µm, with
dust emission (due to strong PAH features at λ ∼ 6 − 9µm) be-
coming very important at bright fluxes (which correspond to low
average redshifts - see Fig. A1(b) in the Appendix). The 8.0 µm
counts thus are predicted to be dominated by dust emission from
quiescently star-forming galaxies, except at the faintest fluxes. The
counts at 5.8 µm show behaviour which is intermediate, with mild
emission effects at bright fluxes and mild extinction at faint fluxes.
Comparing the solid blue and dashed magenta lines, we see that the
predicted number counts in the IRAC bands are almost the same
whether or not we assume a top-heavy IMF in bursts, consistent
with the counts being dominated by quiescent galaxies.
Consider next the results for the MIPS bands, shown in Fig. 2.
We again see remarkably good agreement of the standard model
with the observational data. The agreement is especially good at
faint fluxes (corresponding to higher redshifts). In particular, the
model matches well the observed 24 µm counts at the “bump”
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 9
around fluxes Sν ∼ 0.1 − 1mJy. Accurate modelling of the PAH
emission features is obviously crucial for modelling the 24 µm
number counts, since the PAH features dominate the flux in the
24 µm band as they are redshifted into the band at z & 0.5.
On the other hand, the standard model overpredicts the number
counts at bright fluxes (corresponding to low redshifts) in all three
MIPS bands. The evolution at these wavelengths predicted by our
ΛCDM-based model thus seems to be not quite as strong as indi-
cated by observations.
In the MIPS bands, emission from galaxies is completely
dominated by dust, which is why no dashed blue lines are shown in
Fig. 2. Comparing the red and green curves, we see that quiescent
(but star-forming) galaxies tend to dominate the number counts in
these bands at brighter fluxes, and bursts at fainter fluxes. This re-
flects the increasing dominance of bursts in the mid- and far-IR
luminosity function at higher redshifts. Comparing the solid blue
and dashed magenta curves, we see that our standard model with a
top-heavy IMF in bursts provides a significantly better overall fit to
the observed 24 µm counts than the variant model with a normal
IMF in bursts (although at the brightest fluxes, the variant model
fits better). The faint number counts at 70 µm also favour the top-
heavy IMF model, while the number counts at 160 µm cover a
smaller flux range, and do not usefully distinguish between the two
variants of our model with different burst IMFs.
We can use our model to predict the flux levels at which
sources should become confused in the different Spitzer bands.
We estimate the confusion limit using the source density cri-
terion (e.g. Vaisanen et al. 2001; Dole et al. 2003): if the tele-
scope has an FWHM beamwidth of θFWHM , we define the ef-
fective beam solid angle as ωbeam = (π/(4 ln 2)) θ
FWHM =
1.13θ2FWHM , and then define the confusion limited flux Sconf
to be such that N(> Sconf ) = 1/(Nbeamωbeam), where
N(> S) is the number per solid angle of sources brighter than
flux S. We choose Nbeam = 20 for the number of beams
per source, which gives similar results to more detailed analy-
ses (e.g. Vaisanen et al. 2001; Dole et al. 2004b). We use values
of the beamsize θFWHM = (1.66, 1.72, 1.88, 1.98) arcsec for
the four IRAC bands (Fazio et al. 2004b) and (5.6, 16.7, 35.2)
arcsec for the three MIPS bands (Dole et al. 2003). Our stan-
dard model then predicts confusion-limited fluxes of Sconf =
(0.62, 0.62, 0.69, 0.70)µJy in the (3.6, 4.5, 5.8, 8.0)µm IRAC
bands, and Sconf = (0.072, 2.6, 43)mJy in the (24, 70, 160)µm
MIPS bands. These confusion estimates for the MIPS bands are
similar to those of Dole et al. (2004b), which were based on ex-
trapolating from the observed counts. These values for the confu-
sion limits are indicated in Figs. 1 and 2 by vertical dashed lines.
Our galaxy evolution model does not compute the contribution
of AGN to the IR luminosities of galaxies. On the other hand, the
observed number counts to which we compare include both normal
galaxies, in which the IR emission is powered by stellar popula-
tions, and AGN, in which there is also IR emission from a dust
torus, which is expected to be most prominent in the mid-IR. How-
ever, multi-wavelength studies using optical, IR and X-ray data in-
dicate that even at 24 µm, the fraction of sources dominated at
that wavelength by AGN is only 10-20% (e.g. Franceschini et al.
2005), and the contribution of AGN-dominated sources in the other
Spitzer bands is likely to be smaller. Therefore we should not make
any serious error by comparing our model predictions directly with
the total number counts, as we have done here.
Figure 2. Galaxy differential number counts in the three MIPS bands. The
curves show model predictions while the symbols with error bars show ob-
servational data. The meaning of the different model lines is the same as
in Fig. 1. (a) 24 µm, with observational data from Papovich et al. (2004).
(b) 70 µm, with observational data from Dole et al. (2004a) (filled sym-
bols), Frayer et al. (2006a) (crosses), and Frayer et al. (2006b) (open sym-
bols). (c) 160 µm (bottom panel), with observational data from Dole et al.
(2004a) (filled symbols) and Frayer et al. (2006a) (crosses).
c© 0000 RAS, MNRAS 000, 000–000
10 Lacey et al.
Figure 3. Predicted evolution of the galaxy luminosity function in our standard model (including dust) at rest-frame wavelengths of (a) 3.6 and (b) 8.0 µm for
redshifts z = 0, 0.5, 1, 1.5, 2 and 3, as shown in the key.
4 EVOLUTION OF THE GALAXY LUMINOSITY
FUNCTION
While galaxy number counts provide interesting constraints on the-
oretical models, it is more physically revealing to compare with
galaxy luminosity functions, since these isolate behaviour at partic-
ular redshifts, luminosities and rest-frame wavelengths. In the fol-
lowing subsections, we compare our model predictions with recent
estimates of luminosity function (LF) evolution based on Spitzer
data.
4.1 Evolution of the galaxy luminosity function at 3-8 µm
We consider first the evolution of the luminosity function in the
wavelength range covered by the IRAC bands, i.e. 3.6-8.0 µm.
Fig. 3 shows what our standard model with a top-heavy IMF in
bursts predicts for LF evolution at rest-frame wavelengths of 3.6
and 8.0 µm for redshifts z = 0 − 3 3. We see that at a rest-frame
wavelength of 3.6 µm, the model LF hardly evolves at all over the
whole redshift range z = 0 − 3. This lack of evolution appears to
be somewhat fortuitous. Galaxy luminosities at a rest-frame wave-
length of 3.6 µm are dominated by the emission from moderately
old stars, but the stellar mass function in the model evolves quite
strongly over the range z = 0 − 3 (as we show in §5). The weak
evolution in the 3.6 µm LF results from a cancellation between a
declining luminosity-to-stellar-mass ratio with increasing time and
increasing stellar masses (see Figs. 13(a) and (e)). On the other
hand, at a rest-frame wavelength of 8.0 µm, the model LF becomes
significantly brighter in going from z = 0 to z = 3. Galaxy lu-
minosities at a rest-frame wavelength of 8.0 µm are dominated by
emission from dust heated by young stars, so this evolution reflects
the increase in star formation activity with increasing redshift (see
Fig. 13(b) in §5).
In Fig. 4, we compare the model predictions for evolu-
tion of the LF at 3.6 µm with observational estimates from
3 In this figure, and in Figs. 4, 5, 8, and 10, the luminosities Lν are calcu-
lated through the corresponding Spitzer passbands.
Babbedge et al. (2006) and Franceschini et al. (2006)4. The
model predictions are given for redshifts z = 0, 0.5 and 1. For
the observational data, the mean redshifts for the different redshift
bins used do not exactly coincide with the model redshifts, so we
plot them with the model output closest in redshift 5. The observa-
tional estimates of the 3.6 µm LF rely on the measured redshifts.
In the case of Babbedge et al. (2006), these are mostly photomet-
ric, using optical and NIR (including 3.6 and 4.5 µm) fluxes, while
for the Franceschini et al. sample, about 50% of the redshifts are
spectroscopic and the remainder photometric. In both samples, the
measured 3.6 µm fluxes were k-corrected to estimate the rest-frame
3.6 µm luminosities.
We see from comparing the blue curve with the observational
data in Fig. 4 that the 3.6 µm LF predicted by our standard model
is in very good agreement with the observations. In particular, the
observational data show very little evolution in the 3.6 µm LF over
the redshift range z = 0 − 1. The largest difference seen is at
z = 1, where the Babbedge et al. data show a tail of objects to
very high luminosities, which is not seen in the model predictions.
However, this tail is not seen in the Franceschini et al. data at the
same redshift, and is also not present in the observational data at
the lower redshifts. More spectroscopic redshifts are needed for the
Babbedge et al. sample to clarify whether this high-luminosity tail
is real. Comparing the red, green and blue lines for the standard
model shows that the model luminosity function is dominated by
quiescent galaxies at low luminosity, but the contribution of bursts
becomes comparable to that of quiescent galaxies at high luminosi-
ties. We have not shown model LFs excluding dust extinction in
this figure, since they are almost identical to the predictions includ-
ing dust. The dashed magenta lines show the predicted LFs for the
4 Babbedge et al. (2006) also compared their measured LFs at 3.6, 8.0 and
24 µm with predictions from a preliminary version of the model described
in this paper
5 Specifically, for z = 0, we compare with the z = 0.1 data from
Babbedge et al. , for z = 0.5 we compare with the z = 0.5 data
from Babbedge et al. and z = 0.3 data from Franceschini et al. , and for
z = 1, we compare with the z = 0.75 (open symbols) and z = 1.25
(filled symbols) data from Babbedge et al. and z = 1.15 data from
Franceschini et al.
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 11
Figure 4. Predicted evolution of the galaxy luminosity function at rest-
frame 3.6 µm compared to observational data. The different panels show
redshifts (a) z = 0, (b) z = 0.5 and (c) z = 1. The predictions for our stan-
dard model are shown by the blue line, with the red and green lines showing
the separate contributions from ongoing bursts and quiescent galaxies. The
dashed magenta line shows the prediction for a variant model with a nor-
mal IMF for all stars. The error bars on the model lines indicate the Poisson
uncertainties due to the finite number of galaxies simulated. The black sym-
bols with error bars show observational data from Babbedge et al. (2006)
(open circles and triangles, for z = 0, 0.5 and 1) and Franceschini et al.
(2006) (filled squares, for z = 0.5 and 1).
Figure 5. Predicted evolution of the galaxy luminosity function at rest-
frame 8.0 µm compared to observational data. The different panels show
redshifts (a) z = 0, (b) z = 1 and (c) z = 2. The coloured lines showing
the model predictions have the same meaning as in Fig. 4. The black sym-
bols with error bars show observational data from Babbedge et al. (2006)
(open circles for z = 0 and 0.7, triangles for z = 1.2), Huang et al. (2007)
(filled circles for z = 0) and Caputi et al. (2007) (filled circles for z = 1
and 2). The observed LFs are for normal galaxies and exclude AGN.
c© 0000 RAS, MNRAS 000, 000–000
12 Lacey et al.
variant model with a normal IMF in bursts. We see that these differ
only slightly from our standard model, but are a somewhat poorer
fit to the observational data at higher luminosities.
In Fig. 5 we show a similar comparison for the LF evolution at
a rest-frame wavelength of 8 µm. The model predictions are given
for redshifts z = 0, 1 and 2, and are compared with observational
estimates by Huang et al. (2007) (for z ∼ 0), Babbedge et al.
(2006) (for z ∼ 0 and z ∼ 1) and Caputi et al. (2007) (for z ∼ 1
and z ∼ 2). These papers all classified objects in their samples
as either galaxies or AGN, and then computed separate LFs for
the two types of objects 6. Our model does not make any predic-
tions for AGN, so we compare our model predictions with the ob-
served LFs for objects classified as galaxies only. We see that for
redshifts around z = 1, the observed LFs from Babbedge et al.
and Caputi et al. are in very poor agreement with each other, with
the Caputi et al. LF being around 10 times higher in number den-
sity at the same luminosity. This differerence presumably results
from some combination of: (a) different methods of classifying ob-
jects as galaxies or AGN (Babbedge et al. used only optical and
IR fluxes to do this, while Caputi et al. also used X-ray data); (b)
different photometric redshift estimators; and (c) different meth-
ods for k-correcting luminosities to a rest-frame wavelength of 8
µm. There are smaller differences between the Huang et al. and
Babbedge et al. LFs at z ∼ 0. Futher observational investigation
appears to be necessary to resolve these issues. Our standard model
is in reasonable agreement with the Babbedge et al. observed LF
at z ∼ 0, and with the Caputi et al. observed LFs at z ∼ 1 and
z ∼ 2, but not with the Babbedge et al. observed LF at z ∼ 1. The
comparison with Caputi et al. favours our standard model with a
top-heavy IMF in starbursts over the variant model with a normal
4.2 Evolution of the galaxy luminosity function at 12-24 µm
In this subsection, we consider the evolution of the galaxy lumi-
nosity function at mid-IR wavelengths, and compare with data ob-
tained using mainly the MIPS 24 µm band.
Fig. 6 shows what our standard model with a top-heavy IMF
in bursts predicts for the evolution of the galaxy LF at rest-frame
wavelengths of 15 and 24 µm for redshifts z = 0 − 3 7. At rest-
frame wavelengths of 15 and 24 µm, galaxy luminosities are typi-
cally dominated by the continuum emission from warm dust grains
heated by young stars (although PAH emission is also significant
at some nearby wavelengths). Fig. 6 shows strong evolution in the
model LFs over the redshift range z = 0− 3 at both wavelengths,
reflecting both the increase in star formation activity with increas-
ing redshift (see Fig. 13(b)) and the increasing dominance of the
burst mode of star formation, for which the top-heavy IMF fur-
ther boosts the mid- and far-IR luminosities compared to a normal
IMF. Comparing Fig. 6 with Fig. 3(a), we also see a difference in
the shape of the bright end of the LF: at 3.6 µm, where the LF is
dominated by emission from stars, the bright end cuts off roughly
6 Note that a variety of criteria have been used for classifying observed IR
sources as AGN or normal galaxies, and these do not all give equivalent
results. Even if an object is classified as an AGN, it is also not clear that in
all cases the AGN luminosity dominates over that of the host galaxy in all
Spitzer bands
7 In this figure, and in Figs. 7 and 8, the 24 µm luminosities are calculated
through the corresponding MIPS passband, while the 15 µm luminosities
are calculated through a top-hat filter with a fractional width of 10% in
wavelength.
Figure 8. Predicted evolution of the galaxy luminosity function at rest-
frame wavelength 24 µm compared to observational data from Shupe et al.
(1998) (at z = 0, open symbols) and from Babbedge et al. (2006) (for the
same redshifts as in Fig. 4). The meaning of the curves showing the model
predictions is the same as in Fig. 4. (a) z = 0, (b) z = 0.5 and (c) z = 1.
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 13
Figure 6. Predicted evolution of the galaxy luminosity function in our standard model at rest-frame wavelengths(a) 15 µm (left) and (b) 24 µm (right) for
redshifts z = 0, 0.5, 1, 1.5, 2 and 3, as shown in the key.
exponentially, while at 15 and 24 µm, where the LF is dominated
by emission from warm dust, the bright end declines more gradu-
ally, roughly as a power-law. This difference reflects the difference
in shape of the galaxy stellar mass function (GSMF) and the galaxy
star formation rate distribution (GSFRD). The GSMF shows an
exponential-like cutoff at high masses, while the GSFRD shows a
more gradual cutoff at high SFRs because of starbursts triggered by
galaxy mergers (see Figs. 13(a) and (b) in §5). This difference was
noticed earlier by observers comparing optical and far-IR LFs of
galaxies, but its origin was not understood (Lawrence et al. 1986;
Soifer et al. 1987b).
In Fig. 7, we compare the model LFs at rest-frame wave-
lengths 12 and 15 µm with observational estimates. For z = 0,
we plot the observational estimates from Soifer & Neugebauer
(1991) and Rush et al. (1993), based on IRAS 12 µm data (with
AGN removed). For z = 0.5 − 1 and z = 1.5 − 2.5, we
plot the data of Le Floc’h et al. (2005) and Perez-Gonzalez et al.
(2005) respectively, which were obtained from galaxy samples
selected on Spitzer 24 µm flux. Le Floc’h et al. k-corrected
their measured 24 µm fluxes to 15 µm rest-frame luminosities,
while Perez-Gonzalez et al. k-corrected to 12 µm rest-frame8.
Le Floc’h et al. obtained most of their redshifts from photometric
redshifts based on optical data, while Perez-Gonzalez et al. used a
new photometric redshift technique based on fitting empirical SEDs
to all of the available broad-band data from the far-UV to 24 µm,
and also removed “extreme” AGN from their observed LF. Note
that the redshifts for the observed LFs do not exactly coincide with
model redshifts in all cases, but are close.
We see from comparing the blue line to the observational dat-
apoints in Fig. 7 that our standard model with a top-heavy IMF in
bursts fits the observations remarkably well up to z = 2. In partic-
ular, the model matches the strong evolution in the mid-IR LF seen
8 The exact passband used for the model LF in each panel depends on
which observational data we are comparing with. For z = 0, we use the
IRAS 12 µm passband; at z = 0.5 and z = 1 we use a top-hat passband
centred at 15 µm; and at z = 1.5, 2 and 2.5, we use a top-hat passband
centred at 12 µm (both top-hat passbands having fractional width 10% in
wavelength).
in the observational data. The model falls below the observational
data at z = 2.5, but here both the photometric redshifts and the
k-corrections are probably the most uncertain. The standard model
also does not provide a perfect fit to the z = 0 data, predicting
somewhat too many very bright galaxies and somewhat too few
very faint galaxies (though the latter discrepancy might be affected
by local galaxy clustering in the IRAS data). Comparing the red,
green and blue lines for the standard model in the figure, we see
that the bright end of the 12 or 15 µm LF is dominated by bursts
at all redshifts. The figure also shows by a dashed magenta line the
predictions for the variant model with a normal IMF in bursts. This
latter model predicts much less evolution in the bright end of the LF
than is observed. This comparison thus strongly favours the model
with the top-heavy IMF in bursts.
Finally, in Fig. 8, we carry out a similar comparison of the evo-
lution of predicted and observed LFs at a rest-frame wavelength
of 24 µm over the redshift range z = 0 − 1, in this case com-
paring with observational estimates from Shupe et al. (1998) (for
z = 0), based on IRAS data, and from Babbedge et al. (2006) (for
z = 0 − 1), based on Spitzer data9. The galaxy redshifts for the
Babbedge et al. data were obtained in the same way as for the 3.6
µm LFs shown in Fig. 4, and the luminosities were k-corrected
from observer-frame 24 µm to rest-frame 24 µm. The LF plotted
from Babbedge et al. is that for normal galaxies, with AGN ex-
cluded.
The conclusions from comparing the model with the 24 µm
LFs are similar to those from the comparison with the 12 and 15
µm LFs. The data favour our standard model over the variant with
a normal IMF in bursts (except possibly for z = 0.5), as the latter
predicts too little evolution at the bright end. At z = 0, the model
fits the 24 µm data rather better than for the corresponding com-
parison at 12 µm. On the other hand, at z = 0.5 and z = 1, the
model LF is a somewhat worse fit to the observational data at 24
µm than at 15 µm. These differences between the 12/15 and 24
µm comparisons might result from the different photometric red-
shifts and k-corrections used in the observational samples in the
9 The model luminosities are all computed through the Spitzer 24 µm pass-
band.
c© 0000 RAS, MNRAS 000, 000–000
14 Lacey et al.
Figure 7. Predicted evolution of the galaxy luminosity function at rest-frame wavelength 12 or 15 µm compared to observational data. The different panels
show redshifts: (a) z = 0, (b) z = 0.5, (c) z = 1, (d) z = 1.5, (e) z = 2 and (f) z = 2.5. The meaning of the curves showing the model predictions is the
same as in Fig. 4. In panel (a), the predictions at 12µm are compared to observational determinations from Soifer & Neugebauer (1991) (open symbols) and
Rush et al. (1993) (filled symbols) based on IRAS data. In panels (b) and (c), the predictions at 15µm are compared to observational data from Le Floc’h et al.
(2005). In panels (d), (e) and (f), the predictions at 12µm are compared to observational data from Perez-Gonzalez et al. (2005).
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 15
Figure 9. The predicted galaxy luminosity function at 60 µm compared
to observational data from IRAS. The meaning of the different lines is
the same as in Fig. 4. The black symbols show observational data from
Saunders et al. (1990) (crosses), Soifer & Neugebauer (1991) (open cir-
cles), and Takeuchi et al. (2003) (filled circles).
two cases. Alternatively, they might result from problems in mod-
elling the dust SEDs in the complex mid-IR range.
4.3 Evolution of the galaxy luminosity function at 70-160 µm
We now briefly consider the evolution of the luminosity function
in the far-IR. The far-IR is the wavelength range where most of
the luminosity from dust in normal galaxies is emitted. The lo-
cal 60 µm luminosity function was very well measured by sur-
veys with IRAS, and so is commonly used as a starting point or
benchmark for modelling the evolution of the galaxy population in
the far-IR. We therefore present in Fig. 9 the model prediction for
the 60 µm luminosity function at z = 0, compared with obser-
vational data from Saunders et al. (1990), Soifer & Neugebauer
(1991) and Takeuchi et al. (2003). As discussed in Baugh et al.
(2005), the local 60 µm luminosity function was used as one of the
primary constraints in fixing the parameters of our galaxy forma-
tion model, and the figure shows that our standard model provides
a good match to the data. The variant model with a normal IMF in
bursts underpredicts the abundance of the brightest 60 µm galaxies.
In Fig. 10, we show the model predictions for the evolution of
the luminosity function in the two longer wavelength MIPS bands,
at rest-frame wavelengths of 70 and 160 µm, from z = 0 to z = 3.
At 70 µm, the luminosity function at high luminosities is predicted
to brighten by about a factor 10 going from z = 0 to z = 2. This
is about a factor 2 less than the brightening predicted in the mid-IR
at 15 µm (compare to Fig. 6), but nearly a factor 2 more evolution
than is predicted at 160 µm. These differences between the amount
of evolution seen at different IR wavelengths reflect evolution in
the shapes of the SEDs of the galaxies responsible for the bulk of
the IR emission. No observational estimates of the evolution of the
luminosity function at 70 and 160 µm have yet been published,
but they are expected to be forthcoming from ongoing surveys with
Spitzer.
Figure 11. Predicted evolution of the total mid+far-IR (8-1000 µm) galaxy
luminosity function for our standard model, for redshifts z = 0, 1, 2, 3, 4
and 6, as shown in the key.
4.4 Evolution of the total mid+far-IR luminosity function
The total mid+far IR luminosity of a galaxy, LIR, integrated over
the whole wavelength range 8-1000 µm, is a very good approxi-
mation to the total luminosity emitted by interstellar dust grains in
all galaxies except those with very small dust contents. In galax-
ies with significant star formation, LIR is mostly powered by dust
heated by young stars, and so provides a quantitative indicator of
the amount of dust-obscured star formation which is independent
of the shape of the IR SED (though still subject to uncertainties
about the IMF). The evolution of the luminosity function in LIR
is therefore a very interesting quantity to compare between models
and observations. We show in Fig. 11 what our standard model pre-
dicts for the evolution of the IR LF over the range z = 0 − 6. We
see that the model predicts substantial evolution in this LF, with the
high luminosity end brightening by a factor ∼ 10 from z = 0 to
z = 2, followed by a “plateau” from z = 2 to z = 4, and a decline
from z = 4 to z = 6.
In Fig. 12, we compare our model predictions with existing
observational estimates of the total IR LF for z = 0 − 2. These
observational estimates are only robust for z = 0, where they are
based on IRAS measurements covering the wavelength range 12-
100 µm. At all of the higher redshifts plotted, the observational
estimates are based on measurements of the mid-IR luminosity de-
rived from Spitzer 24 µm fluxes, converted to total IR luminosi-
ties by assuming SED shapes for the mid- to far-IR emission. The
bolometric correction from the observed mid-IR luminosity to the
inferred total IR luminosity is typically a factor ∼ 10, and is sig-
nificantly uncertain. Therefore, the most robust way to compare
the models with the observations is to compare them at the mid-
IR wavelengths where the measurements are actually made, as we
have done in §4.1 and §4.2. Nonetheless, if we take the observa-
tional determinations at face value, then we see that observed evo-
lution of the total IR LF agrees remarkably well with the predic-
tions of our standard model with a top-heavy IMF. On the other
hand, the variant model with a normal IMF predicts far too few
c© 0000 RAS, MNRAS 000, 000–000
16 Lacey et al.
Figure 10. Predicted evolution of the galaxy luminosity function in our standard model (including dust) at rest-frame wavelengths (a) 70 µm and (b) 160 µm,
for redshifts z = 0, 0.5, 1, 1.5, 2 and 3, as shown in the key.
high LIR galaxies at higher z, and is strongly disfavoured by the
existing data.
5 INFERRING STELLAR MASSES AND STAR
FORMATION RATES FROM Spitzer DATA
In this section, we consider what the models imply about how well
we can infer the stellar masses and star formation rates (SFRs) in
galaxies from measurements of rest-frame IR luminosities. The top
two panels of Fig. 13 show the predicted galaxy stellar mass func-
tion (GSMF, left panel) and galaxy star formation rate distribution
(GSFRD, right panel), for redshifts z = 0− 6. We see that the pre-
dicted stellar mass function shows dramatic evolution over this red-
shift range, with a monotonic decline in the number of high-mass
galaxies with increasing redshift. On the other hand, the SFR distri-
bution shows much less dramatic evolution over this redshift range,
with a mild increase in the number of high-SFR objects up to z ∼ 3,
followed by a decline above that. The lower four panels in Fig. 13
show the relation in the models between stellar masses and SFRs
and rest-frame luminosities at different IR wavelengths. (Note that
in all cases, luminosities are measured in units of the bolometric
solar luminosity.) The middle and bottom left panels respectively
show the mean ratio of luminosity in the rest-frame K (2.2µm) or
3.6 µm bands to stellar mass as a function of stellar mass. The
middle and bottom right panels respectively show the mean ratio
of total mid+far-IR (8− 1000µm) or rest-frame 15 µm luminosity
to SFR as a function of SFR. (The mean L/M∗ or L/SFR ratios
plotted are computed by dividing the total luminosity by the total
mass or SFR, in each bin of mass or SFR.)
The near-IR luminosity is often used as a tracer of stellar mass.
The left panels of Fig. 13 show that the L/M∗ ratio varies strongly
with redshift, reflecting the difference in the ages of the stellar pop-
ulations. At higher redshifts it also shows a significant dependence
on stellar mass, presumably reflecting a trend of age with mass.
However, the variation of mean L/M∗ with redshift is seen to be
much smaller at 2.2 µm than at 3.6 µm, implying that the rest-
frame K-band light should provide a more robust estimator of stel-
lar mass than the light at longer wavelengths. The differences be-
tween L/M∗ values at 2.2 µm and 3.6 µm reflect the larger contri-
bution from AGB compared to RGB stars at the longer wavelength.
AGB stars have higher masses and younger ages than RGB stars,
and so are more sensitive to star formation at recent epochs. The
scatter in L/M∗ at a given mass is also found in the models to in-
crease with redshift. In the K-band, it increases from ∼ 40% at
z ∼ 0 to a factor ∼ 3 at z ∼ 6. The large scatter at high redshifts
results in part from having two different IMFs.
The luminosity in the mid- and far-IR is widely used as a
tracer of dust-obscured star formation (although in galaxies with
very low star formation rates, the dust heating can be dominated
by older stars). The total mid+far-IR (rest-frame 8-1000 µm) lumi-
nosity is expected to provide a more robust tracer of star formation
than the luminosity at any single IR wavelength, since the shape of
the SED of dust emission depends on the dust temperature distri-
bution (as well as on the dust grain properties). This is borne out
by our model predictions. The middle right panel of Fig. 13 shows
that the LIR/SFR ratio depends weakly on both SFR and red-
shift. This behaviour results mostly from having different IMFs in
the model in quiescent and bursting galaxies, with the fractional
contribution of the bursts increasing both with SFR and with red-
shift. If we look at quiescent and bursting galaxies separately, we
find roughly constant ratios LIR/SFR ≈ 6 × 10
9h−1L⊙/M⊙
and LIR/SFR = 2× 10
10h−1L⊙/M⊙ respectively, for galaxies
where LIR is powered mostly by young stars. However, there is
also a trend at lower redshift for LIR/SFR to be larger at lower
SFR - this reflects the larger fraction of dust heating from older
stars in galaxies with lower SFRs, which more than compensates
for the lower average dust obscuration in these galaxies. The lower
right panel of Fig. 13 shows that the L/SFR ratio in the mid-IR
(in this case at 15 µm in the rest-frame) shows more variation with
SFR and redshift than the ratio for the total IR luminosity. This re-
flects the variation in the mid- to far-IR SED shapes in the model.
The scatter in the L/SFR ratio is roughly a factor 2 around the
average relation for the total IR luminosity, but is larger for the 15
µm luminosity.
The results of this section illustrate why it is not straightfor-
ward to compare theoretical predictions for the evolution of the
galaxy stellar mass function and star formation rate distribution (or
even the stellar mass and star formation rate densities) with obser-
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 17
Figure 12. Predicted evolution of the total mid+far IR (8-1000µm) galaxy LF compared to observational data. The different panels show redshifts (a) z = 0,
(b) z = 0.5, (c) z = 1 and (d) z = 2. For z = 0, we compare with observational data from Sanders et al. (2003) (filled symbols) and Takeuchi et al. (2003)
(open symbols, converting his 60 µm LF to a total IR LF assuming a constant conversion factor, LIR/νLν(60µm) = 2.5). We compare with data from
Le Floc’h et al. (2005) for z = 0.5 and z = 1 (filled and open symbols), and with Caputi et al. (2007) for z = 1 and z = 2 (crosses).
vational estimates. In addition to assumptions about galaxy star for-
mation histories and metallicities (for stellar mass estimates), and
about the SED shapes for dust emission (for SFR estimates from
IR and sub-mm data), observational estimates all rest on some as-
sumed form for the IMF. If the IMF assumed in the observational
analysis is different from the true IMF, the observational estimates
for stellar masses and SFRs can be wrong by large factors. If the
IMFs differ only below 1M⊙, then one can apply a simple rescal-
ing to relate stellar mass and SFR estimates for different IMFs.
However, if our current galaxy formation model is correct, stars
form with different IMFs in quiescent disks and in merger-driven
bursts, and so no observational estimate based on assuming a sin-
gle IMF can give the correct GSMFs and GSFRDs, nor the correct
stellar mass and SFR densities. A direct comparison of the GSMF
and GSFRD evolution predicted by our model with observational
estimates is therefore not meaningful. Instead, the comparison be-
tween models and observations must be made via directly observ-
able (rather than inferred) quantities, such as the K-band luminosi-
ties to constrain stellar masses, and the total IR luminosities to con-
strain SFRs.
6 CONCLUSIONS
We have computed predictions for the evolution of the galaxy pop-
ulation at infrared wavelengths using a detailed model of hierar-
chical galaxy formation and of the reprocessing of starlight by
dust, and compared these predictions with observational data from
the Spitzer Space Telescope. We calculated galaxy formation in
the framework of the ΛCDM model using the GALFORM semi-
analytical model, which includes physical treatments of the hier-
archical assembly of dark matter halos, shock-heating and cool-
ing of gas, star formation, feedback from supernova explosions
and photo-ionization of the IGM, galaxy mergers and chemical en-
richment. We computed the IR luminosities and SEDs of galaxies
using the GRASIL multi-wavelength spectrophotometric model,
which computes the luminosities of the stellar populations in galax-
c© 0000 RAS, MNRAS 000, 000–000
18 Lacey et al.
Figure 13. Model predictions for properties related to stellar masses (left column) and star formation rates (right column), for redshifts z = 0, 1, 2, 3, 4, and
6: (a) galaxy stellar mass function (GSMF); (b) galaxy star formation rate distribution (GSFRD); (c) mean ratio of rest-frame K-band luminosity to stellar
mass, as a function of stellar mass; (d) mean ratio of total mid+far IR luminosity to SFR, as a function of SFR; (e) mean ratio of rest-frame 3.6 µm luminosity
to stellar mass, as a function of stellar mass; (f) mean ratio of rest-frame 15 µm luminosity to SFR, as a function of SFR. (The 15 µm luminosity is here
calculated through top-hat filter with a fractional wavelength width of 10%.)
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 19
ies, and then the reprocessing of this radiation by dust, including
radiative transfer through a two-phase dust medium, and a self-
consistent calculation of the distribution of grain temperatures in
each galaxy based on a local balance between heating and cool-
ing. The GRASIL model includes a treatment of the emission from
PAH molecules, which is essential for understanding the mid-IR
emission from galaxies.
Our galaxy formation model incorporates two different IMFs:
quiescent star formation in galaxy disks occurs with a normal so-
lar neighbourhood IMF, but star formation in bursts triggered by
galaxy mergers happens with a top-heavy x = 0 IMF. In a previ-
ous paper (Baugh et al. 2005), we found that the top-heavy IMF in
bursts was required in order that the model reproduces the observed
number counts of the faint sub-mm galaxies detected at 850 µm,
which are typically ultra-luminous starbursts at z ∼ 2, with total IR
luminosities LIR ∼ 10
− 1013L⊙. This conclusion was arrived
at following a search of a large grid of model parameters, with the
imposition of a variety of detailed observational constraints. The
parameters in the Baugh et al. (2005) model were chosen before
the publication of any results from Spitzer, without reference to any
IR data apart from the local 60 µm luminosity function and the 850
µm galaxy counts. We have kept the same parameter values in the
present paper, in order to test what the same model predicts at other
wavelengths and other redshifts. By doing this, we hope to address
the criticism made of many semi-analytical models that they have
no predictive power, because their parameters are always adjusted
to match the observational data being analysed at that instant.
We first compared the predictions from our model with the
galaxy number counts measured in all 7 Spitzer bands, from 3.6 to
160 µm. We found broad agreement between the model and the
observations. In the 4 IRAC bands (3.6-8.0 µm), where the counts
are mostly dominated by emission from older stellar populations,
we found that the predicted counts were insensitive to whether we
had a top-heavy or normal IMF in bursts. On the other hand, in
the MIPS bands (24-160 µm), where the counts are dominated by
emission from dust in star-forming galaxies, the predicted counts
are more sensitive to the choice of IMF, and the counts are fit better
by the model with a top-heavy IMF. We next investigated the evo-
lution of the galaxy luminosity function at IR wavelengths, where
several groups have now used Spitzer data to try to measure the
evolution of the galaxy luminosity function over the redshift range
z ∼ 0− 2, at rest-frame wavelengths from 3.6 to 24 µm.
Our model predicts that at mid- and far-IR rest-frame wave-
lengths, the luminosity function evolution is very sensitive to the
choice of IMF in bursts. We found that our standard model with a
top-heavy IMF in bursts fits the measured evolution of the mid-IR
luminosity function remarkably well (when allowance is made for
complexity of predicting dust emission in the mid-IR), without any
adjustment of the parameters. On the other hand, a model with a
normal IMF in bursts predicts far too little evolution in the mid-IR
luminosity function compared to what is observed. We made a sim-
ilar comparison with the evolution of the total IR luminosity func-
tion, where in the case of the observations, the total IR luminosities
at high redshifts have been inferred from the 24 µm fluxes by fit-
ting SEDs, and reached the same conclusion. The evolution of the
galaxy luminosity function in the mid-IR found by Spitzer thus sup-
ports our original conclusion about the need for a top-heavy IMF
in bursts, which was based only on the sub-mm counts. This con-
clusion will be further tested by ongoing Spitzer surveys at longer
wavelengths. To assist this, we have also presented predictions for
the evolution of the luminosity function in the Spitzer 70µm and
160µm bands.
We have also presented predictions for the evolution of the
stellar mass function and star formation rate distribution of galax-
ies. We investigated how the L/M∗ and L/SFR ratios varied with
galaxy mass, SFR and redshift in different IR wavelength ranges,
and considered the implications for observational estimates of stel-
lar masses and SFRs from IR observations. Even in the near-IR, the
predicted variations inL/M∗ with mass and redshift can be surpris-
ingly large. The variations in L/M∗ are much larger at a rest-frame
wavelength of 3.6 µm than at 2.2 µm, implying that the 2.2 µm
luminosity is a more robust tracer of stellar mass.
Finally, we have presented in an Appendix the predictions of
our model for the redshift distributions of galaxies selected at dif-
ferent IR fluxes in the Spitzer bands.
One significant limitation of our model is that it does not in-
clude the effects of AGN. Two effects are relevant here. The first is
feedback from AGN on galaxy formation. In several recent galaxy
formation models, AGN feedback is invoked to prevent the forma-
tion of too many massive galaxies at the present day. In the model
presented here, we instead posit feedback from supernova-driven
galactic superwinds, which perform a similar role to AGN feed-
back in suppressing the formation of very massive galaxies. Both
the superwind and AGN feedback models include free parameters
which are tuned to give a match to the present-day optical galaxy
luminosity function. However, the redshift dependence of the feed-
back will be different between our superwind model and the various
AGN feedback models, so in general they will all predict different
evolution of the galaxy population with redshift. We will investi-
gate galaxy evolution in the IR in a model with AGN feedback in
a future paper. The second effect of AGN which we have not in-
cluded is the emission from AGN and their associated dust tori. In
order to compensate for this, we have wherever possible compared
our model predictions with observations from which the AGN con-
tribution has been subtracted out. This was possible for most of our
comparisons of luminosity function evolution. This was not possi-
ble for the number counts comparisons, but in this case the contri-
bution from AGN is thought (based on observations) to be a small
fraction of the total over the flux range explored by Spitzer, even in
the mid-IR where the dust tori are the most prominent. We therefore
believe that emission from AGN does not seriously affect our con-
clusions about the IR evolution of star-forming galaxies. We hope
to include AGN emission directly into our models in the future.
We have thus shown that Spitzer data provide a stringent test of
galaxy formation theory, by probing galaxy evolution, constraining
star formation rates and the role of dust to z ∼ 2. We find that
an ab initio ΛCDM model gives an acceptable fit to the Spitzer
data provided that ∼ 10% of the stars in galaxies today formed
in bursts of star formation with a top-heavy IMF. Future facilities
like Herschel, SPICA, JWST and ALMA will continue to exploit the
valuable information on galaxy formation contained in the IR part
of the electromagnetic spectrum.
ACKNOWLEDGEMENTS
We thank T. Babbedge, K. Caputi, A. Franceschini, E. Le Floch,
and P. Perez-Gonzalez, for providing us with their observational
data in a convenient form. CMB acknowledges the receipt of a
Royal Society University Research Fellowship. CSF is the recip-
ient of a Royal Society Wolfson Research Merit Award. This work
was also supported by the PPARC rolling grant for extragalactic
astronomy and cosmology at Durham.
c© 0000 RAS, MNRAS 000, 000–000
20 Lacey et al.
REFERENCES
Almeida, C., Baugh, C.M., Lacey, C.G, 2007, MNRAS, 376, 1711
Babbedge T. S. R., et al. , 2006, MNRAS, 676
Barnes, J., 1998, in Galaxies: Interactions and Induced Star For-
mation, Saas-Fee Advanced Course 26. Lecture Notes 1996.
Swiss Society for Astrophysics and Astronomy, XIV, Edited by
R. C. Kennicutt, Jr. F. Schweizer, J. E. Barnes, D. Friedli, L.
Martinet, and D. Pfenniger., Springer-Verlag Berlin/Heidelberg,
p.275
Baugh C.M., Lacey C.G., Frenk C.S., Granato G.L., Silva L.,
Bressan A., Benson A.J., Cole S., 2005, MNRAS, 356, 1191 (Pa-
per I)
Baugh, C.M., 2006, Rep.Prog.Phys., 69, 3101
Benson, A.J., Lacey, C.G., Frenk, C.S., Baugh, C.M., Cole, S.,
2002, MNRAS, 333, 156
Benson, A.J., Bower, R.G., Frenk, C.S., Lacey, C.G., Baugh,
C.M., Cole, S., 2003, ApJ, 599, 38
Birnboim, Y., Dekel, A., 2003, MNRAS, 345, 349.
Blain A.W., Jameson A., Smail I., Longair M.S., Kneib J.P., Ivison
R.J., 1999b, MNRAS, 309, 715.
Bower, R.G., Benson, A.J., Malbon, R., Helly, J.C., Frenk, C.S.,
Baugh, C.M., Cole, S., Lacey, C.G., 2006, MNRAS, 370, 645
Bressan, A., Granato, G.L., Silva, L., 1998, A&A, 332, 135
Bressan, A., Silva, L., Granato, G.L., 2002, A&A, 392, 377
Caputi K. I., et al. , 2006, ApJ, 637, 727
Caputi, K.I., Lagache, G., Yan, L., Dole, H., Bavouzet, N., Le
Floc’h, E. Choi, P.I., Helou, G., Reddy, N., 2007, ApJ, 660, 97
Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., Blaizot, J.,
MNRAS, 370, 165
Chary R., Elbaz D., 2001, ApJ, 556, 562.
Chapman, S. C., Blain, A. W., Smail, Ian, Ivison, R. J., 2005, ApJ,
622, 772
Cole S., Lacey C.G., Baugh C.M., Frenk C.S., 2000, MNRAS,
319, 168.
Croton, D.J., Springel, V., White, S.D.M., De Lucia, G., Frenk,
C.S., Gao, L., Jenkins, A., Kauffmann, G., Navarro, J.F.,
Yoshida, N., 2006, MNRAS, 365, 11
Devriendt J.E.G., Guiderdoni B., 2000, A&A, 363, 851.
Dole, H., Gispert, R., Lagache, G., Puget, J.-L., Bouchet, F. R.,
Cesarsky, C., Ciliegi, P., Clements, D. L., Dennefeld, M., Desert,
F.-X., and 10 coauthors, 2001, A&A, 372, 364
Dole, H., Lagache, G., Puget, J.-L., 2003, ApJ, 585, 617
Dole, H., Le Floc’h, E., Perez-Gonzlez, P. G., Papovich, C.,
Egami, E., Lagache, G., Alonso-Herrero, A., Engelbracht, C. W.,
and 16 coauthors, 2004, ApJS, 154, 87
Dole, H., et al. , 2004b, ApJS, 154, 93
Elbaz, D., Cesarsky, C. J., Fadda, D., Aussel, H., Desert, F. X.,
Franceschini, A., Flores, H., Harwit, M., Puget, J. L., Starck, J.
L., and 4 coauthors, 1999, A&A 351, L37
Elbaz, D., Cesarsky, C. J., Chanial, P., Aussel, H., Franceschini,
A., Fadda, D., Chary, R. R., 2002, A&A, 384, 848
Eyles, Laurence P., Bunker, Andrew J., Stanway, Elizabeth R.,
Lacy, Mark, Ellis, Richard S., Doherty, Michelle, 2005, MN-
RAS, 364, 443
Fardal, M.A., Katz, N., Weinberg, D.H., Dave, R., 2006, submit-
ted to MNRAS (astro-ph/0604534)
Fazio, G. G., Ashby, M. L. N., Barmby, P., Hora, J. L., Huang, J.-
S., Pahre, M. A., Wang, Z., Willner, S. P., and 7 coauthors 2004,
ApJS, 154, 39
Fazio, G. G., et al. 2004, ApJS, 154, 10
Figer, D.F., Kim, S.S., Morris, M., Serabyn, E., Rich, R.M.,
McLean, I.S., 1999, ApJ, 525, 750
Franceschini, A., Aussel, H., Cesarsky, C. J., Elbaz, D., Fadda, D.,
2001, A&A 378, 1
Franceschini et al. , 2005, AJ, 129, 2074
Franceschini et al. , 2006, A&A, 453, 397
Frayer D. T., et al. , 2006, AJ, 131, 250
Frayer, D. T., Huynh, M. T., Chary, R., Dickinson, M., Elbaz, D.,
Fadda, D., Surace, J. A., Teplitz, H. I., Yan, L. Mobasher, B.
2006b, ApJ, 647, L9
Granato, G.L., Lacey C.G., Silva, L., Bressan, A., Baugh C.M.,
Cole S., Frenk C.S., 2000, ApJ, 542, 710
Granato G.L., De Zotti G., Silva L., Bressan A., Danese L., 2004,
ApJ, 600, 580.
Gruppioni, C., Lari, C., Pozzi, F., Zamorani, G., Franceschini, A.,
Oliver, S., Rowan-Robinson, M., Serjeant, S., 2002, MNRAS,
335, 831
Gruppioni, C., Pozzi, F., Lari, C., Oliver, S., Rodighiero, G., 2005,
ApJ, 618, L9
Guiderdoni B., Hivon E., Bouchet F.R., Maffei B., 1998, MNRAS,
295, 877
Harayama, Y., Eisenhauer, F., Martins, F., 2007
(arXiv:0710.2882)
Hauser, M. G., Arendt, R. G., Kelsall, T., Dwek, E., Odegard, N.,
Weiland, J. L., Freudenreich, H. T., Reach, W. T., Silverberg, R.
F., Moseley, S. H., and 8 coauthors, 1998, ApJ, 508, 25
Helly, J.C., Cole, S., Frenk, C.S., Baugh, C.M., Benson, A., Lacey,
C., 2003, MNRAS, 338, 913
Huang, J.-S., et al. 2007, ApJ 664, 840
Hughes D.H., et al., 1998, Nat, 394, 241.
Kauffmann, G., White, S. D. M., Guiderdoni, B., 1993 MNRAS,
264, 201
Kauffmann G., 1996, MNRAS, 281, 487.
Kaviani, A., Haehnelt, M. G., Kauffmann, G., 2003, MNRAS,
340, 739
Kennicutt R.C., 1983, ApJ, 272, 54.
Kogut A., et al. (the WMAP team), 2003, ApJS, 148, 161.
Kroupa, P., 2001, MNRAS, 322, 231
Lacey, C., Cole, S., 1993, MNRAS, 262, 627
Lagache, G., Dole, H., Puget, J.-L., 2003, MNRAS, 338, 555
Larson, R.B., 1998, MNRAS, 301, 569
Larson, R.B., 2005, MNRAS, 359, 211
Lawrence, A., Walker, D., Rowan-Robinson, M., Leech, K. J.,
Penston, M. V., 1986, MNRAS, 219, 687
Le Delliou, M., Lacey, C., Baugh, C.M., Guiderdoni, B., Bacon,
R., Courtois, H., Sousbie, T., Morris, S.L., 2005, MNRAS, 357,
Le Delliou, M., Lacey, C., Baugh, C.M., Morris, S.L., 2006, MN-
RAS, 365, 712
Le Floc’h E., Papovich C., Dole H., Bell E., Lagache G., Rieke
G., Egami E., Perez-Gonzalez P., and 9 co-authors, 2005, ApJ,
632, 169
Li, A., Draine, B.T., 2001, ApJ, 554, 778
McCrady, N., Gilbert, A.M., Graham, J.R., 2003, ApJ, 596, 240
Maness, H., Martins, F., Trippe, S., Genzel, R., Graham, J. R.,
Sheehy, C., Salaris, M., Gillessen, S., Alexander, T., Paumard,
T., and 3 coauthors, 2007, ApJ, 669, 1024
Monaco, P., Fontanot, F., Taffoni, G., 2007, MNRAS, 375, 1189
Mortier, A. M. J., Serjeant, S., Dunlop, J. S., Scott, S. E., Ade, P.,
Alexander, D., Almaini, O., Aretxaga, I., Baugh, C., Benson, A.
J., and 68 coauthors, 2005, MNRAS, 363, 563
Nagashima, M., Lacey, C.G., Baugh, C.M., Frenk, C.S., Cole, S.,
2005, MNRAS, 358, 1247
c© 0000 RAS, MNRAS 000, 000–000
http://arxiv.org/abs/astro-ph/0604534
http://arxiv.org/abs/0710.2882
Galaxy evolution in the IR 21
Nagashima, M., Okamoto, T., Lacey, C.G., Baugh, C.M., Frenk,
C.S., Cole, S., 2005, MNRAS, 363, L31
Papovich, C., Dole, H., Egami, E., Le Floc’h, E., Prez-Gonzlez, P.
G., Alonso-Herrero, A., Bai, L., Beichman, C. A., and 15 coau-
thors, 2004, ApJS, 154, 70
Parra, R., Conway, J.E., Diamond, P.J., Thrall, H., Lonsdale, C.J.,
Lonsdale, C.J., Smith, H.E., 2007, ApJ,659,314
Patris, J., Dennefeld1, M., 2, Lagache, G., Dole, H., 2003, A&A
412, 349
Paumard, T., Genzel, R., Martins, F., Nayakshin, S., Beloborodov,
A. M., Levin, Y., Trippe, S., Eisenhauer, F., Ott, T., Gillessen, S.,
and 4 coauthors, 2006, ApJ, 643, 1011
Pearson C., Rowan-Robinson M., 1996, MNRAS, 283, 174.
Perez-Gonzalez P.G., Rieke G.H., Egami E., Alonso-Herrero A.,
Dole H., Papovich C., Blaylock M., Jones J., and 6 co-authors,
2005, ApJ, 630, 82
Pilbratt, G., 2003, SPIE, 4850, 586
Puget, J.-L., Abergel, A., Bernard, J.-P., Boulanger, F., Burton, W.
B., Desert, F.-X., Hartmann, D., 1996, A&A, 308, L5
Rieke, G.H., Loken, K., Rieke, M.J., Tamblyn, P., 1993, ApJ, 412,
Rowan-Robinson, M., 2001, ApJ, 549, 745
Rush B., Malkan M. A., Spinoglio L., 1993, ApJS, 89, 1
Sanders, D.B., Mirabel, I.F., 1996, ARA&A, 34, 749
Sanders, D.B, Mazzarella, J.M., Kim, D.-C., Surace, J.A., Soifer,
B.T., AJ, 126, 1607
Saunders, W., Rowan-Robinson, M., Lawrence, A., Efstathiou,
G., Kaiser, N., Ellis, R.S., Frenk, C.S., 1990, MNRAS, 242, 318
Shupe D. L., Fang F., Hacking P. B., Huchra J. P., 1998, ApJ, 501,
Silva L., Granato G.L., Bressan A., Danese L., 1998, ApJ, 509,
Silva, L., De Zotti, G., Granato, G. L., Maiolino, R., Danese, L.,
2005, MNRAS, 357, 1295
Smail I., Ivison R.J., Blain A.W., 1997, ApJ, 490, L5.
Soifer B. T., Neugebauer G., 1991, AJ, 101, 354
Soifer, B. T., Houck, J.R., Neugebauer, G., 1987, ARA&A, 25,
Soifer, B. T., Sanders, D. B., Madore, B. F., Neugebauer, G.,
Danielson, G. E., Elias, J. H., Lonsdale, Carol J., Rice, W. L.,
1987, ApJ, 320, 238
Spergel D. N., et al. 2007, ApJS 170, 377
Springel, V., et al. , 2005, Nat, 435, 629.
Stolte, A., Brandner, W., Grebel, E.K., Lenzen, R., Lagrange, A.,
2005, ApJ, 628, L113
Takeuchi, T.T., Yoshikawa, K., Ishii, T.T., ApJ, 587, L89
Van Dokkum, P., 2007, ApJ (in press) (arXiv:0710.0875)
Vaisanen, P., Tollestrup, E.V., Fazio, G.G., 2001, MNRAS, 325,
Vega, O., Silva, L., Panuzzo, P., Bressan, A., Granato, G.L.,
Chavez, M., 2005, MNRAS, 364, 1286.
Werner, M.W., et al. , 2004, ApJS, 154, 1
Wright, G. S., Joseph, R. D., Meikle, W. P. S., 1984, Nature 309,
Xu, C., Hacking, P.B., Fang, F., Shupe, D.L., Lonsdale, C.J., Lu,
N.Y., Helou, G., Stacey, G.J., Ashby, M.L.N., 1998, ApJ, 508,
APPENDIX A: REDSHIFT DISTRIBUTIONS
In this Appendix, we present some predictions from our standard
model for the redshift distributions of galaxies selected at different
fluxes in the Spitzer bands. This is principally for completeness,
to assist in interpreting data from current surveys, and to assist in
planning future surveys based on Spitzer data. The set of redshift
distributions at all observed fluxes in principle contains equivalent
information to that in the luminosity functions at different wave-
lengths and redshifts. However, comparing models with observa-
tions via luminosity functions is more physically transparent than
making the comparison via redshift distributions, which is why we
have presented our results on luminosity functions in the main part
of the paper, and why we make only a limited direct comparison
with observed redshift distributions in this Appendix. In addition,
if one only compares the predicted and observed redshift distribu-
tions for galaxies above a single flux limit (e.g. the flux limit of
a survey), this has less information than comparing the luminosity
functions at different redshifts.
We first show in Fig. A1 how the median redshift, and the 10-
90 percentile range, are predicted to change with flux for galax-
ies selected in one of the four Spitzer bands 3.6, 8.0, 24 or 70
µm. While at most wavelengths the median redshift is predicted
to increase smoothly and monotonically with decreasing flux, this
is not true at 24 µm, where there is a bump around Sν ∼ 100µJy.
The structure seen for the 24 µm band as compared to the other
wavelengths results from different PAH emission features moving
through the band with increasing redshift.
In Fig. A2, we show the predictions from our standard model
for the redshift distributions of galaxies in the four IRAC bands.
For each band, we show the redshift distribution for galaxies se-
lected to be brighter than Sν > 10µJy in that band. The flux limit
Sν > 10µJy has been chosen to match that in the observed deep
sample selected at 3.6µm by Franceschini et al. (2006). In each
panel, the blue curve shows the predicted dN/dz for all galax-
ies, normalized to unit area under the curve, and the red and green
curves show the separate contributions of bursting and quiescent
galaxies to the total. For 3.6µm, the black line shows the ob-
served redshift distribution from Franceschini et al. (2006), which
has also been normalized to unit area under the curve. We see that
the observed redshift distribution peaks at a slightly higher redshift
than in the model. However, the luminosity function evolution de-
rived from this same sample is in reasonable agreement with the
model, as was already shown in Fig. 4. Franceschini et al. (2006)
note that the peak seen in their data at z ∼ 0.8 is partly contributed
by large-scale structures in the CDFS field.
In Fig. A3, we show predicted redshift distributions for galax-
ies selected to be at a set of different fluxes in the four IRAC bands.
The curves for the different fluxes are all normalized to have unit
area as before, but in this figure the galaxies are selected to be at a
particular flux, rather than being brighter than a certain flux. As one
would expect, the typical redshift increases as the flux decreases.
Figs. A4 and A5 show for the three MIPS bands the equiv-
alent of Figs. A2 and A3 for the IRAC bands. In Fig. A4, we
show the predicted redshift distributions for galaxies brighter than
a particular flux, where this flux limit is taken to be 83 µJy at 24
µm, 10 mJy at 70 µm and 100 mJy at 160 µm. The flux limit
at 24 µm has been chosen to match that used in the deep obser-
vational samples of Le Floc’h et al. (2005), Perez-Gonzalez et al.
(2005) and Caputi et al. (2006a), while the flux limits at 70 µm
and 160 µm have been chosen to be roughly 3 times brighter than
the source confusion limits in these bands. We see in Fig. A5 that
c© 0000 RAS, MNRAS 000, 000–000
http://arxiv.org/abs/0710.0875
22 Lacey et al.
Figure A1. Model predictions for the median redshift as a function of flux in four Spitzer bands. (a) 3.6 µm, 8.0 µm, (c) 24 µm, (d) 70 µm. In each panel,
the median redshift for galaxies at each flux is shown by a solid line, and the 10- and 90-percentiles are shown by dashed lines.
the redshift distributions at 24 µm show much more structure than
at other wavelengths. This results from different PAH emission fea-
tures moving through the 24 µm band with changing redshift.
In Fig. A4(a), we compare the predicted redshift dis-
tribution at 24 µm with observational determinations from
Perez-Gonzalez et al. (2005) (dashed black line) and Caputi et al.
(2006a) (solid black line). The observed distributions have been
separately normalized to unit area under the curve, as for the
model distribution. Both observed distributions are based pri-
marily on photometric redshifts, but the photometric redshifts of
Caputi et al. (2006a) are likely to be more accurate than those of
Perez-Gonzalez et al. (2005), since the former are based on deeper
optical and K-band data than the latter. (Perez-Gonzalez et al.
found optical counterparts with BAB . 24.7 or RAB . 23.7 for
∼ 70% of their Sν(24µm) > 83µJy sources, but relied on IRAC
fluxes in deriving photo-z’s for the remaining ∼ 30% of their sam-
ple. On the other hand, Caputi et al. found K-band counterparts
with K(V ega) < 21.5 for 95% of their Sν(24µm) > 80µJy
sample, and derived photo-z’s for essentially all of these sources
using optical and K-band data alone). Both observed distributions
are similar, but the Caputi et al. distribution shows more structure.
This is a combination of the effects of more accurate photometric
redshifts but also a 9 times smaller survey area, which means that
fluctuations due to galaxy clustering are larger. Caputi et al. argue
that the separate peaks at z ∼ 0.7 and 1.1 result from large-scale
structure, but that the bump at z ∼ 1.9 results from PAH emis-
sion features entering the observed 24 µm band. We see that the
model also predicts peaks in the redshift distribution at z ∼ 0.3,
z ∼ 1 and z ∼ 2, which can be explained by different PAH fea-
tures moving through the 24 µm band, although the z ∼ 2 peak is
more prominent than is seen in the observational data. Overall, the
model redshift distribution at this flux limit is too skewed to high
redshift compared to the observations, predicting too few galaxies
at z ∼ 0.5− 1, and too many in the peak at z ∼ 2.
We investigate further this apparent discrepancy in the 24 µm
redshift distribution in Fig. A6, where we show the effects of appar-
ent magnitude limits in the R and K-bands on the predicted redshift
distributions for Sν(24µm) > 83µJy. In this plot, the redshift
distributions are plotted as number per solid angle, without nor-
malizing to unit area under the curve. The left and right panels re-
spectively have the redshift distributions of Perez-Gonzalez et al.
and Caputi et al. overplotted. We concentrate on the comparison
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 23
Figure A2. Predicted galaxy redshift distributions in the four IRAC bands, for galaxies brighter than Sν = 10µJy. (a) 3.6 µm, (b) 4.5 µm, (c) 5.8 µm, and
(d) 8.0 µm. The model curves (which all include the effects of dust) are as follows - blue: total; red: ongoing bursts; green: quiescent galaxies. The curves are
normalized to unit area under the curve for the total counts. The median (z50) and 10- and 90-percentile (z10, z90) redshifts for the total counts in each band
are also given in each panel. For 3.6 µm, the model predictions are compared with observational data from Franceschini et al. (2006) (black dashed line),
normalized to unit area as for the models. The error bars plotted on the observational data include Poisson errors only.
with Caputi et al. , since this has the simpler sample selection and
more accurate redshifts. The model prediction for K < 21.5
(which is the magnitude limit used by Caputi et al. ) is shown
by the short-dashed blue line, while the prediction with no limit
on the K-magnitude is shown by the solid blue line. The model
dN/dz with no limit on the K magnitude is most discrepant with
the Caputi et al. data at z ∼ 2, where it predicts ∼ 2 times too
many galaxies. This is directly related to the fact that the pre-
dicted luminosity function at z = 2 at rest-frame wavelength
8 µm (corresponding to observed wavelength 24 µm) and lumi-
nosity ∼ 1011L⊙ is also ∼ 2 times too high compared to what
Caputi et al. estimate from their data, as shown in Fig. 5(c). When
the effect of the K < 21.5 limit is included, the predicted redshift
distribution is closer to the observational data, but only 58% of the
model galaxies are brighter than this K-band magnitude limit, as
against 95% in the observed sample of Caputi et al. . We conclude
that the main reason for the discrepancy between the predicted and
observed redshift distributions at 24 µm is that the model predicts
a rest-frame 8 µm luminosity function at z ∼ 2 which is somewhat
too high at luminosities ∼ 1011L⊙, even though it reproduces quite
well the general features of the evolution of the mid-IR luminosity
function.
c© 0000 RAS, MNRAS 000, 000–000
24 Lacey et al.
Figure A3. Predicted galaxy redshift distributions in the four IRAC bands, for different fluxes. (a) 3.6 µm, (b) 4.5 µm, (c) 5.8 µm, and (d) 8.0 µm. In this
figure, the redshift distributions are for galaxies at a particular flux. Predictions are shown for fluxes Sν = 0.1, 1, 10, 100 and 1000 µJy, as shown in the key.
In all cases, the model curves are normalized to unit area, and include the effects of dust.
c© 0000 RAS, MNRAS 000, 000–000
Galaxy evolution in the IR 25
Figure A4. Predicted galaxy redshift distributions in the three MIPS bands,
for galaxies brighter than a specified flux. (a) 24 µm, Sν > 83µJy, (b) 70
µm, Sν > 10mJy, and (c) 160 µm, Sν > 100mJy. The model curves
are as follows - blue: total; red: ongoing bursts; green: quiescent galax-
ies. The curves are normalized to unit area under the curve for the total
counts. The median (z50) and 10- and 90-percentile (z10, z90) redshifts for
the total counts in each band are also given in each panel. For 24 µm, the
model predictions are compared with observational data from Caputi et al.
(2006a) (solid black line) and Perez-Gonzalez et al. (2005) (dashed black
line), normalized to unit area as for the models. The error bars plotted on
the observational data include Poisson errors only for Caputi et al. , but also
include errors in photometric redshifts for Perez-Gonzalez et al.
Figure A5. Predicted galaxy redshift distributions in the three MIPS bands,
for different fluxes. (a) 24 µm, (b) 70 µm, and (c) 160 µm. In this figure,
the redshift distributions are for galaxies at a particular flux, as shown in the
key in each panel. In all cases, the model curves are normalized to unit area,
and include the effects of dust.
c© 0000 RAS, MNRAS 000, 000–000
26 Lacey et al.
Figure A6. Predicted redshift distributions at 24µm, showing the effects of optical or near-IR magnitude limits. Model galaxies are selected with Sν > 83µJy
together with the optical/NIR magnitude limits as shown in the key. The fraction of 24 µm sources brighter than each magnitude limit is also given. (a) R-band
magnitude limit. The observed redshift distribution from Perez-Gonzalez et al. (2005) is overplotted in black. Note Le Floc’h et al. (2005) used R < 24
and obtained 54% completeness. (b) K-band magnitude limit. The observed redshift distribution from Caputi et al. (2006a) (with K < 21.5) is overplotted.
Magnitudes are on the Vega system.
c© 0000 RAS, MNRAS 000, 000–000
	Introduction
	Model
	GALFORM galaxy formation model
	GRASIL model for stellar and dust emission
	Choice of parameters in the GALFORM+GRASIL model
	Number counts
	Evolution of the galaxy luminosity function
	Evolution of the galaxy luminosity function at 3-8 m
	Evolution of the galaxy luminosity function at 12-24 m
	Evolution of the galaxy luminosity function at 70-160 m
	Evolution of the total mid+far-IR luminosity function
	Inferring stellar masses and star formation rates from Spitzer data
	Conclusions
	Redshift distributions

  We present predictions for the evolution of the galaxy luminosity function,
number counts and redshift distributions in the IR based on the Lambda-CDM
cosmological model. We use the combined GALFORM semi-analytical galaxy
formation model and GRASIL spectrophotometric code to compute galaxy SEDs
including the reprocessing of radiation by dust. The model, which is the same
as that in Baugh et al (2005), assumes two different IMFs: a normal solar
neighbourhood IMF for quiescent star formation in disks, and a very top-heavy
IMF in starbursts triggered by galaxy mergers. We have shown previously that
the top-heavy IMF seems to be necessary to explain the number counts of faint
sub-mm galaxies. We compare the model with observational data from the SPITZER
Space Telescope, with the model parameters fixed at values chosen before
SPITZER data became available. We find that the model matches the observed
evolution in the IR remarkably well over the whole range of wavelengths probed
by SPITZER. In particular, the SPITZER data show that there is strong evolution
in the mid-IR galaxy luminosity function over the redshift range z ~ 0-2, and
this is reproduced by our model without requiring any adjustment of parameters.
On the other hand, a model with a normal IMF in starbursts predicts far too
little evolution in the mid-IR luminosity function, and is therefore excluded.

Introduction
	Model
	GALFORM galaxy formation model
	GRASIL model for stellar and dust emission
	Choice of parameters in the GALFORM+GRASIL model
	Number counts
	Evolution of the galaxy luminosity function
	Evolution of the galaxy luminosity function at 3-8 m
	Evolution of the galaxy luminosity function at 12-24 m
	Evolution of the galaxy luminosity function at 70-160 m
	Evolution of the total mid+far-IR luminosity function
	Inferring stellar masses and star formation rates from Spitzer data
	Conclusions
	Redshift distributions
