The many lives of AGN: cooling flows, black holes and the luminosities and colours of galaxies

21 Pages·2005·5.81 MB·English

Mon.Not.R.Astron.Soc.000,000–000 (2005) Printed15December 2005 (MNLATEXstylefilev2.2)
The many lives of AGN: cooling flows, black holes and the
luminosities and colours of galaxies
Darren J. Croton1, Volker Springel1. Simon D. M. White1, G. De Lucia1, C. S.
Frenk2, L. Gao1, A. Jenkins2, G. Kauffmann1, J. F. Navarro3, N. Yoshida4
1Max-Planck-Institut fu¨r Astrophysik, D-85740 Garching, Germany.
2Institute for Computational Cosmology, Physics Department, Durham. U.K.
3Department of Physics and Astronomy, University of Victoria, B.C., Canada.
4Department of Physics, Nagoya University,Chikusa-ku, Nagoya 464-8602, Japan
Accepted —.Received—;inoriginalform—
ABSTRACT
We simulate the growth of galaxies and their central supermassive black holes by
implementing a suite of semi-analytic models on the output of the Millennium Run,
a very large simulation of the concordance ΛCDM cosmogony. Our procedures fol-
low the detailed assembly history of each object and are able to track the evolution
of all galaxies more massive than the Small Magellanic Cloud throughout a volume
comparable to that of large modern redshift surveys. In this first paper we supple-
mentprevioustreatmentsofthegrowthandactivityofcentralblackholeswithanew
modelfor‘radio’feedbackfromthoseAGNthatlieatthecentreofaquasistaticX-ray
emitting atmosphere in a galaxy groupor cluster. We show that for energeticallyand
observationallyplausible parameterssucha model cansimultaneously explain:(i) the
low observed mass drop-out rate in cooling flows; (ii) the exponential cut-off at the
bright end of the galaxy luminosity function; and (iii) the fact that the most massive
galaxiestendtobe bulge-dominatedsystemsinclustersandto containsystematically
older stars than lower mass galaxies. This success occurs because static hot atmo-
spheresform only inthe mostmassive structures,and radiofeedback (in contrast,for
example, to supernova or starburst feedback) can suppress further cooling and star
formation without itself requiring star formation. We discuss possible physical mod-
els which might explain the accretion rate scalings requiredfor our phenomenological
‘radio mode’ model to be successful.
Key words: cosmology: theory, galaxies: formation, galaxies: evolution, galaxies:
active, galaxies: cooling flows, black hole physics
1 INTRODUCTION structure (Peacock et al. 2001), the baryon fraction in rich
clusters(Whiteet al.1993) and thetheoryofBig Bangnu-
cleosynthesis (Olive et al. 2000). A working model for the
Theremarkableagreementbetweenrecentmeasurementsof
growth of all structurethusappears well established.
cosmic structure over a wide range of length- and time-
scales has established a standard paradigm for structure Inthiscosmogony,galaxiesformwhengascondensesat
formation, the ΛCDM cosmogony. This model can simulta- the centres of a hierarchically merging population of dark
neously match the microwave background fluctuations seen haloes, as originally proposed by White & Rees (1978). At-
atz 1000(e.g.Spergeletal.2003),thepowerspectrumof temptstounderstandthisprocessindetailhaveconsistently
∼
thelowredshiftgalaxydistribution(e.g.Percivaletal.2002; run into problems stemming from a mismatch in shape be-
Tegmarketal.2004),thenonlinearmassdistributionatlow tweenthepredicteddistributionofdarkhalomassesandthe
redshift as characterised by cosmic shear (e.g. Van Waer- observed distribution of galaxy luminosities. Most stars are
beke et al. 2002) and the structure seen in the z=3 Ly α in galaxies of Milky Way brightness; the galaxy abundance
forest (e.g.Mandelbaum etal. 2003). Italso reproducesthe declinesexponentiallyatbrighterluminositiesandincreases
present acceleration of the cosmic expansion inferred from sufficiently slowly at fainter luminosities that relatively few
supernova observations (Perlmutter et al. 1999; Riess et al. stars are in dwarfs. In contrast, the theory predicts a much
1998), and it is consistent with the mass budget inferred broader halo mass distribution – a constant mass-to-light
for the present universe from the dynamics of large-scale ratio would produce more high and low luminosity galax-
(cid:13)c 2005RAS
2 Croton et al.
ies than are observed while underpredicting the number of present below, radio galaxies prevent significant accretion,
galaxies like the Milky Way. Attempts to solve these prob- thuslimitingthemassofthecentralgalaxiesandpreventing
lemsinitiallyinvokedcoolinginefficienciestoreducegascon- them from forming stars at late times when their mass and
densationinmassivesystems,andsupernovafeedbacktore- morphology can still change through mergers. The result is
ducestarformationefficiencyinlowmasssystems(White& a galaxy luminosity function with a sharper high-mass cut-
Rees 1978; White& Frenk 1991). Formation of dwarfs may off in which the most massive systems are red, dead and
also be suppressed by photoionisation heating (Efstathiou elliptical.
1992).AsThoul&Weinberg(1995)emphasised,coolingef- To make quantitative predictions for the galaxy popu-
fectsalonearetooweaktoproducethebrightendcut-offof lation in aΛCDM universeit is necessary tocarry out sim-
theluminosityfunction,andrecentattemptstofitobserved ulations. Present numerical capabilities allow reliable sim-
luminosity functions havebeen forced to include additional ulation of the coupled nonlinear evolution of dark matter
feedback processes in massive systems (e.g. Benson et al. and diffuse gas, at least on the scales which determine the
2003). In this paper we argue that radio sources may pro- global properties of galaxies. Once gas cools and condenses
videtherequiredfeedbackwhileatthesametimeproviding intohalocores,however,bothitsstructureandtheratesat
a solution to two other long-standing puzzles. whichitturnsintostarsandfeedsblackholesaredetermined
An important unanswered question is why the gas at bysmall-scale ‘interstellar medium’processeswhicharenot
the centre of most galaxy clusters is apparently not con- resolved. These are usually treated through semi-analytic
densing and turning into stars when the observed X-ray recipes,parameterisedformulaewhichencapsulate‘subgrid’
emission implies a cooling time which is much shorter than physicsintermsofstarformationthresholds,Schmidt‘laws’
the age of the system. This cooling flow puzzle was noted for star formation, Bondi models for black hole feeding,
assoon asthefirst X-raymapsofclustersbecameavailable etc. The form and the parameters of these recipes are cho-
(Cowie & Binney 1977; Fabian & Nulsen 1977) and it was sen to reproduce the observed systematics of star forma-
made more acute when X-ray spectroscopy demonstrated tion and AGN activity in galaxies (e.g. Kennicutt 1998).
that very little gas is cooling through temperatures even a Withawell-constructedschemeitispossibletoproducesta-
factor of three below that of the bulk of the gas (Peterson bleandnumericallyconvergedsimulationswhichmimicreal
etal. 2001; Tamuraetal.2001). Acluetothesolution may star-forminggalaxiesremarkablywell(Springel&Hernquist
come from the observation (Burns et al. 1981) that every 2003a). Instrongly star-forming galaxies, massivestars and
cluster with a strong cooling flow also contains a massive supernovae produce winds which redistribute energy, mass
andactivecentralradiogalaxy.Tabor&Binney(1993)sug- andheavyelementsoverlargeregions(Heckmanetal.1990;
gested that radio galaxies might regulate cooling flows, and Martin 1999). Even stronger feedback is possible, in princi-
thisideahasgainedconsiderablerecentsupportfromX-ray ple,fromAGN(Begelmanetal.1991).Inbothcasesthede-
mapswhichshowdirectevidenceforaninteractionbetween termining processes occur on very small scales and so have
radiolobesandtheintraclustergas(Fabianetal.2003;Mc- to be included in simulations through parametrised semi-
Namara et al. 2000, 2005). A number of authors have sug- analytic models. Unfortunately, the properties of simulated
gestedwaysinwhichtheradiosourcemightreplacethether- galaxies turn out to depend strongly on how these unre-
mal energy lost to X-ray emission (Binney & Tabor 1995; solved star-formation and feedback processes are treated.
Churazovetal.2002;Bru¨ggen&Kaiser2002;Ruszkowski& Since the diffuse gas distribution and its cooling onto
Begelman 2002; Kaiser & Binney 2003; Omma et al. 2004). galaxiesarestronglyaffectedbythedescriptionadoptedfor
We do not focus on this aspect of the problem here, but the subgrid physics, every modification of a semi-analytic
rather demonstrate that if the scaling properties of radio model (or of its parameters) requires a simulation to bere-
sourcefeedbackaresetsotheycanplausiblysolvethecool- peated.Thismakesparameterstudiesortestsof,say,theef-
ing flow problem they induce a cut-off at the bright end of fectsofdifferentAGNfeedbackmodelsintoaveryexpensive
the galaxy luminosity function which agrees well with ob- computing exercise. A cost-effective alternative is to repre-
servation. sentthebehaviourofthediffusegasalsobyasemi-analytic
Anotherpuzzlingaspectofthegalaxypopulationisthe recipe. Since the dark matter couples to the baryons only
fact that the most massive galaxies, typically ellipticals in through gravity, its distribution on scales of galaxy haloes
clusters, are made of the oldest stars and so finished their andaboveisonlyweaklyaffectedbythedetailsofgalaxyfor-
star formation earlier than lower mass galaxies (Bender & mation. Its evolution can therefore be simulated once, and
Saglia1999).Confirmingevidenceforthiscomesfromlook- the evolution of the baryonic component can be included
backstudieswhichshowthatbothstar-formationandAGN in post-processing by applying semi-analytic models to the
activity take place more vigorously and in higher mass ob- stored histories of all dark matter objects. Since the sec-
jectsatredshiftsof1to2thaninthepresentUniverse(e.g. ond step is computationally cheap, available resources can
Shaver et al. 1996; Madau et al. 1996). Cowie et al. (1996) be used to carry out the best possible dark matter simula-
termed this phenomenon ‘down-sizing’, and prima facie it tion, and then many parameter studies or tests of alterna-
conflicts with hierarchical growth of structure in a ΛCDM tivemodelscanbecarriedoutinpost-processing.Thissim-
cosmogonywheremassivedarkhaloesassembleatlowerred- ulation approach was first introduced by Kauffmann et al.
shift thanlowermasshaloes (e.g.Lacey&Cole 1993). This (1999) and it is the approach we apply here to the Millen-
puzzle is related to the previous two; the late-forming high nium Run, the largest calculation to date of the evolution
mass haloes in ΛCDMcorrespond togroups and clusters of ofstructureintheconcordanceΛCDMcosmogony(Springel
galaxies,andsimpletheoriespredictthatatlatetimestheir et al. 2005b).
centralgalaxiesshouldgrowtomasseslargerthanthoseob- This paper is organised as follows. In Section 2 we de-
servedthroughaccretionfromcoolingflows.Inthemodelwe scribe the Millennium Run and thepost-processing we car-
(cid:13)c 2005RAS,MNRAS000,000–000
Cooling flows, black holes and the luminosities and colours of galaxies 3
ried out to construct merging history trees for all the dark The initial conditions at z=127 were created by dis-
haloeswithinit.Section3presentsthemodelfortheforma- placing particles from a homogeneous, ‘glass-like’ distribu-
tion and evolution of galaxies and their central supermas- tion (White 1996) using a Gaussian random field with a
sive black holes that we implement on these merging trees. ΛCDM linear power spectrum as given by the Boltzmann
Section 4 describes the main results of our modelling, con- code CMBFAST (Seljak & Zaldarriaga 1996). The simula-
centrating on the influenceof ‘radio mode’ feedback on the tion was thenevolved tothepresent epoch usingaleapfrog
properties of the massive galaxy population. In Section 5 integrationschemewithindividualandadaptivetimesteps,
we discuss physical models for black hole accretion which with up to11000 time steps for individual particles.
may explain the phenomenology required for our model to Thesimulationitselfwascarriedoutwithaspecialver-
besuccessful.Finally,Section6summarisesourconclusions sion of theGADGET-2code (Springelet al. 2001b; Springel
and suggests some possible directions for future investiga- 2005) optimised for very low memory consumption so that
tion. it would fit into the nearly 1 TB of physically distributed
memoryavailableontheparallelIBMp690computer1 used
for the calculation. The computational algorithm used the
‘TreePM’method(Xu1995;Bodeetal.2000;Bagla2002)to
2 THE DARK MATTER SKELETON: THE
evaluategravitationalforces,combiningahierarchicalmulti-
MILLENNIUM RUN
poleexpansion,or‘tree’algorithm(Barnes&Hut1986),and
Our model for the formation and evolution of galaxies and a classical, Fourier transform particle-mesh method (Hock-
their central supermassive black holes is implemented on ney & Eastwood 1981). An explicit force-split in Fourier-
top of the Millennium Run, a very large dark matter simu- space produces a very nearly isotropic force law with neg-
lation of the concordance ΛCDM cosmology with 21603 ligible force errors at the force matching scale. The short-
1.0078 1010 particles in a periodic box of 500h−1Mpc o'n range gravitational force law is softened on comoving scale
×
a side. A full description is given in Springel et al. (2005b); 5h−1kpc (Plummer-equivalent) which may be taken as the
here we summarise the main simulation characteristics and spatial resolution limit of the calculation, thus achieving a
the definition and construction of the dark matter merging dynamic range of 105 in 3D. The calculation, performed in
historytreesweuseinourgalaxyformationmodelling.The parallelon512processors,requiredslightlylessthan350000
dark matter distribution is illustrated in the top panel of processor hoursofCPU time,or28daysof wall-clock time.
Fig. 1 for a 330 280 15h−1Mpc slice cut from the full
× ×
volume. The projection is colour coded by density and lo-
cal velocity dispersion, and illustrates the richness of dark 2.2 Haloes, substructure, and merger tree
matter structure for comparison with structure in the light construction
distribution to which we will come later. Dark matter plots
Ourprimary application of theMillennium Run in this pa-
on a wider range of scales may be found in Springel et al.
perusesfinelyresolvedhierarchicalmergingtreeswhichen-
(2005b).
code the full formation history of tens of millions of haloes
andthesubhaloesthatsurvivewithinthem.Thesemerging
2.1 Simulation characteristics history trees are the backbone of the model that we imple-
mentinpost-processinginordertosimulatethewiderange
We adopt cosmological parameter values consistent with a of baryonic processes that are important during the forma-
combined analysis of the 2dFGRS (Colless et al. 2001) and tionandevolutionofgalaxiesandtheircentralsupermassive
first-year WMAP data (Spergel et al. 2003; Seljak et al. black holes.
2005). They are Ωm = Ωdm+Ωb = 0.25, Ωb = 0.045, h = Westorethefullparticledatabetweenz=20andz=0
0.73, ΩΛ = 0.75, n = 1, and σ8 = 0.9. Here Ωm denotes at 60 output times spaced in expansion factor according to
the total matter density in units of the critical density for theformula
closure, ρcrit = 3H02/(8πG). Similarly, Ωb and ΩΛ denote
thedensitiesofbaryonsanddarkenergyatthepresentday. log(1+zn)=n(n+35)/4200 . (1)
The Hubble constant is given as H0 =100hkms−1Mpc−1, This spacing is ‘locally’ logarithmic but becomes smoothly
while σ8 is the rms linear mass fluctuation within a sphere finer at lower redshift, with a temporal resolution by red-
of radius 8h−1Mpc extrapolated to z=0.
shiftzeroofapproximately300Myr.Additionaloutputsare
The chosen simulation volume is a periodic box of
added at z=30,50,80,127 to produce a total of 64 snap-
size 500h−1Mpc, which implies a particle mass of 8.6
shots in all. We note that each snapshot has a total size in
108h−1M(cid:12). This volume is large enough to include inte×r- excessof300GB,givingarawdatavolumeofnearly20TB.
esting objects of low space density, such as quasars or rich
Together with each particle coordinate dump, the sim-
galaxy clusters, the largest of which contain about 3 mil-
ulation code directly produces a friends-of-friends (FOF)
lion simulation particles at z = 0. At the same time, the
group catalogue on the fly and in parallel. FOF groups are
mass resolution is sufficient that haloes that host galaxies
definedasequivalenceclasseswhereanypairoftwoparticles
as faint as 0.1L are typically resolved with at least 100
? is placed into the same group if their mutual separation is
∼
particles.Notethatalthoughdiscretenessnoisesignificantly
less than 0.2 of the mean particle separation (Davis et al.
affects the merger histories of such low mass objects, the
1985). This criterion combines particles into groups with
galaxies that reside in halos with <100 particles are usu-
ally sufficiently far down the lumin∼osity function that any
uncertainty in their properties has little impact on our re- 1 This computer is operated by the Computing Centre of the
sults or conclusions. Max-PlanckSocietyinGarching,Germany.
(cid:13)c 2005RAS,MNRAS000,000–000
4 Croton et al.
Figure 1.Theredshiftzerodistributionofdarkmatter(top)andofgalaxylight(bottom) forasliceofthickness15h−1Mpc,cutfrom
the Millennium Run. For the dark matter distribution, intensity encodes surface density and colour encodes local velocity dispersion.
For the light distribution, intensity encodes surface brightness and colour encodes mean B−V colour. The linear scale is shown by the
barinthetoppanel.
(cid:13)c 2005RAS,MNRAS000,000–000
Cooling flows, black holes and the luminosities and colours of galaxies 5
a mean overdensity of about 200, corresponding approxi- archicallyinCDMuniversestherecanbeseveralprogenitors
mately to that expected for a virialised group. The group for any given halo, but in general there is only one descen-
catalogue saved to disk for each output only kept groups dant. (Typically the cores of virialised dark matter struc-
with at least 20 particles. tures do not split into two or more objects.) We therefore
High-resolutionsimulationslikethepresentoneexhibit construct merger trees based on defining a unique descen-
a rich substructure of gravitationally bound dark matter dantforeachhaloandsubhalo.Thisis,infact,sufficientto
subhaloesorbitingwithinlargervirialisedhalos(e.g.Ghigna define the entire merger tree, since the progenitor informa-
et al. 1998). The FOF group-finder built into our simula- tion then follows implicitly. Furtherdetails can be found in
tion code identifies thehaloes but not their subhalos. Since Springel et al. (2005b).
we wish to follow the fate of infalling galaxies and their We store the resulting merging histories tree by tree.
halos, and these are typically identifiable for a substantial Since each tree contains the full formation history of some
time as a dark matter subhalo within a FOF halo, we ap- z=0 halo, the physical model for galaxy formation can be
ply in post-processing an improved and extended version computed sequentially tree by tree. It is thus unnecessary
of the SUBFIND algorithm of Springel et al. (2001a). This to load all the history information on the simulation into
computesanadaptivelysmoothed darkmatterdensityfield computer memory at the same time. Actually, this would
withineachhalousingakernel-interpolationtechnique,and be currently impossible, since the combined trees contain a
then exploits the topological connectivity of excursion sets total of around 800 million subhaloes for each of which a
above a density threshold to identify substructure candi- numberof attributes are stored. Note that,although evolv-
dates. Each substructure candidate is subjected to a gravi- ingthegalaxypopulation sequentiallyinthiswayallows us
tational unbinding procedure. If the remaining bound part to consider, in principle, all interactions between galaxies
has more than 20 particles, the subhalo is kept for further that end up in the same present-day FOF halo, it does not
analysis and some ofits basic physicalpropertiesare deter- allowustomodellongerrangeinteractionswhichmighttake
mined(angularmomentum,maximumofitsrotationcurve, place between galaxies that end up in different FOF halos
velocity dispersion, etc.). After all subhaloes are identified at z=0.
they are extracted from the FOF halo so that the remain-
ing featureless ‘background’ halo can also be subjected to
theunbindingprocedure. This technique,however,neglects
thefact that substructuresembedded within a halo help to 3 BUILDING GALAXIES: THE
bind its material, and thus removing them can, in princi- SEMI-ANALYTIC MODEL
pal, unbind some of the FOF halo that may otherwise still
3.1 Overview
be loosely bound. We accept this small effect for techni-
cal reasons related to the robustness of our halo definition In the following sub-sections we describe the baryonic
procedures and the need for unambigous particle-subhalo physics of one particular model for the formation and evo-
assignments in our data structures. We note that the to- lution of galaxies and of their central supermassive black
tal mass in substructures is typically below 10% and often holes. A major advantage of oursimulation strategy isthat
substantially smaller. Thus any bias in the bound mass of the effects of parameter variations within this model (or
the parent halo due to additional unbinding is always very indeed alternative assumptions for some of the processes)
small. can be explored at relatively little computational expense
TocomputeavirialmassestimateforeachFOFhalowe since the model operates on the stored database of merger
use the spherical-overdensity approach, where the centre is trees; the simulation itself and the earlier stages of post-
determinedusingtheminimumofthegravitationalpotential processing do not need to be repeated. We have, in fact,
within the group and we define the boundary at the radius explored a wide model and parameter space to identify our
which encloses a mean overdensity of 200 times the critical current best model. We summarise the main parameters of
value.Thevirialmass,radiusandcircularvelocityofahalo thismodel, their‘best’ values,and theirplausible ranges in
at redshift z are then related by Table 1. These choices produce a galaxy population which
100 V3 matches quiteclosely many observed quantities. In this pa-
Mvir = G H2(z)Rv3ir = 10GvHir(z) (2) per we discuss the field galaxy luminosity-colour distribu-
tion;themean stellar mass–stellar agerelation; theTully-
where H(z) is theHubbleconstant at redshift z. Fisherrelation,coldgasfractionsandgas-phasemetallicities
Atz=0ourproceduresidentify17.7 106 FOFgroups, ofSb/cspirals;thecolour–magnituderelationofellipticals;
downfromamaximumof19.8 106 atz×=1.4whengroups thebulge mass – black hole mass relation; and the volume-
×
weremoreabundantbutoflowermassonaverage.Thez=0 averagedcosmicstar-formationandblackholeaccretionhis-
groups contain a total of 18.2 106 subhaloes, with the tories.InSpringeletal.(2005b)wealsopresentedresultsfor
×
largestFOFgroupcontaining2328ofthem.(Notethatwith galaxycorrelationsasafunctionofabsolutemagnitudeand
ourdefinitions,allFOFgroupscontainatleastonesubhalo, colour, for the baryonic ‘wiggles’ in the large-scale power
the main subhalo which is left over after removal of any spectrumofgalaxies,andfortheabundance,originandfate
substructure and the unbound component. We require this of high redshift supermassive black holes which might cor-
main subhalo to contain at least 20 particles.) respond to the z 6 quasars discovered by the SDSS (Fan
∼
Having found all haloes and subhaloes at all output et al. 2001)
snapshots, we then characterise structural evolution by In our model we aim to motivate each aspect of the
building merging trees that describe in detail how haloes physics of galaxy formation using the best available obser-
growastheuniverseevolves.Becausestructuresmergehier- vations and simulations. Our parameters have been chosen
(cid:13)c 2005RAS,MNRAS000,000–000
6 Croton et al.
Table1.Asummaryofourmainmodelparametersandtheirbestvaluesandplausibleranges,asdescribedinthetext.Onceset,these
valuesarekeptfixedforallresultspresentedinthispaper,inparticularformodelsinwhichAGNfeedbackisswitchedoff.
parameter description bestvalue plausiblerange
fb cosmicbaryonfraction(§3.3) 0.17 fixed
z0,zr redshiftofreionization(§3.3) 8,7 fixed
fBH mergercoldgasBHaccretionfraction(§3.4.1) 0.03 0.02−0.04
κAGN quiescent hotgasBHaccretionrate(M(cid:12)yr−1)(§3.4.2) 6×10−6 (4−8)×10−6
αSF starformationefficiency(§3.5) 0.07 0.05−0.15
(cid:15)disk SNfeedbackdiskreheatingefficiency(§3.6) 3.5 1−5
(cid:15)halo SNfeedbackhaloejectionefficiency(§3.6) 0.35 0.1−0.5
γej ejected gasreincorporationefficiency(§3.6) 0.5 0.1−1.0
Tmerger majormergermassratiothreshold(§3.7) 0.3 0.2−0.4
R instantaneous recycledfractionofSFtothecolddisk(§3.9) 0.3 0.2−0.4
Y yieldofmetalsproducedperunitSF(§3.9) 0.03 0.02−0.04
to reproduce local galaxy properties and are stable in the ature and is added to a quasi-static hot atmosphere that
sensethatnoneofourresultsiscritically dependentonany extendsapproximately tothevirial radiusofthedarkhalo.
single parameter choice; plausible changes in one parame- Gasfrom thecentralregion ofthisatmospheremayaccrete
teror recipecan usually beaccommodated through adjust- ontoacentralobject through acooling flow. Atearly times
mentoftheremainingparameterswithintheirownplausible and in lower mass systems the infalling gas still shocks to
range. The particular model we present is thus not unique. thevirialtemperaturebutitspost-shockcoolingtimeissuf-
Importantly, our model for radio galaxy heating in cooling ficiently short that a quasi-static atmosphere cannot form.
flows, which is the main focus of this paper, is only weakly Rather the shock occurs at much smaller radius and the
influencedbytheremaininggalaxyformationandblackhole shocked gas cools rapidly and settles onto a central object,
growth physics. This is because our radio mode feedback is which we assume to be a cold gas disk. This may in turn
active only in massive objects and at late times and it has besubjecttogravitationalinstability,leadingtoepisodesof
noeffectduringtheprincipalgrowthphaseofmostgalaxies star formation.
and AGN. More specifically, the cooling time of a gas is conven-
The distribution of galaxy light in our ‘best’ model is tionally taken as the ratio of its specific thermal energy to
shown in the bottom panel of Fig. 1 for comparison with thecooling rate per unit volume,
the mass distribution in the top panel. For both the vol-
3 µ¯m kT
ume is a projected 330 280 15h−1Mpc slice cut from t = p . (3)
cool
the full 0.125h−3Gpc3 s×imulati×on box. The plot of surface 2ρg(r)Λ(T,Z)
brightnessiscolour-codedbytheluminosity-weightedmean Here µ¯m is the mean particle mass, k is the Boltzmann
p
B Vcolourofthegalaxies.Onlargescaleslightclearlyfol- constant, ρ (r) is the hot gas density, and Λ(T,Z) is the
− g
lowsmass,butnon-trivialbiasesbecomeevidentonsmaller coolingfunction.Thelatterdependsbothonthemetallicity
scales, especially in ‘void’ regions. The redder colour of Z andthetemperatureofthegas.Inourmodelsweassume
galaxies in high density regions is also veryclear. the post-shock temperature of the infalling gas to be the
virial temperature of the halo, T = 35.9(Vvir/kms−1)2K.
When needed, we assume that the hot gas within a static
3.2 Gas infall and cooling atmosphere has a simple ‘isothermal’ distribution,
m
WefollowthestandardparadigmsetoutbyWhite&Frenk ρ (r)= hot , (4)
(1991) asadaptedforimplementationonhighresolutionN- g 4πRvirr2
body simulations by Springel et al. (2001a) and De Lucia where mhot is the total hot gas mass associated with the
et al. (2004). This assumes that as each dark matter halo halo and is assumed to extendto its virial radius Rvir.
collapses,itsown‘fairshare’ofcosmicbaryonscollapsewith Toestimateaninstantaneouscoolingrateontothecen-
it (but see Section 3.3 below). Thus in our model the mass tralobjectofahalo,givenitscurrenthotgascontent,wede-
fraction in baryons associated with every halo is taken to finethecoolingradius,r ,astheradiusatwhichthelocal
cool
be f = 17%, consistent with the first-year WMAP result cooling time(assuming thestructureof Eq.4) is equalto a
b
(Spergeletal.2003). Initially thesebaryonsareintheform suitablydefinedageforthehalo.White&Frenk(1991)took
of diffuse gas with primordial composition, but later they this age to be the Hubble time tH, while Cole et al. (1994,
includegas in several phasesas well as stars and heavyele- 2000) used the time scale over which the main progenitor
ments.Thefateoftheinfalling gasdependsonredshiftand last doubled its mass. Somerville & Primack (1999) argued
on the depth of the halo potential (Silk 1977; Rees & Os- that the time interval since the last major merger is more
triker1977;Binney1977;White&Rees1978).Atlatetimes appropriate since such mergers redistribute hot gas within
andin massive systemsthegasshocks tothevirial temper- thehalo. Herewe follow Springel et al. (2001a) and DeLu-
(cid:13)c 2005RAS,MNRAS000,000–000
Cooling flows, black holes and the luminosities and colours of galaxies 7
Figure 2.Theratioofthecoolingradiustovirialradiusforarandomselectionofvirialisedsystemsatz=0plottedagainsttheirdark
matter virial mass. Systems identified to be in the ‘rapid cooling’ regime are shown in the left panel, while those that have formed a
static hot halo are shown on the right (Section 3.2). A sharp transition between the two regimes is seen close to that found by Keres
etal.(2004), markedbythesolidverticalline.
cia et al. (2004) and define the cooling radius as the point accretion shock at large radius. The high-resolution spheri-
where the local cooling time is equal to thehalo dynamical calinfall simulationsofForcada-Mir´o &White(1997) show
time, Rvir/Vvir = 0.1H(z)−1. This is an order of magni- that in this situation the accretion shock moves inwards,
tudesmaller than tH and so results in substantially smaller thepost-shocktemperatureincreasesandthemassstoredin
cooling radii and cooling rates (typically by a factor of 3) thepost-shock hot atmosphere decreases, because thepost-
than the assumption of White & Frenk (1991). Our choice shockgasrapidlycoolsontothecentralobject.Ineffect,all
isjustifiedbythetestsofYoshidaetal. (2002) whoverified infalling material is accreted immediately onto the central
explicitlythatitresultsingoodobject-by-objectagreement disk.Inthisrapidcoolingregime wethereforesetthecooling
between theamount of gas predicted to condense in galaxy rateontothecentralobjecttobeequaltotherateatwhich
mass haloes and the amount which actually condensed in new diffusegas is added to thehalo.
their high-resolution SPH simulations of the formation of
a galaxy cluster and its environment. These tests assumed
primordial abundances in the cooling function. When we
3.2.1 Rapid cooling or cold accretion?
implementachemicalenrichmentmodelconsistentwiththe
observed element abundances in intracluster gas (see Sec- Although much simplified, the above model was shown by
tion 3.9) cooling rates in galaxy mass haloes are substan- Yoshidaet al. (2002) andHelly et al. (2003) togivereason-
tially enhanced and we find (as did De Lucia et al. 2004) ably accurate, object-by-object predictions for the cooling
that asmaller coefficient thanused in theoriginal White& and accumulation of gas within the galaxies that formed
Frenk cooling model is required to avoid excessive conden- in their N-body+SPH simulations. These neglected star-
sation of gas. formation and feedback effects in order to test the cooling
Using the above definition, a cooling rate can now be model alone. They also assumed primordial abundances in
determined through a simple continuity equation, the cooling function. In Fig. 2 we show the ratio rcool/Rvir
as a function of virial mass for haloes in the ‘rapid cool-
m˙ =4πρ (r )r2 r˙ . (5)
cool g cool cool cool ing’regime(leftpanel)andinthe‘statichalo’regime(right
Despiteitssimplicity,thisequationisagoodapproximation panel) at z=0 for our ‘best’ galaxy formation model. The
to the rate at which gas is deposited at the centre in the two regimes are distinguished by the dominant gas compo-
Bertschinger (1989) similarity solution for a cooling flow. nent in each halo: when the mass of hot halo gas exceeds
Putting it all together we take the cooling rate within a that of cold disk gas, we say thegalaxy hasformed astatic
halo containing a hot gas atmosphere to be halo, otherwise the system is taken to be in the rapid cool-
ingphase.Manyhaloesclassifiedas‘rapidlycooling’bythis
m˙cool =0.5mhotrcoRol2Vvir . (6) criterionhavercool <Rvir,whichwouldapparentlyindicate
vir a static hot halo. This is misleading, however, as systems
Weassumethisequationtobevalidwhenrcool <Rvir.This where cooling is rapid deposit infalling gas ontothe central
is the criterion which White & Frenk (1991) proposed to galactic disk on a short timescale until they have a low-
separate the static hot halo regime from the rapid cooling mass residual halo which satisfies rcool Rvir. This then
∼
regime (see below). persists until the next infall event. Averaging over several
Inlowmasshaloesorathighredshiftstheformalcool- cycles of this behaviour, one finds that the bulk of the in-
ing radius lies outside the virial radius rcool > Rvir. The falling gas cools rapidly. This is why we choose to classify
post-shock gas then cools in less than one sound crossing systemsbytheirdominantgascomponent.Notealsothata
timeandcannotmaintainthepressureneededtosupportan massive hot halo forms immediately once cooling becomes
(cid:13)c 2005RAS,MNRAS000,000–000
8 Croton et al.
inefficient, just as in the 1-D infall simulations of Forcada- tion heating of the intergalactic medium (IGM) which sup-
Mir´o & White (1997) and Birnboim & Dekel (2003). Our presses the concentration of baryons in shallow potentials
classification is thusquiterobust. (Efstathiou 1992). More recent models identify thepossible
The transition between the ‘rapid cooling’ and ‘static low-redshiftsignatureofsuchheatinginthefaintendofthe
halo’regimesisremarkablywelldefined.Atz=0itoccursat galaxy luminosity function, most notably in the abundance
ahalovirialmassof2–3 1011M(cid:12),andisapproximatelyin- of the dwarf satellite galaxies in the local group (e.g. Tully
×
dependentofredshiftouttoatleastz=6.Thisisclosetothe et al. 2002; Benson et al. 2002).
transition mass found for the same cosmology by Birnboim Gnedin (2000) showed that the effect of photoioniza-
&Dekel(2003)usingsphericallysymmetricsimulations,and tion heating on the gas content of a halo of mass Mvir can
byKeresetal.(2004)usingfully3-Dsimulations.Thisis,in be modelled by defining a characteristic mass scale, the so
fact,acoincidencesincebothsetsofauthorsassumecooling calledfiltering mass,M ,belowwhichthegasfractionf is
F b
functionswithnometals,whereasour‘best’modelincludes reduced relative to theuniversal value:
heavy elements which substantially enhance cooling at the fcosmic
t‘reamppidercaotoulrinesg’rreelegvimanetafroirsetshferotmrantshietimoni.xiEnngroicfhinmfeanlltinignpthrie- fbhalo(z,Mvir)= (1+0.26Mb F(z)/Mvir)3 . (7)
mordialgaswithsupernova-enrichedgasejectedinagalactic Thefilteringmassisafunctionofredshiftandchangesmost
wind (see Section 3.6). If we modify our model to assume significantly around the epoch of reionization. It was esti-
a zero-metal cooling function, we find the transition mass mated by Gnedin using high-resolution SLH-P3M simula-
to shift downwards by about a factor of 2–3, resulting in a tions(butseeHoeftetal.2005).Kravtsovetal.(2004) pro-
lowercoolingrateincomparison totheworkofBirnboim& videdananalyticmodelfortheseresultswhichdistinguishes
Dekel(2003) and Keres et al. (2004). three‘phases’intheevolutionoftheIGM:z>z0,theepoch
Thereasonforthemoreefficient(wewouldargueoverly before the first HII regions overlap; z0<z<zr, the epoch
efficient)coolingappearstobedifferentinthesetwostudies. when multiple HII regions overlap; z<zr, the epoch when
The spherical infall simulations of Forcada-Mir´o & White themedium isalmost fully reionized. They findthat choos-
(1997) showed good agreement with a transition at the ingz0=8andzr=7providesthebestfittothenumerically
point predicted using the original cooling radius definition determined filtering mass. We adopt these parameters and
of White & Frenk (1991) rather than our revised definition keepthemfixedthroughoutourpaper.Thischoiceresultsin
which was checked explicitly against SPH simulations by afilteringmass of 4 109M(cid:12) at z=zr,and 3 1010M(cid:12) by
× ×
Yoshidaetal.(2002).Sphericalmodelsthuspredictmoreef- the present day. See Appendix B of Kravtsov et al. (2004)
ficientcoolingthanoccursintypical3-Dsituations.Thisex- for a full derivation and description of the analytic model.
plains,perhaps,whyBirnboim&Dekel(2003)findahigher
transition mass than we would predict for zero metallicity.
Yoshida et al. (2002) also showed that the density estima- 3.4 Black hole growth, AGN outflows, and
tion in the implementation of SPH used by Keres et al. cooling suppression
(2004) leads toovercooling ofgas in galaxy mass objects as
There is a growing body of evidence that active galactic
compared to theirown entropy and energy conserving SPH
nuclei (AGN) are a critical piece in the galaxy formation
scheme;theeffectislargeenoughtoexplainwhyKeresetal.
puzzle. Our principal interest in this paper is their role in
(2004) find a higher transition mass than we find for their
suppressing cooling flows, thereby modifying the luminosi-
assumed cooling function.
ties, colours, stellar masses and clustering of the galaxies
Both Birnboim & Dekel and Keres et al. refer to the
that populatethebright endof thegalaxy luminosity func-
‘rapidcooling’ regimeas‘cold infall’. Thisis,in fact,amis-
tion. To treat this problem, we first need a physical model
nomer.Inthismodetheaccretionshockoccursclosertothe
for thegrowth of black holes within our galaxies.
centralobject,andsodeeperinitspotentialwellthanwhen
thereisastatichothalo.Asaresult,thepre-shockvelocity
of infalling gas is greater, resulting in a larger post-shock
3.4.1 The ‘quasar mode’
temperature. The two modes do not differ greatly in the
temperaturetowhichinfalling gas isshocked,butratherin Inourmodel(whichisbasedcloselyonthatofKauffmann&
how long (compared to the system crossing time) the gas Haehnelt2000)supermassiveblackholesgrowduringgalaxy
spends at the post-shock temperature before its infall en- mergers both by merging with each other and by accretion
ergy is lost toradiation. Finally, we notethat theexistence of cold disk gas. For simplicity, black hole coalescence is
andimportanceofthesetwomodeswasthemajorinsightof modelled as adirectsum of theprogenitormasses andthus
the original work of Silk (1977), Binney (1977) and Rees & ignores gravitational wave losses (including such losses is
Ostriker(1977) and has been built into all modern theories foundtohavelittleeffectonthepropertiesofthefinalgalaxy
forgalaxyformation.Adetaileddiscussioncanbefound,for population).Weassumethatthegasmassaccretedduringa
example, in White& Frenk(1991). mergerisproportionaltothetotalcoldgasmasspresent,but
with an efficiency which is lower for smaller mass systems
and for unequalmergers. Specifically,
3.3 Reionization
f0 m
Abeccirneetffiiocnienatndtoceoxoplilnagininwhlyowdwmarafssgahlaaxloieess ciosnrteaqinuiaredreltao- ∆mBH,Q= 1+(280BHkmsc−o1ld/Vvir)2 , (8)
tivelysmallfractionofallcondensedbaryons(White&Rees wherewehavechangedtheoriginalparameterisationbytak-
1978).Thisinefficiencymayinpartresultfromphotoionisa- ing
(cid:13)c 2005RAS,MNRAS000,000–000
Cooling flows, black holes and the luminosities and colours of galaxies 9
fB0H =fBH (msat/mcentral) . (9)
HerefBH 0.03isaconstantandischosentoreproducethe
≈
observedlocalmBH mbulgerelation(Magorrianetal.1998;
−
Marconi & Hunt 2003; H¨aring & Rix 2004). In contrast to
Kauffmann & Haehnelt (2000) we allow black hole accre-
tionduringbothmajor and minormergersalthoughtheef-
ficiency in the latter is lower because of the msat/mcentral
term. Thus, any merger-induced perturbation to the gas
disk (which might come from a bar instability or a merger-
induced starburst – see Section 3.7) can drive gas onto the
central black hole. In this way, minor merger growth of the
black hole parallels minor merger growth of the bulge. The
fractionalcontributionofminormergerstobothistypically
quitesmall,sothataccretiondrivenbymajormergersisthe
dominantmodeofblackholegrowthinourmodel.Werefer
to this as the ‘quasar mode’. [Note that a more schematic
treatment of black hole growth would suffice for the pur-
poses of this paper, but in Springel et al. (2005b) and in Figure3.Theblackholeaccretionratedensity,m˙BH,asafunc-
future work we wish to examine the build-up of the black tion of redshift for both the ‘quasar’ and the ‘radio’ modes dis-
holepopulationwithingalaxiesinconsiderablymoredetail.] cussedinSection3.4.Thisfigureshowsthatthegrowthofblack
There is substantial evidence for strong hydrodynamic holesisdominatedbythe‘quasarmode’athighredshiftandfalls
and radiative feedback from optical/UV and X-ray AGN offsharplyatz∼<2.Incontrast,the‘radiomode’becomesimpor-
tantatlowredshiftswhereitsuppressescoolingflows,butisnot
(Arav et al. 2001; de Kool et al. 2001; Reeves et al. 2003;
asignificantcontributor totheoverallblackholemassbudget.
Crenshaw et al. 2003). We have not yet explicitly incorpo-
rated such feedback in our modelling, and it may well turn
outtobeimportant(see,forexample,therecentsimulations the cooling flow described in Section 3.2. We take the me-
ofDiMatteoetal.2005;Springeletal.2005a;Hopkinsetal. chanical heating generated by the black hole accretion of
2005). Weassume ‘quasar mode’ accretion to be closely as- Eq.10 to be
sociated with starbursts, so this feedback channel may be
partially represented in our models by an enhanced effec- LBH =ηm˙BHc2 , (11)
tive feedback efficiency associated with star formation and where η = 0.1 is the standard efficiency with which mass
supernovae(see Section 3.6). isassumed to produceenergy neartheeventhorizon,and c
is the speed of light. This injection of energy compensates
in part for the cooling, giving rise to a modified infall rate
3.4.2 The ‘radio mode’
(Eq.6) of
In our model, low energy ‘radio’ activity is the result of
hot gas accretion onto a central supermassive black hole m˙0cool =m˙cool− 1LVBH2 . (12)
once a static hot halo has formed around the black hole’s 2 vir
host galaxy. We assume this accretion to be continual and Forconsistencyweneverallowm˙0 tofallbelowzero.Itis
cool
qmuoidesecl:entandtobedescribedbyasimplephenomenological awnodrtmh˙hneoattin≡g2tLhBaHt/m˙Vcv2oirol∝∝mfBh3o/Ht2fΛh(oVtVvivri)r1/(2EVqv.2ir12t−H),1/s2o(tEhqat. 6)
mw˙hBeHr,eRm=BκHAGisN(cid:16)th1e0m8bBMlaHc(cid:12)k(cid:17)h(cid:16)of0lhe.o1tm(cid:17)a(cid:16)ss2,0f0hVkovmtiriss−1th(cid:17)e3f,raction(10o)f mm˙˙hcoeaotl ∝ fh1o/t2Λm(BVHvirt)1H1//22Vvir . (13)
the total halo mass in the form of hot gas, Vvir ∝ Tv1i/r2 is (ThesescalingsareexactforanEdSuniverse;wehaveomit-
thevirialvelocityofthehalo,andκAGN isafreeparameter ted weak coefficient variations in other cosmologies.) Thus
with units of M(cid:12)yr−1 with which we control the efficiency in our model the effectiveness of radio AGN in suppressing
of accretion. We find below that κAGN = 6 10−6M(cid:12)yr−1 cooling flows is greatest at late times and for large values
×
accurately reproduces the turnoverat the bright end of the of black hole mass. This turns out to be the qualitative be-
galaxy luminosity function. Note that fhotVv3irtH is propor- haviourneededforthesuppressionofcoolingflowstorepro-
tional to the total mass of hot gas, so that our formula is ducesuccessfully theluminosities, colours and clustering of
simplytheproductofthehotgasandblackholemassesmul- low redshift bright galaxies.
tiplied by a constant efficiency and divided by the Hubble
time tH. In fact, we find fhot to be approximately constant 3.4.3 The growth of supermassive black holes
for Vvir > 150kms−1, so the dependence of m˙BH,R on this
quantity∼haslittleeffect.TheaccretionrategivenbyEq.10 Fig. 3 shows theevolution of themean black hole accretion
is typically orders-of-magnitude below theEddington limit. rate per unit volume averaged over the entire Millennium
InSection5wediscussphysicalaccretionmodelswhichmay Simulation.Weseparatetheaccretion intothe‘quasar’ and
reproduce thisphenomenology. ‘radio’ modes described above (solid and dashed lines re-
Weassumethat‘radiomode’feedbackinjectssufficient spectively). Black hole mass growth in our model is dom-
energyintothesurroundingmedium toreduceorevenstop inated by the merger-driven ‘quasar mode’, which is most
(cid:13)c 2005RAS,MNRAS000,000–000
10 Croton et al.
Figure 5. The star formation rate density of the universe as a
functionofredshift.Thesymbolsshowacompilationofobserva-
tionalresultstakenfromFig.12ofSpringel&Hernquist(2003b).
Thesolidlineshowsour‘best’model,whichpredictsthatgalax-
iesformmuchoftheirmassrelativelyearly.Thedashedline(and
rightaxis)indicatetheincreaseinstellarmasswithredshift.Ap-
Figure4.Theblackhole-bulgemassrelationformodelgalaxies proximately50%ofallstarsformbyz=2.
atthepresentday.ThelocalobservationalresultofHa¨ring&Rix
(2004)isgivenbythesolidline,wherethedashedboxshowsthe
approximaterangeoverwhichtheirfitwasobtained. (1996),thiscriticalsurfacedensityatadistanceRfromthe
galaxy centremay beapproximated by
efficient at redshifts of two to four, dropping by a factor Σcrit(R)=120(cid:16)200Vkvmirs−1(cid:17)(cid:16)kRpc(cid:17)−1M(cid:12)pc−2 , (14)
of five by redshift zero. This behaviour has similar form to
We convert this critical surface density into a critical mass
but is weaker than the observed evolution with redshift of
byassumingthecoldgasmasstobeevenlydistributedover
thebrightquasarpopulation(e.g.Hartwick&Schade1990).
thedisk. The resulting critical cold gas mass is:
(SeealsothediscussioninKauffmann&Haehnelt2000).In
cexonptercatsetd,ofruorm‘rathdeiosmcaolidneg’idsisscigunsisfiecdaanbtoovnel,yaantdlaftoerttihmeehs,igahs mcrit =3.8×109(cid:16)200Vkvmirs−1(cid:17)(cid:16)1r0dkispkc(cid:17)M(cid:12) , (15)
feedbackefficiencyassumedinEq.11itcontributesonly5%
of thefinalblack holemass density.Wewill show, however, whereweassumethediskscalelengthtobers=(λ/√2)Rvir
that the outflows generated by this accretion can have a (Moetal.1998),andsettheouterdiskradiustordisk=3rs,
major impact on the final galaxy properties. Finally, inte- based on the properties of the Milky Way (van den Bergh
grating the accretion rate density over time gives a present 2000).Hereλisthespinparameterofthedarkhaloinwhich
dayblackholemassdensityof3 105M(cid:12)Mpc−3,consistent the galaxy resides (Bullock et al. 2001), measured directly
with recent observational estima×tes (Yu & Tremaine 2002; from the simulation at each time step. When the mass of
Merloni 2004). cold gas in a galaxy is greater than this critical value we
assume thestar formation rate to be
The relationship between black hole mass and bulge
mass is plotted in Fig. 4 for the local galaxy population
m˙∗ =αSF (mcold mcrit) / tdyn,disk , (16)
in our ‘best’ model. In this figure, the solid line shows the −
bestfittotheobservationsgivenbyH¨aring&Rix(2004)for wheretheefficiency αSF istypically set so that5to15 per-
asampleof30nearbygalaxieswithwellmeasuredbulgeand cent of the gas is converted into stars in a disk dynamical
black holemasses. Theirresults only probemasses overthe time tdyn,disk, which we define to be rdisk/Vvir. This star
rangeboundedbythedashedlines.Ourmodelgalaxiespro- formation model produces a global star formation history
duces a good match to these observations with comparable consistent with the observed star formation density of the
scatter in theobserved range (see their Fig. 2). universe out to at least z=2, as shown in Fig. 5 (note that
thisfigurealso includesstarformation throughstarbursts–
see Section 3.7) .
When implemented in our model, Eq. 16 leads to
3.5 Star formation
episodic star formation that self-regulates so as to main-
We use a simple model for star formation similar to those tain the critical surface density of Eq. 14. This naturally
adopted by earlier authors. We assume that all star forma- reproduces the observed spiral galaxy gas fractions with-
tion occurs in cold disk gas, eitherquiescently orin a burst out theneed for additional parameterisation, as we demon-
(see Section 3.7). Based on the observational work of Ken- strate in the top panel of Fig. 6 using model Sb/c galaxies
nicutt (1998), we adopt a threshold surface density for the identified as objects with bulge-to-total luminosity: 1.5 6
cold gas below which no stars form, but above which gas M M 6 2.5 (bulge formation is described in
I,bulge I,total
−
starts to collapse and form stars. According to Kauffmann Section 3.7).
(cid:13)c 2005RAS,MNRAS000,000–000