Off the Beaten Path: A New Approach to Realistically Model The Orbital Decay of Supermassive Black Holes in Galaxy Formation Simulations

0.31 MB·English

by  Michael Tremmel

Mon.Not.R.Astron.Soc.000,1–8(2015) Printed12June2015 (MNLATEXstylefilev2.2)
Off the Beaten Path: A New Approach to Realistically Model
The Orbital Decay of Supermassive Black Holes in Galaxy
Formation Simulations
5 M. Tremmel1⋆, F. Governato1, M. Volonteri2,3, T. R. Quinn1
1
1Astronomy Department, Universityof Washington, Box 351580, Seattle, WA, 98195-1580
0
2Institut dÕAstrophysique, UMR 7095 CNRS, Universitńe Pierre etMarie Curie, 98bis Blvd Arago, 75014 Paris, France
2
3 Department of Astronomy, University of Michigan, Ann Arbor, MI, USA

n
u
J
12June2015

1
1
ABSTRACT

]
A
We introducea sub-gridforcecorrectiontermtobetter modelthe dynamicalfric-
G tion(DF)experiencedbyasupermassiveblackhole(SMBH)asitorbitswithinitshost
. galaxy.This new approachaccurately follows a SMBH’s orbital decay and drastically
h
improves over commonly used ‘advection’ methods. The force correction introduced
p
here naturally scales with the force resolution of the simulation and converges as
-
o resolution is increased. In controlled experiments we show how the orbital decay of
r the SMBHcloselyfollowsanalyticalpredictionswhenparticlemassesaresignificantly
t
s smallerthanthatoftheSMBH.Inacosmologicalsimulationoftheassemblyofasmall
a
galaxy, we show how our method allows for realistic black hole orbits. This approach
[
overcomesthe limitations of the advection scheme, where black holes are rapidly and

2 artificiallypushedtowardthehalocenterandthenforcedtomerge,regardlessoftheir
v orbits. We find that SMBHs from merging dwarf galaxies can spend significant time
9 awayfromthecenteroftheremnantgalaxy.ImprovingthemodelingofSMBHorbital
0 decay will help in making robust predictions of the growth, detectability, and merger
6
rates of SMBHs, especially at low galaxy masses or at high redshift.
7
0 Keywords: Numericalmethods:Supermassiveblackholes:cosmologicalsimulations:
1. dynamics
0
5
1
:
v 1 INTRODUCTION (Feruglio et al. 2010; Alatalo et al. 2011) and are thought
i to play a critical role in the quenching of star formation in
X Supermassive Black Holes (SMBHs) represent a crucial,
massive galaxies and the exponential cutoff of the galaxy
r though still poorly understood, aspect of galaxy evolu-
luminosity function at high masses (Springel et al. 2005a;
a tion theory. SMBHs with as much as 109 M⊙ in mass Bower et al. 2006; Croton 2009; Teyssier et al. 2011).
powerluminousquasars observedat z> 6(Fan et al. 2001;
Mortlock et al. 2011), while in the local Universe, SMBHs
Numerical simulations haveprovento bea critical tool
are ubiquitous in both massive galaxies (e.g. Gehren et al.
for understanding the formation and evolution of galaxies
1984; Kormendy & Richstone 1995; Kormendy& Ho 2013)
and extensive work has already been done to implement
and small, bulge-less disk galaxies (Shields et al. 2008;
SMBHs into these simulations (e.g. Di Matteo et al. 2003;
Filippenko & Ho 2003) as well as dwarfs (Reines et al.
Hopkinset al. 2005;Sijacki et al. 2009).Duetothescaleof
2011;Reines & Deller2012;Reines et al.2013;Moran et al.
these simulations, which often include a full cosmological
2014). Empirical scaling relations between the mass of
environment,theresolutionisnecessarilylimited.Accretion
SMBHsandthestellarmassandvelocitydispersionoftheir
ontoSMBHsandthefollowingfeedbackprocessesarehence
host galaxies are indicative of co-evolution (Häring & Rix
implementedviasub-gridprescriptions,underthebroadas-
2004; Gültekin & et al. 2009; Schramm & Silverman 2013;
sumptionthatconditionsatthesmallestresolvedscaledrive
Kormendy& Ho 2013). Powerful outflows resulting from
theBHevolutionatthescaleofjustafewparsecs.Onema-
feedback from an active nucleus have been observed
jorobstacleisthatitbecomesnecessarytoaccuratelyfollow
thedynamicsofasingleobject(theSMBH),somethingthat
cosmological simulations with force resolution larger thana
⋆ email:[email protected] few pc are inherently not well equipped to do. These dy-
©2015RAS
2 M. Tremmel et al.
namics can have important consequences for the growth of tionsbutaddsanadditionalfreeparameter,thedistancethe
SMBHs,particularly duringgalaxymergers. Detailed simu- SMBHispushedeachstep.Evendisregardingtheirindivid-
lationsofgalaxymergersshowthatthedynamicsofSMBHs ual drawbacks, noneof thesesolutions are ideal, as they all
areintimatelyconnectedwith accretion andviceversa(e.g. require the explicit assumption that SMBHs rapidly decay
Callegari et al. 2011). into and remain stable at the center of their host galaxies ,
Inthispaper,wedescribeourapproachtotakeintoac- which is often not accurate (Bellovary et al. 2010).
countunresolveddynamicalfriction(DF)intheorbitalevo- In strongly interacting systems, where SMBHs are
lutionofSMBHsandshowthatourprescriptionisapromis- thought to go through much of their growth and pos-
ing and realistic alternative to advection that can easily be sibly their formation (Mayer et al. 2007), the inner re-
implementedinexistingcodes.Weshowthatthisapproach, gions of a galaxy are likely to experience strong poten-
when applied to simulations with the typical resolution of tialfluctuations,affecting thedynamicsof theSMBH. This
modern cosmological simulations (∼0.1−0.5 kpc) leads to could be especially evident in dwarf galaxies, where black
more realistic dynamics that match well with analytic ap- holes exist (Reineset al. 2013; Moran et al. 2014) within
proximationsfor DFtimescales. InSection 2we summarize a shallow potential with a cored density profile (Oh et al.
somecurrentmethodsforcorrectingSMBHdynamicsincos- 2011; Teyssier et al. 2013), making perturbations to their
mological simulations and why a new method is needed. In orbits more likely. Additionally, SMBHs accreted during
section 3 we present our method. In section 4 we present both major and minor mergers, especially dry mergers
our results, comparing our method with a popular advec- (Kazantzidis et al. 2005), do not immediately sink to the
tionprescriptioninbothanisolateddarkmatter(DM)halo center of the galaxy, creating a population of small ‘wan-
and a fully cosmological zoomed-in simulation of a dwarf dering’ SMBHs (Islam et al. 2003; Volonteri & Perna 2005;
galaxy. In section 5 we summarize and discuss our results. Bellovary et al. 2010). Advection, by assuming that orbital
decaytimescales arealwaysshort,couldthenartificially in-
crease themass of central SMBHs, through inflated merger
and accretion rates.
2 THE NEED FOR A NEW MODEL
Moreadvancedapproacheshavebeenexplored.Forex-
Dynamical friction, the force exerted by the gravitational ample, Dubois et al. (2013) include a sub-grid prescription
wake caused by a massive object moving in an extended forgasdynamicalfrictionactingontheSMBH.Whilethisis
medium (Chandrasekhar 1943; Binney & Tremaine 2008) apromisingmethod,itrequiresassumptionsaboutthemul-
causes the orbits of SMBHs to decay towards the center of tiphasenatureandequationofstateofthegasfarbelowthe
massive galaxies (Governato et al. 1994; Kazantzidis et al. resolutionlimitofanycosmological simulation(∼1−5pc).
2005). In cosmological simulations, modeling DF is partic- Lupiet al.(2015)increase theresolution aroundSMBHsto
ulary challenging, as the particle representing a SMBH is attainveryaccuratedynamics,buttheirmethodisapplica-
often only a few times the mass of the background par- bleonlyforveryhighresolutionsimulationsmeanttoclosely
ticles. In this case SMBH orbits can become dynamically follow BH-BH mergers, not cosmological simulations.
heated due to the limited mass resolution and the result- The method we propose in thisLetter is a simple solu-
ing spurious collisionality from a noisy gravitational poten- tion to correcting SMBH dynamics in cosmological simula-
tial(Hernquist& Barnes1990).Thisnumericalheatingcan tions.Ourapproachistoestimatetheunresolveddynamical
cause the black holes to gain/lose energy and be unrealis- frictionfeltbytheSMBHandapplytheappropriateforceto
tically perturbed away from the center of its host halo. In the SMBH particle. This is a significant improvement over
order to lessen collisionality and improve performance, N- artificial advection,asit makesnoassumption about where
body simulations employ a gravitational softening length, the SMBHs should be located in a galaxy at a given time
a characteristic scale below which the gravitational force and it has no explicit effect on the SMBH surroundings. In
between two particles becomes damped. This mechanism, addition, it requires minimal assumptions about the state
while necessary numerically, hinders DF by preventing the ofthesimulationbelowtheresolutionlimitanditnaturally
close interactions that are necessary to form a wake in the converges with increasing resolution.
vicinity of the SMBH.
To account for the unrealistic dynamics of SMBHs
in cosmological simulations, many groups employ artifi- 3 METHODOLOGY
cial advection schemes that change the motion of SMBHs.
Each method comes with its own drawbacks. For exam- 3.1 The Test Simulations
ple, placing the SMBH at the position of the lowest po-
To test the dynamics of black holes in a realistic DM halo,
tential gas particle around it (e.g. Springelet al. 2005b;
we generate initial conditions for a collapsing overdensity
Booth & Schaye 2009) causes chaotic motions, especially following the procedure of Evrard (1988)1. We begin the
when relative velocity constraints are put on the gas par-
simulationbyapproximatingthestateofahalobeginningto
ticles (Wurster& Thacker 2013). Utilizing a large ‘tracer
collapseathighredshift.Wethenallow thehalotocollapse
mass’ with which the SMBH gravitationally interacts with before implanting a black hole of mass 106 M⊙. We run a
its surroundings (Debuhret al. 2011) affects the morphol-
suiteofsimulationsatdifferentmassandspatialresolutions
ogy in thegalactic centeras well astheaccretion history of
the SMBH (Wurster& Thacker 2013). Pushing the SMBH
a certain distance in the direction of either the local stellar 1 Created using the package ICInG created by M. Tremmel.
density gradient (Okamoto et al. 2008) or the local center The code is publicly available and can be downloaded at
of mass (Wurster& Thacker 2013) can avoid chaotic mo- https://github.com/mtremmel/ICInG.git
©2015RAS,MNRAS000,1–8
Off The Beaten Path 3
Table 1.Information about theresolution ofeach isolateddark hole. We assume that within ǫg from theblack hole the ve-
matterhalotestsimulation. locity distribution is isotropic, giving Chandrasekhar’s Dy-
namical Friction Formula(Chandrasekhar 1943).
Name Number DarkMatter Softening
ofParticles Mass(M∗⊙ length, ǫg (pc)† v vBH
OvLeorwsamRpesled 81..3095××110066 19..2728××110055 331111 aDF =−4πG2MmalnΛvB3BHH Z0 dvava2f(va) (2)
HighRes 6.71×107 1.53×104 77 This can be furthersimplified by substitutingthe inte-
SmallSoft 6.71×107 1.53×104 10 gral for ρ(<vBH), which is the density of particles moving
slower than theblack hole.
∗Massofthedarkmatter particlesinthesimulation
†splinekernel gravitational softening v
aDF =−4πG2Mρ(<vBH)lnΛv3BH (3)
BH
(seeTable1)withamassive‘blackhole’particleinitially1)
at the center while the halo is still actively collapsing or 2) TheCoulomb logarithm dependson themaximumand
onaneccentricorbitafter thehalohasmostlyrelaxed.The minimum impact parameters, bmax and bmin, such that
resolutionofthe‘LowRes’and‘Oversampled’runsarecom- lnΛ ∼ ln(bbmmainx). Because DF is well resolved at scales
parabletocurrentuniformvolumesimulations(Kereš et al. greaterthanthesofteninglength,wesetbmax =ǫg toavoid
2012; Sijacki et al. 2014; Schayeet al. 2015; Dubois et al. doublecountingfrictional forcesthatarealreadyoccurring.
2014a), while the ‘High Res’ simulation is representative For the minimum impact parameter, we take it to be the
ofhighresolution ‘zoomed-in’simulations(Governato et al. minimum90◦ deflectionradius,withalowerlimitsettothe
2012; Christensen et al. 2014; Hopkinset al. 2014). The Schwarzschild Radius,RSch.
force softening lengths adopted are typical of cosmological
rtuionnss,,wthheilefoirnce‘SrmesaollluStoioftn’,isaovnalryia1n0tpocf,tthyep‘iHcaiglhofRsiems’usliamtiuolnas- bmin =max(b90,RSch);b90 = GvM2BH (4)
ofisolatedbinarymergers(Capelo et al.2015)orveryhigh- BH
res cosmological zoomed-in simulations (e.g. Dubois et al. Forthecalculation,weuse64collisionlessparticles(i.e.
2014b). darkmatterandstarparticles,ifpresent)closesttotheblack
We verified that the halo formed in the collapse has a hole.Wecalculatethevelocityofeachparticlerelativetothe
density profile typical of CDM halos (Navarroet al. 1996). COMvelocityofthose64particles. Weverifiedthatourre-
Thehaloweuseherehas(afteratimeconsistentwithz∼0) sults do not depend strongly on the number of neighbors
Mvir ∼ 2×1011 M⊙, Rvir ∼ 115 kpc, and a concentration used, although using too few particles could result in nu-
c∼4.5.Ananalyticexpressionoftheapproximatetimescale merical noise in the calculation of this force. Since we are
forarigidobjecttosinktothecenterofaDMhalowascal- explicitlyassumingthevelocitydistributionisisotropic,the
culatedbyTaffoni et al. (2003).Giventheseconditionsand following must be true.
themass ofthetest black hole, theestimated DFtimescale
(tτaDncFe)ofofr2ankpecccfernomtricth(ev=ha0l.o1vcceinrct)erorisbitapaptroaxniminaittiealyl d1i.s8- ρ(<vBH)= M(<vBH)ρ (5)
Mtotal
Gyr. This initial orbit was chosen as typical of a SMBH af-
ter its parent satellite has been tidally disrupted. It has a Where ρ, the total density around the black hole, is
non-negligible τDF, but still much less than a Hubbletime. calculated by smoothing over the chosen 64 particles, i.e.
ρ= 6i4miW(rBH −ri,h). M(<vBH) is the total mass of
the cPhosen particles that are moving slower than the black
3.2 The Dynamical Friction Prescription
hole relative to the local center of mass velocity and Mtotal
Dynamical friction occurs due to both large (Colpi et al. is thesummed mass of all 64 particles.
1999) and small scale perturbations tothe black hole’s sur- The resulting acceleration (from eq. 3) is added to the
roundings. We consider perturbations on scales larger than blackhole’scurrentacceleration,tobeintegratedthefollow-
the gravitational softening length, ǫg, to be well resolved. ing time step. As the spatial resolution or black hole mass
On these scales the potential should be smooth, so long as increases (or the velocity of the black hole decreases) bmin
enoughparticlesareusedsothatdynamicalheatingismin- will become greater than bmax, in which case we claim DF
imized. is being fully resolved and therefore the correction is not
The acceleration a black hole of mass MBH feels from needed.
DF due to particles of mass ma ≪ MBH is given by the This method is not accounting for DF from gas, which
following formula: can have important effects for supersonic black holes in re-
gions where gas density dominates stars and dark matter
(Ostriker 1999; Chapon et al. 2013). This should not occur
aDF =−4πG2MBHmalnΛZ d3vaf(va)|vvBBHH−−vvaa|3 (1) oofftleanrgoenrtghaelasxciaelsesarreeldeovmanitnaintetdhebsyesstiamruslaantidonsms.aTllheercgeanlatxer-
Where f() is the velocity distribution function, lnΛ is ies have significant star and dark matter fractions within
the Coulomb logarithm, va is the velocity of the surround- thecentral regions (Oh et al. 2008,2011),so thiseffect will
ing background objects relative to the local center of mass only be a minor correction in most cases. DF may be over-
(COM)velocity,andvBH istherelativevelocityoftheblack estimated within resonant DMcores where DFcan become
©2015RAS,MNRAS000,1–8
4 M. Tremmel et al.
No Added Correction With DF Correction
Low Res
Oversampled
)
c
p2000 High Res
(
Small Soft
r
e
t
n
e
C1500

o
al
H

m
o1000
r
f

e
c
n
a
t 500
s
i
D
0
0 1 2 3 4 5 0 1 2 3 4 5
Time (Gyr)
Figure 1.Effects of the DF correction: Distanceofatest106 M⊙ blackholefromhalocenter asafunctionoftimeatdifferent
resolutions.Dashed-dotanddashedlinesindicateǫg and2ǫg respectivelyforboththeLowResandOversampledmodels.Theblackhole
starts onan eccentric (v=0.1 vcirc)orbit withapocenter of 2kpc. The vertical solidline represents the analytically derived timescale
for orbital decay. When the DF correction is applied, marked improvement is seen for all models except ‘Low Res’, which experiences
too muchdynamical heating due tothe lowermass resolution. The orbits of ‘HighRes’ and‘SmallSoft’ arevery nearly the sameonce
thecorrection isimplemented,indicatingnumericalconvergence whentheDMparticlemassis∼104 M⊙.
much less efficient (Read et al. 2006). This effect is sec- 2kpc,placedinthehaloafterithasfinishedmostofitscol-
ondary and mainly important when the gravitational soft- lapse.Thecenterofthehaloisdefinedateachstepusingthe
eninglength isappreciable compared tothesize ofthecore shrinking spheres method (Power et al. 2003). We verified
structure. Often the orbital differences should be smaller thatthedensitymaximumandpotentialminimumcoincide
than or similar to the resolution limit and therefore unim- withinmuchlessthantheforceresolution.Theverticalline
portant.Additionally,interactionswithclumpygashasbeen representsthedynamicalfriction timescale fortheorbit de-
shown to significantly increase the timescale for the or- rived from the analytic model of Taffoni et al. (2003) and
bitaldecayofSMBHbinariesbelow∼100pc(Fiacconi et al. the horizontal lines represent ǫg and 2ǫg. Without the DF
2013;Roškaret al.2015).However,thiseffectwould notbe correction, only the ‘Small Soft’ model, with 10 pc spatial
well resolved byevenhigh resolution ‘zoomed-in’ cosmolog- resolution and DM particle mass almost 100 times smaller
ical simulations. thantheBH,isabletoshowsubstantialorbitaldecaywithin
2τDF. Implementing our DF correction results in a notice-
ableimprovementfortheorbitaldecay,evenattherelatively
modest resolution of the Oversampled model, where it falls
4 RESULTS
towithin2ǫgofthecenterbefore2τDF.Athigherresolution,
4.1 Isolated Dark Matter Halo thedynamics converge to closely match with the analytical
approximation. Note that even for our highest resolution
We find that, for SMBHs placed in the center of a collaps- simulation,SmallSoft,theDFcorrection causestheSMBH
inghalo,onlytheLowRessimulationexperiencessignificant tosinkalmost1Gyrsooner.Theseareveryencouragingre-
dynamicalheating,causingtheBHtobeunrealisticallyper- sults,astheyindicatethatthiscorrectionresultsinrealistic
turbedawayfarfromhalocenter.Intherestoftherunsthe blackholeorbitalevolutionevenatresolutionsattainablein
BH remains within one softening length from halo center largevolumesimulationsandithasimportantconsequences
for 6 Gyrs and shows no sign of heating. This is due to the even at the highest resolutions tested here.
highermassresolutionpresentinthoseruns.Simulationsin-
cludingeither ourDFcorrection or an advection correction Figure2comparestheperformanceofourDFprescrip-
(see below) have the same results, showing little difference tion to that of a commonly used advection scheme used in
from simulations with no correction at all for this scenario. Sijacki et al.(2007) andvariousothersimulations. Thetest
Figure1showstheorbitalevolutionofablackholeini- is done using the Oversampled run, as this most closely re-
tially on an eccentric (v=0.1 vcirc) orbit with apocenter of semblestheresolution of acosmological volumesimulation.
©2015RAS,MNRAS000,1–8
Off The Beaten Path 5
z
5.04.0 3.0 2.0 1.0 0.5
c)
kp101 Dynamical Friction
er (pc)2000 NDAdoyv nCe. ocFrtrrii.ocnt. alo Center (1100-10
ent m H10-2
alo C1500 ce fro10-3
m H stan10-4
o1000 Di10-5
r 2 4 6 8 10
e f
nc c) 5.04.0 3.0 2.0 1.0 0.5
Dista 500 nter (kp110001 Advection
e
0 o C10-1
0 1 2 3 4 5 6 al
Time (Gyr) m H10-2
o
Figure 2. DF correction vs Advection: Results from e fr10-3
the‘Oversampled/modelwhenimplementingourDFcorrection anc10-4
(blue) compared with a commonly used advection routine (red) st
and no correction to dynamics (cyan). Dashed-dot and dashed Di10-5
2 4 6 8 10
linesindicateǫg and2ǫg respectively.The blackholestarts onan Time (Gyr)
eccentric(v=0.1vcirc)orbitwithapocenterof2kpc.Advection
immediatelypushes theoff-center blackholetothecenter, miss-
ingtheorbitaldecaythatourmethodcaptureswell.Withoutany Figure 3. DF correction in a cosmological dwarf sim-
correction, theorbitdecays fartooslowly,remainingfar(>2ǫg) ulation:The dynamics of four black holes in the cosmological
fromhalocenter evenafter6Gyr. zoomed-indwarfgalaxysimulationwithDF(top)andadvection
(bottom).Thesearetheblackholesthatendupinthemostmas-
sivesystembytheendofthesimulation.Eachcoloredlinetraces
The BH is placed on an eccentric orbit, as in Figure 1. The the distance of a black hole from the center of the most mas-
advection scheme adopted repositions the black hole each sive halo. Black dots mark merger events and the dashed lines
time step to the position of the lowest potential particle markthegravitationalsofteninglengthofthesimulation(87pc).
within its 32 nearest neighbors while keeping the velocity Which of the two black holes emerges from a merger event and
whichis‘eaten’isunimportant.DFisabletosustainalong-lived
unchanged. Not surprisingly, this results in the black hole
dualblackholesystem(blueandred)whiletheadvectionscheme
staying very close to the center of the halo even when ini-
causesthemtoquicklymerge.Thegreenblackholeremainsona
tiallysetonanoff-centerorbit.TheDFcorrectioncaptures
very wideorbit inthe DFrun, but isquicklyand unrealistically
the more gradual orbital decay that the advection scheme
pulledtothecenter withadvection.
completely misses. The run with no dynamical correction
fails to havetheBH sink within 2ǫg evenafter 6 Gyr.
Thisisan importantconclusion becausethesedifferent
galaxy for this test because SMBHs are more likely to be-
orbital evolutions would result in drastically different ac-
comeperturbedaway from galactic center,given theirshal-
cretion histories for the black hole. Off center BHs should
lowgravitationalpotentialandactivelyevolving,coredDM
accrete less due to lower gas densities. Additionally, simu-
profile(Governato et al. 2012;Pontzen & Governato 2012).
lations that utilize advection would have black holes merge
This will guarantee a useful test environment for exploring
much sooner than what is predicted by their orbital decay
thedifferences between out method and advection.
timescale. Our improved method should then have impor-
This test is also topical, as there is a growing sam-
tantimplications forthegrowth andmergerrateofSMBHs
ple of dwarf galaxies with detected SMBHs (Reines et al.
in cosmological simulations of galaxy formation.
2013; Moran et al. 2014). Realistic numerical studies of BH
formation and growth in these small galaxies, focusing on
theiroccupationfractionandhowtheyandtheirhostgalax-
4.2 Cosmological Dwarf Galaxy Simulation
ies evolve toward the correlation with the stellar veloc-
As a first test of the dynamics of SMBHs in a fully cosmo- ity dispersion, would provide vital constraints on BH seed
logical setting,werunahighresolution‘zoomed-in’simula- massesandearlygrowthmechanisms(Volonteri et al.2008;
tion that results in two dwarf galaxies with masses ∼ 1010 Volonteri& Gnedin 2009; Volonteri2010).
M⊙ at z = 0. The simulation has a resolution similar to We use the new N-body + SPH code ChaNGa
our High Res isolated halo model, with dark matter par- (Menon et al. 2015), which includes all of the physicsmod-
ticle mass 1.6×104 M⊙, gas particle mass 3.3×103 M⊙, ules previously implemented in Gasoline (Wadsley et al.
andgravitational softeningofonly 87pc.Weshowed inthe 2004) such as hydrodynamics, gas cooling, a cosmic UV
previoussection that at thisresolution ourDFprescription background,starformationandSNefeedback.The‘zoomed-
gives results that match analytic models. Wechose a dwarf in’approachpreservesthelargescaletidalfieldwhileallow-
©2015RAS,MNRAS000,1–8
6 M. Tremmel et al.
10
Dynamical Friction Advection
5
]
c
p
k 0
[

y
−5
−10
−10 −5 0 5 10−10 −5 0 5 10
x [kpc] x [kpc]
Figure 4. BHs in a Cosmological Dwarf: A snapshot of a zoomed-in cosmological simulation of a forming dwarf galaxy at z
= 0.846. The gas density integrated along the line of sight is shown with darker colors indicating higher densities. In the dynamical
friction simulation, aprevious merger has created acentral binary BH system (red and blue, see Figure 3). The separation of ∼1kpc
iswellresolved bythesimulation,whichhas aforce resolutionof 87pc.A morerecent merger has setathirdBH(ingreen) onawide
orbit (see Figure 3). In the image, the green BH is at its closest approach to galactic center. In the advection simulation, all BHs are
quicklypushedtothecenter, wheretheymerge, causingthesimulationtomissthesemorerealisticBHorbits.
ingustomodelasmallregionathighresolution.Thisissim- enoughtobegravitationally bound.Thedashed blacklines
ilar to the simulation described in Governato et al. (2010) indicate thegravitational softening length of 87 pc.
andinShenet al.(2014).Inthisparticularsimulation,stars Theblackholesbecomeperturbedasthered,cyan,and
and DMdensities dominatethat of gas within theinnerre- blue host galaxies interact between z = 3 and 4. In the
gions by more than a factor of 10. Examining other dwarf advectioncasetheblackholesaredrivenquicklytowardthe
galaxysimulations(e.g.Governato et al.2010),wefindthat center where they all merge, leaving only one black hole
gasdensitiescanoftenreachsimilarvaluestoDMandstars. (labeledasred).WiththeDFcorrection,onlythecyanand
Ineithercasethecontributionfrom gastoDFisonlyami- blueBHsmerge.Aftertheredandbluegalaxy hostsmerge
nor correction, at most a factor of a few, and negligible in with DF,the blueand red BHs remain orbiting around the
our current example. center of the merger remnant. The blue BH comes back to
thecenteronlyafter4GyrandtheredBHremainsorbiting
Byz<1thesesystemsareagoodrepresentationofreal
at around 1 kpc (11 times the force resolution) for another
dwarf galaxies, with an extended stellar disk, no bulge, a
4 Gyrbefore sinkingand merging with theblue BH.
high gas fraction and a cored DM profile. We first run the
The more striking difference between the DF and the
simulation until z = 6, when we insert five black holes of
mass 5×105 M⊙ in the centers of the five most massive advectionruninvolvesthegreenblackhole.Whenthemuch
halos at the time, which all have a mass of 2× 108 M⊙ smaller green host galaxy merges with theblue/red host at
z ∼ 1, it is initially far from halo center (∼ 30 kpc) and
orhigher.Inthesesimulationsblackholesdonotaccreteor
is quickly disrupted by the main galaxy. With DF, the BH
producefeedback,asweareonlyinterestedinfollowingtheir
stayson awide orbit, nevercoming muchcloser than afew
dynamics. From z=6 we run two simulations to z < 1, one
kpc from halo center. In the advection case, however, the
with our DF routine and the other with advection. By the
green black hole is quickly pushed to the center where it
endofthesimulation,thereareonlytwomajorstarforming
galaxies with masses ∼1010 M⊙. misearngesunwrietahlitshtieccreensturlatl,raesdtBhHe D(sFeetFimigeusrceasle3oafnda 45).×T1h0i5s
We then run the Amiga Halo Finder (AHF) M⊙ blackholethatfarawayfromthecenterofsuchasmall
(Knollmann & Knebe2009) forallthesavedsnapshotsand galaxy would belonger than a Hubbletime.
calculate thecenterof themain halo at each step using the With this simulation we can clearly see how the choice
shrinking spheres approach. In Figure 3 we follow the tra- of dynamical correction can affect the ability of SMBHs to
jectoriesoffourblackholeswithrespecttothecenterofthis become perturbed during mergers. In the DF simulation,
halo (which originally just has the blue BH at its center). BHs are able to remain off-center for many Gyr while with
Each color represents a different black hole. Black dots in- advectiontheyarequicklydriventothecenter.Additionally,
dicate a black hole merger, which happens when two BHs theDFsimulation allowsforsustainedwideorbitsresulting
come within two ǫg of one another at relative speeds low from minor mergers. Suchdynamicscan havean important
©2015RAS,MNRAS000,1–8
Off The Beaten Path 7
impactoninterpretingtheconnectionbetweentheinitialoc- all phases of galaxy evolution, including merger events.
cupationprobabilityofSMBHseedsindwarfgalaxyprogen- Implementation of DF routines such as the one presented
itors and the observed occupation at low redshift. Methods here will improve the ability of cosmological simulations
suchastheadvectionschemepresentedherewouldpredicta to accurately model SMBH accretion, growth and energy
more direct connection, while the simulation with DF indi- depositionintheIGMand,therefore,increasetheabilityof
catesthatthenatureofthemergers(i.e.massratioandori- simulations to interpret and predict observational results.
entation) can havean impact on which dwarf galaxies have The implementation of this approach in our future cosmo-
observable SMBH activity and when that activity occurs. logical volume and zoom simulations represents an exciting
The DF correction could also have important implications chance to realistically study the growth and merger rate of
for BH merger rates, allowing them tobecome more decou- SMBHs across cosmic time.
pled from galaxy merger rates than advection simulations.
This will affect predictions of gravitational wave detections
as well as estimates for recoiling BHs, which can have an FG and TQ were funded by NSF grant AST-0908499.
important effect on observability (e.g. Madau & Quataert FG acknowledges support from NSF grant AST-0607819
2004). andNASAATPNNX08AG84G.MVacknowledgessupport
This simulation is a useful illustration of the variety of from SAO award TM1-12007X, from NSF award AST
different,realisticBHorbitsourmethodallowscomparedto 1107675, and from a Marie Curie FP7- Reintegration-
commonly used advection schemes. In future work, we will Grant within the 7th European Community Framework
explorethedynamicsofBHsinavarietyofdifferentmerger Programme (PCIG10-GA-2011-303609). Some results
eventswithin a cosmological context. were obtained using the analysis software pynbody
(Pontzen et al.2013).ChaNGawasdevelopedwithsupport
from NationalScienceFoundationITRgrant PHY-0205413
to the University of Washington, and NSF ITR grant
5 DISCUSSION AND SUMMARY
NSF-0205611 to the University of Illinois. Simulations
We have introduced a sub-grid force correction term for were run using NCSA Bluewaters and NASA Pleiades
SMBH motion based on dynamicalfriction. Thiscorrection computers. Resources supporting this work were provided
allows us to better model the orbital decay of SMBHs in by NSF PRAC Award 1144357 and the NASA High-End
numerical simulations. We have shown using controlled ex- Computing (HEC) Program through the NASA Advanced
periments of isolated DM halos that this addition matches Supercomputing (NAS) Division at Ames Research Cen-
analytic predictions of the orbital decay in DM halos with ter. We thank Andrew Pontzen, Aycin Aykutalp, John
resolutions attainable by large-volume cosmological simu- Wise, Jillian Bellovary, and Alyson Brooks for stimulating
lations. We have also demonstrated that our prescription discussions and a careful reading of themanuscript.
naturally converges with resolution.
This method is a significant improvement over exist-
ing ‘advection’ methods that force a short orbital decay
timescale regardless of the dynamical state of the system.
REFERENCES
When applied to a cosmological dwarf galaxy simulation,
our method results in noticeably different black hole dy- Alatalo K. et al., 2011, ApJ,735, 88
namics compared with theadvection scheme. In particular, Aykutalp A., Wise J. H., Spaans M., Meijerink R., 2014,
our prescription: ApJ,797, 139
Bellovary J. M., Governato F., Quinn T. R., Wadsley J.,
• Models the perturbation and gradual orbital decay ShenS., Volonteri M., 2010, ApJ, 721, L148
of a central BH duringa galaxy merger Binney J., Tremaine S., 2008, Galactic Dynamics: Second
Edition. Princeton UniversityPress
• Allows for long-lived dual BH systems with close Booth C. M., SchayeJ., 2009, MNRAS,398, 53
(<1 kpc) orbits. Bower R.G., Benson A.J., Malbon R., Helly J. C., Frenk
C.S.,BaughC.M.,ColeS.,LaceyC.G.,2006,MNRAS,
• Maintains a stable central BH when appropriate. 370, 645
Callegari S.,KazantzidisS.,MayerL.,ColpiM.,Bellovary
• Allows for sustained wide (>5 kpc)orbit BHs. J. M., Quinn T., Wadsley J., 2011, ApJ, 729, 85
Capelo P. R., Volonteri M., Dotti M., Bellovary J. M.,
Correctlymodelingtherichorbitaldynamicsofablack MayerL., GovernatoF., 2015, MNRAS,447, 2123
holewithinitshostgalaxycanhaveimportantconsequences Chandrasekhar S., 1943, ApJ, 97, 255
for its accretion history, duty cycle and observability that Chapon D., Mayer L., Teyssier R., 2013, MNRAS, 429,
was previously neglected in simplified ‘advection’ schemes. 3114
The dynamically complex and more realistic orbits allowed Christensen C.R.,Brooks A.M., Fisher D.B., Governato
by our method will have crucial implications for the early F.,McCleary J., QuinnT.R.,ShenS.,Wadsley J.,2014,
growth of SMBHs, which takes place in small, rapidly MNRAS,440, L51
growing galaxies at high redshift (Aykutalpet al. 2014). Colpi M., Mayer L., Governato F., 1999, ApJ, 525, 720
Additionally, understanding the relative importance of dif- Croton D. J., 2009, MNRAS,394, 1109
ferent accretion mechanisms throughout a SMBH’s lifetime Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412,
requirestheabilitytoaccuratelymodelitsdynamicsduring 1341
©2015RAS,MNRAS000,1–8
8 M. Tremmel et al.
Di Matteo T., Croft R. A. C., Springel V., Hernquist L., Reines A.E., Deller A.T., 2012, ApJ,750, L24
2003, ApJ, 593, 56 Reines A.E., Greene J. E., Geha M., 2013, ApJ,775, 116
Dubois Y., Gavazzi R., Peirani S., Silk J., 2013, MNRAS, ReinesA.E.,SivakoffG.R.,JohnsonK.E.,Brogan C.L.,
433, 3297 2011, Nature, 470, 66
DuboisY.,VolonteriM.,SilkJ.,2014a,MNRAS,440,1590 Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn
Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., T. R., Wadsley J., 2015, MNRAS,449, 494
2014b, MNRAS,440, 2333 SchayeJ. et al., 2015, MNRAS,446, 521
Evrard A.E., 1988, MNRAS,235, 911 Schramm M., Silverman J. D., 2013, ApJ, 767, 13
Fan X.et al., 2001, AJ, 122, 2833 Shen S., Madau P., Conroy C., Governato F., Mayer L.,
Feruglio C., Maiolino R., Piconcelli E., Menci N., Aussel 2014, ApJ, 792, 99
H., Lamastra A., Fiore F., 2010, A&A,518, L155 ShieldsJ.C.,WalcherC.J.,BökerT.,HoL.C.,RixH.-W.,
Fiacconi D., Mayer L., Roškar R., Colpi M., 2013, ApJ, van der Marel R. P., 2008, in IAU Symposium, Vol. 245,
777, L14 IAU Symposium, Bureau M., Athanassoula E., Barbuy
Filippenko A. V.,Ho L. C., 2003, ApJ, 588, L13 B., eds., pp. 259–260
GehrenT.,FriedJ.,WehingerP.A.,WyckoffS.,1984,ApJ, Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007,
278, 11 MNRAS,380, 877
Governato F. et al., 2010, Nature, 463, 203 Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS,
Governato F., Colpi M., Maraschi L., 1994, MNRAS,271, 400, 100
317 Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey
Governato F. et al., 2012, MNRAS,422, 1231 P., Snyder G., Nelson D., Hernquist L., 2014, ArXiv e-
Gültekin K., et al., 2009, ApJ, 698, 198 prints
Häring N., Rix H.-W., 2004, ApJ, 604, L89 Springel V., Di Matteo T., Hernquist L., 2005a, ApJ, 620,
Hernquist L., Barnes J. E., 1990, ApJ, 349, 562 L79
HopkinsP.F.,HernquistL.,CoxT.J.,DiMatteoT.,Mar- Springel V., Di Matteo T., Hernquist L., 2005b, MNRAS,
tiniP., Robertson B., Springel V.,2005, ApJ, 630, 705 361, 776
Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.- Taffoni G., Mayer L., Colpi M., Governato F., 2003, MN-
A.,QuataertE.,MurrayN.,BullockJ.S.,2014,MNRAS, RAS,341, 434
445, 581 Teyssier R., Moore B., Martizzi D., Dubois Y., Mayer L.,
Islam R.R., Taylor J. E., Silk J., 2003, MNRAS,340, 647 2011, MNRAS,414, 195
Kazantzidis S.et al., 2005, ApJ, 623, L67 Teyssier R.,PontzenA.,DuboisY.,ReadJ. I.,2013, MN-
Kereš D., Vogelsberger M., Sijacki D., Springel V., Hern- RAS,429, 3068
quist L., 2012, MNRAS,425, 2027 Volonteri M., 2010, A&ARev., 18, 279
Knollmann S. R.,KnebeA., 2009, ApJS,182, 608 Volonteri M., Gnedin N.Y., 2009, ApJ,703, 2113
Kormendy J., HoL. C., 2013, ARA&A,51, 511 VolonteriM.,LodatoG.,NatarajanP.,2008,MNRAS,383,
Kormendy J., RichstoneD., 1995, 33, 581 1079
Lupi A.,Haardt F., Dotti M., 2015, MNRAS,446, 1765 Volonteri M., Perna R., 2005, MNRAS,358, 913
Madau P., Quataert E., 2004, ApJ, 606, L17 WadsleyJ.W.,StadelJ.,QuinnT.,2004,NewAstronomy,
Mayer L., Kazantzidis S., Madau P., Colpi M., Quinn T., 9, 137
Wadsley J., 2007, Science, 316, 1874 WursterJ., Thacker R.J., 2013, MNRAS,431, 2513
Menon H., Wesolowski L., Zheng G., Jetley P., Kale L.,
Quinn T., Governato F., 2015, Computational Astro- ThispaperhasbeentypesetfromaTEX/LATEXfileprepared
physicsand Cosmology, 2, 1 bythe author.
Moran E.C.,ShahinyanK.,SugarmanH.R.,VélezD.O.,
Eracleous M., 2014, AJ, 148, 136
Mortlock D.J. et al., 2011, Nature,474, 616
NavarroJ.F.,FrenkC.S.,WhiteS.D.M.,1996,ApJ,462,
563
Oh S., de Blok W. J. G., Walter F., Brinks E., Kennicutt
R.C., 2008, AJ, 136, 2761
Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kenni-
cutt,Jr. R.C., 2011, AJ, 141, 193
Okamoto T., Gao L., Theuns T., 2008, MNRAS,390, 920
Ostriker E. C., 1999, ApJ, 513, 252
Pontzen A., GovernatoF., 2012, ArXiv e-prints
Pontzen A., Roškar R., Stinson G., Woods R., 2013, pyn-
body: N-Body/SPH analysis for python. Astrophysics
SourceCode Library
Power C., Navarro J. F., Jenkins A., Frenk C. S., White
S.D.M.,SpringelV.,StadelJ.,QuinnT.,2003,MNRAS,
338, 14
Read J. I., Goerdt T., Moore B., Pontzen A.P., StadelJ.,
LakeG., 2006, MNRAS,373, 1451
©2015RAS,MNRAS000,1–8