Million-Body Star Cluster Simulations: Comparisons between Monte Carlo and Direct $N$-body

1.7 MB·

by  Carl L. Rodriguez

MNRAS000,1–9(2016) Preprint19January2016 CompiledusingMNRASLATEXstylefilev3.0
Million-Body Star Cluster Simulations: Comparisons between
Monte Carlo and Direct N-body
Carl L. Rodriguez1,2(cid:63), Meagan Morscher1,2, Long Wang3,4, Sourav Chatterjee1,2,
Frederic A. Rasio1,2, Rainer Spurzem3,5,6,7
1CenterforInterdisciplinaryExplorationandResearchinAstrophysics(CIERA),NorthwesternUniversity,Evanston,IL,USA
2DepartmentofPhysicsandAstronomy,NorthwesternUniversity,Evanston,IL,USA
3KavliInstituteforAstronomyandAstrophysics,PekingUniversity,YiheyuanLu5,HaidianQu,100871,Beijing,China
6 4DepartmentofAstronomy,SchoolofPhysics,PekingUniversity,YiheyuanLu5,HaidianQu,100871,Beijing,China
1 5NationalAstronomicalObservatoriesandKeyLaboratoryofComputationalAstrophysics,ChineseAcademyofSciences,
0 20ADatunRd.,ChaoyangDistrict,100012,Beijing,China
2 6KeyLaboratoryofFrontiersinTheoreticalPhysics,InstituteofTheoreticalPhysics,ChineseAcademyofSciences,Beijing,100190,China
7AstronomischesRechen-Institut,ZentrumfürAstronomie,UniversityofHeidelberg,Mönchhofstrasse12-14,69120,Heidelberg,Germany
n
a
J
Lastupdated19January2016
6
1

] ABSTRACT
M We present the first detailed comparison between million-body globular cluster simula-
tions computed with a Hénon-type Monte Carlo code, CMC, and a direct N-body code,
I
. NBODY6++GPU.Bothsimulationsstartfromanidenticalclustermodelwith106 particles,
h
and include all of the relevant physics needed to treat the system in a highly realistic way.
p
Withthetwocodes“frozen"(nofine-tuningofanyfreeparametersorinternalalgorithmsof
-
o the codes) we find excellent agreement in the overall evolution of the two models. Further-
r more, we find that in both models, large numbers of stellar-mass black holes (> 1000) are
t
s retainedfor12Gyr.Thus,theveryaccuratedirect N-bodyapproachconfirmsrecentpredic-
a tionsthatblackholescanberetainedinpresent-day,oldglobularclusters.Wefindonlyminor
[
disagreements between the two models and attribute these to the small-N dynamics driving

1 theevolutionoftheclustercoreforwhichtheMonteCarloassumptionsarelessideal.Based
v ontheoverwhelminggeneralagreementbetweenthetwomodelscomputedusingthesevastly
7 differenttechniques,weconcludethatourMonteCarloapproach,whichismoreapproximate,
2 butdramaticallyfastercomparedtothedirectN-body,iscapableofproducingaveryaccurate
2
description of the long-term evolution of massive globular clusters even when the clusters
4
containlargepopulationsofstellar-massblackholes.
0
. Keywords: binaries:close—globularclusters:general—Gravitationalwaves—Methods:
1
0 numerical—Stars:blackholes—Stars:kinematicsanddynamics
6
1
:
v
i 1 INTRODUCTION enablesthemodelingofGCsystemswithrealisticnumberofstars
X
andalltherelevantphysics.
r In the last few years, our understanding of the dynamical evolu-
a tionofglobularclusters(GC),especially,asaresultofthedetailed Sincetheearlytheoreticalworkpredictingtherapidejection
ofBHsfromclusters,severalgroupshaveperformedevermorere-
dynamicalevolutionandfateofthelargenumbersofblackholes
alisticevolutionarysimulations.Whileseveralofthefirstattempts
(BHs)thatareboundtoformintheselarge-N clusters,hasshifted
confirmedthetheoreticalpredictionforcompleteevaporation,the
significantly.Whileitwasoncethoughtthatpresent-dayoldGCs
most recent models have predicted that at least some BHs, and
shouldhaveatmostacoupleofBHsremaining,morerecentstud-
ieshaveshownthatifBH-formationkicksarenotsufficientlyhigh likelymany,mayremaininoldGCstoday.Theseresultsarecom-
ingatanexcitingtimewhenstellarBHsarebeingdiscoveredata
toejectallBHsfromtheGCs,asignificantfractionoftheformed
rapidpaceinGCsintheMilkyWayandinotherGalaxiesthrough
BHsareretaineduptothetypicaloldages(∼12Gyr)oftheGCs
X-ray and radio surveys (Maccarone et al. 2007; Barnard et al.
(Breen&Heggie2013;Heggie&Giersz2014;Mackeyetal.2008;
2011;Maccaroneetal.2011;Shihetal.2010;Straderetal.2012;
Sippel&Hurley2013;Morscheretal.2015).Thisnewperspective
Chomiuketal.2013).
has been driven by major advances in parallel computing, which
We are also beginning to understand why these new results
areatoddswiththeearlytheoreticalpredictionofrapidBHevap-
(cid:63) Contacte-mail:[email protected] oration. The original argument assumed that the BHs would be-
(cid:13)c 2016TheAuthors
2 C.L.Rodriguezetal.
comeSpitzerunstableduetotheirrapidmass-segregation,dynam- meansitcanbeusedtoexploretheparameterspaceofinitialcon-
icallydecouplingfromthecluster.Oncedecoupled,theBHswould ditions and test the robustness of our results. This feature sets it
evolveasanisolatedcluster,withanevaporationtimescaleof∼ 1 apartfromdirectN-bodysimulations,whicharesignificantlymore
Gyr.Whiletherewereafewtheoreticalhintsfromsimulationsthat computationallyexpensive,andthereforeareusuallyrestrictedto
atleastsomeBHscouldberetainedfor12Gyr,therewasnogood N (cid:46) 105.ThespeedoftheMCtechniquecomesataprice,how-
explanationforhowthiscouldbepossible,andwheretheoldpre- ever,inthattwo-bodyrelaxation,aswellasotherdynamicalpro-
diction breaks down. Using two-component direct N-body mod- cesses,aretreatedinanapproximatewayviaasingleinteraction
els, Breen & Heggie (2013) showed that after the BHs began to perpairofstarspertimestep,wherethetimestepischosentobe
segregate, their rate of energy generation was controlled by the a small fraction of the relaxation timescale. Physics occuring on
rate of energy flow through the cluster as a whole, and so it was muchshortertimescales,suchasthedynamicaltimescale,cannot
setbytherelaxationtimescaleoftheentirecluster,notsimplyby beresolvedaccuratelywiththeMCmethod,butishandlednatu-
thatofthesmallsubpopulationofBHs.Recently,Morscheretal. rallywiththedirectN-bodytechnique.
(2015)(hereafterMOR15)presentedagridof 42detailedMonte OurMCcodehasbeencomparedextensivelytodirectN-body
Carlo models for realistic star clusters containing substantial ini- simulationswheneverpossible,andhasshownexcellentagreement
tialpopulationsofstellarBHs.Thisstudywasthefirsttopresenta (Joshietal.2000,2001;Fregeauetal.2003;Fregeau&Rasio2007;
large number of realistic simulations with a range of initial clus- Chatterjeeetal.2010;Umbreitetal.2012).However,thetypesof
ter masses (from ∼ 105 −106M(cid:12)) as well as variation in other clusterswearenowcapableofsimulatingarequitedifferentthan
importantparameters,suchasthevirialradius,thatalsoincludes thosewehavesimulated(andthustested)inthepast.Inparticular,
alloftherelevantphysicsrequiredtodescribethesesystemsaccu- theinclusionoflargepopulationsofBHswithabroadspectrumof
rately.Startingfromstandardassumptionsfortheinitialconditions masseshasasignificanteffectonclusterevolution,producingclus-
of MW-like clusters and BH formation processes, MOR15 found ters at 12 Gyr with very different properties than models without
thatnearlyallofthemodelsdidindeedretainsignificantfractions suchobjects(cf.Chatterjeeetal.2010,inwhichtheBHmassspec-
oftheirBHsuptotheendofthesimulationat12Gyr.Theretained trumwastruncatedatjustabove2M ).Forexample,thedeepBH-
(cid:12)
BHsheatedthefullcluster,leadingtolargefinalcoresizes;how- driven core oscillations that occur ubiquitously in our new CMC
ever,themostcompactinitialclustermodelswereabletoejectthe modelsareanewphenomenonthatwehavenotseenbeforeinpre-
majorityoftheirBHsby12Gyr,causingtheircorestocontract(the viousmodels.Duringthesecollapses,thephysicsmaybegoverned
so-calledsecondcorecollapse)tosizessimilartothoseobservedin byaverysmallnumberofobjects(∼ 10)thatarepartiallydecou-
MWGCs.Workisinprogresstostudytheeffectofvaryingtheini- pledfromtherestofthecluster.Iftherelaxationtimescaleofthe
tialmassfunctionandtheBHbirthkickdistribution,ontheinitial subsetofstarsismuchsmallerthantherelaxationtimescaleofthe
BHpopulationsinclusters,andultimatelytheirimpactoncluster clusterasawhole,thenitsevolutionisperhapsnotbeingcomputed
dynamicsandlong-termBHretention. accuratelybyourapproximateMCscheme.Furthermore,itispos-
Perhapsmostinterestingly,MOR15hasrevealedanewtheo- siblethatourcrudeprescriptionforformingthree-bodybinariesis
reticalunderstandingofthecomplexdynamicalinterplaybetween notcapableofproperlycapturingthephysicsofthisdynamically
BHsandclusters,andhowitmaybepossibleforclusterstoretain importantprocess.
large numbers of BHs for many Gyr. These MC models showed These limitations of our MC calculations could potentially
that,inrealisticsystensm,theSpitzerinstabilitydoesnot involve have an impact on our results regarding the dynamical evolution
alloftheBHsinthecluster.Rather,theBHspowercoreoscilla- ofBHsindenseclusters,includinglong-termBHretentionandthe
tionsduringwhichonlyasmallsubsetoftheBHssegregatefrom resultingheatingthatproducespuffyclustercores.Fortheserea-
the cluster to form a deep cusp, which then promptly re-expands sons, it would be desirable to re-test the results produced by the
uponformationofthree-bodybinaries.Thus,asaresultoftheen- latestversionofourMCcodewiththebestcurrently-availabledi-
ergy generated by the dynamics of the BHs, the BH interaction rectN-bodysimulations.Forourspecificinterests,therehavepre-
rate, and therefore also the evaporation rate, is kept much lower viouslybeennosuitablelarge-N directN-bodysimulationsavail-
thanpreviouslythought,makingitpossibleforstarclusterstore- abletowhichwecouldcompareourlarge-NMCsimulations.Re-
tainsignificantnumbersofBHsfor12Gyr.Itseemsthatweshould cently,however,Wangetal.2015developedanewdirectN-body
notexpecttheBHstosuccombtotheSpitzerInstability;mostof codethatemployshybridparallelizationmethodstospeedupthe
theBHs,infact,alwaysremainspreadthroughoutthecluster,far NBODY6++ code (Spurzem 1999; Spurzem et al. 2008), which
fromthecentralcusp.ThisnewunderstandingofBHdynamicsin is considered to be the gold standard in direct N-body modeling.
starclustershasprovidedatheoreticalbasisforexplainingthere- Thisnewcode,NBODY6++GPU,iscapableofmodelingclusters
centdiscoveriesofseveralBHcandidatesinoldGCs. with106 starswithinayear(Wangetal.2015;Wangetal.2016).
Asalways,wewouldliketobeabletocompareourresultsto Thisprovidesanexcellentopportunityforustotestourpredictions
simulationsdonewithothercodes,especiallythosethatusediffer- regardingtheevolutionofclusterswithBHsthroughadirectcom-
entmodelingtechniques.ThedirectN-bodytechniqueisthemost parison.
accurate and assumption-free method for modeling the dynamics Herewepresentacomparisonbetweenmodelsgeneratedwith
of star clusters. It resolves physics on the dynamical timescale, our MC code (the MC model) and with NBODY6++GPU (the NB
which means it can inprinciple be used to model clusters of any model). We present this as an honest comparison using what are
size,becauseitisvalidevenforclusterswithrelaxationtimescales believedtobethebestversionsofthetwocodes,andwedonotdo
not much longer than their dynamical timescales. This is in con- anyfine-tuningoffreeparameters.Thisrestofthispaperisorga-
trast to the MC technique, whose assumptions rely on the relax- nizedasfollows.InSection2wedescribethetwocodesinmore
ation timescale being long compared to the dynamical timescale. detailandgivetheinitialconditionsforthemodelthatisthefocus
ForlargeenoughN,however,theMCtechniqueusedinourstudy ofthiscomparison.InSection3wepresentandcomparetheresults
is an extremely powerful tool that is capable of computing many ofthetwosimulations,includingtheglobalstructuralevolutionas
large-Ndynamicalclustermodelsinashortamountoftime,which wellastheevolutionoftheBHpopulations.Wediscusssomeof
MNRAS000,1–9(2016)
Million-BodyComparisonsbetweenMonteCarloandDirect N-body 3
thedifferencesbetweenthemodelsanduncertaintiesandcompare ontheamountoffallbackmaterial,andbothcodeshavealsoup-
tootherstudiesinSection4.FinallyinSection5wesummarizeour datedtheirbirthkickprescriptionstodependonfallbackaccording
resultsandstateourconclusions. toBelczynskietal.(2002).ThekicksaredrawnfromaMaxwellian
velocity distribution with σ = 265 km s−1, then reduced propor-
tionallytothemassofthematerialthatfallsbackontothenewly-
formedcompactobject.
2 NUMERICALSETUP
Whilethree-bodybinaryformationoccursnaturallyindirect
WecomparetheresultsofasmallsetofMCmodelstoasingledi- N-body calculations, the MC technique relies on pair-wise inter-
rectN-bodysimulation,allstartingwithidenticalinitialconditions. actionsbetweentwoneighboringparticles,precludingthedynam-
Thecomparisonmodelisamillion-particledirectN-bodysimula- ical formation of a binary from encounters involving three parti-
tion that has been computed with NBODY6++GPU (Wang et al. cles. Therefore we opt to use a simple analytic prescription that
2015;Wangetal.2016),anewlydevelopedoptimizedversionof reliesonanestimateofthelocalprobabilityofbinaryformation,
NBODY6++withimprovedhybridparallelizationmethods(MPI, asdescribedindetailinMOR15.IntheMCcode,weusetheex-
GPU,OpenMPandAVX/SSE). actphysical assumptionsand parametersas presentedin MOR15
This new code combines the MPI parallelized NBODY6++ (e.g.,thoseassociatedwiththree-bodybinaryformation;seeSec-
(Spurzem1999;Hemsendorf,Khalisi,Omarov&Spurzem2003) tion2.2,Equations1and2fromMOR15).Onedifferenceisthatthe
withtheGPUandAVX/SSElibrariesofNitadori&Aarseth2012 MCcodeisnotabletotreatstabletriplesystemsthatformthrough
to significantly accelerate the direct N-body integratcion. They four-body(binary-binary)encounters.Thesesystemsarethusbro-
havealsomadeimprovementstospeedupthetime-stepschedul- kenuponformation.However,sincetheyformrelativelyrarely,we
ingandstellarevolution,whichhadbecomebottlenecksafterthe wouldnotexpectthesetohaveasignificantimpactontheevolution
hybridparallelizationschemewasimplemented.Withthesemodi- ofthecluster.
fications,forclustermodelsconsistingofsingleobjectsonly,dis- TheinitialconditionsusedinthemodelofWangetal.(2016)
tributedover320CPUcores(across16nodes),with86016GPU wereselectedtobesimilartolargeGCsintheMilkyWay.They
cores(on32GPUs),theyachieveaspeed-upfactorof400-2000, useaKingmodelwithconcentrationWo = 6,N = 106,half-mass
dependingonthenumberofstars.NBODY6++GPUisthefirstdi- radiusRh = 7.5pc,binaryfractionof5%,Galactocentricdistance
rectN-bodycodecapableofsimulatingtheevolutionofamillion- RG = 7.1 kpc, and metallicity Z = 0.00016. Stellar masses were
bodycollisionalstarclusterovermanyGyr.Foradetaileddescrip- chosen in the range 0.1 − 100M(cid:12) according to the initial mass
tionofthishybridcodeseeWangetal.(2015). function(IMF)ofKroupa(2001),whichisabrokenpower-lawof
Wangetal.havesharedtheresultsforthemillion-bodysim- theformdN/dm ∝ m−α,withα = 1.3for0.08 (cid:54) m/M(cid:12) < 0.5,
ulation presented in Wang et al. (2016). Provided with the exact and α = 2.3 for m/M(cid:12) (cid:62) 0.5. Binaries are created by choosing
initialclustermodelusedintheirN-bodysimulation,wehaveper- anexistingsinglestarrandomlytobecomeabinary,anddrawing
formedasimulationwithourownMCcode,CMC,whichhasbeen acompanionmassfromthedistribution0.6(m1/m2)−0.4(withm1(cid:54)
describedingreatdetailinseveralearlierpapers(Joshietal.2000, m2)(Kouwenhovenetal.2007).Thesemi-majoraxisisthendrawn
2001;Fregeauetal.2003;Fregeau&Rasio2007;Chatterjeeetal. from a distribution that is uniform in loga in the range 0.005-50
2010; Umbreit et al. 2012). Briefly, we use a variation of the so- AU,andeccentricityischosenfromthethermaldistribution.The
called“orbit-averagedMonteCarlomethod"developedbyHénon maindifferencebetweentheinitialmodelusedhereandthetypes
(1971)forsolvingtheFokker-Planckequation.Withourrecently ofmodelspresentedinMOR15isthatthemodelpresentedhereis
parallelizedMCcode(Pattabiramanetal.2013),wecancalculate muchmoreextended,withahalf-massradiusof7.56pc,compared
million-bodysimulationssuchastheonepresentedhereinlessthan tothetypicalhalf-massradiiofaround2pcforthemodelsfrom
2days1whendistributedoverjust48CPUcores. ourmainstudy.
While these two codes employ very different techniques for WehaveconvertedtheinitialconditionsusedbyWangetal.
modeling stellar dynamics, they include nearly all of the same into a format that can be used in CMC, which means our model
physicalprocesses,includingtwo-bodyrelaxation(treatedinanap- isanidenticalstar-by-starrepresentationoftheN-bodymodel.Of
proximatewaybytheMCcode),directphysicalcollisions,higher- course,toretainsphericalsymmetryweneededtoconvertthestar
orderstrongbinaryencounters(integrateddirectlyinbothcodes), coordinatesx,y,andztoaradialcoordinater,andthevelocityvec-
andsingleandinteractingbinarystellarevolution(usingSSEand torsvx,vy,andvztotheircorrespondingradialandtransversevec-
BSE, Hurley et al. 2000, 2002). These stellar evolution prescrip- tors,vrandvt.Allotherinitialproperties,includingstellarmasses,
tionshavebeenmodifiedfromtheoriginalpubliclyavailablever- positions,andbinaryproperties,areidentical.
sions. In the original SSE/BSE codes, only very low-mass BHs
were formed. This was pointed out by Belczynski et al. (2002),
who also provided a new prescription for forming more realistic 3 COMPARISONTODIRECTN-BODY
BHmassesincludingfallbackbasedonFryer&Kalogera(2001).
In both CMC and NBODY6++GPU, the remnant masses are se- Inwhatfollows,wemakedirectcomparisonsbetweenthetwosim-
ulationsuptoatimeof12Gyr.InSection3.1,weexploretheevo-
lectedfollowingthemetallicity-dependentprescriptionofBelczyn-
lutionanddynamicsoftheBHsinbothsimulations,whileinSec-
skietal.(2002)whichformsBHsintherangeofabout3−30M for
(cid:12)
tion 3.2 we highlight the agreement (and disagreement) between
typicalGCmetallicities.Also,themagnitudesofthenatalkicksre-
thestructuralpropertiesofbothclustermodels.
ceivedbystellarremnantsthatformviaasupernovashoulddepend
3.1 BlackHoleRetentionandEjection
1 Thisistrueformodelswithrelativelylargeinitialhalf-massradii,asis
thecasehere.However,MCsimulationsstartingfromcompactmodelswith Westartwithourmostexcitingresult:thedirect N-bodysimula-
halfmassradiiofabout1pccantakemorethanamonth. tion confirms the findings of MOR15 that large numbers of BHs
MNRAS000,1–9(2016)
4 C.L.Rodriguezetal.
Black Holes, 100 Myr Black Holes, 12 Gyr
400 400
2 Monte Carlo
350 350
N-body
300 300
1
250 250
r
e
b
m 200 0 200
nu 150 30 35 40 45 50 55 150
100 100
50 50
0 0
10 15 20 25 30 10 15 20 25 30
mass (M ) mass (M )
fl fl
1.0 1.0
on 0.8 0.8
uti
b BHs (MC)
stri 0.6 BHs (NB) 0.6
di
e non-BHs (MC)
v
ati 0.4 non-BHs (NB) 0.4
ul
m
u 0.2 0.2
c
0.0 0.0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
radius (pc) radius (pc)
Figure1. TheBHmass(top)andradialdistributions(bottom)forBHsinbothclustermodelsat100Myr(left)and12Gyr(right).Onthetop,theBHmasses
formodelMCareshownbythesolidblackline,andthoseformodelNBareindicatedbythedashedredline.OverallthenumberofretainedBHsandthemass
distributionsagreeverywell.ThereareasmallnumberofBHswithmassesbetween30M(cid:12)and55M(cid:12),highlightedintheinsert,whichareonlyformedinthe
NBmodel.ThebottompanelsshowtheradialdistributionsfortheBHs(solidcurves)andthenon-BHs(dottedcurves)separately,withmodelMCinblackand
modelNBinred.Eachstariscountedindividually,regardlessofwhetheritisapartofabinary(oranevenhigher-ordersystem,whichispossibleinthedirect
N-bodysimulation).
canberetainedformanyGyrinoldGCs.Attheendofthesimula-
tion,thereare1085BHsremaininginmodelMCand1036inmodel
NB.Figure1showstheinitialandfinaldistributionsofBHsmasses
forthetwomodels,whichareinnearlyperfectagreement.Dynami-
cally,thetwomodelsproduceextremelysimilarresults,withmodel
MCejecting336BHsandmodelNBejecting400BHsfromthepop-
ulation that is retained initially, and, in agreement with MOR15,
ejectingthemostmassiveBHsfirst.WedonotethatmodelNBcon-
tainsahandfulofBHswithmassesabove35M ,whereasmodel
(cid:12)
MChasnone.TheseabnormallymassiveBHsareformedasaresult
ofaminorbugintheN-bodytreatmentofmergedstellarbinaries,
whichwasdiscoveredafterthesimulationhadbegun.
Onthebottompanelsinthesamefigureweshowthecumula-
Figure2. ComparisonofthetotalnumberofBHsretainedineachcluster tiveradialdistributionsoftheBHs(solidcurves)andthenon-BHs
modelasafunctionoftime.ModelMCisshowninblackandmodelNBin (dottedcurves)afterBHformationandattheendofthesimulation
red.AsalreadyseeninFigure3,bothmodelsformsandretainsinitially (solidcurves).Thetwomodelsshowexcellentagreementinthera-
(at20Myr,afterBHformation)aroughlyidenticalnumberofBHs.Over
dialdistributionsofbothstartypes,withonlyslightdisagreement
time,BHsareslowlyejectedinbothmodelsatasimilarrate,withmodel
growingintheoutskirtsoftheclustermodelsforthenon-BHs.The
MCejecting336BHsandmodelNBejecting400BHsby12Gyr.
crudetidaltreatmentintheMCcodeisnotexpectedtoreproduce
themoreaccuratethree-dimensionaltreatmentinthedirectN-body
code. It is clear that the BHs, while more centrally concentrated
MNRAS000,1–9(2016)
Million-BodyComparisonsbetweenMonteCarloandDirect N-body 5
1.05
Monte Carlo
1.00 N-body
0.95
0.90
6
0
10.85
N/
0.80
0.75
0.70
0.65
0 2 4 6 8 10 12
time (Gyr)
0.050
0.049
n
o
ti
c
a0.048
r
y f
r
a
n
bi0.047
0.046
0 2 4 6 8 10 12
time (Gyr)
Figure4.Top:Numberofboundstarsremainingintheclusterasafunction
oftimeformodelMC(solidblack)andmodelNB(dashedred).Eachstaris
countedindividually,regardlessofwhetheritisapartofabinary(oraneven
higher-order system, which is possible in the direct N-body simulation).
Bottom:binaryfractionasafunctionoftime,withthesamecolorscheme
asabove.
thanthenon-BHs,arestillquitespreadoutoverthecluster.About
half of the BHs lie outside of the inner 4 pc, and about 10% lie
beyond10pc.
InFigure2weshowthetotalnumberofBHsretainedineach
modelasafunctionoftime.Theagreementisexcellentthroughout
theformationofnearlyalloftheBHs,withtheMCmodelforming
andretaininganinitialpopulationof1421BHsandtheNBmodel
forming 1436 BHs. By about 20 Myr, after all BHs have formed
andtheoneswithlargenatalkickshavebeenejectedfromtheclus-
ter, the number of BHs in each model levels out to the initially
retainednumbersgivenabove.OvertimethenumberofBHsde-
Figure3.ThepropertiesofthebinaryBHsejectedfromeachsimulation. creasesslowlyineachmodelatasimilarrate,withtheNBmodel
The top plot shows the number of ejected binaries over time for the MC ejecting 64 more BHs than the MC model by 12 Gyr. Part of this
model(solidblack)andtheNBmodel(dashedred).Themiddleplotshows
increase can be attributed to the more-massive BHs (35-55 M )
thecomponentmassesofthebinaries(MCinblackcircles,NBinredcrosses), (cid:12)
foundinthe N-bodysimulationwhichwillbeejectedfasterthan
whilethebottomplotshowstheeccentricityandsemi-majoraxisdistribu-
the∼ 25M BHs.However,the18erroneouslylargeBHscannot
tionsofthebinaries(withsamecolorscheme). (cid:12)
accountforthefulldiscrepancyinthenumberofejectedBHsbe-
tweenthetwomodels.
MNRAS000,1–9(2016)
6 C.L.Rodriguezetal.
Finally,wecomparethepropertiesoftheejectedBHbinaries. els.Intheouterpartsofthecluster(beyondafewtensofpc),model
InFigure3,weshowthenumberofBHbinariesejectedbyeach NBexpandsslightlymorethaninmodelMCearlyon,andthemem-
model over time. As the NB model is ejecting BHs more quickly oryofthisseemstolastthroughtheendofthesimulations.Again,
thantheMCmodel,therateofbinaryBHejectionisalsolargerin wedonotnecessarilyexpectthetwocodestoproduceidenticalre-
theN-bodymodel(54versus41binaries).Here,abouthalfofthese sultsbecauseoftheverydifferenttidalstrippingprescriptionsem-
ejectedbinariescanbeattributedtobinarieswhosecomponentsare ployed.Moreimportantly,the10%and50%Lagrangeradiiagree
abnormally large. Even then, this suggests that the Monte Carlo nearlyperfectly.The1%radiusinmodelMCisnoisieranddisplays
methodmayunderestimatetheproductionrateofbinaryBHs.We slightlydeepercollapsesthanthatofmodelNB,althoughtheouter
alsocomparethecomponentmasses,eccentricities,andsemi-major envelopeofthetwocurvesmatchesverywell.Thebehaviorofthe
axesoftheejectedbinaries(themiddleandbottompanelsofthe 1%radiusissimilartothatofthecoreradiusshowninFigure5.
samefigure).Inthisinstance,wefindverygoodagreementbetween Since the MC approach is not designed to handle the physics of
the NB and MC models. Such good agreement is to be expected, smallnumbersofobjects,itisactuallyquiteimpressivethatthebe-
asbothcodesmodelbinaryhardeningandpartnerexchanges(via havioronthesesmallscalesagreesaswellasitdoeswiththedirect
threeandfour-bodyencounters)usingadirectN-bodyintegration. N-bodysimulation.
InFigure7weagainshowtheLagrangeradii,onlythistime
weseparatetheBHsfromtheothertypesofobjectsinordertobet-
3.2 StructuralClusterProperties
ter understand the behavior of the BHs. With model MC in black
Nextwecomparetheoverallstructuralevolutionofthetwomod- andmodelNBinred,weshow(frombottomtotop)the1%,10%,
els.InFigure4weshowthenumberofboundstarsremainingin 50%and90%radiifortheBHs,andthesameLagrangeradiifor
thecluster(top)andthebinaryfraction(bottom)asafunctionof thenon-BHsinblueformodelMCandcyanformodelNB.Wecan
time,withmodelMCshowninblack(solidline)andmodelNBin seethattheinnermost1%oftheBHmassinbothmodelspartici-
red(dashedline).Herethenumberofboundstarsisdeterminedby patesinoscillationssimilartothoseseeninthecoreradiiinFigure
countingeverystarindividually(i.e.,abinarycountsastwostars). 5(i.e.,deeperoscillationsintheMCmodel,moreshallowonesin
By the final time of 12 Gyr, both models have lost roughly 30% theN-body),buttheseoscillationsarenotseeninthenon-BHpop-
oftheinitial1.05×106 stars.Thetwomodelsagreetowithin5% ulation. This tells us that these collapses are primarily driven by
(withmodelsMCandNBhaving6.95×105and7.29×105stars,re- theBHs.Thisinnermost1%massbinfortheBHsconsistsofonly
spectively).TheN-bodycodeusesaradialtidalcriterionsuchthat ∼ 10objects,soitisunreasonabletoexpecttheMCcodetotreat
anystarthatgoesbeyondtwicetheclustertidalradiusisremoved. thissubsystemperfectly.Inallbuttheinnermost1%BHmass,we
Incontrast,theMCcodereliesonanenergycriterion,whichhas findexcellentagreementbetweenthetwomodelsintheradialdis-
showntoproducebetteragreementwithdirect N-bodythanara- tributionsofbothspecies(BHsandnon-BHs)asafunctionoftime.
dialcriterion,buthasstillbeenobservedtostripstarsonaslightly Theoutermost90%LagrangeradiifortheBHsfallsatabout10pc
shorter timescale (Chatterjee et al. 2010). The binary fraction in forbothmodelsattheendofthesimulations,andthe90%radius
the models shows even better agreement to within about 1%, de- for the non-BHs lies just a bit further out at around 30 pc. Thus
creasing rather slowly with time from f =0.05 initially down to the BHs remain quite spread out in the cluster, well-mixed with
b
f =0.0465inMCand f =0.0462inmodelNB. theotherobjects,allthewaytotheendofthesimulation.Despite
b b
In Figure 5 we show the time evolution of the core radius, the deeper collapses of the central BH subsystem in model MC,
r ,andthehalf-massradius,r ,againwithmodelMCinblackand thefactthattheBHejectionrateagreessowellprovidesencour-
c h
model NB in red. Here we use the standard N-body definition of agingevidencethatBHdynamicsandretentionisnotnecessarily
thecoreradiusfromCasertano&Hut(1985),whichisadensity- governedbythedepthoftheBH-drivencoreoscillations.Theonly
weightedestimate ofthe averageof thestar positions.This theo- noticeablechangeintheBHdistributionsbetween100Myrand12
reticalcoreradiusisofcourseunrelatedtothatwhichanobserver Gyr(Figure1)isatthehighmassend,withtheheavyBHsbeing
wouldcalculate,sinceitdependsonmassratherthanluminosity. amongthefirsttobeejected.Thetwomodelsagreeverywellin
The half-mass radius is the radius that encloses half of the total thisregard.Thus,notonlydothenumbersofejectedBHsagree,
cluster mass. We find nearly perfect agreement between the two butsotoodothemassesoftheejectedBHs.
modelsintheevolutionofr fortheentirespanofthesimulations.
h
The agreement in the core radius is also very good, although at
early times (within a few hundred Myr) the core radius in model
4 DISCUSSION
MC expands slightly more than that of model NB. This is likely
causedbysubtledifferencesinourstellarevolutionroutines(e.g. We have found that our MC approach and the new code
windmassloss),whichhavebothbeenmodifiedfromtheoriginal NBODY6++GPUproducesimilarevolutionforamillion-starclus-
publicly-availableSSEandBSEcodes(Hurleyetal.2000,2002). terover12Gyr,bothintermsofglobalstructuralpropertiesandthe
Thesedifferencescanbedifficulttotrackdown,andthereforewe dynamicalevolutionoftheBHpopulations.Despitethesignificant
leavethistopicforfutureinvestigations.Inanycase,sincethecore agreementinthedynamics,wearestillinterestedinexploringwhy
radiiinbothmodelshavecontracteddowntoaverysimilarvalue theBHsinmodelMCexperiencedeepercollapsesthaninmodelNB,
withinaGyr,andthenremainingoodagreementuptotheend,this asthismayhelpusbetterunderstandthenatureoftheinteractions
slightdisagreementearlyonappearstohaveverylittleinfluenceon betweenBHsandtherestofthecluster.Onepotentialcauseforthe
long-termevolutionofthecoreradius. deepercollapsesintheMCmodelwouldbethepresenceofmore
InFigure6weshowacomparisonoftheLagrangeradiifor massive BHs than in the N-body model. MOR15 found that it is
thetwomodels,withmodelMCinblackandmodelNBinred.From always the most massive BHs that are found near the center and
bottom to top, the pairs of black and red lines correspond to the drivingthedeepcoreoscillations.Asaresult,theyarealsothefirst
radiienclosing1%,10%,50%,90%and99%ofthetotalmassin tobeejected.However,weactuallyfindtheoppositetobetrue:for
thecluster.Overallwefindstrongagreementbetweenthetwomod- thefirstGyr,theN-bodymodelcontainsabouttenBHsinthemass
MNRAS000,1–9(2016)
Million-BodyComparisonsbetweenMonteCarloandDirect N-body 7
[
Figure5.Comparisonofthetimeevolutionofthecoreradius(rc)andhalf-massradius(rh)betweenMC(solidblackline)andNB(dashedredline)uptothe
currenttimeoftheN-bodysimulation.Thetwomodelsshowgoodagreementintheevolutionofthehalf-massradius(towithin∼ 1pc),andfairlygood
agreementinthecoreradiusevolution.Wenotethatthecore-collapsesinmodelMCgodeeperthanthoseseeninNB(seeFigure7foramoredetailedlookat
whatishappeninginthecentralregion).
2
10
MC NB
)
c 1
p
10
(

dii
a
r

e
g
n
a
r 0
g 10
a
L
1
10
0 2 4 6 8 10 12
time (Gyr)
Figure6.ComparisonoftheLagrangeradiiformodelsMC(black)andNB(red).Thepairsofblackandredcurvesshowtheradiienclosing1%,10%,50%,
90%and99%(bottomtotop)ofthetotalclustermassforthetworespectivemodelsovertheentiresimulation.
range 35−55M , whereas the MC model contains none in this retainedinthetwomodelsagreesfairlywellthroughouttheentire
(cid:12)
range.By4Gyr,theseheavyBHshavebeenejectedfrommodel simulation.
NB.Inanycase,itisclearthatthedeepcollapsesintheMCmodel
cannot be explained by the presence of more massive BHs. Fur- Anotherpossibleexplanationforthedifferenceinthecollapse
thermore,wesawinFigure1that,besidesthehandfulofBHswith depthscouldbedifferencesinthepopulationofBHbinariesthat
bassesbetween35−55M ,theoveralldistributionofBHmasses arisefromthesimplifiedthree-bodybinaryformationprescription
(cid:12)
usedinCMC.Inourstandardprescription,weallowonlyBHsto
form binaries through three-body encounters, a choice motivated
MNRAS000,1–9(2016)
8 C.L.Rodriguezetal.
Figure7. LagrangeradiifortheBHsandthenon-BHsseparately.Theblack(MC)andred(NB)curvesshowtheradiienclosing1%,10%,50%and90%ofthe
BHmass.Similarly,theblue(MC)andcyan(NB)curvesshowtheradiienclosingthesamefractionsofthenon-BHmass,Herewecanseethatitisprimarily
theBHsthataredrivingthecoreoscillationsshowninFigure5.The10%,50%and90%BHLagrangeradiiforbothmodelsagreeextremelywell,althoughin
modelMCthe10%radiusisabitnoisier.The1%BHLagrangeradiusinmodelMCcollapsessignificantlydeeperthaninmodelNB,howeverthisbintypically
containsonlyabout10BHs.TheBHsremainverywellmixedwithotherstars(evenat12Gyr,90%oftheBHmassisspreadthroughoutthe∼ 105stars
comprisingtheinner10%ofnon-BHmass).
bypreviousstudiessuggestingthatfornon-compactstars,thehigh Forfulldetailsofthethree-bodybinaryformationprescriptionem-
densitiesrequiredtoformthree-bodybinarieswouldinsteadlead ployedinCMC,seeMOR15.
to physical collisions of these stars, rather than forming binaries Totesttheeffectofmodificationstoourbinaryformationpro-
(Chernoff&Huang1996),andisthusneverimportant.Thismeans
cedure we performed two additional simulations identical to MC,
that in model MC we are suppressing the creation of BH-non-BH
exceptforaslightmodificationtoourthree-bodybinarytreatment.
binariesbythree-bodyinteractions,andourpopulationofBH-non-
We have performed one simulation in which all stars, including
BHbinariesislimitedtothenumberofBHsthatcanexchangeinto
non-BHs, are allowed to participate in three-body binary forma-
existingbinariesthroughbinaryinteractions.Thesmallernumber tion.Wefind,however,thatthischangehasessentiallynoeffecton
of BH binaries (of any type) in the MC model may be responsi-
theBH-drivencollapsesnorontheBHevaporationrate.Nextwe
ble for the deep core collapses which would be overturned much
reducedη to0.5(afactoroftensmallerthanusedinmodelMC).
earlier in the collapse if we had a larger number of BH binaries. min
Inthiscasewesawaslightreductioninthecollapsedepths,butthe
Alternatively, the collapse depths could be sensitive to the hard- numberofretainedBHsremainedunaffected.Itseemsthatwhen
nessofthethree-bodybinariesthatareformed,ortheactualrate
softerbinariesareallowedtoform,whichalsoincreasestheover-
of binary formation, both of which are determined in part by the
allrateofproduction,theBH-drivencollapsescanbeoverturned
parameterη ,theminimumhardnessofthebinariesthatareal-
min slightlysooner,andtheverydeepestcollapsesareavoided.While
lowedtoforminthesimpleprescription(seeMOR15forthefull
thesemodificationstothesimpletreatmentofbinaryformationdo
detailsoftheprocedure).Thisparameteraffectsbinaryformation notseemtohavemuchofaneffectonoverallclusterevolutionor
intwodifferentways:first,ηminenterstherateequation,sotherate ontheretentionofBHs,itmayverywellaffectthedetailsofthe
ofbinaryformationinagiventimestepwiththisminimumhard-
BH-binarypopulations,suchastheproductionofBH-X-raybina-
ness Γ(η (cid:62) η ) ∼ η−3.5, meaning that it is easier to form softer
min min riesandtightBH-BHbinaries.Itwouldbepossibletostudythese
binaries(smallerη).Secondly,ifourproceduredeterminesthata
binarypopulationsingreaterdetail,andperhapsevenusethedirect
binaryshouldform,thenwechoosetheactualhardnessforthebi-
N-bodymodeltocalibratethethree-bodybinaryprocedureusedin
nary from a distribution according to the differential rate, dΓ/dη,
CMC,butweleavethisforafuturestudy.
withlowerlimitη .Thisfunctionfallsoffrapidlyforincreasing
min
valuesofη,soalowerlimitwilltendtoresultinsignificantlysofter IncontrasttotheMCmodelspresentedinMOR15,theinitial
binariesbeingformed.InmodelMC,aswellasinthemodelspre- modelusedinthisstudyismuchmoreextended.Otherthanslight
sentedinMOR15,wehavestrictlyusedahardnesscutoffη =5. differences in the initial binary populations, the rest of the initial
min
model details and the physics implementation remains the same.
The behavior of the core radius and the inner BH Lagrange radii
MNRAS000,1–9(2016)
Million-BodyComparisonsbetweenMonteCarloandDirect N-body 9
in model MC is nonetheless very similar to that seen in the much MM was supported by an NSF GK-12 Fellowship, award DGE-
morecompactMCmodelspresentedinMOR15.Also,theorderin 0948017. W. L. acknowledge the support by the Silk Road
whichBHsareejected,fromheavytolight,agreeswiththefind- Project at the National Astronomical Observatories of China
ingsofMOR15.ThemaindifferenceintheevolutionofmodelMC (NAOC, http://silkroad.bao.ac.cn), the Max-Planck-Institute for
presentedhereisthat,beingmuchmoreextended,theclusterhasa Astrophysics and the Max Planck Computing and Data Facility
significantlylongerrelaxationtimescale,andthereforeitprocesses (MPCDF, http://www.mpcdf.mpg.de/) in Garching, Germany. All
(andthusejects)BHsmuchmoreslowlythanmorecompactmod- directN-bodysimulationswererunontheirHydraGPUcluster.
els with similar N from MOR15, which have typically ejected ∼
1000sofBHsonthesamephysicaltimescale.
REFERENCES
BarnardR.,GarciaM.,LiZ.,PriminiF.,MurrayS.S.,2011,ApJ,734,79
BelczynskiK.,KalogeraV.,BulikT.,2002,ApJ,572,407
5 SUMMARYANDCONCLUSIONS
BreenP.G.,HeggieD.C.,2013,MNRAS,432,2779
Inthisstudywehavecomparedtheresultsofrealisticlarge-Nstar CasertanoS.,HutP.,1985,ApJ,298,80
clustersimulationscreatedusingtwoverydifferenttechniques:an ChatterjeeS.,FregeauJ.M.,UmbreitS.,RasioF.A.,2010,ApJ,719,915
ChernoffD.F.,HuangX.,1996,inHutP.,MakinoJ.,eds,IAUSympo-
orbit-averagedMCapproach(CMC)andnewdirectN-bodycode
(NBODY6++GPU). In terms of overall dynamical structure and siumVol.174,DynamicalEvolutionofStarClusters:Confrontationof
TheoryandObservations.p.263
evolution,aswellasthedynamicsoftheBHsandlong-termBHre-
ChomiukL.,StraderJ.,MaccaroneT.J.,Miller-JonesJ.C.A.,HeinkeC.,
tention,wefindquiteremarkableagreementbetweenthetwotech-
NoyolaE.,SethA.C.,RansomS.,2013,ApJ,777,69
niques.Mostnotably,thedirectN-bodymodelconfirmsthefinding FregeauJ.M.,RasioF.A.,2007,ApJ,658,1047
ofMOR15thatverylargenumbersofBHs(∼1000)canberetained FregeauJ.M.,GürkanM.A.,JoshiK.J.,RasioF.A.,2003,ApJ,593,772
for∼10GyrtimescalesinoldGCs. FryerC.L.,KalogeraV.,2001,ApJ,554,548
Weseeslightdifferencesintheevolutionoftheinnermost1% HarrisW.E.,1996,VizieROnlineDataCatalog,7195,0
BHLagrangeradius:modelMCdisplayscoreoscillationsdrivenby HeggieD.C.,GierszM.,2014,MNRAS,439,2459
theheaviestBHs,whichdonotoccurinmodelNB.However,this HemsendorfM.,KhalisiE.,OmarovC.T.,SpurzemR.,2003,HighPer-
formanceComputinginScienceandEngineering.SpringerVerlag,71,
disagreementdoesnotseemtohaveanysignificantimpactonthe
388
overallevolutionoftheclustermodels.Forexample,thehalf-mass
HénonM.H.,1971,Ap&SS,14,151
radiusevolutionofmodelMCmodelsagreeswiththatoftheN-body
HurleyJ.R.,PolsO.R.,ToutC.A.,2000,MNRAS,315,543
model. Additionally, both models eject BHs at roughly the same
HurleyJ.R.,ToutC.A.,PolsO.R.,2002,MNRAS,329,897
rate,losingjustover300oftheinitiallyretainedBHsovertheentire JoshiK.J.,RasioF.A.,PortegiesZwartS.,2000,ApJ,540,969
simulationandthusstillretainingmorethanathousandBHsatthe JoshiK.J.,NaveC.P.,RasioF.A.,2001,ApJ,550,691
endofthesimulation.Thisverifiesthat,whileourapproximateMC KouwenhovenM.B.N.,BrownA.G.A.,PortegiesZwartS.F.,KaperL.,
approachmaynot becapableoftreating thephysicsofthe small 2007,A&A,474,77
numberofBHsinthecoreperfectly,ourtechniqueisnonetheless KroupaP.,2001,MNRAS,322,231
capable of reproducing the bulk evolution of the cluster as seen MaccaroneT.J.,KunduA.,ZepfS.E.,RhodeK.L.,2007,Nature,445,183
MaccaroneT.J.,KunduA.,ZepfS.E.,RhodeK.L.,2011,MNRAS,410,
intheN-bodymodelratherwell,makingitaviabletechniquefor
1655
modeling the dynamical evolution of realistic clusters containing
MackeyA.D.,WilkinsonM.I.,DaviesM.B.,GilmoreG.F.,2008,MN-
largenumbersofstellar-massBHs.WithourparallelCMCcode,a
RAS,386,65
clustermodelliketheonestudiedherecanbecompleted(evolved
MorscherM.,PattabiramanB.,RodriguezC.,RasioF.A.,UmbreitS.,2015,
upto12Gyr)withinaboutaday.Thesamemodelcomputedwith ApJ,800,9
NBODY6++GPUrequiresmorethan6monthstocomplete. NitadoriK.,AarsethS.J.,2012,MNRAS,424,545
MOR15 argued that the long-term survivability of BHs in Pattabiraman B., Umbreit S., Liao W.-k., Choudhary A., Kalogera V.,
clusters relies on the fact that the BHs mostly avoid the Spitzer MemikG.,RasioF.A.,2013,ApJS,204,15
instability,incontrasttowhathasoftenbeenassumed.Thelackof ShihI.C.,KunduA.,MaccaroneT.J.,ZepfS.E.,JosephT.D.,2010,ApJ,
BH-drivencollapsesintheN-bodymodelpresentedhereprovides 721,323
SippelA.C.,HurleyJ.R.,2013,MNRAS,430,L30
evidence thatthe BHsare even less susceptible tothe Spitzerin-
Spurzem R., 1999, Journal of Computational and Applied Mathematics,
stabilitythanpredictedpreviouslybyMOR15.Thiswouldprovide
109,407
strongsupportforthemainconclusionfromMOR15,namelythat
SpurzemR.,BerentzenI.,BerczikP.,MerrittD.,Amaro-SeoaneP.,Harfst
if large numbers of BHs are retained initially (as they are under
S.,GualandrisA.,2008,inAarsethS.J.,ToutC.A.,MardlingR.A.,
standard assumptions regarding star cluster initial conditions and eds, Lecture Notes in Physics, Berlin Springer Verlag Vol. 760, The
BHformationprocesses),thenmanywillberetainedat12Gyr,es- Cambridge N-Body Lectures. p. 377, doi:10.1007/978-1-4020-8431-
peciallyforamodelwithalargeinitialvirialradius(asstudiedin 7_15
thiswork). StraderJ.,ChomiukL.,MaccaroneT.J.,Miller-JonesJ.C.A.,SethA.C.,
2012,Nature,490,71
UmbreitS.,FregeauJ.M.,ChatterjeeS.,RasioF.A.,2012,ApJ,750,31
WangL.,SpurzemR.,AarsethS.,NitadoriK.,BerczikP.,Kouwenhoven
M.B.N.,NaabT.,2015,MNRAS,450,4070
6 ACKNOWLEDGMENTS
WangL.,SpurzemR.,AarsethS.,GierszM.,AskarA.,BerczikP.,NaabT.,
ThisworkwassupportedbyNSFGrantAST-1312945andNASA SchadowR.,KouwenhovenM.B.N.,2016,MNRAS,submitted
Grant NNX14AP92G. All computations using CMC were per-
formedonNorthwesternUniversity’sHPCclusterQuest.CRwas ThispaperhasbeentypesetfromaTEX/LATEXfilepreparedbytheauthor.
supported by an NSF GRFP Fellowship, award DGE-0824162.
MNRAS000,1–9(2016)