On the Dynamics of Supermassive Black Holes in Gas-Rich, Star-Forming Galaxies: the Case for Nuclear Star Cluster Coevolution

2.4 MB·

by  Pawel Biernacki

MNRAS000,000–000(0000) Preprint20January2017 CompiledusingMNRASLATEXstylefilev3.0
On the Dynamics of Supermassive Black Holes in Gas-Rich,
Star-Forming Galaxies: the Case for Nuclear Star Cluster
Coevolution
Pawel Biernacki(cid:63), Romain Teyssier and Andreas Bleuler
CenterforTheoreticalAstrophysicsandCosmology,InstituteforComputationalScience,UniversityofZurich,
Winterthurerstrasse190,8057Zurich,Switzerland
7
1
0
2
20January2017
n
a
J
ABSTRACT
8
We introduce a new model for the formation and evolution of supermassive black holes
1
(SMBHs) in the ramses code using sink particles, improving over previous work the treat-

ment of gas accretion and dynamical evolution. This new model is tested against a suite of
]
A high-resolutionsimulationsofanisolated,gas-rich,coolinghalo.Westudytheeffectofvari-
G ousfeedbackmodelsontheSMBHgrowthanditsdynamicswithinthegalaxy.
Inrunswithoutanyfeedback,theSMBHistrappedwithinamassivebulgeandisthere-
.
h foreabletogrowquickly,butonlyiftheseedmassischosenlargerthantheminimumJeans
p mass resolved by the simulation. We demonstrate that, in the absence of supernovae (SN)
- feedback,themaximumSMBHmassisreachedwhenActiveGalacticNucleus(AGN)heat-
o
ingbalancesgascoolinginthenuclearregion.
r
t WhenourefficientSNfeedbackisincluded,itcompletelypreventsbulgeformation,so
s
a that massive gas clumps can perturb the SMBH orbit, and reduce the accretion rate signifi-
[ cantly.Toovercomethisissue,weproposeanobservationallymotivatedmodelforthejoint

evolutionoftheSMBHandaparentnuclearstarcluster(NSC),whichallowstheSMBHto
1
remaininthenuclearregion,growfastandresistexternalperturbations.Inthisscenario,how-
v
ever, SN feedback controls the gas supply and the maximum SMBH mass now depends on
0
9 thebalancebetweenAGNheatingandgravity.WeconcludethatSMBH/NSCco-evolutionis
1 crucialforthegrowthofSMBHinhigh-zgalaxies,theprogenitorsofmassiveellipticaltoday.
5
Keywords: methods:numerical-galaxies:evolution-galaxies:active-quasars:supermas-
0
siveblackholes-galaxies:starclusters:general
.
1
0
7
1
: 1 INTRODUCTION 2007; Nandra et al. 2007; Fabian 2012; Yesuf et al. 2014; Che-
v
ungetal.2016).TheformationoftheSMBHthemselvesremains
i Supermassive Black Holes (SMBH) are found in the central re-
X amystery.Twomainscenariosareconsideredleadingtomassive
gionofmassivegalaxiesatallredshifts,mostlyintheformofAc-
r enoughSMBH:1)directcollapseofmassiveclumpsofpristinegas
a tive Galactic Nuclei (AGN). There is accumulating evidence that (Loeb&Rasio1994;Bromm&Loeb2003)or2)mergersofstellar
SMBHaretightlylinkedtotheevolutionoftheirhostgalaxy(Rich-
remnantsindensestellarclusters(Quinlan&Shapiro1990;Porte-
stoneetal.1998;Ferrarese&Merritt2000;Gebhardtetal.2000;
giesZwartetal.1999;Devecchi&Volonteri2009),eachscenario
Marconi&Hunt2003;Ha¨ring&Rix2004;Kormendy&Ho2013),
havingclearstrengthsandweaknesses,asexplainedinthereviews
puttingAGNphysicsatthecentreofourunderstandingofgalaxy
ofBegelmanetal.(2006)andVolonteri(2010).
evolution.ThestrongcorrelationofSMBHmassesandstellarve-
locitydispersion,forexample,suggestsapossibleco-evolutionof Motivatedbytheseobservationalhints,theoreticalmodelsof
thecentralSMBHanditshostgalaxy(Magorrianetal.1998;Laor SMBH growth and their associated feedback (mostly based on
2001; McLure & Dunlop 2002;Ha¨ring & Rix 2004). AGNfeed- complexnumericalsimulations)becameinrecentyearsmoreand
back is also often invoked as one of the possible origins of the moresophisticated,withmixedsuccesseswhencomparedagainst
quenchingofstarformationinellipticalgalaxies(Schawinskietal. observational data (Springel et al. 2005; Di Matteo et al. 2005;
Boweretal.2006;Crotonetal.2006;Hopkins&Hernquist2006;
Ciottietal.2010;Teyssieretal.2011;Duboisetal.2011).AGN
(cid:63) [email protected] feedbackintheoreticalmodelsofgalaxyformationhasprovenvery
(cid:13)c 0000TheAuthors
2 P.Biernackietal.
efficient at regulating the Star Formation Rate (SFR) in massive, automaticallyextractedfromthesimulationusingtherecentlyde-
red and dead galaxies, but the X-ray properties of the intergalac- velopedclumpfinderonboardtheramsescode(Bleuler&Teyssier
tic gas are very difficult to reproduce. One natural explanation to 2014).Wealsoimprovedthedynamicalintegratorofthesinkparti-
the difficulties of these models is the formidable range of scales cle,allowingustoperformdetaileddynamicalstudies.Finally,we
one has to capture, in order to resolve numerically the entire ac- addedtwo newingredientsto themodel, namelyafully momen-
cretionflow fromparsecscalestowards thelaststableorbit (typ- tum conserving drag force and a model for SMBH and Nuclear
ically 10−5 pc). Numerical implementation of SMBH formation, Star Cluster (NSC) co-evolution. Our goal is to apply these var-
theiraccretionflowsandassociatedenergeticoutflows,havetorely ious ingredients to model simulations featuring a cooling, Milky
onstrongapproximations,usuallyreferredtoas“subgridmodels”. Way-sizedhalo(SeeSection3),leadingtotheformationofagas-
Note that the same technique is applied to star formation recipes rich,clumpyandviolentlyturbulentdisc,reminiscentofthehigh-
in galaxy formation simulations, making the whole endeavour of redshift galaxies population detected in deep Hubble Space Tele-
simulatinggalaxiesverychallenging. scopeimages.InSection4,weoutlinethefactthatSMBHdynam-
Astheresolutionofgalaxyformationsimulationsisincreas- ics in this turbulent environment is extremely chaotic, leading to
ing,fromthousandsofpcinlarge-scalecosmologicalsimulations theejectionoftheSMBHfromthecentralregionofthegalaxy,un-
(Duboisetal.2014;Vogelsbergeretal.2014;Schayeetal.2015; lessoneconsidersveryspecificdynamicalmodels.Realisticstellar
Dubois et al. 2016), to hundreds of pc in cosmological zoom-in andAGNfeedbackmodelsmakethesituationevenmorecritical.
simulations of galaxy formation (e.g. Kim et al. 2011; Angle´s- InSection5,wefinallydiscussamodelwhereSMBHsareeither
Alca´zaretal.2014;Duboisetal.2015),ultimatelyreachingafew hostedandprotectedbyaparentNSC,ormassiveenoughtosus-
pcinisolateddiscssimulations(Gabor&Bournaud2013;Hopkins taintheviolentperturbationsfromtheirhostgalaxy.InSection6,
et al. 2014), these subgrid models need to be tuned and adapted wediscussvariousobservationalargumentsinfavourofthisnew
totheincreasinglybetterresolvedinterstellarmedium(ISM)struc- scenario.
ture,withanincreasinglystrongersupersonicturbulence.
The goal of this paper is precisely to study such a model of
SMBHformation,growthandfeedbackinhighlyresolved,turbu-
lent and clumpy galactic discs, typical of high redshift, gas-rich
galaxies(Elmegreenetal.2008a;Dekeletal.2009;Bournaudetal.
2 ANEWMODELFORSMBHFORMATIONAND
2012).ThisenvironmentisparticularlyrelevanttoSMBHphysics,
EVOLUTION
astheseclumpydiscsarebelievedtobetheprogenitorsofthegiant
ellipticalshostingthemostmassiveSMBHsinourpresentepoch The first generation of SMBH models was developed in the con-
(Kormendy&Ho2013;McConnell&Ma2013). textofcosmologicalsimulations,withresolutionaround1kpcor
NumericalmodelsofSMBHformationandevolutionareall more(Bellovaryetal.2010;Vogelsbergeretal.2014;Schayeetal.
basedontheso-called“sinkparticle”technique.TheSMBHisrep- 2015;Duboisetal.2016)orforrelativelysmoothgalaxymodels,
resented by a point mass, moving through the fluid and interact- using either a pressurised ISM equation of state (Truelove et al.
ingwithitthroughaccretionandejectionofmass,energyandmo- 1997;vandeVoortetal.2011)oralowgasfractionrelevantfor
mentum. Sink particles were first implemented in simulations of low-redshiftgalaxyevolution.Thesinkparticlewasnotallowedto
star-forming turbulent molecular clouds (Bate et al. 1995), using moveawayfromthegalaxycentre,byeitherforcingittoremain
aSmoothedParticleHydrodynamics(SPH)code.Krumholzetal. close to the gravitational potential minimum, or by using various
(2004)wasthefirstonetoproposeasinkparticleimplementation drag forces (Springel et al. 2005; Okamoto et al. 2008; Gabor &
for grid-based codes, using Adaptive Mesh Refinement (AMR). Bournaud2013).ThenextgenerationofSMBHmodelsneedtobe
The sink particle technique was then adapted to the SMBH for- abletoresolvetheSMBHdynamicswithinthegalaxy,andmore
mationandevolution,hereagainfirstinSPHcodes(Springeletal. importantly,tofollowitsevolutionwithinhighlyturbulent,gas-rich
2005;DiMatteoetal.2005)andthenlaterinAMRcodes(Dubois environments typical of galaxy evolution at high redshift. In this
et al. 2010; Kim et al. 2011). The key ingredients in our SMBH section,wepresentthenew-generationSMBHmodelimplemented
formationandevolutionmodelsarethefollowings:a)theforma- intheramsescode.Itisheavilybasedontheoldmodelpresentedin
tionoftheSMBHparticleandinparticularthechoicetotheinitial Duboisetal.(2010)andTeyssieretal.(2011),andcapitalisesover
seedmass(e.g.Begelmanetal.2006;Volonteri2010),b)thedy- the new sink particle implementation we have developed within
namicsoftheSMBHparticle,withthepossibleinclusionofadrag the context of star-forming molecular clouds (Bleuler & Teyssier
force(seetherecentworkofTremmeletal.2015),c)thegrowthof 2014).
theSMBHparticlemassasafunctionoftime,withtwofundamen- AlthoughwemodelSMBHsascollisionlessparticles,wedo
talingredientsbeingtheBondi-Hoyle-Lyttleton(Hoyle&Lyttleton notusetheParticleMeshsolverdesignedforthedarkmattercom-
1939;Bondi&Hoyle1944;Bondi1952)accretionrate,limitedto ponent.Instead,weplacearoundeachsinkasphericaluniformdis-
theEddingtonaccretionrate(forobservationalconstrainsseee.g. tributionoftestparticles(wecallthem“cloudparticles”)ofradius
Kollmeieretal.2006;Steinhardt&Elvis2010),andfinally,d)the r = 4∆x , where ∆x is the size of a cell at the highest
sink min min
feedbackfromtheSMBHparticlethataffectsthesurroundinggas refinement level. These cloud particles are evenly spaced within
(Ostrikeretal.2010;Angle´s-Alca´zaretal.2013;Choietal.2014, thesphere(withroughly8cloudparticlespergridcell)andfollow
2015; Nayakshin 2014; Costa et al. 2014), and therefore couples thesinkparticleasarigidbody.Thesecloudparticlesareusedto
backtoallthepreviousingredientsofthemodel. probethegasdistributionaroundthesinkandtodistributetheac-
Inthiswork,wepresentanewimplementationoftheSMBH cretionandtheejectionofmass,momentumandenergy.Notethat
formationandevolutionmodelintheramsesAMRcode(Teyssier the value for the sink sphere radius can be modified by the user,
2002),inheritedfromtheearlierworkofDuboisetal.(2010)and withrecommendedvaluesrangingfrom1to4∆xmin.
Teyssieretal.(2011),butsignificantlyimprovedinmanyaspects Inthefollowingsubsectionswegivemoredetailsontheim-
(seeSection2).Forexample,oursinkparticleformationsitesare provementsofourSMBHsinkparticleimplementation.
MNRAS000,000–000(0000)
SMBHsdynamicsingas-richgalaxies 3
2.1 SMBHformation follows
GM
ThelifeoftheSMBHinoursimulationsbeginswiththeformation r = sink, (4)
ofthesinkparticle.Itisaproblemwhichdeservesitsowncareful Bondi v2
Bondi
consideration,butherewereduceittotheidentificationofapossi- (cid:113)
v = c2+v2 , (5)
bleformationsiteandtothechoiceoftheinitialmass Mseed.The Bondi s rel
twomainscenariosforSMBHformationare1)directgascollapse
wherecs isthelocalsoundspeedofthegasandvrel istherelative
or2)formationthroughstellarremnantscollisionsinadensestellar
velocitybetweenthesinkvelocityv andthegasaveragevelocity
sink
system.Inbothcase,SMBHformationisassociatedtoexception-
withinthesinkspherev¯
allydenseregions,probablyatveryhighredshift,withproperties
leadingfirsttotheformationofanintermediatemassblackhole, vrel=vsink−v¯ (6)
whichwillaccretegasandgrowevenmoreintotheSMBHregime.
Onecandefinethefree-fallvelocityontothesinkparticleas
Modellingtheseprocessesisclearlyoutofthescopeofthis
(cid:115)
paper,asitwouldrequiremuchhigherresolutionandtheaddition
GM
ofphysicalingredientsthatareabsentfromoursimulations,orthat vff,sink= r sink (7)
arenotevenreallyunderstoodtoday.Wethereforedirectlycreate sink
ourfirstandonlySMBHwhenthefirstdenseclumpofgasforms. Thedimensionlessradiuscanbewrittenas
Thisallowsthesinkparticletoevolveinadenseenvironment,mim-
ickingtheearlyphaseofSMBHgrowth.Forthis,weuseourbuilt- xsink=v2Bondi/v2ff,sink (8)
inclumpfinderphew(Bleuleretal.2015)andformthesinkparticle
andobviouslyindicateswhethertheaccretionflowaroundthesink
inthemostmassivegasclumpatachosentime(seeSection3).It
issupersonicforx <1orsubsonicforx >1.
sink sink
isworthemphasisingthatinthisformationscenarioseedSMBHis
Inthestrongsupersonicregimewhere x (cid:28)1,thedimen-
sink
trappedinnucleargasclump;iftheSNfeedbackisincluded,then
sionlessdensityprofileoftheBondisolutionasymptotestoα(x)(cid:39)
theinitialhostclumpisquicklydestroyed. x−3/2.Onecanre-writetheaccretionrateinthestrongsupersonic
Thevalueoftheinitialseedmassisratherarbitrary.Atypical
limitas
dvaylnuaemoifcaMlsseimedu=la1ti0o5nsM(e(cid:12).,g.isBuosoutahll&ySadchoapyteed2i0n09la)r.gDe-irseccatlecohlyladprsoe- M˙ (cid:39)4πρ¯r3/2 (cid:112)GM =3 Mgas (9)
scenariosofSMBHformationdopredictseedmassesofthismag- Bondi sink sink tff,sink
nitude(e.g.Begelmanetal.2006).Inthispaper,weprefertoadopt wherethesinkfree-falltimeisdefinedas
amorepragmaticapproachandconsidertheseedmassasafreepa- (cid:115)
rameter.TheBondiaccretionmodelwedescribeinthenextsection r3 r
isbasedonthestrongassumptionthatthesinkparticlegravityfield tff,sink= GMsinskink = vffs,isninkk (10)
dominatesoverthegasself-gravity.Aminimumseedmassequalto
andtheavailablegasmasswithinthesinksphereradiusis
thesimulationminimumJeansmassappearstobetherightchoice,
asournumericalexperimentsinSection4indicate. 4π
Mgas= 3 ρ¯rs3ink (11)
Oneconcludesthatinthestrongsupersoniclimit,theaccretionrate
does not depends on the gas properties anymore, but only on the
2.2 SMBHaccretion
availablegasmassandthesinkfree-falltime.Thiscorrespondsto
OncetheSMBHhasformed,itcontinuestogrowinmassviaaccre- amaximumphysicallymotivatedaccretionrateontothesink.
tionofgasfromitssurroundings.Themostpopularandphysically Inthestrongsubsoniclimit,wherex (cid:29)1,onehasα(x)(cid:39)1,
sink
motivatedapproachtocomputetheaccretionrateontotheSMBH andtheaccretionratecanbewrittenas
particleistousetheBondi-Hoyle-Lyttletonformulae(laterBondi
for short; Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi M˙Bondi(cid:39)4πρ¯rB2ondivBondi (12)
1952):
Thisisthisformulathatisusedinmostsinkparticleimplementa-
M˙Bondi=4πρ∞rB2ondivBondi, (1) ttihoins,laasntdfworemwuolaueldisliokneltyovsatrleidssi,nasthiensKurbusmonhioclzreegtiaml.e(,2w00h4e)r,eththaet
where Bondi radius is much smaller than the sink radius. Manipulating
slightly the previous equation, one can rewrite the accretion rate
ρ¯
ρ∞= (2) formulaeas
α(x )
sink M˙ (cid:39)3 Mgas 1 (13)
αisthedimensionlessdensityprofileoftheBondiself-similarso- Bondi tff,sink x3/2
lution(seee.g.Chapter6ofShu1992),ρ¯istheaveragegasdensity sink
withinthesinksphere,and This shows explicitly that the subsonic accretion rate is much
smaller than the supersonic one. This also shows that if the sim-
xsink=rsink/rBondi (3) ulationtimestepiscontrolledbythesinkCourantcondition
isthedimensionlessradiusevaluatedatthesinksphereradius.This ∆t(cid:54) tff,sink (14)
functionα,firstintroducedbyKrumholzetal.(2004)intabulated 3
form,isacrucialingredienttodescribetheaccretionflow,andisof- onecannotremovemorethantheavailablegasmasswithinthesink
tenmissinginmanysinkparticlealgorithmimplementations.The sphereinonetimestep.
Bondi radius r and the Bondi velocity v are defined as It has been proposed by Springel et al. (2005) and Booth &
Bondi Bondi
MNRAS000,000–000(0000)
4 P.Biernackietal.
Schaye(2009)toboostthepreviousaccretionrateformula,toac- oftheaccretedgaswithinthesinkspheredoesnotcoincidewiththe
countforunresolveddensityandtemperaturefluctuationsatscales centreofthespherex .
sink
lowerthanthecellsize.Inthispaper,wefollowthesameidea,al-
lowingthesoundspeedofthegastobereduced,owingtosmaller
unresolvedtemperaturefluctuations.Thisboilsdowntoreplacing
inthepreviousformulaethesoundspeedby
2.3 SMBHdynamics
cs→cs/βboost(ρ¯) (15)
The next fundamental requirement of our sink particle algorithm
wheretheboostfactorisdefinedas is to model properly the dynamics of the SMBH. The sink parti-
cletrajectoryfollowsfromthedynamicalevolutionofapointmass
βboost(ρ)=max[(ρ/ρ∗)2/3,1.0], (16) particle,subjecttothegravitationalforceofthegas,starsanddark
Incaseofzerorelativevelocity,thisformulacorrespondsexactlyto matterparticles,andalsosubjecttoadragforceduetoatightcou-
themodelproposedbyBoothandSchaye.Wewouldliketostress plingbetweentheaccretedgasandthesink.
thattheonlyeffectofthisboostistochangethetransitionfromsu- Notethatthelatterhasbeenofteninvokedintheliteratureto
personictosubsonicaccretion,butthestrongsupersonicaccretion justify why one could artificially locked the sink particle coordi-
rate will not be modified from its maximally physically allowed natestotheminimumofthepotentialwell(e.g.Sijackietal.2007;
valuederivedabove.Wewouldlikealsotostressthatonecannot Costaetal.2014),orartificiallypushedinthedirectionofthehalo
modifiedtherelativevelocityv fromphysicalgrounds.Sinkpar- centre(Gabor&Bournaud2013).Thereisnophysicalmotivation
rel
ticles with very high relative velocities are therefore likely to ac- forthesemodels.LowermassSMBHcanbeexpectedtogetscat-
creteverylittlemass,astheyshould.Reducingtherelativevelocity teredbymassivegasclumps(e.g.Gabor&Bournaud2013).Other
artificiallyhasbeenalsousedinthepasttoboosttheaccretionrate, physicallymotivatedmodelsdoexistintheliterature,thatcanhelp
withoutanyphysicalmotivation. preventingthesinkparticlefromwanderingaroundthegalaxy.For
An important ingredient specific to SMBH accretion is the example, Tremmel et al. (2015) proposed to estimate the amount
maximal allowed accretion rate onto the black hole, namely the ofdynamicalfrictionthatismissingduetopoorresolution,which
Eddingtonrate, consistsinasub-gridmodelforadragforcebetweenthesinkand
the collisionless component. Similar sub-grid model can be con-
M˙ = 4πGMsinkmp = Msink (17) structedforthepotentiallymissingdragforcebetweenthesinkand
Edd (cid:15)σ c t
r T S thesurroundinggasmedium(Chandrasekhar1943;Ostriker1999;
wherempistheprotonmass,σTistheThomsoncrosssectionand Chaponetal.2013).Wewillproposehereanotherphysicallymo-
(cid:15)ristheShakura&Sunyaev(1973)radiativeefficiencyforaSMBH tivatedmodelbasedontheEddingtonlimitedaccretion.
accretion;(cid:15)r=0.1.TheseconstantsarecombinedintotheSalpeter First, the gravitational interaction between the sink and the
time,astS(cid:39)45Myr.Finally,theaccretionrateontotheSMBHis matterdistribution,aswellasbetweenthesinkandpossibleother
computedusing sinksinthecomputationalbox,arebothtreatedusingadirectsum-
M˙acc=min(M˙Bondi,M˙Edd). (18) mationmethodofasoftened1/r2Newtonianacceleration.Wepre-
fer this new approach than using the Particle Mesh method, as it
We would like to stress that the Eddington rate comes from the givesmoreaccuratetrajectories,especiallyiftheSMBHmassdom-
followingpicture:gasisaccretedusingtheBondiratetowardsthe inates the local potential. The softening radius used in the force
SMBH accretion disc, and the accretion energy is converted into calculationsissetto2∆x ,asinBleuler&Teyssier(2014).
min
accretion luminosity, which in turn will remove the fully ionised Whenthesinkaccretesgasfromwithinthesinksphere,italso
gasinthevicinityoftheSMBH,ifitexceedstheEddingtonlumi- accretesthecorrespondingmomentum,whichtranslatesintoanef-
nosity.Sinceouraccretionmodelisappliedtoverylargegalactic fectivedragforcebetweenthegasandthesink.Whentheaccretion
ofISMscales(saybetween10pcto1000pc),wedonotresolve rateontotheSMBHisEddington-limited,thesituationishowever
theregionwhereradiationpressurewillremovethegasandcontrol morecomplicated.Asdescribedbefore,theEddingtonlimitforthe
the accretion onto the SMBH. Eddington limited accretion there- radiationisenforcedinthevicinityoftheSMBH,wherethegasis
foremeansthatgasisaccretedattheBondirate,andthendecreted fullyionisedandhasreachedtheSMBHaccretiondisc.Wecon-
ataslightlysmallerrate,thenetbudgetbeingthe(small)Eddington sider in this paper that the gas accretion rate towards the SMBH
rate.Thispictureisquitedifferentfromwhatisconsideredusually accretiondiscissetbytheBondiformula,andcorrespondstothe
and will be used later in the paper to introduce an additional gas large scale flow, while the gas accretion rate onto the SMBH is
dragforceonthesinkparticle. setbytheEddingtonlimit.Thedifferencebetweenthetworates,
We discuss finally one important technical detail: once we Bondi minus Eddington, corresponds to gas being decreted from
knowthesinkparticle’scurrentaccretionrate,weremovegasfrom the accretion disc region and redistributed on large scale, in our
thesinkspherebyintegratingthepreviousaccretionrateoverthe casewithinthesinksphere.
timestep.
∆Mgas=−M˙acc∆t (19) M˙dec=M˙Bondi−M˙acc (21)
In order to avoid emptying very low density gas cells in the sink This process of accretion and ejection will lead to an additional
sphere,weremovefromeachcell(labelledi)thefollowingmass- exchange of momentum between the gas and the sink, hence an
weightedcontribution, additionaldragforce.
∆M Wemodelthisadditionaldragforcebyrequiringthatthecen-
∆ρi=−ρi acc (20) treofmassofthejointgas+sinksystemremainfixedduringthe
M
gas accretion,andthatitstotalmomentumisconserved.Ifwenotethe
Animportantconsequenceofthisstrategyisthatthecentreofmass gascentreofmasswithinthesinksphereasxgas,thistranslatesinto
MNRAS000,000–000(0000)
SMBHsdynamicsingas-richgalaxies 5
ashiftinthesinkcoordinatesgivenby First,wehavethecoldaccretionregime,forwhichcoolingdom-
inates over heating. The Bondi accretion rate is so high that we
Mgasdxdgtas = M˙decxsink−M˙Bondixgas, considertheaccretiontobeEddingtonlimited,
Msinkdxdstink = M˙Bondixgas−M˙decxsink, (22) M˙acc= MtsSink. (27)
andasimilarmomentumtransferbetweenthesinkandthegas(in WeconsiderforthecoolingfunctiononlyBremsstrahlungsothat
otherwordsadragforce)givenby: Λ(T)=Λ T1/2 (28)
0
Mgasdvdgtas = M˙decvsink−M˙Bondivgas, wmhaetiroenΛfo0r(cid:39)h1ig.2h×te1m0p−e27raeturgres−an1dcmlo3wKm−e0.t5a.llTichiitsyigsaas.gWooedcaopnpcrlouxdie-
Msinkdvdsitnk = M˙Bondivgas−M˙decvsink, (23) ismphmereed,iactoeolylinthgawt,ilfloarlawagyivsewnianvoevraegrehegaatsindge,nasnitdytwheitshiinnkthweilslirnek-
Theseequationsaresolvedforeachtimestep,andareusedtomod- maininthecoldaccretionregime,unlesstheSMBHmassbecomes
ifythesinkpositionandvelocity,butalsothegasdensity,momen- largeenough,sothat
tum and total energy within the sink sphere. More details on the t
M >n2Λ T1/2 S V . (29)
numericalimplementationaregivenintheAppendix.Notethatin sink H 0 (cid:15) (cid:15)c2 sink
c r
case of zero decreted mass (pure unlimited Bondi accretion), the
Becausethesinkisnowmassiveenough,heatingdominatesover
momentum transfer only comes from the accreted gas mass onto
cooling,andthesinksphereentersthesecondphase,namelythehot
thesink,asitshould.Intheoppositecase,whentheaccretionrate
accretionregime.Forthis,wenowassumethatthegastemperature
isstronglyEddingtonlimited,themassdecretionrateismaximal
isalwayslargeenoughthattheaccretionrateisequaltotheBondi
and almost equal to the Bondi rate. This results in a strong drag
rate.WealsoconsidertheSMBHtobeatrestinthecentreofthe
forcebetweenthesinkandthegas.
galaxy.Wethenobtainfortheaccretionrate
(GM )2
2.4 SMBHfeedback M˙acc(cid:39)4πρ c3s(itn)k (30)
s
Inthispaper,weonlyconsideramodelforwhichthermalenergy Wecannowsolvetheenergyequation,ignoringthecoolingterm,
isinjectedwithinthesinksphere,usingfortheSMBHluminosity and obtain the time evolution of the sound speed within the sink
thefollowingformula sphere
where (cid:15)r =0.1 is the aLccArGeNtio=n(cid:15)dciMs˙cacrca(cid:15)drcia2t,ive efficiency and (cid:15)(c24is) cs(t)=125(cid:15)c(cid:15)rc2(cid:32)GRMsisniknk(cid:33)2Rstink1/5 (31)
afreeparameterrepresentingthecouplingefficiencybetweenthe
Obviously, the temperature in the sink region will not grow in-
blastwaveenergyatsmallscaleandtheresultingthermalenergy
depositedatlargescale.Basedonpreviousworkusingtheramses definitely. As soon as it reaches a high enough value, the gas
in the vicinity of the SMBH will expand and cool adiabatically.
code(Teyssieretal.2011;Duboisetal.2012),wefixeditsvalue
to(cid:15)c=0.15,whichisquitetypicalofthecorrespondingliterature, We consider that we have reached the maximum temperature af-
with values ranging from 0.05 (Springel et al. 2005; Wurster & terone soundcrossingtimeof thesinksphere,namely tcross(t)=
Thacker2013)to0.15(Booth&Schaye2009;Gabor&Bournaud Rsink/cs(t)=t.Combiningthiswiththepreviousequationgivesus
themaximumpossiblesoundspeedinthehotaccretionphase
2013).
impleAmnenimtaptioorntanistitmhaptrowveemdeenptosciotmnpoawretdhetormthael perneevrigoyusatraemvesreys cs,max=125(cid:15)c(cid:15)rc2(cid:32)GRMsink(cid:33)21/6 (32)
finetimestep(i.e.thetimestepofthemaximumlevelofrefinement sink
(cid:96)max), and not only at main coarse time steps as before. We also Itcanbecomparedtothegalaxyescapevelocitytoassessthepos-
donotconsideraminimuminjectiontemperature,asinBooth& sibilityfortheSMBHtounbindthegasfromthenuclearregion.
Schaye(2009)orTeyssieretal.(2011).Moreover,thethermalen- Besides various constants that we set to our fiducial values
ergyisdistributedineverygascellwithinthesinkspherepropor- ((cid:15)c=0.15and(cid:15)r=0.1),weseethattheonlyvariablesenteringthe-
tionallytothegasdensity.Thismass-weighteddepositionscheme ses various formulae are the SMBH mass, Msink, the sink sphere
preventstheapparitionofunrealisticallylargegastemperature,as radius Rsink, and finally the average gas density within the sink
opposedtothevolume-weighteddepositionscheme. spherenH.Insertingtypicalvaluesforourpresentsimulation,we
Theseimportantchangesnowallowustomodelthecompe- can compute first the critical SMBH mass beyond which heating
titionbetweenheatingandcoolingwithinthesinksphere.Indeed, dominates over cooling, so that the sink sphere can exit the cold
onecanwriteanenergyequationfortheaveragegasspecificinter- accretionregimeandactuallyheatsthegasaroundtheSMBH
nalenergywithinthesinksphereas (cid:32) n (cid:33)2(cid:32) R (cid:33)3
ρD(cid:15) = LAGN−n2Λ(T) (25) Msink,crit(cid:39)8×104 M(cid:12) 100HH/cc 10s0inpkc , (33)
Dt V H
sink whereweassumedthegastemperaturetobefixedat106 Kinthe
wherethespecificinternalenergyisrelatedtothetemperatureand
coolingfunction.Ifthisisthecase,thenthetemperaturewithinthe
thesoundspeedby
sinkspherewillsteadilyincreaseaccordingtoEq.(32)andreach
k T themaximumsoundspeed
(cid:15)(cid:39) B (cid:39)c2(t) (26)
µmH s (cid:32) M (cid:33)1/3(cid:32) R (cid:33)−1/3
Wewantnowtodistinguishtworegimesofaccretiononthesink. cs,max(cid:39)750km/s 108siMnk(cid:12) 10s0inpkc (34)
MNRAS000,000–000(0000)
6 P.Biernackietal.
Insummary,ifenoughgasmakesitintothesinksphere,thedensity
Table1.SummaryoffiducialparametersrelatedtoSMBHsinkparticlesin
willbehighandcoolingwilldominate,maintainingthegastemper- ramsessimulations.
aturetorelativelylowvaluesandtheaccretionratetotheEddington
limit.If,ontheotherhand,thegasdensitywithinthesinksphereis Parameter Fiducial Description
toolow,orifthesinkmassistoolarge,weenterthehot,adiabatic value
regimeforwhichthegastemperatureisquicklyrisingtoitsmax-
imumvalue.Unfortunately,allthesequantitiesdependsensitively Mseed 106M(cid:12) Sinkseedmass
ontheadoptedresolution.Abetterspatialresolution,resultingin Mclump 108M(cid:12) Massoftheclumpinwhichweseedthe
asmallersinkradius,canreducethecriticalSMBHmass,butcan sink
alsoincreaseitbyallowingforlargergasdensities.Betterspatial
resolutioncanalsoincreasethegastemperatureinthehotaccretion Direct yes ThedirectN-bodysolverusedtoevolve
solver thetrajectoryofasink
regimesignificantly.Ifonewantstoalleviatethissensitivitytores-
olution,onecanrelatethegasdensityinthecriticalmasstotheav- Drag yes Gasdragforcefromaccretion
eragedensityoftypicalgascloudsthatarebombardingtheSMBH
inthenuclearregion,andonecanchoosetodepositthefeedback αboost Eq.(16) BoostfactorfortheBondivelocity
energywithinafixedradius,invokingotherphysicalprocessesto
setthisenergydepositionscale.Inwhatfollows,wewillusethis
simpleanalyticalargumentstointerpretournumericalresults.
length
3WeNuUseMthEeRAICMARLcSoEdTeUrPamses (Teyssier 2002) and its second- λJ=cs(cid:114)Gπρ = (cid:115)Γmπ2HkGBTn∗∗ (cid:39)332pc(cid:39)4∆xmin (36)
order, unsplit Godunov scheme to solve the Euler equations. The andintheminimumJeansmass
evolutionofdarkmatterandstarsisperformedwiththeAdaptive
4π (cid:18)λ (cid:19)3
Particle-Meshsolverwithcloud-in-cellinterpolation.Thedynam- MJ= n∗mH J (cid:39)4×106M(cid:12) (37)
ical evolution ofthe sink particle isperformed the direct gravita- 3 2
tionalacceleration(seeSubsection2.3). Themeshrefinementstrategywehaveadoptedforalloursimula-
Ourinitialconditionsfeatureanisolated,gas-rich,slowlyro- tionsisaquasi-Lagrangianapproach,wherecellsarerefinedonce
tating (spin parameter of 0.04) dark matter halo of 2×1012 M(cid:12) their mass exceed 8×mres, where our mass resolution is set to
sampled using one million dark matter particles. The halo has a mres (cid:39)1.5×105M(cid:12), so that our minimum Jeans mass is always
truncatedNFW(Navarroetal.1997)profilewithaconcentration sampledbyatleast32resolutionelements.Inallsimulations,star
parameterc=10andwiththecircularvelocityV200=160kms−1, formation is modelled with a Schmidt law with a rather low ef-
whichresultsintheradiusR200=230kpc,whilethehaloistrun- ficiency (cid:15)∗ =0.01 coming from observations of local molecular
catedat514kpc.Initially,thegaseoushaloisinhydrostaticequi- clouds(Krumholz&Tan2007).Collisionlessstarparticlesoffixed
libriumandhastheuniversalgasfractionof fgas=15%.Theini- mass1.3×105 M(cid:12) arespawnedstochasticallywithaPoissondis-
tialisationfollowsthesetupofTeyssieretal.(2013).Ourfiducial tributionifthegasdensityinthecellislargerthann∗mH(Rasera&
runhasaspatialresolutionof∆xmin=78pc. Teyssier2006).Feedbackfromsupernovae,ifconsidered,ismod-
Using an isolated cooling halo is dictated by a compromise elled with a non-thermal energy injection with efficiency of 10%
betweenrealisticbutexpensivecosmologicalsimulationsandide- and yield of 10% (1 M(cid:12) of metals for each 10 M(cid:12) of ejected
alisedbuthighlyresolvedisolateddiscsimulations.Sinceweare material). The non-thermal energy dissipation timescale is set to
usingarealisticinitialangularmomentumprofileinspiredfromthe 10Myr.Weboosttheefficiencyofoursupernovaefeedbackrecipe
average angular momentum distribution from N body simulation by grouping stochastically multiple star particles into one single
(Bullock et al. 2001), gas will be continuously accreted from the starclusterofmass108M(cid:12).
halo into the disc, with the right amount of angular momentum, AsitwasalreadymentionedinSection2.1,weallowonlyone
giving us the possibility to feed the nuclear region, and possibly sink to form in our galaxy. We form the sink at around 200 Myr
thecentralSMBH. after the start of the simulation. We use the phew clump finder
We use the Sutherland & Dopita (1993) model for radiative (Bleuleretal.2015)toidentifythemostmassivegasclumpasthe
coolingofgasforH,Heandmetallinesforgashotterthan104K formation site for the SMBH, and let the sink evolve from there.
andfrommetalfine-structurecoolingprocessesatlowertempera- AllfiducialparametersofourSMBHmodelarelistedinTable1.
tures.Weadvectthemetallicityintheformofapassivescalarand
we choosethe initialmetallicity tobe Zini=0.05 Z(cid:12). Apressure
floorisintroducedathighdensityandlowtemperature,toprevent
theuncontrolledfragmentationofgasbeyondthespatialresolution,
possiblyleadingtotheformationofnumericalsingularities(espe- 4 RESULTS
cially because we are using a low star formation efficiency). The
We now present our simulation results, including each important
temperaturecorrespondingtothepressurefloorissetto
processonebyone,inordertocomparethem,andgaugetheirrel-
(cid:32)n (cid:33)Γ−1 ative importance. These processes are listed in Table 2. For each
Tfloor=T∗ nH∗ (35) feedbackprocess,weusetheparametersdescribedintheprevious
section.WehoweverconsidertheSMBHseedmassasafreepa-
withacriticalgasnumberdensityn∗=9cm−3,acriticaltempera- rameter, and we explore values ranging from 105 to 109 M(cid:12), as
tureT∗=2×103K,andΓ=2.ThisresultsintheminimumJeans listedinTable2.
MNRAS000,000–000(0000)
SMBHsdynamicsingas-richgalaxies 7
enoughtohaveasustained,largerthanEddington,Bondiaccretion
Table2.Summaryofsimulationrunsandparametersusedinthisstudy.Pa-
rametersvariedwithrespecttothefiducialrunarehighlightedinboldprint. rate.
Columns:(1)subsectioninwhichthesimulationsareanalysed(withexcep- Wearguethatthecriticalmassforthesinkparticletoaccrete
tionforfiducialrun);(2)maximumallowedrefinementlevel;(3)fractionof fastenoughistheminimumJeansmassassociatedtoouradopted
SNenergydepositedinthegas;(4)dragforcemodelled(orinclusionofa mesh resolution. Indeed, assuming that v =0, we can re-write
rel
nuclearstarcluster);(5)initialseedmassinlog10M(cid:12);(6)AGNfeedback. theparameterthatcontrolswhetherBondiaccretionissubsonicor
supersonic(seeEq.8)as
Se(c1ti)on l(m2a)x (cid:15)(S3N) d(r4a)g m(s5e)ed AG(N6)fbk. x = c2sλJ (cid:39) MJ (38)
sink GM M
sink sink
4.1 14 0.0 yes 5 no where we used the fact that r =4∆x =λ . For our fiducial
sink min J
44..11 1144 00..00 yyeess 76 nnoo resolution,theJeansmass MJ is4×106 M(cid:12).Forthelowestseed
mass,whichisbelowtheJeansmass,accretionfollowstheBondi
4.1 14 0.0 yes 8 no
rate,andisratherlow,becausetheaccretionissubsonic.Notethat
4.1 14 0.0 yes 9 no
inthisregime,becausetheaccretionrateislow,thedynamicalcou-
4.2 14 0.0 yes 5 yes plingbetweenthesinkandthegasisweak,makingthesinkvery
4.2 14 0.0 yes 6 yes
sensitive to external perturbations. One can see in Figure 1 that
4.2 14 0.0 yes 7 yes
the trajectory of the sink particle is quite perturbed, with visible
4.2 14 0.0 yes 8 yes
oscillationsaroundthecentreofthegalaxy.Theseoscillationsin-
4.2 14 0.0 yes 9 yes
crease the relative velocity between the gas and the sink, further
4.3 14 0.1 yes 5 no contributing to the low accretion rate. Once the sink mass grows
4.3 14 0.1 yes 6 no beyond106M(cid:12),about800Myrafterthestartofthesimulationfor
4.3 14 0.1 yes 7 no thesmallseedmassorimmediatelyaftercreationfortheotherseed
4.3 14 0.1 yes 8 no
masses,theBondiaccretionevolvesfromsubsonictosupersonic.A
4.3 14 0.1 yes 9 no
muchmorerapid,Eddington-limitedexponentialgrowthfollows.
4.4 14 0.1 yes 5 yes Afterthisphase,theSMBHmassseemstosaturate,andgrows
4.4 14 0.1 yes 6 yes onlymildly,mostlybecauseoftheslowaccretionoffreshgasinto
4.4 14 0.1 yes 7 yes the nuclear region. Indeed, since we didn’t include any feedback
4.4 14 0.1 yes 8 yes processesinthisfirstexperiment,thefinalSMBHmassisregulated
4.4 14 0.1 yes 9 yes
bytheavailablegasmasswithinthenuclearregion.Thisregime,
4.5 14 0.1 NSC 5 yes calledhereaccretion-limitedgrowth,wasfirstdiscussedinBour-
4.5 14 0.1 NSC 6 yes naud et al. (2011). The late accretion phase is controlled by an-
4.5 14 0.1 NSC 7 yes gularmomentumtransferinthegalacticdisk,triggeredbyvarious
4.5 14 0.1 NSC 8 yes instabilitiesandslowlyfeedingtheSMBHwithfreshgas.Inthis
4.5 14 0.1 NSC 9 yes case,theSMBHtrajectoryremainswellwithinthenuclearregion,
adenseandmassivestellarbulgethatprovidesaverystableenvi-
4.6 15 0.0 yes 6 yes
4.6 15 0.1 yes 6 yes ronmentfortheSMBH.Asaresult,thesinkparticleneverleaves
4.6 15 0.1 NSC 6 yes thenuclearregion.
4.2 AGNfeedback-limitedgrowth
Wehaverepeatedthesamesimulationsasintheprevioussection,
4.1 Accretion-limitedgrowth
but this time with AGN feedback. The resulting dynamical and
Ourfirstsuiteofsimulationshasbeenperformedwithoutanyfeed- mass evolutions of the SMBH are shown in Figure 2a and Fig-
backprocessesandwithonlyonesinkparticleseededinthefirst, ure2b.Theonlydifferencewiththeprevioussetupisthefinalmass
massiveenough,nucleargasclump,growingviaEddington-limited ofthesink,whichisnowregulatedbyAGNfeedback.
Bondiaccretion.Becauseoftherelativelylowangularmomentum TheinitialgrowthoftheSMBHinoursimulationswithAGN
inourcoolinghalo,mimickingwhatweexpectfromcosmological feedbackisverysimilartotherunswithoutfeedback.Duetothe
simulations,thesesimulationswithoutfeedbackleadtotheforma- largeamountsofgasinthenuclearregion,feedbackheatingdoes
tion of a gas-rich, clumpy and bulge-dominated galaxy (see also no affect the gas surrounding the sink, as cooling dominates. As
Teyssieretal.2013;Duboisetal.2016)thatresemblesmanyob- soonastheaccretionrateishighenough,heatingdominatesover
served high-z galaxies, in particular the so-called “blue nuggets” coolingandtheSMBHquicklyreachesitsmaximummass,which
(seee.g.Damjanovetal.2009). inourcaseisaround2×108M(cid:12).
ThetrajectoryandthemassgrowthoftheSMBHareshownin Themaximum,self-regulatedmassisrelatedtotheheating-
Figure1.Forallouradoptedseedmasses,theSMBHremainswell coolingbalancewehavediscussedinSubsection2.4(seeEq.(32)
withinthenuclearregion(definedhereasthecentralkiloparsec),in and(34)).Usingthesimulationwith Mseed=106M(cid:12) asanexam-
whichtheywereborn.Interestingly,thelowestseedmass105M(cid:12) ple,weseethatataround420Myr,thesink’sgrowthisterminated.
showsaverydifferentbehaviourthantheother,largerseedmasses. Initially, the gas density in the sink sphere is quite large, around
Itsgrowthisveryslow,foralmost1Gyr,andonlywhenitreaches ∼600H/cm3,sothatclearlycoolingdominatestheenergybudget
106 M(cid:12) does it have a high enough accretion rate and grows ex- inthesinksphere.Gradually,astheSMBHmassgrows,feedback
ponentially.Theotherseedmassesstartgrowingexponentiallyim- isabletoheatthegasmoreandmoreinthesinksphere,untilthe
mediatelyaftertheircreation,whichmeansthattheyaremassive SMBHmassreachesthecriticalvalueforwhichheatingdominates.
MNRAS000,000–000(0000)
8 P.Biernackietal.
c]6
p
1e5 1e8
k
e [5 1e6 1e9 1010
r 1e7
t
n 109
e4
c ]
o fl
l M 108
ha3 [
m BH 107
o2 M 1e5
r
f 1e6
nce 1 106 11ee78
a
st 105 1e9
Di0
0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 1200 1400
Time [Myr] Time [Myr]
(a)Distancebetweenthesinkandthecentreofthehalo (b)SMBHmassgrowth
Figure1.EvolutionofdistancetothecentreofhaloandsinkmassfortherunswithoutneitherSNandAGNfeedbacksforfivedifferentseedmasses:105M(cid:12)
-red,106M(cid:12)-blue,107M(cid:12)-green,108M(cid:12)-purple,and109M(cid:12)-orange.Thesinkparticleoccupiespositioninthecentreofthehaloanditsgrowthis
limitedfirstbyEddingtonrateandlaterbyangularmomentumlossinthegas.
c]6
p
1e5 1e8
k
e [5 1e6 1e9 1010
r 1e7
t
n 109
e4
c ]
o fl
l M 108
ha3 [
m BH 107
o2 M 1e5
r
f 1e6
nce 1 106 11ee78
a
st 105 1e9
Di0
0 200 400 600 800 1000 1200 1400 0 200 400 600 800 1000 1200 1400
Time [Myr] Time [Myr]
(a)Distancebetweenthesinkandthecentreofthehalo (b)SMBHmassgrowth
Figure2.EvolutionofdistancetothecentreofhaloandsinkmassfortherunswithoutSNfeedbackbutwithAGNfeedbackforfivedifferentseedmasses:
105M(cid:12)-red,106M(cid:12)-blue,107M(cid:12)-green,108M(cid:12)-purple,and109M(cid:12)-orange.Thesinkparticleresidesinthecentreofthehalotravellingwithmost
massiveclumpanditsgrowthislimitedfirstbyEddingtonrateandlaterterminatedatself-regulationscaleduetoitsfeedbackheating.
Veryquicklythegassoundspeedwithinthesinksphererises,until right),afterwhichtheaveragesoundspeedquicklyexceedsvesc,
itexceedstheescapevelocityofthehalo.Whenthishappens,gas whichthenmarkstheendofthecoldaccretionregime(bottomleft)
isremovedfromthenuclearregionbyablastwave,whichreduces and the beginning of the hot mode of accretion. For comparison,
theaveragegasdensitydowntoorevenbelowρ¯(cid:39)10H/cm3,and wehaveplottedinFigure3oursimpleanalyticalpredictionsfrom
makesfeedbackevenmoreefficient.Feedbackisabletomaintain Subsection 2.4. We can predict quite nicely the onset of efficient
thesoundspeedtoahighvalue(cs(cid:39)400km/s),stronglyreducing heating, when the SMBH mass reaches its critical value, as well
theaccretionrate(seeEq.(5)andbottomleftpanelofFigure3). as the end of the mass growth, when the maximum sound speed
This is only when feedback processes are able to accelerate reachestheescapevelocityofthehalo.
gastotheescapevelocitythatthegrowthofSMBHishalted(see
also Silk & Rees 1998; Fabian 1999). This can be seen on Fig-
ure3,whereweplotvariousaveragequantitiesmeasuredinthesink The case with Mseed =109M(cid:12) is very different than all the
sphere.Thegasdensity(topleft)dropsbytwoordersofmagnitude othercases.Here,theinitialseedmassisalreadyabovethemaxi-
assoonasthemaximumsoundspeedsignificantlyexceedsthehalo mum,self-regulatedmass.AGNfeedbackimmediatelyblowsaway
escapevelocity(vesc,topright).ThecriticalSMBHmassMsink,crit, thegasfromthenuclearregion.Asaresult,thegasinthevicinity
forwhichheatingbalancescooling,isreachedat420Myr(bottom ofthesinkremainsveryhotandtheaccretionrateverylow.
MNRAS000,000–000(0000)
SMBHsdynamicsingas-richgalaxies 9
104 2500
cs
cs,max
cs,theo 2000
103
c] vesc /s]
c m
H/ 1500k
[ [
y 102 y
t t
si 1000ci
n o
e l
e
D
V
101
500
ngas
100 0
103 10 17
Bondi
r] 102 Eddington 10 18
y
/ ]
)
fl 10 19 3
M 101 m
[ c
e 10 20 s
at 100 g/(
n r 10 21 er
o [
i 10 1 e
et 10 22 at
cr R
c 10 2
A heating 10 23
cooling
10 3 10 24
300 350 400 450 500 300 350 400 450 500
Time [Myr] Time [Myr]
Figure3.Timeevolutionof1)averagegasdensitywithinthesinksphere(topleft);2)average,mass-weighted,soundspeed(blue)andmaximumsoundspeed
(red)withinthesinksphere(topright).Wehavealsorepresentedoursimpletheoreticalmodel(Eq.32and34)(green)comparedtothehaloescapevelocity
(orangedashed);3)Bondi(red)andEddington(blue)accretionrates(bottomleft)and4)averageheating(red)andcooling(blue)rateswithinthesinksphere
(bottomright)forsimulationwithAGNfeedbackonlyandMseed=106M(cid:12).
4.3 Supernovaefeedback-limitedgrowth proceedsmuchmorerapidlyandthesinkmasscangrowuptoits
accretion-limitedvalue,asinSection4.1.
We now remove AGN feedback from the picture, but include in-
In order to estimate the mass of the typical clumps that
stead supernova feedback from dying massive stars. We use the
will perturb the trajectory of the SMBH, we use the classical
samesimulationsuitethanbefore,withseedmassesfromM =
seed Toomre analysis of gas fragmentation in an idealised razor thin
105M(cid:12)toMseed=109M(cid:12).OnFigure4,weagainplotthetimeevo- disc(Toomre1964).ThelargestunstablewavelengthistheToomre
lutionofthedistanceofthesinkparticletothehalocentreandofits
length
mass.Hereagain,wecanseetwodifferentregimes.Lowandinter-
mediateseedmassesarequicklyremovedfromthecentralkilopar-
sec.There,supernovaefeedbackisefficientenoughtodestroythe GΣ f
gas gas
parent clump, and the sink particles are perturbed by interaction λ (cid:39) (cid:39) R . (39)
T κ2 π gal
withnearbyclumps.Asaconsequence,thetrajectoryofthesinks
becomes more complicated and eccentric, and the relative veloc-
itybetweenthesinkandthegaswithinthesinkspheregrowssig- In this approximate formula G is gravitational constant,
nificantly,reducingtheaccretionrateandthecorrespondingdrag Σgas=Mgas/πR2gal is the gas surface density, κ(cid:39)Vgal/Rgal is the
forceaccordingly.ForseedmasseslargerthanMseed=108M(cid:12),the epicyclicfrequency,Vgal= (cid:112)GMtot/Rgalisthegalaxycircularve-
sinktrajectoryappearsasmuchlessperturbedandthesinkmanage locityand fgasisthegas-to-totalmassfractioninthedisc.Inorder
to remain within the nuclear region. As a consequence, accretion forthiswavelengthtobetrulyunstable,theToomreparametermust
MNRAS000,000–000(0000)
10 P.Biernackietal.
satisfy In a previous section, we have seen that a seed mass lower
Q= πGcΣsκ <1. (40) WthaenseMeJneoanws trheastulatsseiendamnaasrstilfiocwiaelrlythalonwM, sTurbessounltiscianctchreetsioinnkrpaater-.
gas
ticlebeingscatteredoutofthenuclearregionbylargegasclumps.
wherecscanbetakenaseitherthesoundspeedorthevelocitydis- Similarly, a large seed mass with a dynamical friction time scale
persionofthegas.Undersuchconditions,onecanthenestimatethe comparableto(orshorterthan)theorbitaltimet (cid:39)200Myrwill
orb
massofthemostmassiveclumpsastheToomremassMT defined helpmaintainingthesinkparticlewithinthenuclearregion.
by Insummary,largeinitialseedmasses(108and109M(cid:12))have
a larger accretion rate, as M˙ ∝M2 , so they can grow fast,
MT=Σgasπ(cid:18)λ2T(cid:19)2(cid:39) Mt4oπtf2g3as. (41) atottohrebiirtaEldpdeirntugrtobnat-iloimnsi.teFdurrathtBeeo,rnmadniodreb,eMcBoHmiescqoumicpkalyralbelsestsoenMsitive,
T seed
ForaMilkyWay-likegalaxy,onehas Mtot(cid:39)1011 M(cid:12) inthedisc thussinkparticlesdonotsufferfromencounterswithlargermass
(not to be confused with the total mass in the halo, which is one perturbers.Also,itsdynamicalfrictiontimescaleisrelativelyshort,
orderofmagnitudelarger).Atlowredshift,ingalaxiessimilarto helpingtheSMBHtoremaininthecentre.
ourownMilkyWay,onefinds fgas(cid:39)0.1,whichresultsinatypi-
calclumpmassof MT (cid:39)2.5×106 M(cid:12).Athighredshift,however, 4.4 AGNfeedback-limitedgrowthwithsupernovaefeedback
likethecoolinghaloset-upweareadoptinginthispaper,thegas
fraction is much higher, fgas(cid:39)0.5, for a similar total mass. This WenowcombinesupernovaandAGNfeedback,repeatingthesame
leads to much bigger clumps, with MT (cid:39)3×108 M(cid:12). This value numerical experiments. As before, the low and intermediate seed
istypicalformassiveandgasrichgalaxies(see.e.g.Genzeletal. masses Mseed=105M(cid:12), Mseed=106M(cid:12) and Mseed=107M(cid:12) do
2008,2011;Guoetal.2012;Tacconietal.2013;Tamburelloetal. not really grow, as can be seen in Figure 5 for the red, blue and
2015,forin-depthdiscussion).Inconclusion,asinkparticlewith greenlines,andasitwasalreadythecaseforoursupernova-only
massM (cid:54)M willhaveitstrajectoryeasilydisruptedbyclumps feedback model. The large seed mass, on the other hand, are al-
sink T
inthedisc.Largersinkmasses,ontheotherhand,willresultina readytoocloseorevenlargerthantheirmaximum,self-regulated
muchmorestableorbitalevolution(seebelow). SMBHmass,asitwasalreadythecaseforourAGN-onlyfeedback
In order for the sink particle to reach (or remain in) the nu- model.Soeventheselargeseedmassesdonotfavourafastgrowth
clearregionofthegalaxy,weneedtoestimatethedynamicalfric- ofthesinkparticles,whicharecontinuouslyperturbedbyclumps
tion timescale as introduced by Chandrasekhar (1943). Although withmasscomparableorsmallerthantheToomremass.Moreover,
theoriginalformulawasderivedforacollisionlessfluid(darkmat- sincethesinkmassisnotgrowingmuchbeyond109 M(cid:12),thedy-
terandstars),averysimilarformulacanbeusedtocomputethe namical friction time scale remains longer or comparable to the
dynamical friction on the gas (Ostriker 1999). For the gas drag, orbitaltimescaleandthesinkparticleskeepmovingaroundwith
a correction factor must be introduced, compared to the original eccentricorbitsandlargepericentreradii(seeFigure5withviolet
collisionless case, but only for transonic relative velocities. For a andorangelines;alsoFigure10,leftcolumn).
SMBHwithatypicalorbitalvelocityof200km/s,thedragforceis SNandAGNfeedbacksworkhandinhandtocompletelypre-
likelytobeinthestrongsupersonicregime,forwhichnocorrection vent SMBH growth in this gas rich, highly turbulent and clumpy
isrequired. environment. We argue that only SMBH already as massive as
Using Chandrasekhar’s formula, we compute the dynamical 1010 M(cid:12) can survive in the nuclear region of such a galactic en-
frictiontimescalet (e.g.Eq.8.12ofBinney&Tremaine2008) vironment,becausetheyresisttheperturbationsfromclumpsand
df
because they have a short-enough dynamical friction time scale.
t = 1.65 R2orbσ (42) This conclusion is of course valid only if one considers that our
df lnΛGM twofeedbackmodelsarerealisticenough,whichisofcoursehighly
BH
speculative,sincetheyrelyonsub-gridphysics.Thesemodelsare
wheretheCoulomblogarithmisgivenby
neverthelessquitestate-of-the-art,andarerequiredtoexplainthe
lnΛ=lnRGgaMlVo2rb (43) lfoowrmsattairofnoqrmueanticohninegffiincimenacsysi(vfeorgSalNaxfieeesd(bfoacrkA)GanNdfteoedebxapclaki)n.star
BH
ThefactthatSMBHcannotgrowatall(excepttheextremely
Rorb and Vorb are the orbital radius and orbital velocity of the massiveones)ifonecombinesthetwosourcesoffeedbackenergy
SMBH.Assumingthatthevelocitydispersionofthecollisionless isthereforeafundamentalprobleminthetheoryofSMBHgrowth
componentsσ(cid:39)Vgal,theorbitalradiusandvelocityoftheSMBH andco-evolutionwithgalaxies.Thisalsoexplainswhymanyau-
to be of the order of the galaxy radius Rgal(cid:39)5 kpc and circular thorshavetorelyonartificialtrickstomaintaintheSMBHwithin
velocityVgal(cid:39)200kms−1,andfinallyusinglnΛ(cid:39)6.9asatypical thenuclearregionsofgalaxies,especiallywhenperforminghigh-
valueforourpurposes,wefindthedynamicalfrictiontimescaleto resolutionsimulations.
be
108M(cid:12) 4.5 GrowthwithinaNuclearStarCluster
t (cid:39)2.7Gyr . (44)
df M
BH Oneofthekeydifferencebetweenthesimulationwithsupernova-
OnlySMBHwithmassesgreaterthan109M(cid:12)willbeabletodecay feedback and the simulations without supernova-feedback is the
quickly enough to the centre of the galaxy, as they will have an presenceofamassivebulge,orinotherwords,amassivenuclear
orbitaldecayratecomparableorfasterthantheirrotationrate.Itis concentrationofstars(seeFigure8).Indeed,inthenosupernova
interestingtoseethattheToomremassandthecriticaldynamical feedbackcases,wedoformmassiveclumpsofgasandstarswith
frictionmassarebothcomparableto109M(cid:12)inhigh-redshiftMilky massesoftheorderof(orsmallerthan)theToomremass,thatap-
Wayanalogues(seealsoBournaudetal.2014). pear as bumps in the stellar surface density profile in Figure 8.
MNRAS000,000–000(0000)