Efficient Merger of Binary Supermassive Black Holes in Non-Axisymmetric Galaxies

0.26 MB·English

by  Peter Berczik

DRAFTVERSION5THFEBRUARY2008
PreprinttypesetusingLATEXstyleemulateapjv.6/22/04
EFFICIENTMERGEROFBINARYSUPERMASSIVEBLACKHOLESINNON-AXISYMMETRICGALAXIES
PETERBERCZIK1,2,3,DAVIDMERRITT1,RAINERSPURZEM2,HANS-PETERBISCHOF4
Draftversion5thFebruary2008
Abstract
Binarysupermassiveblackholesformnaturallyingalaxymergers,buttheirlong-termevolutionisuncertain.
In sphericalgalaxies, N-bodysimulationsshowthat binaryevolutionstalls at separationsmuchtoo largefor
significantemissionofgravitationalwaves(the“finalparsecproblem”). Here,wefollowthelong-termevolu-
tionofamassivebinaryinmorerealistic,triaxialandrotatinggalaxymodels. Wefindthatthebinarydoesnot
stall. ThebinaryhardeningratesthatweobservearesufficienttoallowcompletecoalescenceofbinarySBHs
6
in10Gyrorless,evenintheabsenceofcollisionalloss-conerefillingorgas-dynamicaltorques,thusproviding
0
0 apotentialsolutiontothefinalparsecproblem.
2 Subjectheadings:

n
a
J 1. INTRODUCTION nearly spherical galaxy models. But it has been suggested
0 When two galaxies containing supermassive black holes (Merritt&Poon2004) thatbinaryhardeningmightbemuch
moreefficientinnon-axisymmetricgalaxiesduetothequali-
3 (SBHs) merge, a binarySBH formsat the center of the new
tativelydifferentcharacterofthestellarorbits. Here,wetest
galaxy. The two SBHs can eventually coalesce, but only
thatsuggestionbycarryingoutthefirstN-bodysimulationsof
1 after stellar- or gas-dynamical processes bring them close
v enough together (<10 2 pc) that gravitational radiation is massive binaries in strongly non-axisymmetric galaxy mod-
8 emitted. There is∼stron−g circumstantial evidence that rapid els. We find that the hardeningrate is independentof N for
9 coalescence is the norm. For instance, no binary SBH has particle numbersup to at least 0.4 106. To the extent that
×
6 ever been unambiguouslyobserved(Merritt&Milosavljevic ourgalaxymodelsaresimilartorealmergerremnants,these
1 2005). Furthermore, in a galaxy containing an uncoalesced resultsimplythatbinarySBHscanefficientlyhardenthrough
0 binary, mergers would eventually bring a third SBH into purely stellar-dynamical interactions in many galaxies, thus
6 the nucleus, precipitating a gravitational slingshot interac- providingaplausiblesolutiontothefinalparsecproblem.
0
tion that would eject one or more of the SBHs from the nu-
/ 2. METHOD
h cleus(Mikkola&Valtonen1990;Volonteri,Haardt&Madau
p 2003). This could produce off-center SBHs, and could Our N-body models were generated from the phase-space
- also weakenthe tightcorrelationsthatare observedbetween distributionfunction
o
str SGBraHhammaesstaaln.d20g0a1la;xMyaprcroopneir&tieHsu(Fnter2r0a0re3s)e. &Merritt 2000; f(E,Lz)=const×(cid:16)e−b E−1(cid:17)e−bW 0Lz (1)
:a rapUidnllyesbsrtihnegsbtihnearsymmalalsesrrSaBtiHoiisnteoxtaredmiseta,ndcyenamGicµa/lsfr2icfrtioomn (Lagoute&Longaretti1996). Here E =v2/2+F is the en-
v ∼ ergyperunitmassofastar,F (v ,z)isthegravitationalpoten-
thelargerSBH,whereµ M M /(M +M )isthebinaryre-
Xi duced mass and s is the≡1D1vel2ocity1dispe2rsion of the stars. tialinthemeridionalplane,andLz istheangularmomentum
per unit mass in the direction of the symmetry (z) axis; the
At this separation – of order 1 pc – the two SBHs begin to
r
potentialissettozeroattheradiuswherethedensityfallsto
a actlikea“hard”binary,ejectingpassingstarswithvelocities
zero. The quantity in parentheses on the right hand side of
largeenoughtoremovethemfromthenucleus. N-bodysim-
equation(1)istheenergy-dependentKing(1966)distribution
ulations (Makino&Funato 2004; Szell,Merritt&Mikkola
function. Theadditional,angular-momentum-dependentfac-
2005; Berczik,Merritt&Spurzem 2005) show that contin-
torhastheeffectofflatteningthemodelsandsimultaneously
uedhardeningofthebinarytakesplaceataratethatdepends
givingthemanetrotationaboutthezaxis. Thedegreeofflat-
stronglyonthenumberNof“star”particlesusedinthesimu-
teningcanbespecifiedviathedimensionlessrotationparam-
lation. AsN increases,thehardeningratefalls,asexpectedif
thebinary’slossconeisrepopulatedbystar-stargravitational eterw 0 9/(4p Gr 0)W 0, with r 0 thecentralmassdensity
encounters (Yu 2002; Milosavljevic&Merritt 2003). When ofthega≡laxpy.Theparameterb determinesthecentralconcen-
extrapolatedto the much larger N of real galaxies, these re- trationofthemodel;itsvaluewaschosensuchthatthespher-
sults suggestthat binaryevolutionwouldgenerallystall (the icalisotropicmodelgeneratedfromequation(1)hadadimen-
“finalparsecproblem”). sionlesscentralconcentrationW0=6(King1966). Hereand
To date, N-bodysimulationsof thelong-termevolutionof belowweadoptstandardN-bodyunits, i.e. the gravitational
binary SBHs have only been carried out using spherical or constantandtotalmassofthegalaxyareone,andthegalaxy’s
energyis 1/4.
−
1 Department ofPhysics, Rochester Institute ofTechnology, Rochester, ApairofmassiveparticlesrepresentingthetwoSBHswere
NY14623,USA introduced into the models at time t =0. The two particles
2AstronomischesRechen-Institut,ZentrumfürAstronomie,Monchhofs- weregivenequalmasses,M =M M /2,andwereplaced
trasse12-14,69120Heidelberg,Germany oncoplanar,circularorbitsa1tdista2n≡ces •0.3fromthegalaxy
3 Main Astronomical Observatory, National Academy of Sciences of ±
Ukraine,ZabolotnohoStreet,27,Kiev,Ukraine,03680 centerintheequatorialplane. Inmostofthesimulationsde-
4DepartmentofComputerScience,RochesterInstituteofTechnology,102 scribedbelow,M =0.04. Thisisratherlargerthanthetyp-
LombMemorialDrive,Rochester,NY14623 ical ratio, 1 •10 3 (Merritt&Ferrarese 2001), observed
−
∼ ×
2 Bercziketal.
between SBH mass and galaxy mass; such a large mass for
the SBH particles was chosen in order to minimize the rate
of relaxation-driven loss cone refilling, which occurs more
rapidly for smaller M (Berczik,Merritt&Spurzem 2005),
and to come as close•as possible to the “empty loss cone”
regimethatcharacterizesreal(axisymmetric)galaxies. Inor-
dertoestimatethedependenceofthebinarydecayrateonM ,
wecarriedoutalimitedsetofadditionalsimulationswithdi•f-
ferentvaluesofM ,asdescribedbelow.
Integrations of•the particle equations of motion were
carried outusing a high-accuracy,direct-summationN-body
code (Berczik,Merritt&Spurzem 2005) on two paral-
lel supercomputers incorporating special-pupose GRAPE
(Fukushige,Makino&Kawai 2005) accelerator boards:
gravitySimulator5 and GRACE.6 Integration parameters
weresimilartothoseadoptedinBerczik,Merritt&Spurzem
(2005)andwereferthereadertothatpaperfordetailsabout
the performance of the code. Integrations were carried out
for various N in the range 0.025 N 0.4 106 and for
various values of the galaxy rotat≤ion p≤arame×ter w in the
0
range0 w 1.8. In addition, each modelwas integrated
0
≤ ≤
with two choices for the orientation of the binary’s angular
momentum,eitherparalleltothatofthegalaxy(“prograde”)
orcountertoit(“retrograde”).
3. RESULTS
After the two SBH particles come close enough together
to form a bound pair, the parameters of their relative Ke- FIG. 1.— Evolution ofthe binary inverse semi-major axis, 1/a, inN-
plerian orbit can be computed. Figure 1a shows the evo- bodysimulations withvarious N. (a)Spherical, nonrotating galaxy model
lution of 1/a, the binary inverse semi-major axis, in a set (w 0=0). (b)Flattened, rotatinggalaxymodel(w 0=1.8). Att 10,this
of simulations with w = 0 and various N. These spheri- modelformsatriaxialbar(cf.Figure2). ≈
0
cal, non-rotatingmodelsare verysimilar to the modelscon-
sidered in Berczik,Merritt&Spurzem (2005), and the bi-
naryevolutionfoundhereexhibitsthesamestrongN depen-
dence that was observed in that study: the hardening rate,
s(t) (d/dt)(1/a),isapproximatelyconstantwithtimeand
≡
decreasesroughlyasN 1. Thisbehaviorhasbeendescribed
−
quantitatively(Milosavljevic&Merritt2003) onthebasisof
loss-conetheory:starsejectedbythebinaryarereplacedina
timethatscalesasthetwo-bodyrelaxationtime,andthelatter
increasesroughlyasN inagalaxyoffixedmassandsize.
When the rotationparameterw is increasedto 0.6, the
0
∼
initiallyaxisymmetricmodelsbecomeunstabletotheforma-
tion of a bar, yielding a slowly-tumbling, triaxial spheroid.
Figure 2 illustrates the instability via snapshots of the w =
0
1.8 model integration. This model is moderately flattened
initially, with mean short-to-long axis ratio of 0.46, and
∼
stronglyrotating,withroughly40%ofthetotalkineticenergy
intheformofstreamingmotion.Moviesbasedonthesimula-
tion7revealthatthetwoSBHparticlesinitiallycometogether
byfallinginwardalongthebarbeforeformingaboundpair.
The long-term behavior of the binary (Figure 1b) is strik-
ingly different in this rotating model than in the spherical FIG.2.— Snapshotsatfourtimes(t=0,5,10,30)oftheparticlepositions,
model: not only is the hardeningrate high, but more signif- projectedontothe(x,y)plane,inanN-bodyintegrationwithw 0=1.8and
N=200k. TheSBHparticles areindicatedingreen. Theevolutionofthe
icantly, it shows no systematic dependenceon particle num- binarysemi-majoraxisinthisintegrationisshowninFigure1b.
ber. Infact,thetwosimulationsofthew =1.8modelwith
0
largestN (200k and 400k)exhibitalmostidenticalevolution
ofthebinary.
tegrationswere carriedoutforvarious(w ,N) andforM =
Todeterminethedependenceofthebinary’sevolutionrate 0
0.04. Axisratiosofthegalaxywerecomputedusingitsm•o-
on the propertiesof the galaxy model, a suite of N-body in-
ment of inertia tensor, as described in Dubinski&Carlberg
5http://www.cs.rit.edu/grapecluster/clusterInfo/grapeClusterInfo.shtml (1991). The results are summarized in Figure 3. After a
6http://www.ari.uni-heidelberg.de/grace strongbar formsatt 10 in the unstablemodels, itevolves
≈
7http://www.cs.rit.edu/grapecluster/BinaryEvolution gradually toward rounder shapes, but the system maintains
EfficientMergerofBinaryBlackHoles 3
FIG. 3.—(a)Evolutionoftheintermediate-to-long axisratiob/aofthegalaxyinN-bodyintegrations withN=200kandvarious values oftherotation
parameterw 0;M1=M2=0.02. Forw 0>0.6,thegalaxyisunstabletonon-axisymmetricdeformationsandformsaslowly-tumblingtriaxialspheroid. The
lengthcoftheshortaxisisnearlyconstantwithtimeandisdeterminedbytheinitialflatteningofthemodel,i.e.byw 0.Thethinblacklineshowstheresultofan
integrationthatomittedthetwoSBHparticles. (b)Hardeningrateofthebinarymeasuredatt=150. Thevalueofw 0isindicatedviathekeyonthelowerleft;
blacksymbolsare“prograde”integrations,i.e. thebinaryrevolvesinthesamesenseasthegalaxy,andredsymbolsare“retrograde”integrations. Thedashed
linesshowstheN 1dependencethatischaracteristicofacollisionally-resuppliedlosscone.
−
a significant triaxiality until the end of the simulation. This evolutionwouldbe substantiallymorerapidthanimpliedby
slow evolution appears similar to that in the simulations Figure3binthecaseM 10 3M .
− gal
of Theis&Spurzem(1999), where two-bodyrelaxationwas •≈
identifiedasthedrivingmechanism. Thepresenceofamas-
4. DYNAMICALINTERPRETATION
sive binary in our simulations might also tend to destroy
the triaxiality (Athanassoula,Lambert&Dehnen 2005), al- A hardeningrate that is independentof N implies a colli-
thoughan integrationexcludingthe binary showeda similar sionless, i.e. relaxation-independent, mode of loss-cone re-
degreeofevolution;seeFigure3a. filling. Just such a mode is expected in triaxial galaxies:
Thesteep, N 1 dependenceofthebinaryhardeningrate the lack of an axis of symmetry implies that stellar orbits
−
inthespherica∼lmodelchangestoanessentiallyconstanthard- neednotconserveanycomponentoftheangularmomentum,
ening rate for w 1.2 (Figure 3b). Even the w = 0.6 hencetheycanpassarbitrarilyclosetothecenterafterafinite
0 0
model– which, Fig≥ure3a suggests, is close to marginalsta- timeandinteractwithacentralobject(Norman&Silk1983;
bility–yieldedahardeningratethatwassubstantiallylarger Gerhard&Binney1985).
than in the spherical models, consistent with suggestions A fullderivationof the expectedrate of supply of stars to
(Merritt&Poon 2004) that even slight departures from ax- thebinaryinthesemodelswouldbeverydifficult,butwecan
isymmetry could significantly influence the binary’s evolu- do an approximate calculation. The rate per unit of orbital
tion.Thehardeningratewasfoundnottodependontheinitial energyatwhichcentrophilicorbitssupplymass(i.e. stars)to
sense(progradevs.retrograde)ofthebinaryorbit(Figure3b). aregionofradiusrt atthecenterofagalaxyis
Since M /M =0.04 is considerably larger than the ra-
tio 1 1•0 3goablservedinrealgalaxies(Merritt&Ferrarese M˙(E)dE=rtA(E)Mc(E)dE (2)
−
∼ ×
2001; Marconi&Hunt 2003), we carried out an additional
whereA(E)distherateatwhichasinglestaronacentrophilic
setofsimulationsinordertoevaluatetheM -dependenceof
(e.g.boxorchaotic)orbitofenergyEexperiencesnear-center
thebinaryhardeningrate. Theseintegrations•usedw =1.8,
0 passages with pericenter distances d, and M (E)dE the
N =0.1 106, and 0.01 M 0.08. We found that the ≤ c
mass in stars on centrophilic orbits with energies from E to
hardening×rate increased w≤ith•d≤ecreasing M . At t = 100,
E+dE (Merritt&Poon2004). Settingr =Ka,withK 1,
the hardening rates were s=(20.0,13.0,8.2•,3.4) for M = t ≈
givesthemass fluxinto thebinary’ssphereofinfluence; the
(0.01,0.02,0.04,0.08). These results should be interpre•ted
impliedhardeningrateis(Berczik,Merritt&Spurzem2005)
with caution since we did not vary N and therefore can not
state with certainty whether the s values are N-independent.
d 1 2 C
Thelinearsizeofthebinary’slossconeisproportionaltoM s h i M˙(E)dE. (3)
makingiteasierfortwo-bodyscatteringtoaffecttheharden•- ≡ dt(cid:18)a(cid:19)≈ aM Z
ing rate as M is decreased. But the M =0.04 integrations •
areclearlyno•tinthecollisionally-repopu•latedregime(Figure Here, C 1.25 is the average value of the dimensionless
1b)andsothedifferencesthatweobservebetweenthehard- energyhciha≈nge during a single star-binary encounter, C
ening rates for M =0.04 and 0.08 are likely to be robust. [M /2m ](D E/E). ≡
⋆
Basedontheseres•ults,itisreasonabletoconcludethatbinary O•urN-bodymodelshavedensityr r 2 beyondthecore
−
∼
4 Bercziketal.
radiusr 0.25.Inar(cid:181) r 2galaxy, is a conservativeinterpretationsince (1) forreasonablescal-
c −
≈
ings of our galaxy model to real galaxies (e.g. total mass
A(E)≈rs 2e−(E−Eh)/s 2, Mc(E)= fc(E)×2√96rGhe(E−Eh)/2s 2 =25010in11MN-⊙b,odhyaluf-nmitasscsorraredsiupson=ds1t0o3p<c)1, aGnyre;la(p2s)etdhetibmineaoryf
h
(4) is continuing to harden at the final ti∼me-step in our simula-
withr =GM /s 2,E =F (r ),and f (E)thefractionofthe tions(Figure1b);(3)ourexperimentswithdifferentM found
h h h c
orbitsatenergy• Ethatarecentrophilic(Merritt&Poon2004). s M 1,implyingsubstantiallymorerapidhardening• inthe
−
Theimpliedbinaryhardeningrateis ca∼seM• /M 10 3;(4)thebinaryhadnonzeroeccentricity
gal −
inours•imulati≈ons.Inaddition,gasisasignificantcomponent
4√6 C Kf s
s≈ 9 h sir2 cZ e−(E−Eh)/2s 2dE≈2.5fcr2. (5) oofftdhieskbignaalrayxwieosualnddb,einacmcealenryatmederbgyergsa,st-hdeynfianmalichaalrtdoerqnuinegs
h h
(Escalaetal.2005;Dottietal.2005).
Here fc isanenergy-weighted,meanfractionofcentrophilic Oursimulationsofbinaryevolutionaresubstantiallymore
orbits,andthelowerintegrationlimitwassettoEh;thelatter realistic than existing ones based on spherical or nearly-
canonlybeapproximatesincethetruedensityofourgalaxy spherical galaxy models. Even more realistic simulations,
modelsdepartsfromr−2atr<rc rh. Substitutings 0.47 which follow both the early and late stages of a merger
≈ ≈
and rh 0.18 from the galaxy models gives s 40fc. By between two galaxies, are probably beyond the capabili-
≈ ≈
comparison,thehardeningratesintheN-bodymodelsreach ties of current algorithms and hardware due to the need to
a peak value at t 20 of s 16, consistent with the de- accurately treat both large ( 10 kpc) and small ( 0.01
≈ ≈ ∼ ∼
rived expression if f 1/2. The gradual drop observed in pc) spatial scales. However, our galaxy models (slowly-
c
≈
the hardening rate at later times, s(t) 16. 5.2ln(t/20), tumbling triaxial spheroids) are similar to those produced
≈ −
20 t 250, suggests that the numberof centrophilicstars in full merger simulations (Bournaud,Jog&Combes 2005;
≤ ≤
isbecomingsmaller,duetodepletionbythebinaryandtothe Naab,Khochfar&Burkert2006), suggestingthatourresults
gradualchangeinthegalaxy’sshape. for the long-termevolution of the binary are probablyfairly
Takenatfacevalue,equation(5)impliess(cid:181) M 2;however genericinspiteoftheratherartificialinitialconditions.
−
forsmallM ,r r andtheassumptionthatr •r 2, r>r Uncertaintiesabouttheresolutionofthe“finalparsecprob-
h c − h
breaksdow•n. In≪any case, the observed depend∼enceof s on lem” have been a major impediment to predicting the fre-
M isslightlyweaker,s M 1 (§3). quencyofSBHmergersingalacticnuclei,andhencetocom-
−
• ∼ • puting event rates for proposed gravitational wave interfer-
5. IMPLICATIONS ometerslike LISA.8 Ifbinarycoalescenceratesare assumed
Thetimescaleforgravitationalwaveemissionbyabinary tobesimilartogalaxymergerrates,gravitationalwaveevents
blackholeis(Peters1964) integratedover the observable universe could be as frequent
as102 yr 1 (Haehnelt1994;Sesanaetal.2004). Ourresults,
−
5 Gµ3c5 a 4 combinedwiththeindirectevidencethatbinarySBHcoales-
t = . (6)
gr 164F(e) s 8M2 (cid:18)a (cid:19) cenceisefficient,suggestthatsuchhigheventratesshouldbe
h
• takenseriously.
Here M =M +M , µ=M M /M2 is the reducedmass of
1 2 1 2
the bina•ry, s is the 1D central velo•city dispersion of stars
in the nucleus, and ah Gµ/4s 2 is the semi-major axis of WethankM.MilosavljevicandS.Harfstforcommentson
≡
the binary when it first becomes “hard,” i.e. tightly bound; the manuscript. This work was supported by grants AST-
the factor F(e) depends on the binary’s orbital eccentricity 0206031, AST-0420920 and AST-0437519 from the NSF,
andF(0)=1. Inorderthatgravitational-wave-drivencoales- grantNNG04GJ48Gfrom NASA, grantHST-AR-09519.01-
cencetakeplaceinlessthan1010 yr,anequal-mass,circular- A from STScI, grant I/80 041 GRACE from the Volkswa-
orbit binary with M = 108M must first reach a separa- genFoundation,bySFB439ofDeutscheForschungsgemein-
tion a<0.05ah (Mer•ritt&Mil⊙osavljevic 2005). This is just schaft,andbyINTASgrantIA-03-59-11. WethanktheCen-
achieve∼d in our simulations: a 1.1 10 2, and the final ter for the Advancementof the Study of Cyberinfrastructure
h −
value of a in the bar-unstable m≈odels×is 6 10 4. This atRITfortheirsupport.
−
∼ ×
8http://lisa.jpl.nasa.gov/
REFERENCES
Athanassoula,E.,Lambert,J.C.,&Dehnen,W.1005,MNRAS,363,496 Iwasawa,M.,Funato,Y.,&Makino,J.2005,ArXivAstrophysicse-prints,
Begelman,M.Blandford,R.D.&Rees,M.J.1980,Nature,287,307 arXiv:astro-ph/0511391
Berczik,P.,Merritt,D.,&Spurzem,R.,ApJ,633,680 King,I.R.1966,AJ,71,64
Bournaud,F.,Jog,C.J.,&Combes,F.2005,AAp,437,69 Kormendy,J.,&Richstone,D.1995,ARAA,33,581
Dotti, M., Colpi, M., & Haardt, F. 2005, ArXiv Astrophysics e-prints, Lagoute,C.,&Longaretti,P.-Y.1996,AAp,308,441
arXiv:astro-ph/0509813 Makino,J.,&Funato,Y.2004,ApJ,602,93
Dubinski,J.,&Carlberg,R.G.1991,ApJ,378,496 Marconi,A.&Hunt,L.K.2003,AJ,589,L21
Escala,A.,Larson,R.B.,Coppi,P.S.,&Mardones,D.2005,ApJ,630,152 Merritt,D.,&Ferrarese,L.2001,MNRAS,320,L30
Ferrarese,L.&Merritt,D.,ApJ,539,L9 Merritt,D.&Milosavljevic,M.2005,Massiveblackholebinaryevolution,
Fukushige,T.,Makino,J.,&Kawai,A.2005,PASJ,57,1009 LivingReviewsinRelativity,http://relativity.livingreviews.org
Gebhardt,K.,etal.2000,ApJ,539,L13 Merritt,D.,&Poon,M.Y.2004,ApJ,606,788
Gerhard,O.E.,&Binney,J.1985,MNRAS,216,467 Mikkola,S.,&Valtonen,M.J.1990,ApJ,348,412
Graham,A.W.,Erwin,P.,Caon,N.,&Trujillo,I.2001,ApJ,563,L11 Milosavljevic´,M.,&Merritt,D.2003,ApJ,596,860
Haehnelt,M.G.1994,MNRAS,269,199 Naab,T.,Khochfar,S.,&Burkert,A.2006,ApJ,636,L81
Norman,C.,&Silk,J.1983,ApJ,266,502
EfficientMergerofBinaryBlackHoles 5
Peters,P.C.1964,Phys.Rev.,136,1224 Volonteri,M.,Haardt,F.,&Madau,P.,ApJ,582,559
Sesana,A.,Haardt,F.,Madau,P.,&Volonteri,M.2004,ApJ,611,623 Yu,Q.2002,Mon.Not.R.Astron.Soc.,331,935
Szell,A.,Merritt,D.&Mikkola,S.2004,AnnNYAcadSci,1045,225
Theis,C.,&Spurzem,R.1999,A&A,341,361
Tremaine,S.,etal.2002,ApJ,574,740