Massive black hole binary plane reorientation in rotating stellar systems
aa r X i v : . [ a s t r o - ph . GA ] N ov Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 9 November 2018 (MN L A TEX style file v2.2)
Massive black hole binary plane reorientation in rotatingstellar systems
Alessia Gualandris ⋆ , Massimo Dotti , and Alberto Sesana Max-Planck Institut f¨ur Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany Dipartimento di Fisica G. Occhialini, Universit`a degli Studi di Milano Bicocca, Piazza della Scienza 3, 20126 Milano, Italy Albert Einstein Institut, Am M¨uhlenberg 1, Golm, D-14476, Germany.
ABSTRACT
We study the evolution of the orientation of the orbital plane of massive black holebinaries (BHBs) in rotating stellar systems in which the total angular momentum ofthe stellar cusp is misaligned with respect to that of the binary. We compare resultsfrom direct summation N -body simulations with predictions from a simple theoreti-cal model. We find that the same encounters between cusp stars and the BHB thatare responsible for the hardening and eccentricity evolution of the binary, lead to areorientation of the binary orbital plane. In particular, binaries whose angular mo-mentum is initially misaligned with respect to that of the stellar cusp tend to realigntheir orbital planes with the angular momentum of the cusp on a timescale of a fewhardening times. This is due to angular momentum exchange between stars and theBHB during close encounters, and may have important implications for the relativeorientation of host galaxies and radio jets. Key words: black hole physics – methods: numerical – stellar dynamics
In the standard ΛCDM cosmological scenario, massive blackhole binaries (BHBs) are the natural outcome of mergersof massive galaxies (Begelman et al. 1980). The orbital de-cay of BHBs has been studied in detail under the assump-tion of purely stellar systems (e.g. Makino & Funato 2004;Berczik et al. 2005; Baumgardt et al. 2006; Merritt et al.2007; Sesana et al. 2007; Sesana 2010) as well as gaseousenvironments (e.g. Armitage & Natarajan 2002, 2005;Escala et al. 2004, 2005; Dotti et al. 2006, 2007, 2009;Cuadra et al. 2009; Lodato et al. 2009; Roedig et al. 2011).Most studies of the evolution of BHBs in stellar envi-ronments consider non rotating, spherically symmetric nu-clei, despite the natural occurrence of rotation in mergerremnants (e.g. Gualandris & Merritt 2011). Simulations ofMBHs in spherical stellar systems predict a stalling of thebinary evolution at parsec scale separations, the “final par-sec problem” (Milosavljevi´c & Merritt 2003), which may ormay not be avoided in the presence of gas dissipation. ABHB can, however, reach coalescence even in a completelygas poor environment, if embedded in a non-sphericallysymmetric stellar distribution. Purely dissipationless simu-lations of non-spherical galaxy models (Berczik et al. 2006;Berentzen et al. 2009) or galaxy mergers (Khan et al. 2011; ⋆ E-mail: [email protected]
Preto et al. 2011; Gualandris & Merritt 2011) find continu-ing hardening of the binary down to separations where en-ergy loss due to emission of gravity waves becomes efficient.Simulations by Gualandris & Merritt (2011) show that effi-cient hardening is achieved even in systems close to axisym-metry. This departure from spherical symmetry is due tothe rotation introduced by the merger, which results in anenhanced rate of stars interacting with the binary.A comprehensive and systematic study of the evo-lution of BHBs in rotating stellar cusps is still missing.Sesana et al. (2011) (hereafter paper I) studied the evolu-tion of the eccentricity of unequal mass BHBs in rotatingsystems by means of an hybrid analytical/3-body scatteringformalism as well as full N -body simulations. They found astrong dependence of the eccentricity evolution of the binaryon the degree of co-rotation in the stellar cusp, with binariesincreasing their eccentricity in cusps containing a significantfraction of counter-rotating stars. This is relevant to the caseof galactic minor mergers, in which counter-rotating systemscan be produced as a result of rotation in the larger galaxy,thereby forming very eccentric binaries.Another important and potentially observable param-eter is the orientation of the binary orbital plane. Its evo-lution has been studied in spherical and isotropic systems(Merritt 2002; Gualandris & Merritt 2007). Due to the lackof any preferential direction in these systems, the orbitalplane of the binary can only undergo a random walk, result- c (cid:13) Gualandris et al. ing in small changes in the orbital plane on long timescales.Here we study the evolution of the orbital plane of BHBsin rotating stellar systems, where, as will be discussed, asignificantly larger reorientation is expected.The letter is organized as follows. In section 2 we de-scribe the evolution of the binary orbital plane in sphericallysymmetric models. In section 3 we describe the case of ro-tating models, comparing results from N -body simulationswith a simple analytical model. Final remarks are presentedin Section 4. The orientation of the orbital plane of BHBs embedded innon rotating, isotropic stellar nuclei has been studied forthe first time by Merritt (2002). In such systems, repeatedencounters with passing stars cause the orientation of thebinary to undergo a random walk (a “rotational Brownianmotion”, in Merritt (2002)). From simple analytical argu-ments, the expected change in orientation on one hardeningtime T h = | a / ˙a | of the binary (which is of the order of120 initial binary orbital periods) is (Gualandris & Merritt2007) ∆ θ ∝ q − / (cid:18) m ∗ M BHB (cid:19) / (cid:0) − e (cid:1) − / (1)where q = M /M is the binary mass ratio, M BHB the totalmass, and e the orbital eccentricity. Numerical scattering ex-periments have constrained the normalization in Eq. 1 onlyfor circular binaries. Neither the dependence on the eccen-tricity nor on the mass ratio has been tested numerically. Inthe case of a binary with M BHB = 10 m ∗ and q = 10 − , theexpected reorientation is ∆ θ ≈ ◦ .Here we check the validity of Eq. 1 by performing N -body simulations of binaries with different eccentricities(in the range 0 − .
99) embedded in non-rotating stellarcusps. We use the direct summation N -body code φ GRAPE(Harfst et al. 2007) in combination with the
Sapporo li-brary (Gaburov et al. 2009) to accelerate the computationson GPU hardware. In our integrations we adopted M =10 M ⊙ and q = 1 /
81 and set the initial semi-major axis to a i = 0 .
06 pc. The MBHs are embedded into a stellar cuspfollowing a Bahcall-Wolf ρ ( r ) ∼ r − / density profile at dis-tances smaller than 1 pc, with total mass M c ∼ . × M ⊙ and a mass enclosed in the binary orbit equal to 2 M . Thecusp is modelled with N = 32k equal-mass particles, result-ing in a black hole to star mass ratio of m ∗ /M = 7 . × − .Fig 1 shows the degree of re-alignment of the binary planeas a function of the initial eccentricity. The dashed line in-dicates the dependence expected from the analytical model,while the points show the results of our N -body integrations.The agreement is remarkable.For a given binary, the re-orientation in a non-rotatingsystem does not depend only on the mass density and ve-locity field of the stars, but also on the mass ratio betweenthe binary and a single star in the cusp . A binary embed- Note that this ratio is properly defined only for a single massstellar cusp. In the presence of a mass spectrum, the mean massof the interacting stars is a good proxy for m ∗ . Figure 1.
Change over one hardening time in the orientationof the angular momentum of the binary, as a function of theinitial eccentricity of the BHB. The dashed line represents thetheoretical dependence given in Eq. 1, normalized to match thedata at e = 0 .
25. The amplitude of the reorientation in our runs ishowever comparable to the analytical estimates in Merritt (2002). ded in a stellar cusp would, on average, experience a greaterrealignment as a result of few encounters with massive starsthan as a result of many encounters with low-mass stars.This feature is typical of any random walk process, in whichthere is no preferential direction for the change of the binaryplane in a single encounter and changes due to different en-counters tend to cancel each other out. The same randomwalk behaviour is observable in rotating cusps, if the angu-lar momentum of the binary is aligned with that of the cusp(Amaro-Seoane et al. 2010). This behaviour is illustrated inthe lower panel of Fig 2, which shows the dependence ofthe change in θ on the number of particles used to modelthe cusp, i.e., on the m ∗ /M BHB mass ratio. The agreementbetween the numerical results (full circles) and Eq. 1 (solidline) is outstanding.
Here we study in detail the re-orientation of a massive BHBembedded in a cusp with net rotation. To this aim, we per-form high accuracy N -body simulations similar to those pre-sented in the previous section, the only difference being thedegree of rotation that we enforce in the cusp. This is ob-tained, as in Paper I, by reversing the sign of all velocitycomponents for a random subset of cusp stars, at the timewhere the BHB is added. In particular, we generated mod-els with a fraction F = 0 .
875 of co-rotating stars and variedthe initial angle θ between the total angular momentumof the cusp and the angular momentum of the binary. Weconsidered cases with θ = 0 ◦ , ◦ , ◦ , ◦ , ◦ , ◦ . The evolution of θ as a function of time is shown inFig. 3 (solid lines) for all runs. Binaries whose angular mo-mentum is initially mis-aligned with respect to that of thecusp tend to realign it on a time-scale of a few hardeningtimes. The re-alignment can be quite significant, with anaverage change in θ of the order of 50% of the initial value,and up to 70% in some cases. For example, the case with θ = 150 ◦ reaches a final value of θ ∼ ◦ , with a changeof more than 100 ◦ over one hardening time. This is about100 times larger than the reorientation measured for a sim-ilar isotropic model (see Fig. 1). As will be discussed be-low, such large reorientation of the orbital plane may have c (cid:13)000
875 of co-rotating stars and variedthe initial angle θ between the total angular momentumof the cusp and the angular momentum of the binary. Weconsidered cases with θ = 0 ◦ , ◦ , ◦ , ◦ , ◦ , ◦ . The evolution of θ as a function of time is shown inFig. 3 (solid lines) for all runs. Binaries whose angular mo-mentum is initially mis-aligned with respect to that of thecusp tend to realign it on a time-scale of a few hardeningtimes. The re-alignment can be quite significant, with anaverage change in θ of the order of 50% of the initial value,and up to 70% in some cases. For example, the case with θ = 150 ◦ reaches a final value of θ ∼ ◦ , with a changeof more than 100 ◦ over one hardening time. This is about100 times larger than the reorientation measured for a sim-ilar isotropic model (see Fig. 1). As will be discussed be-low, such large reorientation of the orbital plane may have c (cid:13)000 , 000–000 assive black hole binary plane reorientation Figure 2.
Change in the binary plane orientation per hardeningtime versus number of particles in the cusp (in the range N = 2kto N = 128k) or, equivalently, ratio between star mass and MBHmass. In all the runs e = 0 .
5. The lower panel refers to isotropicmodels, while the upper panel refers to rotating models with F =0 .
875 and an initial angle between the binary angular momentumand that of the stellar cusp θ = 30 ◦ , ◦ , ◦ . Figure 3.
Evolution of the angle θ between the angular momen-tum vector of the binary and that of the stellar cusp. Time isexpressed in units of the initial binary orbital period (bottom)and of the hardening time (top). Different solid lines are for mod-els with different initial values of θ , from θ = 0 ◦ to θ = 150 ◦ .Dotted lines represent the predictions from the theoretical modelwhile dashed lines indicate the predictions modified to take intoaccount the finite reservoir of angular momentum in the stellarcusp. Figure 4.
Solid lines: Distribution of the angle θ cusp betweenthe angular momentum vector of a single star and the initialangular momentum vector of the cusp, for the run with θ = 150 ◦ .Only unbound stars are considered. The different curves refer todifferent times during the evolution, from 1 T h to 8 T h . Stars areejected roughly isotropically from early times. Dotted line: Initialdistribution of θ cusp for all cusp stars, normalized by a factor of10, for an easy comparison. observational consequences. However, we note that the re-alignment is not complete but stops at a finite value of θ .The evolution of θ can be explained in terms of a simpleanalytical model. Let l B = L B /µ be the angular momentumper unit mass of the binary, with µ = M M / ( M + M ) thereduced mass, and l ⋆ = L ⋆ /m ⋆ the angular momentum perunit mass of a star interacting with the binary. Conservationof angular momentum during an encounter (assuming m ⋆ ≪ M + M ) gives ∆ l B = m ⋆ µ ∆ l ⋆ . (2)Considering that, once normalized per unit mass, theaverage change per encounter in the stellar angular momen-tum is of the order of the binary angular momentum (Merritt2002) , each interaction contributes, on average, a factor < ∆ l B > ∼ m ⋆ µ l B (3)to the binary angular momentum. The rate at which l B varies with time isd l B d t = d n enc d t < ∆ l B > , (4)where d n enc / d t is the rate of binary/star interactions. The This allows us to ignore the evolution of the magnitude of thebinary angular momentum due to hardening when modelling theevolution of the plane orientation. The change in θ would be thesame even if we considered the decay of the binary orbit.c (cid:13) , 000–000 Gualandris et al. L cusp M M L B,1 θ θ L B,0 L cusp L B,0 ∆ L θ θ L B,1 L B,1 L cusp θ θ L B,0 ∆ L ∆ L (a) (b) (c) Figure 5.
Sketch of the model explaining the re-alignment of thebinary plane. L B , is the initial angular momentum of the binary(perpendicular to the orbital plane, shown for clarity in the leftpanel). L B , is the angular momentum of the binary after re-alignment. ∆ L is the amount of angular momentum transferredat any interaction with a single star (dashed line). Note that weassume ∆ L parallel to the total angular momentum L cusp ofthe cusp. θ and θ indicate, respectively, the initial and finalangles between L B and L cusp . Panel (a) shows a case with L B , partially aligned with L cusp ( θ < ◦ ). Panels (b) and (c) referto the opposite case ( θ > ◦ ). Vectors are not on a physicalscale. rate can be estimated asd n enc d t = π n ∗ R , v rel ≈ π n ∗ (cid:18) GM v (cid:19) v rel ≈ π n ∗ (cid:18) M M (cid:19) p GM a , (5)where n ∗ is the number density of stars in the cusp, R i , isthe influence radius of the secondary hole, and v rel is themean relative velocity between the stars and M . Note thatwe have implicitly assumed that: i) only stars efficiently in-teracting with the secondary modify the angular momentumof the binary; ii) the dynamics outside the influence radiusof the secondary is completely dominated by the primary .In order to apply Eqs. 2-5, we need to specify the aver-age direction of < ∆ l B > . This is opposite to the directionof the variation of the angular momentum of the star duringthe interaction. In a reference frame in which the z directionis aligned with the total angular momentum of the cusp, theaverage x and y components of the angular momentum ofthe stars before an encounter are zero. In this simple modelwe assume that the final angular momentum of a star thatexperienced an interaction is isotropically distributed . Thisis justified by the observation that, at any time, stars afteran interaction with the binary are much more isotropicallydistributed than before, as shown in Fig. 4 for the popu-lation of unbound stars. As a consequence, the variation < ∆ l ∗ > = < l ∗ , − l ∗ , > points in the direction oppositeto the angular momentum of the cusp l cusp , and < ∆ l B > ,averaged over repeated encounters, points in the z direction,parallel to l cusp . The fast alignment that we observe in the N -body simulations is a direct consequence of such prefer- ential direction of ∆ l B . Fig. 5 shows a schematic view of ourmodel, for different initial orientations of the binary.The evolution of θ predicted by our model is shown withdotted lines in Fig. 3. Naturally, the case with l B alignedwith l cusp does not show any evolution. Despite the simpli-fying assumptions made, the agreement between the initialevolution of θ in the runs and in the model is very good.The predictions of the model and the results of the N -body runs disagree at late times. In the simulationsthe binary stops re-orienting well before perfect alignmentis reached, while the model predicts θ → l max that can be transferred from the stars to the binary .Since in our model we assumed the transferred angular mo-mentum to depend on the properties of the cusp and noton those of the binary (the stars are ejected isotropically), l max must be the same in every run. We computed the valueof l max that needs to be transferred to the binary to matchexactly the evolution of the θ = 120 ◦ run, and we checked ifsuch a limiting angular momentum can account for the otherruns’ behaviour. The results of this test are shown in Fig. 3,with dashed lines. The agreement with the N -body results isgreatly improved by the introduction of a maximum angularmomentum that can be transferred to the binary.As a final note, we emphasize that the re-alignment inrotating cusps does not depend on the m ∗ /M BHB ratio, asopposed to the non-rotating case. Since the rate of encoun-ters under every assumption depends on the number densityof stars, and the amount of angular momentum transferredper encounter is proportional to m ∗ , d l B / d t depends on themass density ρ of the cusp, but does not show any otherdependence on m ∗ /M BHB . This is a consequence of the factthat the alignment in each encounter has a preferential direc-tion (i.e. towards the angular momentum of the cusp). Theupper panel in Fig. 2 shows the alignment ∆ θ experiencedover an hardening time by binaries in rotating cusps, as afunction of the number of stars in the cusp. Different linesare for different initial angles between l B and l cusp . Whilein the non-rotating case the ratio between the re-alignmentsin the lowest and highest mass resolution runs is ∼
30, inthis case the re-alignment shows only statistical fluctuations,with no systematic trends.
We studied the evolution of the orientation of the orbitalplane of unequal mass BHBs in isotropic and rotating stel-lar systems. The evolution in isotropic systems is character-ized by a random walk on long time-scales, as predicted byMerritt (2002). The evolution in systems in which a frac-tion of stars is counter-rotating with respect to the binarywas briefly discussed in Sesana et al. (2011). If the fractionof counter-rotating stars is small, the change in the orbital c (cid:13) , 000–000 assive black hole binary plane reorientation plane is consistent with the change expected for isotropiccusps. If, on the other hand, a large fraction of stars iscounter-rotating, the plane of the binary evolves consider-ably and the binary angular momentum even reverses in thecase of a fully counter-rotating cusp.Here, we investigated for the first time the evolution ofthe binary plane in rotating systems in which the binaryangular momentum is initially misaligned with respect tothat of the stellar cusp. We find that mis-aligned binaries inrotating cusps tend to align their angular momentum withthat of the cusp on a time-scale of a few hardening times.The alignment process is linear with time, hence, unlike therandom walk observed in fully isotropic systems, does notdepend on the mass ratio between the MBHs and the stars.The realignment of the binary plane that we observe inthe simulations may have significant implications. The direc-tion of the spin axis of the single MBH that results from themerger of the holes due to emission of gravitational wavesis affected by the orientation of the binary plane before co-alescence. The spin axis, in turn, determines the orientationof the accretion disk around the remnant black hole via theBardeen-Peterson effect (Bardeen & Petterson 1975) and, inradio-loud AGNs, the direction of the radio jet.After the coalescence of a BHB, the spin of the rem-nant j r is not necessarily aligned with the spins of the twoprogenitors, j and j . If one of the two MBHs in the bi-nary is radio–loud (and the remnant stays radio–loud), thechange in spin direction would result in a reorientation ofthe relativistic jet (Merritt & Ekers 2002). The orientationof the orbital plane of the binary has implications for the di-rection of j r . The interaction with a rotating cusp tends toforce the binary to co-rotate with the central cusp. At coa-lescence, the binary orbital angular momentum contributionto j r will result in a spin (and possibly a jet) preferentiallyaligned with the cusp angular momentum.We do not expect full alignment, however, because: (i)the alignment we find in the simulations is not complete, andcould halt in the absence of further refilling of stars or in thepresence of isotropic refilling, leaving misalignment angles ofup to 60 ◦ , (ii) in a very unequal-mass binary (e.g. q = 1 / j r depends strongly on j . Using theresults of Rezzolla et al. (2008), if j ≈
1, the angle between j and j r is of only a few degrees, while j ≈ . ∼ > ◦ .Such a preferential alignment is in fact observed. Inweak radio-loud AGNs, Browne & Battye (2010) found atendency for the axis of the radio emission to align with theminor axis of the starlight of the host, an oblate rotationallysupported elliptical. By contrast, they found no preferredradio-optical alignment among the radio-louder objects, (seealso Saripalli & Subrahmanyan 2009), possibly hosted in tri-axial non-rotating ellipticals (e.g. Kormendy et al. 2009). This is valid in gas poor mergers. In gas rich environments, gasaccretion causes the spins of the MBHs to align with the angu-lar momentum of the binary (Bogdanovi´c et al. 2007; Dotti et al.2010; Volonteri et al. 2010; Kesden et al. 2010), preventing anydrastic change in the spins direction at coalescence.
ACKNOWLEDGMENTS
We thank Pau Amaro-Seoane and David Merritt for inter-esting discussions and the anonymous referee for useful com-ments on the manuscript.
REFERENCES
Amaro-Seoane P., Eichhorn C., Porter E. K., Spurzem R.,2010, MNRAS, 401, 2268Armitage P. J., Natarajan P., 2002, ApJL, 567, L9Armitage P. J., Natarajan P., 2005, ApJ, 634, 921Bardeen J. M., Petterson J. A., 1975, ApJL, 195, L65+Baumgardt H., Gualandris A., Portegies Zwart S., 2006,MNRAS, 372, 174Begelman M. C., Blandford R. D., Rees M. J., 1980, Na-ture, 287, 307Berczik P., Merritt D., Spurzem R., 2005, ApJ, 633, 680Berczik P., Merritt D., Spurzem R., Bischof H.-P., 2006,ApJL, 642, L21Berentzen I., Preto M., Berczik P., Merritt D., Spurzem R.,2009, ApJ, 695, 455Bogdanovi´c T., Reynolds C. S., Miller M. C., 2007, ApJL,661, L147Browne I. W. A., Battye R. A., 2010, in L. Maraschi,G. Ghisellini, R. Della Ceca, & F. Tavecchio ed., Accretionand Ejection in AGN: a Global View Vol. 427 of Astro-nomical Society of the Pacific Conference Series, A Di-chotomy in Radio Jet Orientations in Elliptical Galaxies.pp 365–+Cuadra J., Armitage P. J., Alexander R. D., BegelmanM. C., 2009, MNRAS, 393, 1423Dotti M., Colpi M., Haardt F., 2006, MNRAS, 367, 103Dotti M., Colpi M., Haardt F., Mayer L., 2007, MNRAS,379, 956Dotti M., Ruszkowski M., Paredi L., Colpi M., VolonteriM., Haardt F., 2009, MNRAS, 396, 1640Dotti M., Volonteri M., Perego A., Colpi M., RuszkowskiM., Haardt F., 2010, MNRAS, 402, 682Escala A., Larson R. B., Coppi P. S., Mardones D., 2004,ApJ, 607, 765Escala A., Larson R. B., Coppi P. S., Mardones D., 2005,ApJ, 630, 152Gaburov E., Harfst S., Portegies Zwart S., 2009, New A,14, 630Gualandris A., Merritt D., 2007, ArXiv e-printsGualandris A., Merritt D., 2011, ArXiv e-printsHarfst S., Gualandris A., Merritt D., Spurzem R., PortegiesZwart S., Berczik P., 2007, New Astronomy, 12, 357Kesden M., Sperhake U., Berti E., 2010, ApJ, 715, 1006Khan F. M., Just A., Merritt D., 2011, ApJ, 732, 89Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009,ApJS, 182, 216Lodato G., Nayakshin S., King A. R., Pringle J. E., 2009,MNRAS, 398, 1392Makino J., Funato Y., 2004, ApJ, 602, 93Merritt D., 2002, ApJ, 568, 998Merritt D., 2006, ApJ, 648, 976Merritt D., Ekers R. D., 2002, Science, 297, 1310Merritt D., Mikkola S., Szell A., 2007, ApJ, 671, 53Milosavljevi´c M., Merritt D., 2001, ApJ, 563, 34 c (cid:13) , 000–000 Gualandris et al.
Milosavljevi´c M., Merritt D., 2003, 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, L26+Rezzolla L., Barausse E., Dorband E. N., Pollney D., Reis-swig C., Seiler J., Husa S., 2008, Phys. Rev. D, 78, 044002Roedig C., Dotti M., Sesana A., Cuadra J., Colpi M., 2011,MNRAS, pp 979–+Saripalli L., Subrahmanyan R., 2009, ApJ, 695, 156Sesana A., 2010, ApJ, 719, 851Sesana A., Gualandris A., Dotti M., 2011, MNRAS, 415,L35Sesana A., Haardt F., Madau P., 2007, ApJ, 660, 546Volonteri M., G¨ultekin K., Dotti M., 2010, MNRAS, 404,2143 c (cid:13)000