The final-parsec problem in non-spherical galaxies revisited
aa r X i v : . [ a s t r o - ph . GA ] A p r Draft version July 20, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
THE FINAL-PARSEC PROBLEM IN NONSPHERICAL GALAXIES REVISITED
Eugene Vasiliev , Fabio Antonini , and David Merritt Lebedev Physical Institute, Moscow, Russia School of Physics and Astronomy and Center for Computational Relativity and Gravitation,Rochester Institute of Technology, Rochester, NY 14623, USA and Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario, Canada
Draft version July 20, 2018
ABSTRACTWe consider the evolution of supermassive black hole binaries at the center of spherical, axisym-metric, and triaxial galaxies, using direct N -body integrations as well as analytic estimates. We findthat the rates of binary hardening exhibit a significant N -dependence in all the models, at least for N in the investigated range of 10 ≤ N ≤ . Binary hardening rates are also substantially lower thanwould be expected if the binary “loss cone” remained “full,” as it would be if the orbits supplyingstars to the binary were being efficiently replenished. The difference in binary hardening rates betweenthe spherical and nonspherical models is less than a factor of two even in the simulations with thelargest N . By studying the orbital populations of our models, we conclude that the rate of supply ofstars to the binary via draining of centrophilic orbits is indeed expected to be much lower than thefull-loss-cone rate, consistent with our simulations. We argue that the binary’s evolution in the sim-ulations is driven in roughly equal amounts by collisional and collisionless effects, even at the highest N -values currently accessible. While binary hardening rates would probably reach a limiting valuefor large N , our results suggest that we cannot approach that rate with currently available algorithmsand computing hardware. The extrapolation of results from N -body simulations to real galaxies istherefore not straightforward, casting doubt on recent claims that triaxiality or axisymmetry aloneare capable of solving the final-parsec problem in gas-free galaxies. Subject headings: galaxies: elliptical and lenticular, cD – galaxies: evolution – galaxies: kinematicsand dynamics – galaxies: nuclei INTRODUCTION
According to large-scale simulations of the clusteringof dark matter in the universe, the mean time betweenmajor mergers of dark halos varies between ∼ . z = 10 and ∼ yr at z = 1, witha weak dependence on halo mass (Fakouri et al. 2010).Mergers between halo-sized objects are expected to bringthe baryonic components together in a time compara-ble to the halo coalescence time (White & Rees 1978;Barnes 2001). By the same reasoning, if each merg-ing galaxy contained a central supermassive black hole(SMBH), the two SMBHs would form a bound systemin the merged galaxy – a binary SMBH – shortly af-ter the merger was complete (Begelman et al. 1980; Roos1981). This idea has received considerable attention be-cause the ultimate coalescence of such a binary wouldgenerate an observable outburst of gravitational waves(GWs; Thorne & Braginskii 1976).Begelman et al. (1980) pointed out a potential bottle-neck in the evolution of binary SMBHs toward coales-cence. The binary interacts with nearby stars, eject-ing them with velocities comparable to the binary’sorbital velocity. This is the “gravitational slingshot”(Saslaw et al. 1974). But the process is self-limiting,and it is not a priori clear that the orbits will be re-populated in a time shorter than the age of the uni-verse. This has been called the “final-parsec problem”(Milosavljevic & Merritt 2003a); the name derives from [email protected]@[email protected] the fact that the natural separation of a massive binaryat the center of a galaxy is roughly a parsec.Just as in the case of a single SMBH at the cen-ter of a galaxy, a binary SMBH can continue in-teracting with stars only if the relevant orbits –the “loss-cone” orbits – are repopulated. Repopula-tion of loss-cone orbits by gravitational encounters –i.e., collisional relaxation – is well understood in thecontext of single SMBHs in spherical galaxies (e.g.,Lightman & Shapiro 1977; Cohn & Kulsrud 1978). Col-lisional effects are less well understood in the ax-isymmetric case (Magorrian & Tremaine 1999; Yu 2002;Vasiliev & Merritt 2013), and no such treatment existsyet for triaxial systems. Additional complications in thecase of a binary SMBH include the likely anisotropy ofthe initial orbital distribution (just after formation of thehard binary), and the fact that the size of the loss region(= binary semimajor axis a ) changes with time.Steady-state loss-cone theory (Merritt 2013, chapter6) distinguishes between two regimes: the empty-loss-cone regime, in which the rate of repopulation of loss-cone orbits is slow enough that such orbits are fully de-pleted; and the full-loss-cone regime, which is reachedwhen the encounter rate is so high that it is no longera limiting factor. In the former case, the rate of orbitalrepopulation is inversely proportional to the relaxationtime, which itself scales as ∼ N/ ln N with the numberof stars in the galactic nucleus (or particles in the N -body simulation), while in the latter case it saturates ata value that is essentially N -independent. The harden-ing rate of the binary, s ≡ d (1 /a ) /dt , is proportional tothe rate of repopulation of loss-cone orbits. Early N -body simulations of galaxies containing binary SMBHsadopted rather small N -values (Quinlan & Hernquist1997; Milosavljevi´c & Merritt 2001; Hemsendorf et al.2002) and were essentially in the full-loss-cone regime,showing little or no dependence of s on N . More re-cent studies have verified that s drops with N for suf-ficiently high N (Makino & Funato 2004; Berczik et al.2005; Merritt et al. 2007), as expected when approach-ing the empty-loss-cone regime. Approximate Fokker–Planck models of binary evolution in spherical galaxies(Milosavljevi´c & Merritt 2003b; Merritt et al. 2007) alsopredict that the hardening rate of the binary should scaleapproximately as N − for large N .Even in the absence of gravitational encounters,stars can continue to be supplied to the central bi-nary if their orbital angular momenta are modifiedby torques from the nonspherical galaxy. This “colli-sionless” mode of loss-cone repopulation is essentiallyindependent of the number of stars (in a galaxy ofgiven size and mass) and can in principle provide starsto the central binary at high enough rates to ensurecoalescence in a Hubble time (Merritt & Poon 2004;Holley-Bockelmann & Sigurdsson 2006). Some recent N -body simulations of galaxy mergers suggest in fact thatthe rates of binary evolution depend on N weakly, ifat all, a result that the authors have attributed to thenonspherical shapes of the merged galaxies (Khan et al.2011; Preto et al. 2011). However, other interpretationsare possible for these intriguing results, due to the com-plex interplay between collisional and collisionless mech-anisms. For instance, the size of the loss region, fromwhich orbits can be driven into the loss cone via collision-less effects, is much larger in nonspherical systems thanin spherical ones, hence, it can be more readily repopu-lated by collisional relaxation. Considerations like thesesuggest that it might be difficult to design an N -bodysimulation in which collisional effects would be truly neg-ligible, as they are expected to be in many real galaxies.In this paper we carry out direct N -body integrationsof binary evolution in spherical, axisymmetric, and triax-ial galaxies, constructed initially as equilibrium models.We do not simulate the galaxy merger process, and in thissense, our initial conditions can be considered less realis-tic than in some of the studies cited above. On the otherhand, our method lends itself to a more rigorous com-parison between binary evolution rates in galaxies withdifferent geometries. We also carry out a much more de-tailed analysis of the orbital families in our models, andof the connection between orbital types and the binaryevolution rates seen in the simulations. NUMERICAL SIMULATIONS
Model construction
We considered three series of galaxy models havingthe same radial density profile but different degrees ofasymmetry: spherical (S), oblate axisymmetric (A), andtriaxial (T). The mass distributions were given in eachcase by a generalization of the Hernquist (1990) broken- power-law model: ρ ( r ) = M total πabc r (1 + ˜ r ) , (1)˜ r ≡ (cid:16) xa (cid:17) + (cid:16) yb (cid:17) + (cid:16) zc (cid:17) . Axis ratios for the nonspherical models were a : b : c =1 : 1 : 0.8 (A) or 1 : 0.9 : 0.8 (T). We adopt units suchthat abc = M total = 1. The potential of a central pointmass, − GM bin /r , representing the massive binary wasadded to the self-consistent potential of the stars; weconsider a single value for M bin = 10 − , which is typicalof previous studies. This is somewhat larger than thetypical ratio ∼ − of SMBH mass to galaxy mass (e.g.Marconi & Hunt 2003); however, since most stars thatinteract with the binary come from radii smaller than thescale radius of the model, restricting our consideration tothe central parts of a galaxy is not likely to strongly af-fect the results. The distribution function for the spher-ical models was created using Eddington’s inversion for-mula (Merritt 2013, Equation (3.47)), which gives theunique, isotropic f ( E ) corresponding to a specified po-tential and density. The nonspherical { A, T } modelswere constructed by Schwarzschild (1979) orbit super-position method, using the publically available SMILE software (Vasiliev 2013) and 10 orbits. The nonspheri-cal models are intrinsically anisotropic in velocity space,but we imposed “maximal isotropy” by requiring the ve-locity dispersion in the radial direction to equal one-halfthe sum of the velocity dispersions in the two transversedirections (no additional constraint was placed on thelatter).Monte Carlo realizations of each model were con-structed for values of N in the range 8 × ≤ N ≤ .We created several, independent realizations for each N (four for N ≤ N = 250K and 500K, onefor N = 1000K) because simulations with small N ex-hibit considerable scatter in the rate of evolution (e.g.,Merritt et al. 2007). Unless otherwise specified, dataplotted in the figures consist of values averaged over themultiple realizations.We then replaced the single massive particle by twoSMBH particles each of mass M bin / × − lo-cated symmetrically about the origin at x = ± . r m , defined as the ra-dius containing a stellar mass equal to 2 M bin ; initially r m ≈ .
15 and r m increases to about 0.2 as the centraldensity drops. The initial velocities of the SMBH parti-cles were set to 0 .
31 in model units, corresponding to acircular orbit in the x − y plane. In the course of the evo-lution, the eccentricity of the binary was found to remainlow for N > , although it became somewhat larger forsmaller N , with a large scatter between realizations; theaverage eccentricity was ≈ . N/ ) − / .We might also have placed just one SMBH particle atthe origin and allowed the other to spiral in, or placedthe two SMBHs symmetrically into an equilibrium modelcreated without a central mass. We experimented withthese and other configurations but found that the choiceof initial placement affected only the initial stage of evo-lution, not the behavior at the hard binary stage.We verified that the models so constructed were inequilibrium by following their evolution for 20 time units(with a single central point mass) using the N -body codedescribed below. No discernible evolution of the densityprofile or the model shapes was observed. Parameters of the N -body integrations We used the direct N -body integrator φ GRAPEch(Harfst et al. 2008) to follow the evolution of the mas-sive binary. This code combines hardware-acceleratedcomputation of pairwise interparticle forces (using the
Sapporo library (Gaburov et al. 2009), which emulatesthe GRAPE interface utilizing GPU boards) with a high-accuracy chain regularization algorithm to follow the dy-namical interactions of field stars with the two SMBHparticles. The chain radius was set to 4 × − lengthunits. In the present implementation, there could beonly one chain which necessarily includes the first SMBHparticle. Hence, in the early stages of evolution close ap-proaches of field stars with the second SMBH particlewere not regularized, which in principle might have ledto the accumulation of errors; nevertheless, the relativeerror in total energy was typically ∼ − for the accu-racy parameter η = 0 .
01 and even much smaller in thelater stages of evolution ( t ≥ −
40 time units), whenboth SMBHs were included in the chain.We used zero softening for interactions in the chainand set a very small softening length ǫ = 10 − outsidethe chain to prevent energy errors at the early integra-tion stages. This is a much smaller softening length thanthe values typically used in other studies, and also muchsmaller than the distance of strong deflection for encoun-ters between field stars; hence, we are guaranteed not tochange the effective value of the Coulomb logarithm orthe rate of relaxation. We checked this by repeating somesimulations with ǫ = 0 at the late stages and verifyingthat there was no substantial difference in the binaryhardening rate. Our experiments indicated that a largervalue of softening length ( ǫ ≥ − ) decreases the hard-ening rate and reduces the difference between simulationswith different N , as well as between S, A, and T models. Results
The models were evolved for 100 time units, with thefinal value of the binary semimajor axis, a , reaching a . (1 − . × − length units depending on N and onmodel type (Figure 1). The elapsed time until formationof a “hard binary,” a . a hard , was roughly t = 20; wedefine a hard in the standard way as a hard ≡ µM bin r m , µ ≡ M M M + M (2)(Merritt 2013, Equation (8.71)), with µ the reduced massof the binary. For our models, a hard = r m / ≈ − .Formation of the hard binary was accompanied by asubstantial change in the density profile and in the dis-tribution of particle energies in the models. Roughlyspeaking, the original ρ ∼ r − cusp was replaced bya shallower, ρ ∼ r − / density profile inside r m , corre-sponding to a mass deficit (Milosavljevi´c et al. 2002) oforder M bin . The change in the energy distribution wasmore dramatic; almost no particles remained for energies | E | & | Φ | ≈ .
8, where Φ is the depth of the potential well due to the stars alone. In other words, there are al-most no particles left that would be bound to the massivebinary. In the course of subsequent evolution, the densityprofile and distribution function changed only modestly,the changes being spread over a much wider range ofradii. Figure 2 shows that the axis ratios of the mod-els at the hard binary stage, determined iteratively fromthe tensor of inertia of N -body snapshots (Katz 1991),were quite close to their values at t = 0. In other words,the binary has not destroyed the large-scale flattening ortriaxiality, despite introducing some changes within theinfluence radius.For t &
30, the binary hardening rate s ≡ d (1 /a ) /dt was observed to be almost unchanging in each simula-tion. The hardening rate was computed as the slopeof the inverse semimajor axis plotted as a function oftime, calculated for several overlapping intervals of 20time units on the interval 30 ≤ t ≤
100 and aver-aged to get an estimate of scatter. We also averagedthe results over the different realizations with the same N , adding the scatter to the error bars. Figure 3, toppanel, shows that for low N ( . × ) there is lit-tle difference between the three models and a weak de-pendence on N , consistent with earlier studies that alsoused small particle numbers (Quinlan & Hernquist 1997;Milosavljevi´c & Merritt 2001; Chatterjee et al. 2003).For N & the spherical model demonstrates a cleardependence of hardening rate on N , approximately as s ∝ N − . , again in agreement with earlier studies(Makino & Funato 2004; Berczik et al. 2005), and con-sistent with theoretical models of collisional loss-cone re-filling (Merritt et al. 2007).Hardening rates in the nonspherical models A and Tare somewhat larger than in the spherical model, andquite close to each other, the difference appearing onlyfor the N = 10 model. But – contrary to our initialexpectations, based on the galaxy merger studies citedin Section 1 – the hardening rates still exhibit a clear N -dependence at large N , suggesting a continued role forcollisional orbital repopulation in these models. ANALYSIS OF HARDENING RATES
Predictions from scattering experiments
To understand the binary evolution found in the simu-lations, we begin by considering the hardening rate dueto the interaction between the binary and incoming starsin the N -body models. Following Hills (1983), we definethe dimensionless coefficient C describing the energy ex-change in one interaction between a “field” star of mass m ⋆ and a massive binary: C ≡ M bin m ⋆ ∆ E bin E bin , E bin ≡ − GM bin a . (3)The strength of an interaction can also be expressedin terms of the velocity at infinity, v , and impact param-eter, b , of the incoming star. In the hard-binary limit, a ≪ a hard , the scattering outcome depends on a singledimensionless quantity χ ≡ LL bin , L ≡ bv, L bin ≡ p GM bin a . (4) χ = 1 corresponds to a distance of closest approach ofthe field star to the binary’s center of mass equal to a , / a t t t Fig. 1.—
Evolution of binary hardness, 1 /a , as a function of time, for three series of models: spherical (left), axisymmetric (middle), andtriaxial (right). Different curves are for models with N varying from 8 × to 10 (from top to bottom), averaged over several realizationswith the same N , as discussed in the text. a x i s r a t i o r t=0t=30t=100 Fig. 2.—
Axis ratios of the triaxial model with N = 10 as afunction of radius: initially (red), after formation of a hard binary( t = 30, green), and still further into the evolution ( t = 100, blue).The model remains triaxial and close to its original shape. if the binary were replaced by a single point mass. Weadopt the dependence of C on χ displayed in Figure 1 ofSesana et al. (2006), which we find to be reasonably wellapproximated (for a circular, equal-mass binary) by C ( χ ) ≈ (cid:26) . − . χ + 21 . χ − χ , χ < . . χ − . − χ − . , χ ≥ . (cid:18) a (cid:19) = 2 m ⋆ M bin a C ( χ ) ≡ H ( χ ) . (6)Consider for the moment a model in which field starsare drawn from a homogeneous background with isotrop-ically distributed velocities having a single magnitude v . Then the hardening rate is s ≡ ddt (cid:18) a (cid:19) = ρvm ⋆ Z ∞ πb db H ( χ ) = Gρv H , (7a) H ≡ Z ∞ π C ( χ ) χ dχ . (7b)Integrating Equation (7b) using the expression (5) for C ( χ ) yields H ≈ .
5, almost the same as the valuegiven by Quinlan (1996) in the hard-binary limit.Next, we derive similar expressions for the hardeningrate in the more realistic case where the distributionfunction of unbound stars has the form f ( E, L ). We takeinto account that stars interact with the binary once perradial period T rad , and that the number density of starsin { E, L } space is related to the phase-space mass densityby dN = 8 π m − ⋆ T rad ( E, L ) L f ( E, L ) dE dL . Then s = Z dE Z L circ ( E )0 dL π T rad L f ( E, L ) m ⋆ m ⋆ C ( L/L bin ) M bin a T rad = 4 πG Z dE Z L circ /L bin f ( E, χL bin ) 8 πχC ( χ ) dχ. (8)Here Φ is the lowest energy of an orbit unbound tothe central object (equal to the depth of potential wellproduced by the stars), and L circ ( E ) is the angular mo-mentum of a circular orbit of energy E , which is muchlarger than L bin in the limit of a hard binary. If thedistribution function is isotropic, the hardening rate isgiven simply by (e.g., Merritt 2006, Equation 11) s iso = 4 πG H Z f ( E ) dE . (9)In the spherical case, this is the rate expected for the”full-loss-cone” regime, in which the initially isotropicdistribution of orbits in angular momentum remains fixed(i.e. the loss cone is repopulated efficiently enough thatwe can neglect its depletion). In the nonspherical case,the rate would be roughly the same if we keep the or-bit population fixed, even though any individual particlemay precess into and out of the loss cone due to regu-lar changes of angular momentum in addition to random s ( ha r den i ng r a t e ) N N -0.5 Full loss conesphericalaxisymmetrictriaxial 0 1 2 3 4 5 6 7 30 40 50 60 70 80 90 100 s ( ha r den i ng r a t e ) time spherical, simulationaxisymmetric, simulationtriaxial, simulationspherical, loss-cone estimateaxisymmetric, loss-cone estimatetriaxial, loss-cone estimate Fig. 3.—
Top panel: hardening rates s ≡ ( d/dt ) (1 /a ) as afunction of N for three series of models (red, spherical; green, ax-isymmetric; blue, triaxial), measured on the interval 30 ≤ t ≤ s = 18 ± Bottom panel: hardening rates as a function of time for N = 10 models (colors are the same as in the top panel). Filled sym-bols with horizontal error bars show the estimate of the slope of1 /a ( t ) on the corresponding interval of time; dashed lines with opensymbols show the estimates from the loss-cone population (Equa-tion 10) at corresponding moments of time. The two estimatesagree quite well and have a rather moderate variation over time,with a weak tendency to decline. perturbations. Alternatively, we could estimate the full-loss-cone hardening rate by taking the expression for thehardening rate in a homogeneous isothermal background, s = GρH/σ , with H ≈
15 (Quinlan 1996; Sesana et al.2006), and substituting the values of σ and ρ computedat the binary’s radius of influence (say), which yields asimilar number. For our models during the late stages ofevolution ( t & s ≈ ±
2, almostindependent of N and geometry. This value is consistentwith the hardening rates observed in our lowest- N sim-ulations (Figure 3), but is several times higher than that N ( E ) / N t o t a l Er circ Original model at t=0Isotropic model at t=30SphericalAxisymmetricTriaxial
Fig. 4.—
Population of stars in the loss cone (with angularmomenta
L < L bin ≡ √ GM bin a ) at t = 30 ( a = 1 / particles; cyan dot-dashed line: population in an isotropic modelwith the same density profile (i.e., full loss cone); dashed magentaline: same in the original model (before the formation of a hardbinary). of the simulations with N = 10 , for all three geometries:additional evidence that even our triaxial models are farfrom being in the full-loss-cone regime.Another way to justify the conclusion just reached isby calculating s directly from the N -body discrete distri-bution function, by summing the contributions to energyexchange (6) for each particle with angular momentum L i during its orbital period T rad ,i : s N − body = N X i =1 m i M bin a T rad ,i C ( L i /L bin ) . (10)Computed in this way, the predicted hardening ratesfor all three series of models agree quite well with thosemeasured in the simulations (Figure 3, bottom panel).We then artificially randomized the directions of thevelocities of all the stars, leaving their magnitudes un-changed, thus creating an isotropic stellar system, whichwe integrated forward in time. (This was only donefor a spherical system, since it would break the self-consistency of the nonspherical ones). The measuredhardening rate was found to briefly jump to the full-loss-cone value for 1–2 time units before returning to theprevious value.The degree of loss-cone depletion can be quantified byplotting the actual population of particles on loss-coneorbits as a function of energy and comparing it with theisotropic case. Figure 4 shows that indeed the populationof low angular momentum stars is much depleted com-pared with that of an isotropic model. For the S model,however, there is an excess of stars with low binding en-ergies in the loss cone. These are the stars that have in-teracted once with the binary but that remained boundby the galactic potential, making them candidates forthe secondary slingshot process (Milosavljevi´c & Merritt2003b). Interestingly, those stars are absent in A and Tmodels, presumably because angular momentum is notconserved and they precess away from the loss cone. Theoretical estimates
The analysis presented above showed that there is onlya moderate difference in the degree of loss-cone depletionbetween the S, A and T models. The goal of this sectionis to explain that result, and to show that for the presentsimulations, we would not expect to be able to reliablydistinguish between collisional and collisionless loss-conerefilling processes.For this purpose, we develop a simplified model of col-lisionless draining of the “loss region,” defined as theregion of phase space from which orbits of stars in anonspherical potential are able to reach low enough an-gular momenta that they can interact with the binary.(By contrast, the “loss cone” is the subset of trajectoriesin the loss region that are destined to pass near the bi-nary in a single radial period or less.) This treatment issimilar to the orbital draining models considered by Yu(2002), Merritt & Poon (2004), Merritt & Wang (2005),and Sesana et al. (2007) for a binary black hole, or byVasiliev & Merritt (2013) for a single black hole in anaxisymmetric galaxy, but is more realistic in that we:(1) use the actual orbit population of the N -body simu-lations; (2) take into account the movement of orbits intoand out of the loss cone due to torquing by the mean-field potential; and (3) adopt a time-dependent size ofthe loss region around the central object.First, we take a snapshot from the actual simulationat t = 30, when the initial stage of rapid evolution ofthe density profile has finished and a constant harden-ing rate of the binary has set in. We study the orbitalpopulation of the model by extracting a sample of 10 particles from the snapshot and evolving them in a fixedbackground potential corresponding to the same snap-shot (plus one central point of mass M bin ), but repre-sented as a combination of smooth functions of radiustimes spherical harmonic functions of the angles ( θ, φ )(up to l max = 6 in angular harmonics, keeping only tri-axial or axisymmetric terms as appropriate; see Vasiliev(2013) for details). We then follow each orbit for 200 T rad ( E ), which is long enough to build a meaningful dis-tribution of values of the angular momentum at times ofpericenter passage (see the Appendix for details). For agiven orbit i , the probability of having L at pericenterbelow a certain value X is found to be well described bya linear function with slope S − i : P ( L < X ) ≈ X − L min2 ,i S i . (11)A zero value of L min ,i indicates a truly centrophilic or-bit, which can only exist in a triaxial potential; however,all orbits with L min ,i . L bin are “useful” for loss-conerepopulation. The combined mass of these orbits is 6%(1.5%) of the entire model for the T (A) case, i.e., sub-stantially higher than the mass of the binary. We assumethat values of L at subsequent pericenter passages areuncorrelated, which is reasonable given that most orbitsof interest appear to be chaotic. Next, we consider a time-dependent model for binaryevolution that includes a depletion of orbits in the lossregion. At each time step, we compute the instanta-neous hardening rate according to Equation (10), with C ( L i /L bin ) for each particle being averaged over all pos-sible values of angular momentum at pericenter, weightedwith the probability distribution (11): s ( t ) ≡ ddt (cid:18) a (cid:19) = N X i =1 m i ( t ) M bin a ( t ) T rad ,i Z d P C ( χ ( P )) == N X i =1 G m i ( t ) π S i T rad ,i H i ( t ) , (12a) H i ≡ Z ∞ χ min π C ( χ ) χ dχ , χ min ≡ L min ,i L bin ( t ) (12b)The quantity H i equals the constant H (7b) for acentrophilic orbit ( χ min = 0), and tends to zero for orbitswith L min ,i & L bin .Next we need to account for the decrease in the mass m i associated with each orbit, since a star once scatteredis assumed not to interact again with the binary. Toaccount for this, we relate the evolution of a to the massejection rate, as described by a dimensionless coefficient J ≡ aM bin dM ej da (13)(Quinlan 1996), and write dM ej dt = − N X i =1 dm i dt = J M bin a s. (14)Identifying each term in this sum with the correspond-ing term in Equation (12) allows us to write the equationfor the time evolution of the m i : dm i dt = − m i GM bin a ( t ) JH i ( t ) π S i T rad ,i (15)This simplified treatment does not account for the vari-ance in the outcome of scattering events with a givenimpact parameter, but is sufficient for the purposes ofour estimate. The value of J is taken from Sesana et al.(2006); it has little dependence on a and lies in the range0.5–1 and we took the former value, which was also foundin the N -body simulations of Milosavljevi´c & Merritt(2001). A higher value of J would decrease the late-timehardening rate, as the orbits would be depleted faster.Equations (12) and (15), together with the initial con-ditions a (0) = a init , m i (0) = 1 /N and the coefficients L min ,i , S i , T rad ,i derived from the orbit analysis, describethe time-dependent evolution of binary hardness in thepresence of orbital draining due to non-conservation ofangular momentum in a nonspherical background poten-tial.Figure 5 shows the predicted evolution of 1 /a in mod-els T and A starting at t = 30 and 1 /a = 250, comparedwith the actual evolution of 1 /a in the simulations with N = 10 . We also show, for the spherical model, theevolution rate expected for the full-loss-cone case andfor the actual population of loss-cone orbits. Two re-sults are apparent: (1) binary hardening rates predicted / a t F u ll l o ss - c one A c t u a l o r b i t p o p u l a t i o n Simulation, sphericalSimulation, axisymmetricSimulation, triaxialDraining model, axisymmetricDraining model, triaxial
Fig. 5.—
Evolution of binary hardness as a function of time.Solid lines: actual simulations with N = 10 , from top to bottom:T, A and S models; dotted lines: predictions from collisionlessdraining models, obtained by integrating Equations (12) and (15)forward in time, starting at t = 30 with 1 /a ≃ by our (collisionless) machinery in the A and T modelsfall below the actual rates, but only by modest factors.In other words, collisional replenishment of the orbits isstill contributing to the evolution of the binary in thesemodels. (2) The hardening rates in the A and T mod-els – both predicted and actual – fall substantially belowthe full-loss-cone rate. In other words, replenishment oforbits by torquing due to the nonspherical potential isnot efficient enough to keep the loss-cone orbits fully oc-cupied.Figure 5 suggests that – even for this large number ofparticles – we cannot reliably discern the relative con-tribution of collisional and collisionless processes to thehardening rate. A consequence of this result is that it isdifficult to extrapolate our results to the much higher N values relevant to real galaxies.At the same time, our simplified model suggests thatorbital draining can sustain a hardening rate that is sev-eral times below that of a full loss cone for a fairly longtime in the triaxial case, as the total mass of stars oncentrophilic orbits is a few times larger than M bin . Thesituation is less clear for the axisymmetric models: on theone hand, there are no “genuinely centrophilic” orbits inthis case; on the other hand, the reservoir of orbits thathave low enough L min to be able to interact with thebinary is still much larger than the volume of the losscone in the spherical case. The simplified calculationabove predicts that the draining rate of this loss regionis merely a factor of two lower than in the triaxial case,but that it may depend strongly on time. Discussion
The rate of repopulation of loss-cone orbits is deter-mined by a combination of collisional effects (due togravitational encounters) and collisionless effects (due totorquing by a nonspherical potential). The former arewell studied in the context of standard loss-cone theory.Approximate Fokker–Planck models for the binary evolu-tion in the spherical geometry (Milosavljevi´c & Merritt2003b; Merritt et al. 2007) predict that the hardeningrate should scale as s ∝ N − for large N (when thesystem is fully in the empty-loss-cone regime), and thenumerical simulations cited above show substantial N -dependence (although less steep than the prediction) for N & .The second mechanism of loss-cone repopulation,which can torque orbits from a much larger “loss region”into the loss cone, is collisionless and therefore does notdepend on N . Still, in our N -body integrations, we haveobserved a substantial N -dependence in both the axisym-metric and triaxial cases. A simple model for drainingof the loss region found that this mechanism could ac-count for about one-half of the hardening rate seen in thesimulations; the remainder would be attributed to N -dependent collisional effects. To reliably drive the latterbelow the expected rate of collisionless loss-cone repop-ulation would require a still much larger value of N –difficult to achieve with existing algorithms and hard-ware.An additional complication arises from the complex in-terplay between collisional and collisionless factors: thesize of the loss region, from which the orbits can betorqued into the loss cone due to collisionless effects, ismuch larger in nonspherical than in spherical systems,and hence can be more readily repopulated by collisional relaxation. For instance, Vasiliev & Merritt (2013) haveshown that in a steady state, the rate of loss-cone repop-ulation for a single SMBH in an axisymmetric galaxy isa few times higher than in the corresponding sphericalsystem. This makes it even more difficult to design an N -body simulation in which the collisional effects would benegligible, as they are expected to be in many real galax-ies. While we certainly expect there to exist an effectivelower limit, as a function of N , on the binary hardeningrates in nonspherical galaxies, we are unable to make aprecise statement concerning how low that rate might be.Our results suggest only that we cannot approach thatrate with currently available algorithms and computinghardware.The simple draining model considered in this pa-per does not account for a number of other pro-cesses that may be effective in the N -body sim-ulations or in real galaxies: Brownian motion ofthe binary (Merritt 2001), the “secondary slingshot”(Milosavljevi´c & Merritt 2003b), time-dependent per-turbations to stellar orbits even far from the bi-nary (Kandrup et al. 2003), changes in the stel-lar density profile due to continuous ejection ofstars (Milosavljevi´c & Merritt 2001), possible long-termchanges in degree of triaxiality (Merritt & Quinlan1998), among others. Clearly, a still more elaboratemodel, combining both collisionless and collisional pro-cesses, is desired to better understand the dynamics ofbinary SMBHs in galactic nuclei.Even in the context of models like ours, binary harden-ing rates will vary as a function of the degree of velocityanisotropy (i.e., the detailed orbital population), the ra-dial density profile, the shape of the nuclear isodensitycontours, etc. While exploring the full range of such vari-ation is beyond the scope of this paper, we can make somegeneral remarks. Observed galaxies exhibit a variety ofnuclear density profiles, from nearly flat cores to nuclearstar clusters having ρ ∼ r − . At larger radii the densityis almost always well fit by a Sersic or Einasto function.Varying the nuclear density profile would certainly affectthe early hardening rate of a binary. As discussed above(Section 2.3), during the formation of a “hard” binary,the initial density cusp is converted into a shallower pro-file as stars bound to the binary are ejected. The samewill be qualiatively true for any initial nuclear densityprofile, although one expects some dependence of the fi-nal profile on the initial profile (Milosavljevi´c & Merritt2001; Khan et al. 2012). A velocity distribution that ismore or less anisotropic would also lead to higher orlower binary hardening rates, at least initially. With re-gard to changes in the shape of the model, the resultsobtained here suggest that in nonspherical geometries,even fairly radical shape changes (axisymmetric → triax-ial) have only modest consequences for binary hardeningrates, and we expect the same to be true in nonsphericalmodels with axis ratios different from those consideredhere.Our results present an interesting contrast to thoseof other recent studies based on similar techniques.Khan et al. (2013) integrated spherical and axisymmet-ric models created initially in equilibrium, with densityprofiles and black hole masses essentially the same asin our models. While their binary hardening rate fora spherical model with N = 10 is comparable to ours,their flattened models have a much higher hardening ratethan ours, even exceeding our estimate for the full-loss-cone case (9). We also observed much less differencebetween the spherical and nonspherical models up to N = 10 . The reasons for these differences are unknownto us; however, at face value, our results call into ques-tion the robustness of the conclusion reached by thoseauthors that the final-parsec problem is “solved” in ax-isymmetric galaxies.On the other hand, triaxial models formed by a barinstability in a rotating galaxy (Berczik et al. 2006), orby mergers of two galaxies (Khan et al. 2011; Preto et al.2011), have shown essentially no dependence of the hard-ening rate on N . We speculate that this difference withour results might be due to the rotation of those models,to non-stationary clumpy structures in the case of themerger remnants, or to some other factor. More detailedstudy of the orbit populations could shed light on thismystery.Merritt & Poon (2004) (MP04) analyzed the orbitalpopulations in self-consistent, triaxial models of nucleicontaining central SMBHs. In the context of scale-freemodels with a steep, ρ ∼ r − density profile, they arguedthat collisionless feeding rates might be high enough toensure coalescence of massive binaries in less than 10Gyr, even in galaxies where the fraction of centrophilicorbits was small. The models of MP04 were both extreme(a steep density profile, maximal triaxiality) and ideal-ized (scale-free, fixed potential). We attempted to ver-ify their conclusions by creating a non-scale-free model with similar central properties (a steep density cusp andstrong triaxiality) and an SMBH mass of 10 − M gal ; inthis model the region extending to a few influence radiiwas still well inside the break radius. This model had ∼
15% of its mass on chaotic/centrophilic orbits, some-what more than in the models presented above. Thetime-dependent draining rate was then computed as de-scribed above, assuming an initial binary separation of0 . a hard . Our results for a ( t ) were not precisely the sameas those in MP04; in particular, the hardening rate wasfound to drop more rapidly with time, ( a hard /a ) ∼ t . .But the value of ( a/a hard ) − after a time correspondingto several Gyr was ∼ , i.e., more than enough to en-sure GW coalescence on Gyr timescales.This comparison highlights an important point: even amodest rate of binary evolution (due to stellar-dynamicalinteractions) can result in a separation small enough forGW emission to induce coalescence in less than a Hubbletime. The time for a circular binary’s orbit to evolve,from a = a to a = 0, due to GW emission is t GW ≡ t ( a = 0) − t ( a = a ) = 5256 c a G M M M bin (16) ≈ . × (1 + q ) q (cid:18) a − pc (cid:19) (cid:18) M bin M ⊙ (cid:19) − yrwith q ≡ M /M ≤ (cid:18) a a hard (cid:19) . (cid:0) . × yr (cid:1) (1 + q ) q G M c r m (17)i.e. if a a hard . .
03 (1 + q ) / q / (cid:18) M bin M ⊙ (cid:19) / (cid:18) r m
10 pc (cid:19) − (18)or in the units of our N -body models ( a hard ≈ . q =1) a . − (cid:18) M bin M ⊙ (cid:19) / (cid:18) r m
10 pc (cid:19) − . (19)Assuming that we have nearly reached the asymptotic(large- N ) limit in our simulations, Figure 1 suggests thatthe binaries in many galaxies would indeed be able toreach coalescence in a Hubble time. CONCLUSIONS
We carried out direct N -body integrations of binarysupermassive black holes in spherical, axisymmetric, andtriaxial models of galaxies, constructed initially as equi-librium models. Our integrations with particle numbersup to N = 10 demonstrated that in all three geometriesconsidered, the binary hardening rate s ≡ d (1 /a ) /dt ( a =binary semimajor axis) at late times does depend on N and is several times below the rate computed assuming afull loss cone. The difference in hardening rates betweenthe three models was quite modest – within a factor oftwo even for the simulations with largest N – and onlyin the largest- N case was there a noticeable difference in s between the axisymmetric and triaxial geometries.To assist in understanding these results, we computedthe expected hardening rates based on known resultsfrom three-body scattering experiments, together withthe distribution of particles in energy and angular mo-mentum in the N -body models. These predictions werefound to agree well with the hardening rates obtainedin the actual simulations. We also estimated, using asimple model for collisionless draining of orbits in the“loss region” (the collection of orbits that are able toreach the binary’s interaction sphere), the contributionof nonspherical torques to the rate of loss-cone repopu-lation, and we found it to be below or comparable to thecontribution from collisional effects, even for the highest-resolution simulation of our set.Based on these results, we argued that in order toreach a regime that is characteristic of massive galax-ies (in which collisional effects are believed to be negli-gible), substantially higher values of N might be neededin the simulations. Until this is done, it is premature tostate that the final-parsec problem in gas-free galaxies is “solved” by assuming nonspherical geometries.This work was supported by the National ScienceFoundation under grant no. AST 1211602 and by theNational Aeronautics and Space Administration undergrant No. NNX13AG92G. Part of this work was done atthe Alajar meeting and in the Henri Poincar´e institutein Paris. Computations were performed on the GPUsof the CITA Sunnyvale cluster, as well as the ARC su-percomputer at the SciNET HPC Consortium. SciNetis funded by the Canada Foundation for Innovation un-der the auspices of Compute Canada, the Governmentof Ontario, Ontario Research Fund Research Excellence,and the University of Toronto. We thank John Dubinskifor assistance with the cluster. REFERENCESBarnes, J. 2001, in ASP conf. series, 245, Astrophysical Ages andTimes Scales, , ed. T. von Hippel, C. Simpson, & N. Manset,(San Francisco, CA: ASP), 382Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature,287, 307Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680Berczik, P., Merritt, D., Spurzem, R., & Bischof, H. 2006, ApJL,642, L21Chatterjee, P., Hernquist, L., & Loeb, A. 2003, ApJ, 592, 32Cohn, H., & Kulsrud, R. 1978, ApJ, 226, 1087Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS,406, 2267Gaburov, E., Harfst, S., & Portegies Zwart, S. 2009, NewA, 14,630Harfst, S., Gualandris, A., Merritt, D., & Mikkola, S. 2008,MNRAS, 389, 2Hemsendorf, M., Sigurdsson, S., & Spurzem, R. 2002, ApJ, 581,1256Hernquist, L. 1990, ApJ, 356, 359Hills, J. G. 1983, AJ, 88, 1269Holley-Bockelmann, K., & Sigurdsson, S. 2006,arXiv:astro-ph/0601520Kandrup, H., Sideris, I., Terzi´c, B., & Bohn, C. 2003, ApJ, 597,111Katz, N. 1991, ApJ, 368, 325Khan, F. M., Holley-Bockelmann, K., Berczik, P., & Just, A.2013, ApJ, 773, 100Khan, F. M., Just, A., & Merritt, D. 2011, ApJ, 732, 89Khan, F. M., Preto, M., Berczik, P., et al. 2012, ApJ, 749, 147Lightman, A., & Shapiro, S. 1977, ApJ, 211, 244Magorrian, J., & Tremaine, S. 1999, MNRAS, 309, 447Makino, J., & Funato, Y. 2004, ApJ, 602, 93 Marconi, A., & Hunt, L. K. 2003, ApJL, 589, 21Merritt, D. 2001, ApJ, 556, 245Merritt, D. 2006, ApJ, 648, 976Merritt, D., 2013, Dynamics and Evolution of Galactic Nuclei(Princeton, NJ: Princeton Univ. Press)Merritt, D., Mikkola, S., & Szell, A. 2007, ApJ, 671, 53Merritt, D., & Poon, M.-Y. 2004, ApJ, 606, 788Merritt, D., & Quinlan, G. D. 1998, ApJ, 498, 625Merritt, D., & Wang, J. 2005, ApJL, 621, L101Milosavljevi´c, M., & Merritt, D. 2001, ApJ, 563, 34Milosavljevi´c, M., & Merritt, D. 2003a, in AIP Conf. Proc., 686,The Astrophysics of Gravitational Wave Sources, ed.J. Centrella and S. Barnes (Melville, NY: AIP), 201Milosavljevi´c, M., & Merritt, D. 2003b, ApJ, 596, 860Milosavljevi´c, M., Merritt, D., Rest, A., & van den Bosch, F. C.2002, MNRAS, 331, L51Preto, M., Berentzen, I., Berczik, P., & Spurzem, R. 2011, ApJL,732, L26Quinlan, G. 1996, NewA, 1, 35Quinlan, G., & Hernquist, L. 1997, NewA, 2, 533Roos, N. 1981, A&A, 104, 218Saslaw, W. C., Valtonen, M. J., & Aarseth, S. J. 1974, ApJ, 190,253Schwarzschild, M. 1979, ApJ, 232, 236Sesana, A., Haardt, F., & Madau, P. 2006, ApJ, 651, 392Sesana, A., Haardt, F., & Madau, P. 2007, ApJ, 660, 546Thorne, K. S., & Braginskii, V. B. 1976, ApJL, 204, 1Vasiliev, E. 2013, MNRAS, 434, 3174Vasiliev, E., & Merritt, D. 2013, ApJ, 774, 87White, S. D. M., & Rees, M. 1978, MNRAS, 183, 341Yu, Q. 2002, MNRAS, 331, 935APPENDIX
The analysis of draining rates in Section 3.2 necessitated an estimate of the minimum angular momentum L min attained by a given orbit in the smooth potential. Of course, no orbit can reach zero angular momentum on any finitetime interval (unless it is specially arranged to do so), but one can nevertheless estimate whether there is a positivelower limit on L min , or whether it is compatible with being zero, by the following procedure.We record the values of the squared angular momentum at pericenter passages L ,k ( k = 1 ..N peri ) and sort themin ascending order. As discussed in the text, it happens that the distribution usually follows a linear trend at low L .We therefore fit a linear regression with and without a constant term: L ,k = L + s kN peri + δ k = s ′ kN peri + δ ′ k , k = 1 ..N fit , N fit = 0 . N peri . (1)Here, s, L and s ′ are the coefficients of two- and one-parameter fits, and δ k and δ ′ k are the corresponding residuals.We assign the intrinsic dispersion of the values L ,k from the condition that χ per degree of freedom is unity in http://members.aei.mpg.de/amaro-seoane/ALM13 σ ≡ (cid:16)P N fit k =1 δ k (cid:17) / ( N fit − χ between the one-parameter and two-parameter fits is given by∆ χ = ( N fit − P N fit k =1 δ ′ k P N fit k =1 δ k − ! . (2)Of course, the residuals in the one-parameter fit are always greater than in the two-parameter fit, but if they are “nottoo much” greater then we accept the hypothesis that L = 0. More quantitatively, we accept the one-parameterfit if it is less than 3 σ away from the two-parameter fit, i.e., if ∆ χ < ∆ χ σ ≡ .
8, the latter value being the 3 σ deviation for a χ distribution with two degrees of freedom. Usually if this hypothesis is rejected (i.e., the orbit islabeled centrophobic), it is at the level of significance of many hundreds or thousands of σ . Finally, we take thevalues of L and s (or s ′ ) from the adopted regression to estimate the probability of having a given value of L2peri