Cosmological simulations of black hole growth: AGN luminosities and downsizing

24 Pages·2014·3.38 MB·English

Mon.Not.R.Astron.Soc.000,1–24(2002) Printed3June2014 (MNLATEXstylefilev2.2)
Cosmological simulations of black hole growth: AGN
luminosities and downsizing
4 Michaela Hirschmann1⋆, Klaus Dolag2,3, Alexandro Saro2, Lisa Karin Bachmann2,
1
0 Stefano Borgani1,4,5, Andreas Burkert2,6
2 1INAF - Astronomical Observatory of Trieste,viaG.B. Tiepolo 11, I-34143 Trieste,Italy
n 2Universita¨tsSternwarte Mu¨nchen, Scheinerstr.1, D-81679 Mu¨nchen, Germany
u 3Max-Plank-Institute fu¨r Astrophysik, Karl-Schwarzschild Strasse 1, D-85740 Garching, Germany
J 4Astronomy Unit,Department of Physics, Universityof Trieste,via G.B. Tiepolo 11, I-34131 Trieste, Italy
1 5INFN-National Institute for Nuclear Physics, Via Valerio2, I-34127 Trieste,Italy
6Max-Planck-Institut fu¨r extraterrestrische Physik (MPE), Giessenbachstr. 1, 85748 Garching, Germany

]
O
C Accepted ???.Received???inoriginalform???
.
h
p ABSTRACT
-
In this study, we present a detailed, statistical analysis of black hole growth and the
o
r evolution of active galactic nuclei (AGN) using cosmological hydrodynamic simula-
t tions run down to z = 0. The simulations self-consistently follow radiative cooling,
s
a star formation, metal enrichment, black hole growth and associated feedback pro-
[ cesses from both supernovae typeII/Ia and AGN. We consider two simulation runs,
one with a large co-moving volume of (500 Mpc)3 and one with a smaller volume of
3
(68 Mpc)3 but with a by a factor of almost 20 higher mass resolution. We compare
v
3 the predicted statistical properties of AGN with results from large observational sur-
3 veys. Consistently with previous results, our simulations can widely match observed
3 black hole properties of the local Universe.Furthermore, our simulations can success-
0 fully reproducethe evolutionof the bolometric AGN luminosity function for both the
. low-luminosity and the high-luminosity end up to z = 3.0, only at z = 1.5−2.5,
8
0 the low luminosity end is over-estimatedby up to 1 dex. In addition, the smaller but
3 higher resolution run is able to match the observational data of the low bolometric
1 luminosity end at higher redshifts z = 3−4. We also perform a direct comparison
: withthe observedsoftandhardX-rayluminosityfunctions ofAGN, including anem-
v
pirical correction for a torus-level obscuration, and find a similarly good agreement.
i
X These results nicely demonstrate that the observed “anti-hierarchical” trend in the
r AGN number density evolution (i.e. the number densities of luminous AGN peak at
a higher redshifts than those of faint AGN) is self-consistently predicted by our simu-
lations. Implications of this downsizing behaviour on active black holes, their masses
and Eddington-ratios are discussed. Overall, the downsizing behaviour in the AGN
number density as a function of redshift can be mainly attributed to the evolution of
the gas density in the resolved vicinity of a (massive) black hole (which is depleted
with evolving time as a consequence of star formation and AGN feedback).
Keywords: galaxies:active,galaxies:evolution,galaxies:formation,galaxies:nuclei,
quasars: general, methods: numerical
1 INTRODUCTION Gebhardt et al. 2000; Tremaine et al. 2002; H¨aring & Rix
2004; Graham & Driver 2007; Gu¨ltekin et al. 2009;
Itisgenerallyacceptedthatpresent-dayspheroidalgalaxies
Burkert & Tremaine 2010; McConnell & Ma 2013) as the
host supermassive black holes (BHs) at their centres
bulge mass, the stellar velocity dispersion, the Sersic
(Magorrian et al.1998;Genzel & Eckart1999).Inaddition,
index or the number of globular clusters. This if often
strong correlations are found between BH masses and
interpreted as an evidence for a co-evolution between the
properties of their host galaxies (Ferrarese & Merritt 2000;
spheroidal component of host galaxies and their BHs (see
Kormendy& Ho(2013)forarecentreview).However,these
⋆ E-mail:[email protected] relations were also shown to possibly have a “non-causal”
2 Hirschmann et al.
origin and thus, may be just a consequence of statistical lution of their host galaxies, a large amount of studies
merging processes (Peng 2007; Hirschmann et al. 2010; were published based on semi-analytic models of galaxy
Jahnke& Maccio` 2011). In addition, some recent observa- formation within which mechanisms for BH formation
tionsofBHgrowthandstarformation inindividualobjects and evolution were included (e.g. Kauffmann & Haehnelt
appear to contradict the picture of a simple one-to-one 2000;Volonteri et al.2003;Granato et al.2004;Menci et al.
co-evolutionovertime(e.g.Mullaney et al.2012;Bongiorno 2004; Bromley et al. 2004; Monaco & Fontanot 2005;
2012;Rosario2012;Hickoxet al.2014).Thus,thedetailsof Croton 2006; Bower et al. 2006; Marulli et al. 2008;
a connected growth of BHs and their host galaxies remain Somerville et al. 2008; Bonoli et al. 2009; Fanidakis et al.
an unresolved puzzle. 2012; Hirschmann et al. 2012; Neistein & Netzer 2013;
Duringtheirlifetime, BHsareassumedtoundergosev- Menci et al. 2013). It was shown by some of the above
eral episodes of significant gas accretion, during which this studies that an anti-hierarchical behaviour can be ex-
accretion powers luminous quasars or active galactic nuclei plained within a hierarchical structure formation scenario
(AGN)(Salpeter1964;Zel’Dovich1964;Lynden-Bell1969). (Marulli et al. 2008; Bonoli et al. 2009; Fanidakis et al.
By estimating the total energy radiated by AGN overtheir 2012;Hirschmann et al.2012;Menci et al.2013),whencon-
wholelifetime,itcanbeshownthatnearlyallthemassseen sideringdifferentmodelsforBHaccretionandAGNtrigger-
in dormant BHs today can be accumulated during the pe- ing and/or accounting for dust obscuration. Unfortunately,
riods of observed bright AGN activity (Soltan 1982). This these different models have not been able to draw a com-
implies that there is not much room left for “dark” or ob- pletely uniform picture for the origin of the downsizing be-
scured accretion. haviouras they differ in many details.
Recent progress in detecting AGN (particularly faint In addition to these semi-analytic models, also a
and obscured ones) was achieved by analysing data from large number of other studies investigated the predic-
X-ray surveys in the hard and soft X-ray range (XMM- tions of the cosmological ΛCDM model for the formation
Newton,Chandra,ROSAT,ASCA,eRositae.g.Miyaji et al. and evolution of super-massive BHs (SMBH) and AGN.
2000; La Franca et al. 2002; Cowie et al. 2003; Fiore et al. Some works made predictions based on nearly purely an-
2003; Barger et al. 2003; Ueda et al. 2003; Hasinger et al. alytic models (Efstathiou & Rees 1988; Haehnelt & Rees
2005; Barger & Cowie 2005; Sazonov & Revnivtsev 2004; 1993; Haiman & Loeb 1998) and on semi-empirical models
Nandraet al. 2005; Ebrero et al. 2009; Aird et al. 2010; (Shankaret al. 2009, 2010, 2012).
Fiore et al. 2012; Kolodzig et al. 2013). Interestingly these Several years ago, numerical hydrodynamic simula-
studies revealed a possibly at first sight puzzling behaviour tions of isolated galaxy mergers started to include pre-
of AGN concerning the evolution of the co-moving number scriptions for BH growth and AGN feedback using “sub-
density of AGN: the number density of successively more grid” recipes, where BH accretion according to the Bondi-
luminous AGN peaks at higher redshifts than the one of Hoyle-Littletonformulawasassumed(Springel et al.2005a;
less luminous AGN,with thelowest luminosity AGN show- Hopkinset al.2006;DiMatteo et al.2005;Robertson et al.
ing an almost constant co-moving number density. Making 2006; Li et al. 2007; Sijacki et al. 2007; Johansson et al.
the simplified (and naive) assumption that AGN luminos- 2009; Fabjan et al. 2010;Barai et al. 2013;Choi et al. 2013
ity (i.e. BH accretion) is proportional to BH mass (which and references therein). Feedback from AGN, however was,
we would expect if BHs are accreting at the Eddington implemented in different ways, ranging from a purely ther-
rate, L M•) would imply that very massive BHs seem mal energy injection (e.g. Springel et al. 2005a) and a fur-
∝
tobealreadyinplaceatveryearlytimes,whereaslessmas- therdistinctionbetweenradio-andquasarmode(byinject-
sive BHs seem to evolve predominantly at lower redshifts. ing bubbles,e.g. Sijacki et al. 2007, or just upscaling of the
Thisbehaviouriscalled“downsizing”or“anti-hierarchical” feedback efficiencies, e.g. Fabjan et al. 2010) to mechanical
growthofBHs.Thedownsizingtrendisalsoseenintheopti- momentum-drivenwinds (e.g. Choi et al. 2013).
cal(Cristiani et al.2004;Croom et al.2004;Fan et al.2004; To assess statistical properties of BHs and AGN in
Huntet al.2004;Richardset al.2006;Wolf et al.2003)and a full cosmological context and in a more self-consistent
thenear-infrared (NIR) (e.g. Matute et al. 2006). way than it can be achieved by e.g. semi-analytical
Thisobservationalresultseemstobeinconflictwiththe or semi-empirical models, large cosmological, hydrody-
“naive”expectationsarisingfromthecurrentlyfavouredhi- namic simulations were performed by several recent
erarchical structure formation paradigm based on the Cold studies (DiMatteo et al. 2008; McCarthy et al. 2010;
Dark Matter (CDM) model (Peebles 1965; White & Rees Degraf et al. 2010; McCarthy et al. 2011; Degraf et al.
1978; Blumenthal et al. 1985). In this framework, low mass 2011; Booth & Schaye 2011; DiMatteo et al. 2012;
halos are expected to form first and more massive halos to DeGraf et al. 2012; Haas et al. 2013; Rosas-Guevara et al.
grow later over time via subsequent merging and smooth 2013; Khandai et al. 2014) including BH growth and
accretion. However, it is now well known that the evolu- feedback from AGN. Compared to semi-analytical models,
tionary history of observable galaxies also follows an anti- they havethe advantage that the dynamics of thebaryonic
hierarchical or downsizing behaviour, with several indepen- component (gas and stars) and the interaction between
dent observational indicators suggesting that more massive baryonic matter can be followed directly in a cosmological
galaxies formed their stars and had their star formation context,even if the spatial and mass resolution, at present,
quenched earlier than low-mass galaxies, which continue isnothigh enoughtoaccurately simulatesmall scalephysi-
forming stars tothepresentday (seeFontanot et al. (2009) cal processes like e.g. the formation of stars and BHs with
for an overview). theassociated feedback.Instead,theyhavetobecomputed
Totheoretically exploretheformation and evolution of inasimplifiedmannerwithsub-grid/sub-resolutionmodels.
BHs and AGN and its (possible) connection to the evo- Nevertheless, they provide a powerful tool to assess the
Cosmological simulations of black hole growth 3
cosmic evolution of statistical AGN and BH properties and which are compared to observational data. We study the
thus, go clearly beyond the hydrodynamic isolated merger evolution of the AGN luminosity function in Section 4 at
simulations and/or semi-analytic models. (z = 0 5), considering bolometric, soft and hard X-ray
−
Di Matteo et al.(2008)andBooth & Schaye(2011),for and (for z = 0 only) also radio luminosities, and compare
example, investigated the global history of black hole mass ourresultstoobservationaldata.InSection5,weexplicitly
assembly and the evolution of the BH-bulge mass relations showtheevolution oftheAGNnumberdensityand discuss
in cosmological simulations with a maximum box size of the anti-hierarchical trend in BH growth, its consequences
(50 Mpc)3. Such simulations have also been successful on the connection between BH masses and AGN luminosi-
∼
inreproducingthelow-luminosityendoftheAGNluminos- ties and its numerical origin in our simulations. In theend,
ityfunction(Degraf et al.2010)orBHclusteringproperties in Section 6,we summarise and discuss our main results.
(Degraf et al.2011)buthavebeenlimitedbyaboxsizebe-
ing not large enough to follow the evolution of luminous
AGN with Lbol > 1045 erg/s (i.e. AGN luminosities above 2 THE SIMULATIONS
the exponential cut-off of the luminosity function). In a re-
centstudybyKhandai et al.(2014)theyhaveusedalarger In this paper, we analyse a subset of cosmological boxes
cosmological volume of (100 Mpc/h)3 (where the luminous from the Magneticum Pathfinder simulation set (Dolag
endisalsoaccessible),buttheytypicallyunderestimatedlu- et al., in prep.). From this set, we selected a large size
minousAGNatz >0.5andoverestimatedthematz<0.5. cosmological simulation with a co-moving box size of
Focusingonlyonhighredshifts,Di Matteo et al.(2012)and (500Mpc)3, simulated with an initial particle number of
DeGraf et al. (2012) investigated the growth of the first 2 1,5643 (500Mpc/hr,whichisourlarge“high-resolution”
×
quasarsappearing in theUniversebyperforming cosmolog- runthroughoutthisstudy)andonecosmological simulation
icalsimulationswithalargeboxsizeof(500Mpc/h)3,how- with a smaller co-moving box size of (68Mpc)3 (only run
ever, only run down to redshift z =5 findinga rapid, early downtoz=1),withaparticlenumberof2 5763 andthus,
×
growth of BHs drivenby cold, infalling gas flows. an increased (“ultra-high”) resolution (68Mpc/uhr). Table
Inthiswork,weextendtheirstudiesofthecosmicevo- 1 gives an overview of the two simulation runs. We have
lution of statistical properties of BH growth and AGN evo- adopted a ΛCDM model with parameters chosen according
lution at lower redshifts, from z = 5 down to z = 0. We to the seven-year Wilkinson Microwave Anisotropy Probe
analyse a subset of cosmological, hydrodynamical simula- WMAP7 data (Komatsu et al. 2011) with Ωm = 0.272,
tionsfromtheMagneticumPathfinder simulationset(Dolag Ωb = 0.0456 ΩΛ = 0.728 and h = 0.704. The initial power
etal.,inprep),whichincludesa very large cosmological box spectrumfollowsn=0.963andisnormalisedtoσ8 =0.809.
providingavolumeof(500Mpc)3.Thesimulationsincludea
self-consistentmodelforBHgrowthusingtheGadget3code
2.1 The numerical code
containingseveralmodificationsformodellingthegrowthof
BHs with respect to previous studies (e.g. DiMatteo et al. Our simulations are based on the parallel cosmologi-
2008,2012; Khandaiet al. 2014). cal TreePM-SPH code P-GADGET3 (Springel 2005). The
In the course of this study, we will particularly focus code uses an entropy-conserving formulation of SPH
on the evolution of the AGN luminosity function and the (Springel& Hernquist 2002). We are using a higher order
connectedanti-hierarchical/downsizing trendin AGNnum- kernel based on the bias-corrected, sixth-order Wendland
ber density evolution with the aim to understand its origin kernel(Dehnen& Aly2012) with 295 neighbours which to-
within a hierarchical structure formation scenario. In this gether with the usage of a low-viscosity SPH scheme al-
context, we test the effect of empirically motivated dust lows us to properly track turbulence within galaxy clus-
obscuration models and discuss implications of the down- ters (Dolag et al. 2005; Donnert et al. 2013). We also allow
sizing trend on the interplay between AGN luminosities, for isotropic thermal conduction with 1/20 of the classical
BH mass and Eddington-ratios. Overall, in this study, we Spitzer value (Dolag et al. 2004). The simulation code in-
go beyond previous papers of e.g. DiMatteo et al. (2008), cludesa treatment of radiative cooling, heating from a uni-
Degraf et al. (2010) or Khandaiet al. (2014) as we exam- form time-dependent ultraviolet background and star for-
ine the evolution of the AGN luminosity function in a by mationwiththeassociated feedbackprocesses.Thelatteris
a factor of 45 larger cosmological volume with compara- basedonasub-resolutionmodelforthemultiphasestructure
∼
blygood(evenifslightlylowerthaninKhandaiet al.2014) of theinterstellar medium (Springel & Hernquist 2003).
resolution allowing us to probe particularly luminous AGN Radiative cooling rates are computed by following the
with Lbol > 1045 erg/s. In addition, the simulation code is sameprocedurepresentedbyWiersma et al.(2009).Weac-
improved in various details, regarding the SPH implemen- countfor thepresence ofthecosmic microwave background
tation used as well as the richness of physical processes in- (CMB)andofultraviolet(UV)/X-raybackgroundradiation
cluded,includingdetailsofthehandlingoftheBHsinkpar- fromquasarsandgalaxies,ascomputedbyHaardt & Madau
ticles.Herethemostnoticeabledifferenceisreflectedinour (2001).Thecontributionstocooling from eachoneof11el-
ability to follow BHs in galaxies which are inside of galaxy ements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) have been
clusters. pre-computedusingthepubliclyavailableCLOUDYphoto-
This study is organised as follows. Section 2 gives ionisationcode(Ferland et al.1998)foranopticallythingas
an overview of the simulation code and the corresponding in (photo-)ionisation equilibrium.
model for black hole growth and AGN feedback adopted in In the multiphase model for star-formation
our simulations. In Section 3 we present some basic prop- (Springel& Hernquist 2003), the ISM is treated as a
erties of present-day and high-redshift BHs and galaxies, two-phase medium where clouds of cold gas form from
4 Hirschmann et al.
Table 1.Overviewofthetwosimulationrunswhichareanalysedinthisstudy.
Name Box size Resolution Initial particle m(dm) m(gas) m(stars) Softening length
level number (dm,gas,stars)
[Mpc/h] [M⊙/h] [M⊙/h] [M⊙/h] [kpc/h]
500Mpc/hr 352 hr 2×1,5643 6.9×108 1.4×108 3.5×107 3.75,3.75,2.0
68Mpc/uhr 48 uhr 2×5763 3.6×107 7.3×106 1.8×106 1.4,1.4,0.7
2.2 The BH growth model
Most importantly, our simulations also include a pre-
scription for BH growth and for a feedback from active
galactic nuclei (AGN) based on the model presented in
Springel et al. (2005b) and Di Matteo et al. (2005) includ-
ing the same modifications as in the study of Fabjan et al.
(2010) and some new, minor changes for BH seeding and
BH “pinning” which are explained in later in this section.
As for star formation, the accretion onto BHs and the
associated feedbackadoptsasub-resolutionmodel.BHsare
represented by collision-less “sink particles” that can grow
in mass by accreting gas from their environments, or by
merging with other BHs.
The gas accretion rate M˙• is estimated by using the
Bondi-Hoyle-Lyttleton approximation (Hoyle& Lyttleton
1939; Bondi & Hoyle1944; Bondi 1952):
4πG2M2αρ
Figure 1.Histogramofthedistances oftheblackholes totheir M˙• = (c2+v2•)3/2, (1)
halo potential minimum in the 500Mpc/hr simulation. All BHs s
aremaximal2kpcawayfromthepotentialminiumwhichisrea- where ρ and c are the density and the sound speed of the
s
sonablegivenasofteninglengthof5.2kpcinthisrun. surrounding(ISM) gas, respectively,v is thevelocity of the
black hole relative to the surrounding gas and α is a boost
factor for the density and the sound speed which typically
is set to 100 as in most related works (unless a more de-
tailed description as introduced in Booth & Schaye (2009)
is used) and accounts for the fact that in cosmological sim-
ulationswecannotresolvetheintra-clustermedium(ICM)
properties within the vicinity of the BH. The BH accretion
cooling of hot gas and are embedded in the hot gas phase is always limited to the Eddington rate (maximum possi-
assuming pressure equilibrium whenever gas particles ble accretion for balance between inwards directed grav-
are above a given threshold density. The hot gas within itational force and outwards directed radiation pressure):
the multiphase model is heated by supernovae and can M˙• =min(M˙•,M˙edd).Notethatthedetailedaccretionflows
evaporate the cold clouds. A certain fraction of massive ontotheBHsareunresolved,wecanonlycaptureBHgrowth
stars (10 per cent) is assumed to explode as supernovae dueto thelarger scale gas distribution, which is resolved.
type II (SNII). The released energy by SNII (1051 erg) is Oncetheaccretionrateiscomputedforeachblackhole
modelled totrigger galactic windswith a mass loading rate particle the mass continuously grows. To model the loss of
being proportional to the star formation rate (SFR) to thisaccretedgasfromthegasparticles,astochasticcriterion
obtain a resulting wind velocity of vwind =350 km/s. isusedtoselectthesurroundinggasparticlestobeaccreted.
Oursimulations also includeadetailed modelof chem- Unlike in Springelet al. (2005b), in which a selected gas
ical evolution according to Tornatore et al. (2007). Met- particlecontributestoaccretionwithall itsmass,weinclude
als are produced by SNII, by supernovae type Ia (SNIa) the possibility for a gas particle to accrete only with a slice
and by intermediate and low-mass stars in the asymp- of its mass, which corresponds to 1/4 of its original mass.
totic giant branch (AGB). Metals and energy are released This way, each gas particle can contribute with up to four
by stars of different mass to properly account for mass- generations of BH accretion events, thus providing a more
dependent life-times (with a lifetime function according continuousdescription of theaccretion process.
to Padovani & Matteucci 1993), the metallicity-dependent Thetotalreleased energyE˙ isrelatedtotheBHaccre-
stellar yields by Woosley & Weaver (1995) for SNII, the tion rate by
ystiaelrdssanbdy tvhaenydieelndsHboyekT&hieGlermoeannenweegteanl.((12909073)) ffoorr SANGIaB. E˙ =ǫrM˙•c2, (2)
Starsof differentmass are initially distributed according to where ǫr is the radiative efficiency, for which we adopt a
a Chabrier initial mass function (IMF; Chabrier 2003). fixed value of 0.2. Here we are using a slightly larger value
Cosmological simulations of black hole growth 5
500 Mpc
z = 3 z = 2
z = 1 Box2/hr z = 0
Figure 2. Shown is a 500 Mpc wide, 70 Mpc thick slice through the cosmological baryonic mass distribution (stellar and gaseous
density) of the Box2/hr simulation at different redshift steps (z = 3,2,1 and 0). This is the result of a ray tracing visualisation using
SPLOTCH(Dolagetal.2008).White,blueandredcirclesindicatethe20BHswithinthisslicewhichhavethehighestmasses,thehighest
Eddington-ratios and the highest accretion rates, i.e. AGN luminosities,respectively. The sizes of the circles arescaled logarithmically
withthedifferentvalues,normalisedtothemaximumvalueofeachquantity.Azoomontoaregionwherethemostmassiveclusterforms
isshown,together withallgalaxieshavingastellarmasslargerthan1010M⊙ illustratedbywhitecrosses,whereasthewhitediamonds
showallBHs.
than the one standardly assumed (=0.1 ) for a radiatively which follow the metal depending cooling function, see for
efficient accretion disk onto a non-rapidly spinning BH ac- example Booth & Schaye 2011). The energy is distributed
cordingtoShakura& Sunyaev1973,(seealsoSpringel2005; kernelweighted to the surrounding gas particles in an SPH
Di Matteo et al.2005).Thischoiceismotivatedbyobserva- like manner.
tions of Davis & Laor (2011) who find higher ǫ for higher
r
BH masses and by ournumerical resolution. Additionally,weincorporatedthefeedbackprescription
Weassumethatafractionǫf ofthisenergyisthermally according to Fabjan et al. (2010): we account for a tran-
coupled tothe surroundinggas so that sition from a quasar- to a radio-mode feedback (see also
E˙f =ǫrǫfM˙•c2 (3) SanijaEckdideintgatlo.n2-0r0a7ti)owohfefneedvder:=theM˙arc/cMr˙eetdiodn<rat1e0−fa2l.lsDbuerlionwg
is the rate of the energy feedback. ǫf is a free parameter this radio-mode feedback we assume a 4 times larger feed-
and typically set to 0.15 (as usually done in simulations backefficiencythaninthequasarmode.Thisway,wewant
6 Hirschmann et al.
2.5 Mpc
Figure 3.Shown isthe25Mpcwidezoom onto the galaxycluster cluster atz=0,whereanalogous tofigure2allgalaxies withstellar
mass larger than 1010M⊙ are shown as white crosses, and all BHs are shown as white diamonds. The right panel visualises a further
zoom into the cluster. The region shown is 2.5 Mpc wide and correspond to roughly one third of the virial size of the cluster. In the
raytracingvisualisation,thewhitecoloursreflectthestellarcomponent, whilethelightbluecolourscorrespondtothehotphaseofthe
ICM.BlackdiamondsmarkalltheBHsinthesimulation.
toaccountformassiveBHs,whichareradiativelyinefficient smaller values due to the better underlying resolution (BH
(havinglowaccretionrates),butwhichareefficientinheat- seed masses of 8 104M⊙/h in galaxies with a minimum
ing the ICM by inflating hot bubbles in correspondence of stellar mass of 2.5× 109M⊙/h). While the BH then grows
×
the termination of AGN jets. The total efficiency in thera- very fast until it reaches the stellar-BH-mass relation, this
diomodeisveryclosetothevalueof0.1(=0.15 0.2 4). recovers the BH feedback within the galaxies which would
× ×
Thisisthecanonicalvalue,whichChurazov et al.(2005)es- havebeenpresentifresolutionhadallowedtoseedBHsear-
timated to beneeded to balance cooling by AGNfeedback. lier.Thisalsoavoidstoimprintanystellar-BH-massrelation
Notethatwealso,incontrasttoSpringel et al.(2005b), fromthebeginning.Finally,wechoosetheseededBHsatthe
modifythemassgrowthoftheblackholebytakingintoac- position of thestarparticle with thelargest bindingenergy
countthefeedback,e.g.∆M• =(1 ǫr)M˙•∆t.Furthermore, withintheFoFgroup,insteadofatthedarkmatterparticle
−
we introduced some additional, technical modifications of with themaximum density,as originally implemented.
theoriginal implementation which we will now summarise:
(II) In the original implementation by Springel et al.
(I) Onedifference with respect to the original implementa-
(2005b), black holes are forced to remain within the host
tion by Springel et al. (2005b) concerns the seeding of BH
galaxybypinningthemtothepositionoftheparticlefound
particles. In the implementation by Springel et al. (2005b),
having the minimum value of the potential among all the
BHparticlesareseededinahalowheneveritfirstreachesa
particles lying within the SPH smoothing length computed
minimum (total) friends-of-friends (FoF) halo mass, where
at theBH position. Within acosmological context an aside
the FoF is performed on the dark matter particles only. In
effectofthiscriterionisthat,duetotherelativelylargeval-
order to guarantee that BHs are seeded only in halos rep-
ues of SPH smoothing lengths, a BH can be removed from
resenting clearly resolved galaxies, where sufficient star for-
thehostgalaxywheneveritbecomesasatellite,andisspuri-
mationtookplace,ourimplementationperformsaFoFalgo-
ouslymergedintotheBHhostedbythecentralhalogalaxy.
rithmonstarparticles,groupingthemwithalinkinglength
Wehaverelaxedthiscriterionanddonotapplyanypinning
of 0.05 times themean separation of theDM particles1.
of the BH particles to the minimum potential within the
In the “hr” simulation presented here, a total stel-
smoothing length.
lar mass of roughly 1010M⊙/h is needed (corresponding to
a couple of hundreds of star particles) for a halo to be ToavoidthattheBHparticlesarewanderingawayfrom
seeded with a BH particle (starting with a seed mass of the centre of galaxies by numerical effects, we take several
3.2 105M⊙/h).Inthe“uhr”simulationweareusingslightly measures, in addition to theoriginal implementation of the
× BH treatment: first, we enforce a more strict momentum
conservation within the implementation of gas accretion by
1 Note that this linking length is thus much smaller than that, forcing momentum conservation for the smooth accretion
0.15−0.20,originallyused,toidentifyvirialisedhalos. of the gas and then do not model any momentum trans-
Cosmological simulations of black hole growth 7
ferwhenswallowing gas2.Additionallyweimplementedthe arescaled logarithmically withthedifferentvaluesandnor-
conservationofmomentumandcentreofmasswhentwoBH malised to the maximum valueof each quantity.
particles are merging3. This visualises that AGN luminosity does not directly
(III) Moreover, in contrast to theoriginal implementation, trace the mass of a BH as AGN when selected by their
wehaveexplicitlyincludedadynamicalfrictionforce,which BH mass (white) seem to be a better tracer of the un-
is switched on unless the underlying simulation has a high derlying matter distribution, whereas when selected by the
enough resolution so that the cosmological simulations can Eddington-ratio(blue)orbytheirluminosity(red)theyare
numerically resolve dynamical friction reasonably well. To morelocatedinlessdenseenvironments,whatalreadyindi-
estimate thetypical friction force induced onto a BH parti- catesthepresenceofa“downsizing”trendintheirevolution.
cle we are using the following approximation of the Chan- In particular, one can see that at z = 3 many most
drasekhar formula (Chandrasekhar 1943): massiveBHsalsohavethehighestaccretionratesastheac-
cretion rate/luminosity is at that time still largely related
2
Fdf = 4π GM• ρln(Λ) err(x) 2xe−x2 ~v, (4) with BH mass (evenif with somescatter, as discussed later
− (cid:18) v (cid:19) (cid:18) − √π (cid:19)v in Fig. 12). With evolving time, this is not the case any-
where G is the gravitational constant and M• is the mass more and massive BHs accrete at a very broad range of
of the BH. The local density ρ in the vicinity of the black Eddington-ratiossothatatthoselatertimes,massiveblack
hole as well as for therelative velocity~v is calculated using holeshardlycoincidewithhighEddington-ratiosandaccre-
onlythestellarandthedarkmattercomponentsaroundthe tion rates. In addition, at low redshifts, there are several
black hole. The Coulomb logarithm is calculated as objects with high Eddington-ratios located close to voids,
these are low mass BHs (not necessarily satellites), which
Rv are growing by smooth gas accretion (as they do not have
ln(Λ)=ln (5)
(cid:18)GM•(cid:19) experienced much feedback).
Theinlayshowsazoomontotheregionwherethemost
and x=v (2)/σ, where we estimate σ as one third of the
massiveclusterformstoday.Over-plottedhereareallgalax-
maximumpcircularvelocityofthehostingsub-haloandforR ieswithstellarmassesabove1010M⊙ (whitecrosses)andall
(astypicalsizeofthesystem)weusedthehalfmassradiusof
BHs within this slice of the simulation (white diamonds).
thesub-halo,hostingtheBH.Theparametersofthehosting
Fig. 3 shows a zoom onto theredshift zero slice (left panel)
sub-halofor each BH particle are updatedevery time when
an a zoom onto the central, 2.5 Mpc wide region of the
SubFind(Dolag et al. 2009) is executedon-the-fly.
cluster.Here,thestellar populationsofthegalaxies (white)
ThiswayaBHparticleremainswithinthehostgalaxy,
within our ray tracing visualisation can be nicely seen con-
evenifit becomesasatellite ofalarger haloand,compared
trasted to the light blue ICM. Note that not all galaxies
to the original scheme, we are able to track BHs also in
containaBHastheinitial seedingisalso relatedtotheon-
satellite galaxies in cluster environments. When the BHs
going SF in a galaxy. The position of all BHs within this
are not placed artificially on the minimum of thepotential,
simulation are marked as black diamonds and nicely reflect
of course, there is no guarantee (due to numerical noise, 2
the ability of our implementation to keep the BHs at the
body scattering or when two BHs are merging) that black
centreof thesatellite galaxies (see also Fig. 1).
hole particles are stayingalways exactly at thelocal poten-
tial minimum. But due to the above handling of the dy-
namical friction, with evolving time during the simulation,
BHs sink towards the minimum potential so that typical
3 FUNDAMENTAL GALAXY AND BH
displacementsfrom thetruepotentialminimumaresmaller
PROPERTIES AND THEIR EVOLUTION
than the effective gravitational softening. Therefore, they
areordersofmagnitudesmaller thanthetypicalsmoothing Inthissection,wewilldiscusssomefundamentalproperties
radius used for estimating the parameters in the accretion of the simulated galaxies and BHs in the 500Mpc/hr run,
model or used for distributing the feedback energy. Fig. 1 as the stellar mass function, the BH mass function and the
shows explicitly a histogram for the distances of BHs from BH-stellarmassrelation inthepresent-dayUniverseandat
the potential minimum in the 500Mpc/hr run. The distri- higherredshiftsuptoz=4.Toshowtheeffectofresolution
bution is peaking at 0.7 kpc and the vast majority of BHs we additionally present the results for the 68Mpc/uhr run
aremaximal2kpcawayfromthepotentialminimumwhich down to z =1.
is reasonable for theused smoothing length of 5.2 kpc.
Fig.2showsavisualisation(ofthegasandstellarmass
density) of the 500Mpc/hr run at different redshift steps 3.1 The stellar mass function
(z = 3 0, different panels) focusing on a thin (70 Mpc
− Fig. 4 shows the stellar mass function at different redshift
thick)slice.Wehaveadditionallyindicatedthe20BHswith
steps as indicated in the legend. Simulation results (red
highestmasses (whitecircles),thehighest Eddington-ratios
solidlinesforthe500Mpc/hrrunandgreendashedlinesfor
(fedd =L/Ledd,bluecircles)andthehighestaccretionrates
the 68Mpc/uhr run) are compared with observational data
(and thus, luminosities, red circles). The sizes of the circles
from different studies(blacksymbols: P´erez-Gonz´alez et al.
2008; Bundyet al. 2005; Drory et al. 2004; Fontana et al.
2 Notethatotherwiseonewouldstatisticallyaccountforthemo- 2006; Marchesini et al. 2007;Ilbert et al. 2010; bluedashed
lines: Ilbert 2013). At z = 4, the amount of galaxies in the
mentumtransferofaccretedgastwice.
3 Note that the original scheme the merged BH have had the 500Mpc/hr run isslightly under-estimated,which isa reso-
positionandvelocityoftheBHwiththesmallerparticleID. lutioneffect:forthe68Mpc/uhrsimulationthelow-massend
8 Hirschmann et al.
Figure4.Evolutionofthestellarmassfunctioninthe500Mpc/hr(redlines)andinthe68Mpc/uhrrun(greendashedlines)comparedto
observationaldata(blacksymbols,P´erez-Gonz´alezetal.2008;Bundyetal.2005;Droryetal.2004;Fontanaetal.2006;Marchesinietal.
2007; Ilbertetal. 2010 and blue dashed lines, Ilbert2013). Atz <1the massiveend is over-estimated inthe simulations due to atoo
inefficientlyworkingradio-modefeedback.
of the stellar mass function at these high redshifts is con- massgalaxies.Thisismostlikelyaconsequenceofthether-
sistent with the observational data. Instead, down to z =1 mal energy injection scheme in the “radio-mode” adopted
the500Mpc/hrsimulationresultsprovideareasonablygood in our model (see also Puchwein & Springel 2013 who also
match with the observational data, while the low-mass end over-estimate the massive end of the stellar mass func-
in the68Mpc/uhr run is over-estimated by up to 1 dex. tion).Herewemayspeculatethatamechanical-momentum
The over-estimation of low-mass galaxies is a well- input from an AGN coupling to the ambient gas via a
known problem and most likely a consequence of our bipolar wind would be more efficient in limiting the infall
adopted model for stellar winds assuming a constant wind and accretion onto the central BH and also star forma-
velocityfortheejectedgas(e.g.Oppenheimer& Dav´e2006; tion and thus, could help making elliptical galaxies ”red
Dav´eet al. 2011). It was repeatedly shown in literature and dead” by suppressing late star formation (Choi et al.
that energy- or momentum-driven wind models can sig- 2014). Such mechanisms are for example investigated by
nificantly reduce the baryon conversion efficiencies (e.g. Choi et al.(2012);Debuhret al.(2012);Barai et al.(2013);
Hirschmann et al. 2013 and references therein) and thus, Duboiset al.(2012)employingsimulationsofisolatedgalax-
also the low-mass end of the stellar mass function result- ies or galaxy mergers.
ing in an improved match with the observational data (e.g.
Dav´eet al.2013;Puchwein & Springel2013).Inthesestud-
3.2 The BH mass function
ies, energy-driven winds, for example, are shown to reduce
the low mass end of the stellar mass function by at least In Fig. 5, the present-day BH mass function of the
1dexcomparedtoaconstantwindmodelwhichcanaccount 500Mpc/hr run (red line) is compared to observations
forthediscrepancybetweenour68Mpc/uhrsimulationpre- from Marconi et al. (2004); Shankaret al. (2004) and
dictionsandtheobservationalconstraints.Inaddition,mod- Shankaret al. (2009) (black lines and symbols and grey
els including “early” stellar feedback (Stinson et al. 2013; shaded areas) and is found to be in reasonably good agree-
Kannan et al. 2014) or radiation pressure (Hopkinset al. ment for a BH mass range of 7.5 < log(M•/M⊙) < 9.5.
2013)seemtobeparticularlyefficientindelayingstarforma- Below, the simulation is under-predicting the amount of
tion in low mass halos towards later times (i.e. in breaking low massiveBHsbyalmost oneorderof magnitude.This is
the hierarchical formation of galaxies) and thus, predicting mainlyrelatedtoacombinationoftoolowresolution(main
low baryon conversion efficiencies in these halos down to effect, see dashed coloured lines indicating the 68Mpc/uhr
z=0 – consistent with observational constraints. predictions) and BH seed masses.
Turning towards lower redshifts (z < 1), the massive At the high mass end, the amount of massive BHs
end of the stellar mass function (log(Mstellar/M⊙) > 11) log(M•/M⊙) > 9.5 is significantly over-estimated by up
is significantly over-estimated in the 500Mpc/hr run, at to 2 dex. The high mass end is mostly influenced by the
z = 0 by more than one order of magnitude for 1012M⊙- choice for the parameter regulating the strength for AGN
Cosmological simulations of black hole growth 9
Figure 5.EvolutionoftheBHmassfunctioninthe500Mpc/hr
simulation (solid coloured lines) and in the 68Mpc/uhr simula-
tion (dashed coloured lines). The present-day BH mass function
(redline)iscomparedtoobservationsfromMarconietal.(2004);
Shankaretal. (2004) and Shankaretal. (2009) (black lines and
symbolswiththegreyshadedareas).Wefindareasonableagree-
ment between observations and simulations for BH masses be-
tween5×107M⊙<M•<3×109M⊙,whileabove, theamount
ofmassiveBHs(>3×109M⊙)isover-estimatedbyupto1dex
(asaconsequence oftooinefficientradio-modefeedback).
feedback,buthardlydependenton theefficiency parameter
for regulating the radio-mode feedback. Although the lat-
Figure 6.Upper panel: present-dayrelationbetween blackhole
teris supposedtoregulatethelate timestarformation and andstellarmass.Colour-codedcontoursshowthesimulatednum-
BH growth, in the implementation of the thermal injection ber density, the red lineillustrates the mean for the 500Mpc/hr
scheme in the “radio-mode” adopted in our model, is still run. The simulated relation perfectly matches the observational
not efficient enough to lower the high mass end of the BH data(symbolsillustratetheirdatasetofH¨aring&Rix(2004),the
mass function (as discussed for the stellar mass function). blacksolidlineshowsthefittothedataandtheblackdashedlines
Therefore,AGNfeedbackathighredshiftismuchmoreeffi- the corresponding 1-σ scatter). Bottom panel: Evolution of the
blackhole-stellarmassrelationinthe500Mpc/hr(solidlines)and
cientinsuppressingstarformationalsoatlatetimesbypre-
the68Mpc/uhrrun(dashedlines).Forcomparison,theobserved
heatingtheenvironment.Insteadofapurelythermalenergy
present-day relation is also illustrated as in the upper panel. In
injection, we speculate that a mechanical-momentum input
thesimulations,theBH-stellarmassrelationisinplaceatz=3
fromanAGNcouplingtotheambientgasviaabipolarwind
andhardlyevolvesafterwards.
would be more efficient in limiting the infall and accretion
ontothecentralBH(e.g.Choi et al.2012;Barai et al.2013)
andthus,reducingthemassiveendoftheBHmassfunction.
Turning towards higher redshifts (see coloured solid withtheincreasedlow-massendofthestellarmassfunction
lines in Fig. 5 indicating the 500Mpc/hr run), the num- in the 68Mpc/uhr run due to inefficient stellar feedback.
ber density at a given BH mass shows a very rapid evo- Above log(M•/M⊙) > 8, however, the 68Mpc/uhr predic-
lution until z = 1. At that redshift, the high mass end tions roughly converge against the lower resolution results
is mostly in place and only the amount of BHs below of the500Mpc/hr run.
log(M•/M⊙) < 8 further increases until z = 0. This im- Compared to observational constraints for the evolu-
plies that BHs grow most strongly untilz=1 which is also tion of the BH mass function (derived from integrating the
in agreement when considering the cosmic evolution of the continuity equation, e.g. Merloni & Heinz 2008), our pre-
total BH accretion rates, which peak around z =2 and fall dictions are broadly consistent: also in theobservations the
below 0.1 M⊙ yr−1 Mpc−3 after z=1 (see also Fig. 7). massive end of the BH mass function is hardly evolving
To demonstrate the effect of resolution, we also show after z = 1, while the low-mass end is still increasing.
theBH mass functions at z=1 4 for the 68Mpc/uhr run This trend already reflects the anti-hierarchical behaviour
−
(dashed coloured lines). The low mass end (log(M•/M⊙)< we will discuss in the course of this study in more detail.
8) is significantly increased byupto1dexcompared tothe In addition, at z = 2, the observed BH mass function is
500Mpc/hr runasaconsequenceofafastermass growth of peaking at log(M•/M⊙) 8 and has a number density of
∼
low mass BHs in the higher resolution run. This is in line log(Φ) 3 for log(M•/M⊙)=7 and a number density of
∼−
10 Hirschmann et al.
Between z = 1 and z = 3, the BH-to-stellar mass ra-
tioisslightlyincreasingforagivenstellarmass,i.e.BHsare
growingslightlyfasterthanthecorrespondinggalaxystellar
masses. In contrast, afterwards (z <1), thisratio decreases
untiltoday,i.e.therelationisjustshiftedtowardsmoremas-
sive galaxy masses, what is most likely due to the strong
lategrowth of massivegalaxies. Overall, thisismainly con-
sistent with previous results from cosmological simulation
(DiMatteo et al. 2008).
3.4 The SFR and BH accretion rate densities
Observationsrevealthatstarformation rate(SFR)andBH
accretion rate densities (ρ˙stellar and ρ˙•, respectively) trace
eachotherovercosmictimewithρ˙stellar 2 103 ρ˙• (e.g.
∼ × ×
Zhenget al.2009andreferencestherein).InFig.7weshow
the predicted cosmic evolution of the SFR (orange lines)
andBHaccretionratedensities(redlines,multipliedwitha
factorof2 103)forthe500Mpc/hr(solidlines)andforthe
×
68Mpc/uhr run (dashed lines). At high redshifts z > 1.5,
Figure7.Cosmicevolutionofmeanstarformation(orangelines)
both the SFR and BH accretion rate densities are larger
andBHaccretionratedensitiesmultipliedbyafactorof2×103
in the 68Mpc/uhr run than in the 500Mpc/hr one. This
(red lines) in the 500Mpc/hr (solid lines) and the 68Mpc/uhr
explains the increased low mass end of the stellar and BH
run (dashed lines). Both, star formation and BH accretion rate
densities peak between z ≈ 1−2 followed by a decline towards mass function at high redshifts in the 68Mpc/hr run (see
lowerredshifts.Comparedtoobservedcosmicstarformationrate Figs. 4and5),wherestarformation and gasaccretion onto
densities derived from different wavebands (coloured symbols, BHs is more efficient.
Hopkins&Beacom 2006), simulations predict too high star for- In both runs, SFR and BH accretion rates densities
mation rates at z <1resulting intoo massivegalaxies (see Fig. peak between z = 1 2 followed by a decline below
4). 0.1 M⊙ yr−1 Mpc−3 aft−er z = 1 (for the 500Mpc/hr run).
This is qualitatively consistent with the observational com-
pilation for the star formation rate densities derived from
different wavebands (coloured squares: Hopkins& Beacom
log(Φ) 5for log(M•/M⊙)=9.Thisisalso predictedby
∼− 2006). At high redshifts z > 4, however, simulated BH ac-
the68Mpc/uhr run.
cretion ratedensities aretoo low compared totheobserved
starformation rate densities,which is aconsequenceofres-
olution (see improvement for the 68Mpc/uhr run shown by
3.3 The BH-bulge mass relation
the dashed red lines). Between z = 2 4, both the BH
−
ThetoppanelofFig.6illustratesanotherfundamentalprop- accretion and SFR are slightly higher compared to the ob-
ertyofBHs,theBH-stellarmassrelation inthepresent-day servational data, particularly in the higher resolution run.
Universe (colour-coded 2-dimensional histogram, the mean Below z = 1, the BH accretion rate densities are in good
relationisindicatedbytheredline).Forsimplicity,wehave agreement with the observed data (for the SFR), whereas
used the total stellar mass and not only the bulge mass, thesimulated SFRdensitiesaretoohigh(byafactorof3-4
which is considered in observations. However, in our simu- at z =0) – a consequence of the too inefficient radio-mode
lations,galaxieswithmassesabovelog(Mstellar/M⊙)>10.5 feedback implementation resulting in too massive present-
arelargelysystemsonlyconsistingofasphericalcomponent day galaxies, as discussed earlier in Fig. 4.
as the resolution is not high enough to provide sufficient
morphologicalinformationofthegalaxies.Comparedtoob-
servations from H¨aring & Rix (2004) (black lines and sym- 4 THE AGN LUMINOSITY FUNCTION
bols), we find an excellent agreement with the simulations.
Thisisadirectconsequenceofthechoiceofthefeedbackeffi- Inthissection,wefocusonbothradiativelyefficientandin-
ciency,whichincounter-playtothecooling,setstheselfreg- efficientAGNinthe500Mpc/hrandthe68Mpc/uhrrunby
ulated state of the late time evolution (see Churazov et al. comparingtheirbolometricluminositiesatredshiftsbetween
2005). z = 0 5 as well as their present-day radio luminosities to
−
In the bottom panel of Fig. 6 we show the evolution of observationaldata.Foranimprovedcomparisonwithobser-
theBH-stellarmassrelationatz=0 4forthe500Mpc/hr vations,wealso derivethesoft andhard X-rayluminosities
(coloured solid lines) and the 68M−pc/uhr run (coloured andadoptadustobscurationmodel(onatoruslevel)being
dashed lines). We find no significant effect of resolution on dependenton both redshift and AGN luminosity.
the relation and thus, in both runs, the relation is in place
at z=3 (see turquoise lines). Below z=3, theslope of the
4.1 Bolometric luminosities
relationissimilartotheoneofthepresent-dayrelation,but
the BH masses are mainly found to be under-massive for a The bolometric luminosity Lbol (=radiative luminosity) is
given stellar mass. calculated following the study of Churazov et al. (2005),